mqc_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 MQCLinesearch : public Linesearch<Functional_, Filter_>
    {
      public:
        typedef Filter_ FilterType;
        typedef typename Functional_::VectorTypeR VectorType;
        typedef typename Functional_::DataType DataType;
        typedef Linesearch<Functional_, Filter_> BaseClass;
 
      protected:
        DataType _alpha_hard_max;
        DataType _alpha_hard_min;
        DataType _alpha_soft_max;
        DataType _alpha_soft_min;
 
      public:
        explicit MQCLinesearch(Functional_& functional, Filter_& filter, bool keep_iterates = false) :
          BaseClass("MQC-LS", functional, filter, keep_iterates),
          _alpha_hard_max(DataType(0)),
          _alpha_hard_min(DataType(0)),
          _alpha_soft_max(DataType(0)),
          _alpha_soft_min(DataType(0))
          {
          }
 
        explicit MQCLinesearch(const String& section_name, const PropertyMap* section,
        Functional_& functional, Filter_& filter) :
          BaseClass("MQC-LS", section_name, section, functional, filter),
          _alpha_hard_max(DataType(0)),
          _alpha_hard_min(DataType(0)),
          _alpha_soft_max(DataType(0)),
          _alpha_soft_min(DataType(0))
          {
          }
 
        virtual ~MQCLinesearch()
        {
        }
 
        virtual String name() const override
        {
          return "MQCLinesearch";
        }
 
        virtual void reset() override
        {
          BaseClass::reset();
          _alpha_hard_max = DataType(0);
          _alpha_hard_min = DataType(0);
          _alpha_soft_max = DataType(0);
          _alpha_soft_min = DataType(0);
        }
 
        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()));
 
          static constexpr DataType extrapolation_width = DataType(4);
          Status status(Status::progress);
 
          // 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);
          //alpha = this->_alpha_0;
          if(st != Status::progress)
          {
            return st;
          }
 
          // Scaling if we are to use step sizes wrt. to the non-normalized search direction
          // This appears to be the right thing theoretically, but in practice not using the search direction norm as
          // initial guess for the step size works better. Stupid reality...
          //if(this->_dir_scaling)
          //{
          //  alpha = this->_norm_dir;
          //}
 
          // Set hard limits to default values if they have not been set
          _alpha_hard_min = DataType(0);
          // This bogus value should be around 1e50 for double precision. It is chosen to make comparing results with
          // ALGLIB easier
          if(_alpha_hard_max < Math::eps<DataType>())
          {
            _alpha_hard_max = Math::pow(Math::huge<DataType>(), DataType(0.1622));
          }
 
          // Set the intervall of uncertainty
          _alpha_soft_min = DataType(0);
          _alpha_soft_max = this->_alpha_0 + extrapolation_width*this->_alpha_0;
          _alpha_soft_max = Math::min(_alpha_soft_max, _alpha_hard_max);
 
          //std::cout << "Linesearch initial alpha " << this->_alpha_0 << "\n";
 
          DataType alpha_lo(0);
          // It is critical that _f_0 was set from the outside!
          DataType fval_lo(this->_fval_0);
          DataType df_lo(this->_delta_0);
 
          DataType alpha_hi(0);
          DataType fval_hi(this->_fval_0);
          DataType df_hi(this->_delta_0);
 
          // This is the width of the search interval
          DataType width(Math::abs(_alpha_hard_max - _alpha_hard_min));
          DataType width_old(DataType(2)*width);
 
          // Do we know the interval of uncertainty?
          bool interval_known(false);
          // Is the minimum in the interval of uncertainty?
          bool min_in_interval(false);
          // Does the minimum lie outside the current search interval, requiring us to drive the step to its
          // boundary?
          bool drive_to_bndry(false);
 
          // start iterating
          while(status == Status::progress)
          {
 
            IterationStats stat(*this);
            ++(this->_num_iter);
 
            // If we know the minimum is in the search interval, the interval of uncertainty is the search interval
            if(min_in_interval)
            {
              _alpha_soft_min = Math::min(alpha_lo, alpha_hi);
              _alpha_soft_max = Math::max(alpha_lo, alpha_hi);
            }
            // Enlarge the interval of uncertainty
            else
            {
              _alpha_soft_min = alpha_lo;
              _alpha_soft_max = alpha + extrapolation_width*Math::abs(alpha - alpha_lo);
            }
            //std::cout << "Set alpha_smin " << _alpha_soft_min << " alpha_smax " << _alpha_soft_max << "\n";
 
            // Update solution: sol <- initial_sol + _alpha*dir
            vec_sol.copy(this->_vec_initial_sol);
            vec_sol.axpy(this->_vec_pn, alpha);
            //std::cout << "Linesearch alpha " << alpha << "\n";
            //std::cout << "initial_sol " << *(this->_vec_initial_sol) << "\n";
            //std::cout << "dir " << *this->_vec_pn << "\n";
            //std::cout << "sol " << *vec_sol << "\n";
 
            // Prepare and evaluate
            this->_functional.prepare(vec_sol, this->_filter);
 
            // Compute and filter the gradient
            this->_functional.eval_fval_grad(fval, this->_vec_grad);
            this->trim_func_grad(fval);
            this->_filter.filter_def(this->_vec_grad);
 
            // New directional derivative and new defect
            df = this->_vec_pn.dot(this->_vec_grad);
            status = this->_check_convergence(fval, df, alpha);
 
            // If success is reported, check if it is a real success or if something fishy is going on
            if(status == Status::success)
            {
              this->_vec_tmp.copy(this->_vec_initial_sol);
              this->_vec_tmp.axpy(vec_sol, -DataType(1));
              if(fval >= this->_fval_0 || this->_vec_tmp.norm2() == DataType(0))
              {
                status = Status::stagnated;
              }
            }
            // If we are not successful, check if the interval of uncertainty has become too small
            else if(min_in_interval && (_alpha_soft_max - _alpha_soft_min) <= this->_tol_step*_alpha_soft_max)
            {
              //std::cout << "interval width " << _alpha_soft_max - _alpha_soft_min << " : " << this->_tol_step*_alpha_soft_max << "\n";
              status = Status::stagnated;
            }
            // Stagnation due to rounding errors
            //if(min_in_interval && (alpha <= _alpha_soft_min || alpha >= _alpha_soft_max))
            //{
            //  std::cout << "Rounding errors\n";
            //  status = Status::stagnated;
            //}
            // This is not used at the moment because it is only relevant if there are constraints limiting the
            // step length
            //if( alpha == _alpha_hard_max)
            //if( alpha == _alpha_hard_min)
 
            if(status != Status::progress)
            {
              break;
            }
 
            if(!interval_known
                && (fval < this->_fval_0 + this->_tol_decrease*alpha*this->_delta_0)
                && (df >= Math::min(this->_tol_decrease, this->_tol_curvature)))
                {
                  interval_known = true;
                }
 
            //std::cout << "interval known " << interval_known << "\n";
            // If we do not know that the minimum was in the previous interval of uncertainty, we need to compute a
            // new step size to expand the interval of uncertainty at the start of the next iteration.
            if(!interval_known && (fval <= fval_lo)
                && fval > this->_fval_0 + alpha*this->_tol_decrease*this->_delta_0)
            {
              DataType fval_m(fval - alpha*this->_tol_decrease*this->_delta_0);
              DataType df_m(df - this->_tol_decrease*this->_delta_0);
              DataType fval_lo_m(fval_lo - alpha_lo*this->_tol_decrease*this->_delta_0);
              DataType df_lo_m(df_lo - this->_tol_decrease*this->_delta_0);
              DataType fval_hi_m(fval_hi - alpha_hi*this->_tol_decrease*this->_delta_0);
              DataType df_hi_m(df_hi - this->_tol_decrease*this->_delta_0);
 
              // Note that the expansion step might already give us the information that the minimum is the the
              // new search interval
              _polynomial_fit(
                alpha, fval_m, df_m, alpha_lo, fval_lo_m, df_lo_m,
                alpha_hi, fval_hi_m, df_hi_m, min_in_interval, drive_to_bndry);
 
              fval_lo = fval_lo_m + alpha_lo*this->_tol_decrease*this->_delta_0;
              df_lo = df_lo_m + this->_tol_decrease*this->_delta_0;
              fval_hi = fval_hi_m * alpha_hi*this->_tol_decrease*this->_delta_0;
              df_hi = df_hi_m + this->_tol_decrease*this->_delta_0;
 
            }
            else
            {
              _polynomial_fit(alpha, fval, df, alpha_lo, fval_lo, df_lo, alpha_hi, fval_hi, df_hi, min_in_interval,
              drive_to_bndry);
 
              //std::cout << "width " << width << " width_old " << width_old << " min_in_interval " << min_in_interval <<"\n";
 
              if(min_in_interval)
              {
                if(Math::abs(alpha_hi - alpha_lo) >= DataType(0.66)*width_old)
                {
                  //std::cout << "Forcing " << alpha << " to ";
                  alpha = alpha_lo + DataType(0.5)*(alpha_hi - alpha_lo);
                  //std::cout << alpha << "\n";
                }
                width_old = width;
                width = Math::abs(alpha_lo - alpha_hi);
              }
            }
 
          } // while(status == Status:progress)
 
          this->_alpha_min = alpha;
          this->_fval_min = fval;
 
          // If we are successful, we save the last step length as the new initial step length
          if(status == Status::success)
          {
            this->_alpha_0 = this->_alpha_min;
          }
          // If we are not successful, we update the best step length and need to re-evaluate everything for that
          // step
          else
          {
            this->_alpha_min = alpha_lo;
            this->_fval_min = fval_lo;
            //std::cout << "Unusual termination alpha_lo " << alpha_lo << "fval_lo " << fval_lo << "\n";
            vec_sol.copy(this->_vec_initial_sol);
            vec_sol.axpy(this->_vec_pn, alpha_lo);
 
            // Prepare and evaluate
            this->_functional.prepare(vec_sol, this->_filter);
            this->_functional.eval_fval_grad(fval, this->_vec_grad);
            this->trim_func_grad(fval);
            this->_filter.filter_def(this->_vec_grad);
          }
 
          Statistics::add_solver_expression(std::make_shared<ExpressionEndSolve>(this->name(), status, this->get_num_iter()));
          return status;
        }
 
        Status _polynomial_fit(
          DataType& alpha, DataType& fval, DataType& df,
          DataType& alpha_lo, DataType& fval_lo, DataType& df_lo,
          DataType& alpha_hi, DataType& fval_hi, DataType& df_hi,
          bool& min_in_interval, bool& drive_to_bndry)
        {
          DataType alpha_c(0);
          DataType alpha_q(0);
          DataType alpha_new(0);
 
          // Case 1: The function value increased.
          if( fval > fval_lo)
          {
            // Computation of cubic step
            alpha_c  = _argmin_cubic(alpha_lo, fval_lo, df_lo, alpha, fval, df);
            // The first derivative was negative, but the function value increases, so the minimum has to be in this
            // interval (or have been in the original interval, which this interval is a subset of)
            min_in_interval = true;
            drive_to_bndry = true;
            // Use 2nd order approximation NOT using the derivative df because it might be garbage (i.e. if we hit a
            // singularity) because the function value increased
            alpha_q = _argmin_quadratic(alpha_lo, fval_lo, df_lo, alpha, fval, df, false);
 
            // If the cubic minimum is closer to the old best step length, take it. Otherwise, take the mean
            if(Math::abs(alpha_c - alpha_lo) < Math::abs(alpha_q - alpha_lo))
              alpha_new = alpha_c;
            else
              alpha_new = (alpha_q + alpha_c)/DataType(2);
            //std::cout << "Case 1: alpha " << alpha_new << " q " << alpha_q << " c " << alpha_c << "\n";
 
          }
          // Case 2: The first derivative changes sign
          // We also know that the function value increased in the last trial step
          else if( df_lo*df < DataType(0) )
          {
            alpha_c  = _argmin_cubic(alpha, fval, df, alpha_lo, fval_lo, df_lo);
            // Because the derivative changed sign, it has to be zero somewhere in between, so the minimum is in the
            // current interval
            min_in_interval = true;
            drive_to_bndry = false;
            // Default: Quadratic interpolant using fval_lo, df_lo and df
            alpha_q = _argmin_quadratic(alpha_lo, fval_lo, df_lo, alpha, fval, df, true);
            // Take the step closer to the new step alpha
            Math::abs(alpha - alpha_c) > Math::abs(alpha - alpha_q) ? alpha_new = alpha_c : alpha_new = alpha_q;
            //std::cout << "Case 2: alpha " << alpha_new << " q " << alpha_q << " c " << alpha_c << "\n";
          }
          // Case 3: The absolute value of the derivative increases
          // We also know that the function value increased in the last trial step and that the derivative did not
          // change sign.
          // This means we have to further search in the direction of alpha.
          else if(Math::abs(df) < Math::abs(df_lo))
          {
            alpha_c  = _argmin_cubic_case_3(alpha, fval, df, alpha_lo, fval_lo, df_lo);
            drive_to_bndry = true;
 
            // Quadratic interpolant using fval_lo, df_lo and df
            alpha_q = _argmin_quadratic(alpha, fval, df, alpha_lo, fval_lo, df_lo, true);
 
            // If the cubic step is closer to the trial step
            if(Math::abs(alpha - alpha_c) < Math::abs(alpha - alpha_q))
            {
              min_in_interval ? alpha_new = alpha_c : alpha_new = alpha_q;
            }
            else
            {
              min_in_interval ? alpha_new = alpha_q : alpha_new = alpha_c;
            }
            //std::cout << "Case 3: alpha " << alpha_new << " q " << alpha_q << " c " << alpha_c << "\n";
          }
          // Case 4: The absolute value of the derivative did not increase
          // We also know that the function value increased in the last trial step and that the derivative did not
          // change sign.
          else
          {
            drive_to_bndry = false;
            // Note that the arguments for argmin_cubic are different in this case
            // Here we can just take the cubic step because the cubic interpolation sets it to the alpha_soft_max
            // if it recognises that the minimum is outside the current interval.
            if(min_in_interval)
            {
              alpha_c  = _argmin_cubic(alpha, fval , df, alpha_hi, fval_hi, df_hi);
              alpha_new = alpha_c;
            }
            else
            {
              alpha > alpha_lo ? alpha_new = _alpha_soft_max : alpha_new = _alpha_soft_min;
            }
            //std::cout << "Case 4: alpha " << alpha_new << "\n";
          }
 
            // Update the inverval of uncertainty. Has to happen befor we clamp the step to the admissible interval.
            // Check which step(s) we need to replace
            if(fval > fval_lo)
            {
              alpha_hi = alpha;
              fval_hi = fval;
              df_hi = df;
            }
            else
            {
              if(df * df_lo < DataType(0))
              {
                alpha_hi = alpha_lo;
                fval_hi = fval_lo;
                df_hi = df_lo;
              }
              alpha_lo = alpha;
              fval_lo = fval;
              df_lo = df;
            }
 
            _clamp_step(alpha_new);
 
            // If we know the minimum is in the interval, we can drive alpha to the corresponding end more quickly
            if(min_in_interval && drive_to_bndry)
            {
              //std::cout << "Bracketing adjust " << alpha_lo << " " << alpha_hi << " " << alpha_new << "\n";
              //std::cout << "Bracketing adjust alpha_new from " << alpha_new;
              if(alpha_lo < alpha_hi)
                alpha_new = Math::min(alpha_lo + DataType(0.66)*(alpha_hi - alpha_lo), alpha_new);
              else
                alpha_new = Math::max(alpha_lo + DataType(0.66)*(alpha_hi - alpha_lo), alpha_new);
              //std::cout << " to " << alpha_new << "\n";
            }
 
            alpha = alpha_new;
 
            //std::cout << "Polynomial fit: new " << alpha << " " << fval << " " << df <<"\n";
            //std::cout << "Polynomial fit: lo  " << alpha_lo << " " << fval_lo << " " << df_lo <<"\n";
            //std::cout << "Polynomial fit: hi  " << alpha_hi << " " << fval_hi << " " << df_hi <<"\n";
 
          return Status::success;
        }
 
        void _clamp_step(DataType& alpha_new) const
        {
          alpha_new = Math::max(alpha_new, _alpha_hard_min);
          alpha_new = Math::min(alpha_new, _alpha_hard_max);
 
          alpha_new = Math::max(alpha_new, _alpha_soft_min);
          alpha_new = Math::min(alpha_new, _alpha_soft_max);
        }
 
        DataType _argmin_quadratic(const DataType alpha_lo, const DataType fval_lo, const DataType df_lo, const DataType alpha_hi, const DataType fval_hi, const DataType df_hi, const bool interpolate_derivative) const
        {
          DataType alpha(alpha_lo);
 
          // Quadratic interpolation using fval_lo, df_lo, df_hi
          if(interpolate_derivative)
            alpha += df_lo /(df_lo - df_hi) * (alpha_hi - alpha_lo);
          // Quadratic interpolation using fval_lo, fval_hi, df_lo
          else
            alpha += df_lo/( (fval_hi - fval_lo)/( alpha_hi - alpha_lo) - df_lo)/DataType(2)*(alpha_lo - alpha_hi);
 
          return alpha;
        }
 
        DataType _argmin_cubic(DataType alpha_lo, DataType fval_lo, DataType df_lo,
        DataType alpha_hi, DataType fval_hi, DataType df_hi) const
        {
          DataType alpha(alpha_lo);
 
          DataType d1 = DataType(3)*(fval_lo - fval_hi)/(alpha_hi - alpha_lo) + df_lo + df_hi;
          // Scale the computation of r for better numerical stability
          DataType scale = Math::max(Math::abs(df_lo), Math::abs(df_hi));
          scale = Math::max(scale, Math::abs(d1));
 
          DataType r = Math::sqr(d1/scale) - (df_lo/scale) * (df_hi/scale);
 
          //std::cout << "fval_hi " << fval_hi << " fval_lo "<< fval_lo << "\n";
          //std::cout << "alpha_hi " << alpha_hi << " alpha_lo "<< alpha_lo << "\n";
          //std::cout << "df_hi " << df_hi << " df_lo "<< df_lo << "\n";
          //std::cout << "d1 " << d1 << "\n";
          //std::cout << "scale " << scale << "\n";
 
          DataType d2(0);
 
          if(r > DataType(0))
          {
            d2 = Math::signum(alpha_hi - alpha_lo) * scale * Math::sqrt(r);
 
            DataType p(d2 - df_lo + d1);
 
            DataType q(0);
            // This sorting is done to avoid annihilation
            df_hi*df_lo > DataType(0) ?  q = d2 +(df_hi-df_lo) + d2 : q = d2 - df_lo + d2 + df_hi;
 
            //std::cout << "d2 " << d2 << "\n";
            //std::cout << "p " << p << " q " << q << "\n";
 
            alpha += (alpha_hi - alpha_lo)*(p/q);
          }
          // r <= 0 means that the minimum is not in the interior, so it has to lie across the endpoint with the
          // lower function value
          else
          {
            //d2 = Math::signum(alpha_hi - alpha_lo) * scale *
            //  Math::sqrt(Math::max(r, DataType(0)));
            //DataType q(0);
            //if(df_hi*df_lo > 0)
            //  q = d2 +(df_hi-df_lo) + d2;
            //else
            //  q = d2 - df_lo + d2 +df_hi;
            //std::cout << "d2 " << d2 << "\n";
            //std::cout << "p " << (d2 - df_lo + d1) << " q " << q << "\n";
 
            (alpha_lo < alpha_hi && df_lo > DataType(0) )? alpha = _alpha_soft_min: alpha = _alpha_soft_max;
          }
 
          return alpha;
        }
 
        DataType _argmin_cubic_case_3(
          const DataType alpha_lo, const DataType fval_lo, const DataType df_lo,
          const DataType alpha_hi, const DataType fval_hi, DataType df_hi) const
        {
          DataType alpha_c(alpha_lo);
 
          if(alpha_lo == alpha_hi)
          {
            return alpha_lo;
          }
 
          DataType d1 = DataType(3)*(fval_hi - fval_lo)/(alpha_lo - alpha_hi) + df_hi + df_lo;
          DataType scale = Math::max(Math::abs(df_lo), Math::abs(df_hi));
          scale = Math::max(scale, Math::abs(d1));
 
          DataType r = Math::sqr(d1/scale) - (df_lo/scale) * (df_hi/scale);
 
          DataType d2 = Math::signum(alpha_hi - alpha_lo) * scale * Math::sqrt(Math::max(r, DataType(0)));
 
          DataType p(d2 - df_lo + d1);
          DataType q(d2 +(df_hi - df_lo) + d2);
 
          //std::cout << "fval_hi " << fval_hi << " fval_lo "<< fval_lo << "\n";
          //std::cout << "alpha_hi " << alpha_hi << " alpha_lo "<< alpha_lo << "\n";
          //std::cout << "df_hi " << df_hi << " df_lo "<< df_lo << "\n";
          //std::cout << "d1 " << d1 << "\n";
          //std::cout << "scale " << scale << "\n";
          //std::cout << "d2 " << d2 << "\n";
          //std::cout << "p " << (d2 - df_lo + d1) << " q " << q << "\n";
 
          if( p/q < DataType(0) && d2 != DataType(0))
            alpha_c += p/q*(alpha_hi - alpha_lo);
          else
          {
            //std::cout << "Weird ";
            alpha_lo > alpha_hi ? alpha_c = _alpha_soft_max : alpha_c = _alpha_soft_min;
          }
 
          return alpha_c;
        }
 
    }; // class MQCLinesearch
 
    template<typename Functional_, typename Filter_>
    inline std::shared_ptr<MQCLinesearch<Functional_, Filter_>> new_mqc_linesearch(
      Functional_& functional, Filter_& filter, bool keep_iterates = false)
      {
        return std::make_shared<MQCLinesearch<Functional_, Filter_>>(functional, filter, keep_iterates);
      }
 
    template<typename Functional_, typename Filter_>
    inline std::shared_ptr<MQCLinesearch<Functional_, Filter_>> new_mqc_linesearch(
      const String& section_name, const PropertyMap* section, Functional_& functional, Filter_& filter)
      {
        return std::make_shared<MQCLinesearch<Functional_, Filter_>> (section_name, section, functional, filter);
      }
 
  } // namespace Solver
} // namespace FEAT