ROL_LineSearchStep.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//               Rapid Optimization Library (ROL) Package
//
// Copyright 2014 NTESS and the ROL contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#ifndef ROL_LINESEARCHSTEP_H
#define ROL_LINESEARCHSTEP_H

#include "ROL_Types.hpp"
#include "ROL_Step.hpp"
#include "ROL_LineSearch.hpp"

// Unconstrained Methods
#include "ROL_GradientStep.hpp"
#include "ROL_NonlinearCGStep.hpp"
#include "ROL_SecantStep.hpp"
#include "ROL_NewtonStep.hpp"
#include "ROL_NewtonKrylovStep.hpp"

// Bound Constrained Methods
#include "ROL_ProjectedSecantStep.hpp"
#include "ROL_ProjectedNewtonStep.hpp"
#include "ROL_ProjectedNewtonKrylovStep.hpp"

#include <sstream>
#include <iomanip>

namespace ROL {

template <class Real>
class LineSearchStep : public Step<Real> {
private:

  ROL::Ptr<Step<Real> >        desc_;       
  ROL::Ptr<Secant<Real> >      secant_;     
  ROL::Ptr<Krylov<Real> >      krylov_;     
  ROL::Ptr<NonlinearCG<Real> > nlcg_;       
  ROL::Ptr<LineSearch<Real> >  lineSearch_; 

  ROL::Ptr<Vector<Real> > d_;

  ELineSearch         els_;   
  ECurvatureCondition econd_; 

  bool acceptLastAlpha_;  

  bool usePreviousAlpha_; 

  int verbosity_;
  bool computeObj_;
  Real fval_;

  ROL::ParameterList parlist_;

  std::string lineSearchName_;  

  Real GradDotStep(const Vector<Real> &g, const Vector<Real> &s,
                   const Vector<Real> &x,
                   BoundConstraint<Real> &bnd, Real eps = 0) {
    Real gs(0), one(1);
    if (!bnd.isActivated()) {
      gs = s.dot(g.dual());
    }
    else {
      d_->set(s);
      bnd.pruneActive(*d_,g,x,eps);
      gs = d_->dot(g.dual());
      d_->set(x);
      d_->axpy(-one,g.dual());
      bnd.project(*d_);
      d_->scale(-one);
      d_->plus(x);
      bnd.pruneInactive(*d_,g,x,eps);
      gs -= d_->dot(g.dual());
    }
    return gs;
  }

public:

  using Step<Real>::initialize;
  using Step<Real>::compute;
  using Step<Real>::update;

  LineSearchStep( ROL::ParameterList &parlist,
                  const ROL::Ptr<LineSearch<Real> > &lineSearch = ROL::nullPtr,
                  const ROL::Ptr<Secant<Real> > &secant = ROL::nullPtr,
                  const ROL::Ptr<Krylov<Real> > &krylov = ROL::nullPtr,
                  const ROL::Ptr<NonlinearCG<Real> > &nlcg = ROL::nullPtr )
    : Step<Real>(), desc_(ROL::nullPtr), secant_(secant),
      krylov_(krylov), nlcg_(nlcg), lineSearch_(lineSearch),
      els_(LINESEARCH_USERDEFINED), econd_(CURVATURECONDITION_WOLFE),
      verbosity_(0), computeObj_(true), fval_(0), parlist_(parlist) {
    // Parse parameter list
    ROL::ParameterList& Llist = parlist.sublist("Step").sublist("Line Search");
    ROL::ParameterList& Glist = parlist.sublist("General");
    econd_ = StringToECurvatureCondition(Llist.sublist("Curvature Condition").get("Type","Strong Wolfe Conditions") );
    acceptLastAlpha_ = Llist.get("Accept Last Alpha", false); 
    verbosity_ = Glist.get("Print Verbosity",0);
    computeObj_ = Glist.get("Recompute Objective Function",false);
    // Initialize Line Search
    if (lineSearch_ == ROL::nullPtr) {
      lineSearchName_ = Llist.sublist("Line-Search Method").get("Type","Cubic Interpolation"); 
      els_ = StringToELineSearch(lineSearchName_);
      lineSearch_ = LineSearchFactory<Real>(parlist);
    } 
    else { // User-defined linesearch provided
      lineSearchName_ = Llist.sublist("Line-Search Method").get("User Defined Line-Search Name",
                                                                "Unspecified User Defined Line-Search");
    }

  }

  void initialize( Vector<Real> &x, const Vector<Real> &s, const Vector<Real> &g, 
                   Objective<Real> &obj, BoundConstraint<Real> &bnd, 
                   AlgorithmState<Real> &algo_state ) {
    d_ = x.clone();

    // Initialize unglobalized step
    ROL::ParameterList& list
      = parlist_.sublist("Step").sublist("Line Search").sublist("Descent Method");
    EDescent edesc = StringToEDescent(list.get("Type","Quasi-Newton Method"));
    if (bnd.isActivated()) {
      switch(edesc) {
        case DESCENT_STEEPEST: {
          desc_ = ROL::makePtr<GradientStep<Real>>(parlist_,computeObj_);
          break;
        }
        case DESCENT_NONLINEARCG: {
          desc_ = ROL::makePtr<NonlinearCGStep<Real>>(parlist_,nlcg_,computeObj_);
          break;
        }
        case DESCENT_SECANT: {
          desc_ = ROL::makePtr<ProjectedSecantStep<Real>>(parlist_,secant_,computeObj_);
          break;
        }
        case DESCENT_NEWTON: {
          desc_ = ROL::makePtr<ProjectedNewtonStep<Real>>(parlist_,computeObj_);
          break;
        }
        case DESCENT_NEWTONKRYLOV: {
          desc_ = ROL::makePtr<ProjectedNewtonKrylovStep<Real>>(parlist_,krylov_,secant_,computeObj_);
          break;
        }
        default:
          ROL_TEST_FOR_EXCEPTION(true,std::invalid_argument,
            ">>> (LineSearchStep::Initialize): Undefined descent type!");
      }
    }
    else {
      switch(edesc) {
        case DESCENT_STEEPEST: {
          desc_ = ROL::makePtr<GradientStep<Real>>(parlist_,computeObj_);
          break;
        }
        case DESCENT_NONLINEARCG: {
          desc_ = ROL::makePtr<NonlinearCGStep<Real>>(parlist_,nlcg_,computeObj_);
          break;
        }
        case DESCENT_SECANT: {
          desc_ = ROL::makePtr<SecantStep<Real>>(parlist_,secant_,computeObj_);
          break;
        }
        case DESCENT_NEWTON: {
          desc_ = ROL::makePtr<NewtonStep<Real>>(parlist_,computeObj_);
          break;
        }
        case DESCENT_NEWTONKRYLOV: {
          desc_ = ROL::makePtr<NewtonKrylovStep<Real>>(parlist_,krylov_,secant_,computeObj_);
          break;
        }
        default:
          ROL_TEST_FOR_EXCEPTION(true,std::invalid_argument,
            ">>> (LineSearchStep::Initialize): Undefined descent type!");
      }
    }
    desc_->initialize(x,s,g,obj,bnd,algo_state);

    // Initialize line search
    lineSearch_->initialize(x,s,g,obj,bnd);
    //const ROL::Ptr<const StepState<Real> > desc_state = desc_->getStepState();
    //lineSearch_->initialize(x,s,*(desc_state->gradientVec),obj,bnd);
  }

  void compute( Vector<Real> &s, const Vector<Real> &x,
                Objective<Real> &obj, BoundConstraint<Real> &bnd, 
                AlgorithmState<Real> &algo_state ) {
    Real zero(0), one(1);
    // Compute unglobalized step
    desc_->compute(s,x,obj,bnd,algo_state);

    // Ensure that s is a descent direction
    // ---> If not, then default to steepest descent
    const ROL::Ptr<const StepState<Real> > desc_state = desc_->getStepState();
    Real gs = GradDotStep(*(desc_state->gradientVec),s,x,bnd,algo_state.gnorm);
    if (gs >= zero) {
      s.set((desc_state->gradientVec)->dual());
      s.scale(-one);
      gs = GradDotStep(*(desc_state->gradientVec),s,x,bnd,algo_state.gnorm);
    }

    // Perform line search
    ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState();
    fval_ = algo_state.value;
    step_state->nfval = 0; step_state->ngrad = 0;
    lineSearch_->setData(algo_state.gnorm,*(desc_state->gradientVec));
    lineSearch_->run(step_state->searchSize,fval_,step_state->nfval,step_state->ngrad,gs,s,x,obj,bnd);

    // Make correction if maximum function evaluations reached
    if(!acceptLastAlpha_) {
      lineSearch_->setMaxitUpdate(step_state->searchSize,fval_,algo_state.value);
    }

    // Compute scaled descent direction
    s.scale(step_state->searchSize);
    if ( bnd.isActivated() ) {
      s.plus(x);
      bnd.project(s);
      s.axpy(static_cast<Real>(-1),x);
    }
  }

  void update( Vector<Real> &x, const Vector<Real> &s,
               Objective<Real> &obj, BoundConstraint<Real> &bnd,
               AlgorithmState<Real> &algo_state ) {
    ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState();
    algo_state.nfval += step_state->nfval;
    algo_state.ngrad += step_state->ngrad;
    desc_->update(x,s,obj,bnd,algo_state);
    step_state->flag = desc_->getStepState()->flag;
    step_state->SPiter = desc_->getStepState()->SPiter;
    step_state->SPflag = desc_->getStepState()->SPflag;
    if ( !computeObj_ ) {
      algo_state.value = fval_;
    }
  }

  std::string printHeader( void ) const {
    std::string head = desc_->printHeader();
    head.erase(std::remove(head.end()-3,head.end(),'\n'), head.end());
    std::stringstream hist;
    hist.write(head.c_str(),head.length());
    hist << std::setw(10) << std::left << "ls_#fval";
    hist << std::setw(10) << std::left << "ls_#grad";
    hist << "\n";
    return hist.str();
  }
  
  std::string printName( void ) const {
    std::string name = desc_->printName();
    std::stringstream hist;
    hist << name;
    hist << "Line Search: " << lineSearchName_;
    hist << " satisfying " << ECurvatureConditionToString(econd_) << "\n";
    return hist.str();
  }

  std::string print( AlgorithmState<Real> & algo_state, bool print_header = false ) const  {
    const ROL::Ptr<const StepState<Real> > step_state = Step<Real>::getStepState();
    std::string desc = desc_->print(algo_state,false);
    desc.erase(std::remove(desc.end()-3,desc.end(),'\n'), desc.end());
    std::string name = desc_->printName();
    size_t pos = desc.find(name);
    if ( pos != std::string::npos ) {
      desc.erase(pos, name.length());
    }

    std::stringstream hist;
    if ( algo_state.iter == 0 ) {
      hist << printName();
    }
    if ( print_header ) {
      hist << printHeader();
    }
    hist << desc;
    if ( algo_state.iter == 0 ) {
      hist << "\n";
    }
    else {
      hist << std::setw(10) << std::left << step_state->nfval;              
      hist << std::setw(10) << std::left << step_state->ngrad;
      hist << "\n";
    }
    return hist.str();
  }
}; // class LineSearchStep

} // namespace ROL

#endif