Newton.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_NEWTON_HH
#define OPTIMIST_ROOTFINDER_NEWTON_HH
 
#include "Optimist/RootFinder.hh"
 
namespace Optimist {
  namespace RootFinder {
 
    /*\
     |   _   _               _
     |  | \ | | _____      _| |_ ___  _ __
     |  |  \| |/ _ \ \ /\ / / __/ _ \| '_ \
     |  | |\  |  __/\ V  V /| || (_) | | | |
     |  |_| \_|\___| \_/\_/  \__\___/|_| |_|
     |
    \*/

    template <typename Vector>
      requires TypeTrait<Vector>::IsEigen &&
               (!TypeTrait<Vector>::IsFixed || TypeTrait<Vector>::Dimension > 0)
    class Newton : public RootFinder<Vector, Newton<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 RootFinder<Vector, Newton<Vector>>::FirstDerivative;
      using Factorization =
          std::conditional_t<VectorTrait::IsSparse,
                             Eigen::SparseLU<FirstDerivative>,
                             Eigen::FullPivLU<FirstDerivative>>;
 
      OPTIMIST_BASIC_CONSTANTS(Scalar)
 
     private:
      Factorization m_lu; 
 
     public:
      Newton() {}

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

      template <typename FunctionLambda, typename JacobianLambda>
      bool solve_impl(FunctionLambda &&function,
                      JacobianLambda &&jacobian,
                      const Vector &x_ini,
                      Vector &x_sol) {
#define CMD "Optimist::RootFinder::Newton::solve(...): "
 
        // Reset internal variables
        this->reset_counters();
 
        // Print header
        if (this->m_verbose) {
          this->header();
        }
 
        // Initialize variables
        bool damped, success;
        Scalar residuals, step_norm;
        Vector x_old, x_new, function_old, function_new, step_old, step_new;
        FirstDerivative jacobian_old;
 
        // Set initial iteration
        x_old = x_ini;
        success =
            this->evaluate_function(std::forward<FunctionLambda>(function),
                                    x_old,
                                    function_old);
        OPTIMIST_ASSERT(success,
                        CMD "function evaluation failed at the initial point.");
 
        // Algorithm iterations
        Scalar tolerance_residuals{this->m_tolerance};
        Scalar tolerance_step_norm{this->m_tolerance * this->m_tolerance};
        for (this->m_iterations = 1;
             this->m_iterations < this->m_max_iterations;
             ++this->m_iterations) {
          // Calculate step
          success = this->evaluate_jacobian(jacobian, x_old, jacobian_old);
          OPTIMIST_ASSERT(success,
                          CMD "jacobian evaluation failed at iteration "
                              << this->m_iterations << ".");
          this->m_lu.compute(jacobian_old);
          if constexpr (VectorTrait::IsSparse) {
            OPTIMIST_ASSERT_WARNING(this->m_lu.info() == Eigen::Success,
                                    CMD "jacobian factorization failed.");
            step_old = -this->m_lu.solve(function_old).eval();
          } else {
            OPTIMIST_ASSERT(this->m_lu.isInvertible(),
                            CMD "singular Jacobian detected.");
            step_old = -this->m_lu.solve(function_old);
          }
 
          // Check convergence
          residuals = function_old.norm();
          step_norm = step_old.norm();
          if (this->m_verbose) {
            this->info(residuals);
          }
          if (residuals < tolerance_residuals ||
              step_norm < tolerance_step_norm) {
            this->m_converged = true;
            break;
          }
 
          if (this->m_damped) {
            // Relax the iteration process
            damped = this->damp(std::forward<FunctionLambda>(function),
                                x_old,
                                function_old,
                                step_old,
                                x_new,
                                function_new,
                                step_new);
            OPTIMIST_ASSERT_WARNING(
                damped,
                "Optimist::RootFinder::Newton::solve(...): damping failed.");
          } else {
            // Update point
            x_new = x_old + step_old;
            success =
                this->evaluate_function(std::forward<FunctionLambda>(function),
                                        x_new,
                                        function_new);
            OPTIMIST_ASSERT(success,
                            CMD "function evaluation failed at iteration "
                                << this->m_iterations << ".");
          }
 
          // Update internal variables
          x_old        = x_new;
          function_old = function_new;
          step_old     = step_new;
        }
 
        // Print bottom
        if (this->m_verbose) {
          this->bottom();
        }
 
        // Convergence data
        x_sol = x_old;
        return this->m_converged;
 
#undef CMD
      }
 
    };  // class Newton
 
  }  // namespace RootFinder
 
}  // namespace Optimist
 
#endif  // OPTIMIST_ROOTFINDER_NEWTON_HH