Broyden.hh

Go to the documentation of this file.

/* * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *\
 * Copyright (c) 2025, Davide Stocco, Mattia Piazza and Enrico Bertolazzi.                       *
 *                                                                                               *
 * The Optimist project is distributed under the BSD 2-Clause License.                           *
 *                                                                                               *
 * Davide Stocco                          Mattia Piazza                        Enrico Bertolazzi *
 * University of Trento               University of Trento                  University of Trento *
 * davide.stocco@unitn.it            mattia.piazza@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 Real, Integer N>
    class Broyden : public QuasiNewton<Real, N, Broyden<Real, N>>
    {
    public:
      static constexpr bool requires_function{true};
      static constexpr bool requires_first_derivative{true};
      static constexpr bool requires_second_derivative{false};
 
      OPTIMIST_BASIC_CONSTANTS(Real) 
 
      using Method = enum class Method : Integer {GOOD = 0, BAD = 1, COMBINED = 2}; 
      using Vector = typename QuasiNewton<Real, N, Broyden<Real, N>>::Vector;
      using Matrix = typename QuasiNewton<Real, N, Broyden<Real, N>>::Matrix;
      using FunctionWrapper = typename QuasiNewton<Real, N, Broyden<Real, N>>::FunctionWrapper;
      using JacobianWrapper = typename QuasiNewton<Real, N, Broyden<Real, N>>::JacobianWrapper;
      using QuasiNewton<Real, N, Broyden<Real, N>>::solve;
 
    private:
      Method m_method{Method::COMBINED}; 
 
    public:
      Broyden() {}

      std::string name_impl() const
      {
        std::ostringstream os;
        os << "Broyden";
        if (this->m_method == Method::GOOD) {
          os << "Good";
        } else if (this->m_method == Method::BAD) {
          os << "Bad";
        } else if (this->m_method == Method::COMBINED) {
          os << "Combined";
        }
        return os.str();
      }

      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(
        Vector const & delta_x_old, Vector const & delta_function_old, Matrix const & jacobian_old,
        Vector const & delta_x_new, Vector const & delta_function_new, Vector const & /*function_new*/,
        Matrix       & jacobian_new
      ) {
        Vector tmp_1(jacobian_old * delta_function_new);
        Real 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.transpose() * delta_x_old) / std::abs(delta_x_new.transpose() * tmp_1)
            < std::abs(delta_function_new.transpose() * 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.transpose() * 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