Broyden.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
 * Copyright (c) 2025, Davide Stocco and Enrico Bertolazzi.                  *
 *                                                                           *
 * The Optimist project is distributed under the BSD 2-Clause License.       *
 *                                                                           *
 * Davide Stocco                                           Enrico Bertolazzi *
 * University of Trento                                 University of Trento *
 * davide.stocco@unitn.it                         enrico.bertolazzi@unitn.it *
\* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * */
 
#pragma once
 
#ifndef OPTIMIST_ROOTFINDER_BROYDEN_HH
#define OPTIMIST_ROOTFINDER_BROYDEN_HH
 
#include "Optimist/RootFinder/QuasiNewton.hh"
 
namespace Optimist {
  namespace RootFinder {
 
    /*\
     |   ____                      _
     |  | __ ) _ __ ___  _   _  __| | ___ _ __
     |  |  _ \| '__/ _ \| | | |/ _` |/ _ \ '_ \
     |  | |_) | | | (_) | |_| | (_| |  __/ | | |
     |  |____/|_|  \___/ \__, |\__,_|\___|_| |_|
     |                   |___/
    \*/

    template <typename Vector>
      requires TypeTrait<Vector>::IsEigen &&
               (!TypeTrait<Vector>::IsFixed || TypeTrait<Vector>::Dimension > 0)
    class Broyden : public QuasiNewton<Vector, Broyden<Vector>> {
     public:
      static constexpr bool RequiresFunction{true};
      static constexpr bool RequiresFirstDerivative{true};
      static constexpr bool RequiresSecondDerivative{false};
 
      using VectorTrait = TypeTrait<Vector>;
      using Scalar      = typename TypeTrait<Vector>::Scalar;
      using typename QuasiNewton<Vector, Broyden<Vector>>::FirstDerivative;
 
      OPTIMIST_BASIC_CONSTANTS(Scalar)
 
      
      using Method = enum class Method : Integer {
        GOOD     = 0,
        BAD      = 1,
        COMBINED = 2
      };
 
     private:
      Method m_method{Method::COMBINED}; 
 
     public:
      Broyden() {}

      constexpr std::string name_impl() const {
        return "Broyden";
      }

      Method method() const {
        return this->m_method;
      }

      void method(Method t_method) {
        this->m_method = t_method;
      }

      void enable_good_method() {
        this->m_method = Method::GOOD;
      }

      void enable_bad_method() {
        this->m_method = Method::BAD;
      }

      void enable_combined_method() {
        this->m_method = Method::COMBINED;
      }

      void set_method(Method t_method) {
        this->m_method = t_method;
      }

      void update_impl(const Vector &delta_x_old,
                       const Vector &delta_function_old,
                       const FirstDerivative &jacobian_old,
                       const Vector &delta_x_new,
                       const Vector &delta_function_new,
                       const Vector & /*function_new*/,
                       FirstDerivative &jacobian_new) {
        Vector tmp_1(jacobian_old * delta_function_new);
        Scalar tmp_2{delta_function_new.squaredNorm()};
        // Selection criteria: |(dx_new'*dx_old) / (dx_new'*J_old*dF_new)| <
        // |(dF_new'*dF_old) / (dF_new'*dF_new)|
        if (this->m_method == Method::COMBINED ||
            this->m_method == Method::GOOD || this->iterations() < Integer(2) ||
            std::abs(delta_x_new.dot(delta_x_old)) /
                    std::abs(delta_x_new.dot(tmp_1)) <
                std::abs(delta_function_new.dot(delta_function_old)) / tmp_2) {
          // Broyden's Good solver: J_new = J_old - (J_old*dF_new-dx_new) /
          // (C'*dF_new)*C', where C = J_old'*dx_new;
          Vector C_g(jacobian_old.transpose() * delta_x_new);
          jacobian_new = jacobian_old - (tmp_1 - delta_x_new) /
                                            (C_g.dot(delta_function_new)) *
                                            C_g.transpose();
        } else {
          // Broyden's Bad solver: J_new = J_old - (J_old*dF_old-dx_new) /
          // (C'*dF_old)*C', where C = dF_old;
          jacobian_new = jacobian_old - (tmp_1 - delta_x_new) / tmp_2 *
                                            delta_function_old.transpose();
        }
      }
 
    };  // class Broyden
 
  }  // namespace RootFinder
 
}  // namespace Optimist
 
#endif  // OPTIMIST_ROOTFINDER_BROYDEN_HH