secant_linesearch.hpp

// FEAT3: Finite Element Analysis Toolbox, Version 3
// Copyright (C) 2010 by Stefan Turek & the FEAT group
// FEAT3 is released under the GNU General Public License version 3,
// see the file 'copyright.txt' in the top level directory for details.
 
#pragma once
#include <kernel/base_header.hpp>
#include <kernel/solver/linesearch.hpp>
 
namespace FEAT
{
  namespace Solver
  {
    template<typename Functional_, typename Filter_>
    class SecantLinesearch : public Linesearch<Functional_, Filter_>
    {
      public:
        typedef Filter_ FilterType;
        typedef typename Functional_::VectorTypeR VectorType;
        typedef typename Functional_::DataType DataType;
        typedef Linesearch<Functional_, Filter_> BaseClass;
        static constexpr DataType secant_step_default = DataType(1e-2);
 
      protected:
        DataType _secant_step;
 
      public:
        explicit SecantLinesearch(
          Functional_& functional, Filter_& filter,
          const DataType secant_step = secant_step_default,
          const bool keep_iterates = false) :
          BaseClass("S-LS", functional, filter, keep_iterates),
          _secant_step(secant_step)
          {
            this->_alpha_0 = _secant_step;
          }
 
        explicit SecantLinesearch(const String& section_name, const PropertyMap* section,
        Functional_& functional, Filter_& filter) :
          BaseClass("S-LS", section_name, section, functional, filter),
          _secant_step(secant_step_default)
          {
            auto secant_step_p = section->query("secant_step");
            if(secant_step_p.second)
            {
              set_secant_step(DataType(std::stod(secant_step_p.first)));
            }
 
            this->_alpha_0 = _secant_step;
          }
 
        virtual ~SecantLinesearch()
        {
        }
 
        virtual String name() const override
        {
          return "SecantLinesearch";
        }
 
        virtual void reset() override
        {
          BaseClass::reset();
          this->_alpha_0 = _secant_step;
        }
 
        void set_secant_step(DataType secant_step)
        {
          XASSERT(secant_step > DataType(0));
 
          _secant_step = secant_step;
        }
 
        virtual Status apply(VectorType& vec_cor, const VectorType& vec_dir) override
        {
          // clear solution vector
          vec_cor.format();
          this->_functional.prepare(vec_cor, this->_filter);
 
          // apply
          this->_status = _apply_intern(vec_cor, vec_dir);
          this->plot_summary();
          return this->_status;
        }
 
        virtual Status correct(VectorType& vec_sol, const VectorType& vec_dir) override
        {
          this->_functional.prepare(vec_sol, this->_filter);
 
          // apply
          this->_status = _apply_intern(vec_sol, vec_dir);
          this->plot_summary();
          return this->_status;
        }
 
      protected:
        virtual Status _apply_intern(VectorType& vec_sol, const VectorType& vec_dir)
        {
          Statistics::add_solver_expression(std::make_shared<ExpressionStartSolve>(this->name()));
 
          // The step length wrt. to the NORMALISED search direction
          DataType alpha(0);
          // The functional value
          DataType fval(0);
          DataType df(0);
 
          // Perform initializations and checks
          Status st = this->_startup(alpha, fval, df, vec_sol, vec_dir);
          // Use the additional information about the preconditioned search direction's length?
          if(this->_dir_scaling)
          {
            alpha *= this->_norm_dir;
          }
 
          DataType alpha_update(alpha);
 
          // The second secant point in the first iteration is x + _secant_step * vec_dir
          // The first "other" point for the secant
          vec_sol.copy(this->_vec_initial_sol);
          vec_sol.axpy(this->_vec_pn, alpha);
 
          // start iterating
          while(st == Status::progress)
          {
            IterationStats stat(*this);
 
            // Increase iteration count
            ++this->_num_iter;
 
            this->_functional.prepare(vec_sol, this->_filter);
            this->_functional.eval_fval_grad(fval, this->_vec_grad);
            // DO NOT call trim() here, as this simple linesearch accepts bad steps too, which would then be trimmed
            // and used.
            //this->trim_func_grad(fval);
            this->_filter.filter_def(this->_vec_grad);
 
            if(fval < this->_fval_min)
            {
              this->_fval_min = fval;
              this->_alpha_min = alpha;
            }
 
            // Set new defect, do convergence checks and compute the new _alpha
            DataType df_prev = df;
            df = this->_vec_pn.dot(this->_vec_grad);
 
            st = this->_check_convergence(fval, df, alpha);
 
            // Check if the diffence in etas is too small, thus leading to a huge update of relative size > sqrt(eps)
            if(Math::abs(df - df_prev) < Math::abs(alpha_update*df)*Math::sqrt(Math::eps<DataType>()))
            {
              st = Status::stagnated;
            }
 
            if(st != Status::progress)
            {
              break;
            }
 
            // Update alpha according to secant formula
            alpha_update *= df/(df_prev - df);
            alpha += alpha_update;
 
            // Update the solution
            vec_sol.copy(this->_vec_initial_sol);
            vec_sol.axpy(this->_vec_pn, alpha);
          }
 
          // If we are successful, we could save the last step length as the new initial step length. This is
          // disabled by default, as it can lead to stagnation e.g. for minimising the Rosenbrock function.
          //if(st == Status::success)
          //{
          //  this->_alpha_0 = this->_alpha_min/this->_norm_dir;
          //}
          // If we are not successful, we update the best step length and need to re-evaluate everything for that
          // step
          if(st != Status::success)
          {
            vec_sol.copy(this->_vec_initial_sol);
            vec_sol.axpy(this->_vec_pn, this->_alpha_min);
 
            // Prepare and evaluate
            this->_functional.prepare(vec_sol, this->_filter);
            this->_functional.eval_fval_grad(fval, this->_vec_grad);
            // DO NOT call trim() here, as this simple linesearch accepts bad steps too, which would then be trimmed
            // and used.
            //this->trim_func_grad(fval);
            this->_filter.filter_def(this->_vec_grad);
          }
          Statistics::add_solver_expression(std::make_shared<ExpressionEndSolve>(this->name(), Status::undefined, this->get_num_iter()));
          return st;
        }
 
    }; // class SecantLinesearch
 
    template<typename Functional_, typename Filter_>
    inline std::shared_ptr<SecantLinesearch<Functional_, Filter_>> new_secant_linesearch(
      Functional_& functional, Filter_& filter,
      typename Functional_::DataType secant_step = SecantLinesearch<Functional_, Filter_>::secant_step_default,
      bool keep_iterates = false)
      {
        return std::make_shared<SecantLinesearch<Functional_, Filter_>>(functional, filter, secant_step,
        keep_iterates);
      }
 
    template<typename Functional_, typename Filter_>
    inline std::shared_ptr<SecantLinesearch<Functional_, Filter_>> new_secant_linesearch(
      const String& section_name, const PropertyMap* section,
      Functional_& functional, Filter_& filter)
      {
        return std::make_shared<SecantLinesearch<Functional_, Filter_>>(section_name, section, functional, filter);
      }
 
  } // namespace Solver
} // namespace FEAT