ROL_TrustRegionStep.hpp

Go to the documentation of this file.

// @HEADER
// ************************************************************************
//
//               Rapid Optimization Library (ROL) Package
//                 Copyright (2014) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact lead developers:
//              Drew Kouri   (dpkouri@sandia.gov) and
//              Denis Ridzal (dridzal@sandia.gov)
//
// ************************************************************************
// @HEADER

#ifndef ROL_TRUSTREGIONSTEP_H
#define ROL_TRUSTREGIONSTEP_H

#include "ROL_Step.hpp"
#include "ROL_Types.hpp"
#include "ROL_Secant.hpp"
#include "ROL_TrustRegion.hpp"
#include <sstream>
#include <iomanip>

namespace ROL {

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

  Teuchos::RCP<Secant<Real> >       secant_;      
  Teuchos::RCP<TrustRegion<Real> >  trustRegion_; 

  Teuchos::RCP<Vector<Real> >       xnew_;        
  Teuchos::RCP<Vector<Real> >       xold_;        
  Teuchos::RCP<Vector<Real> >       gp_;          

  ETrustRegion      etr_;        
  ESecant           esec_;       

  bool useSecantHessVec_;        
  bool useSecantPrecond_;        

  bool useProjectedGrad_;        

  std::vector<bool> useInexact_; 
  int               TRflag_  ;   
  int               TR_nfval_;   
  int               TR_ngrad_;   
  int               CGflag_;     
  int               CGiter_;     

  Real              delMax_;     

  Real              alpha_init_; 
  int               max_fval_;   

  Real              scale0_; 
  Real              scale1_; 

  bool              softUp_;
  Real              scaleEps_;

  void updateGradient( Vector<Real> &x, Objective<Real> &obj, BoundConstraint<Real> &con, 
                       AlgorithmState<Real> &algo_state ) {
    Teuchos::RCP<StepState<Real> > state = Step<Real>::getState();
    if ( useInexact_[1] ) {
      Real c = scale0_*std::max(1.e-2,std::min(1.0,1.e4*algo_state.gnorm));
      Real gtol1  = c*(state->searchSize);
      Real gtol0  = scale1_*gtol1 + 1.0;
      while ( gtol0 > gtol1*scale1_ ) {
        obj.gradient(*(state->gradientVec),x,gtol1);
        algo_state.gnorm = computeCriticalityMeasure(*(state->gradientVec),x,con);
        gtol0 = gtol1;
        c = scale0_*std::max(1.e-2,std::min(1.0,1.e4*algo_state.gnorm));
        gtol1 = c*std::min(algo_state.gnorm,state->searchSize);
      }
      algo_state.ngrad++;
    }
    else {
      Real gtol = std::sqrt(ROL_EPSILON);
      obj.gradient(*(state->gradientVec),x,gtol);
      algo_state.ngrad++;
      algo_state.gnorm = computeCriticalityMeasure(*(state->gradientVec),x,con);
    }
  }

  Real computeCriticalityMeasure( const Vector<Real> &g, const Vector<Real> &x, BoundConstraint<Real> &con ) {
    if ( con.isActivated() ) {
      if ( useProjectedGrad_ ) {
        gp_->set(g);
        con.computeProjectedGradient( *gp_, x );
        return gp_->norm();
      }
      else {
        xnew_->set(x);
        xnew_->axpy(-1.0,g.dual());
        con.project(*xnew_);
        xnew_->axpy(-1.0,x);
        return xnew_->norm();
      }
    }
    else {
      return g.norm();
    }
  }

public:

  virtual ~TrustRegionStep() {}

  TrustRegionStep( Teuchos::ParameterList & parlist )
    : Step<Real>(),
      secant_(Teuchos::null), trustRegion_(Teuchos::null),
      xnew_(Teuchos::null), xold_(Teuchos::null), gp_(Teuchos::null),
      etr_(TRUSTREGION_DOGLEG), esec_(SECANT_LBFGS),
      useSecantHessVec_(false), useSecantPrecond_(false),
      useProjectedGrad_(false),
      TRflag_(0), TR_nfval_(0), TR_ngrad_(0),
      CGflag_(0), CGiter_(0),
      delMax_(1.e4),
      alpha_init_(1.), max_fval_(20),
      scale0_(1.), scale1_(1.),
      softUp_(false), scaleEps_(1.) {
    Teuchos::RCP<StepState<Real> > step_state = Step<Real>::getState();
    // Trust-Region Parameters
    step_state->searchSize = parlist.sublist("Step").sublist("Trust Region").get("Initial Radius", -1.0);
    delMax_                = parlist.sublist("Step").sublist("Trust Region").get("Maximum Radius", 1000.0);
    // Inexactness Information
    useInexact_.clear();
    useInexact_.push_back(parlist.sublist("General").get("Inexact Objective Function", false));
    useInexact_.push_back(parlist.sublist("General").get("Inexact Gradient", false));
    useInexact_.push_back(parlist.sublist("General").get("Inexact Hessian-Times-A-Vector", false));
    // Trust-Region Inexactness Parameters
    scale0_ = parlist.sublist("Step").sublist("Trust Region").sublist("Inexact").sublist("Gradient").get("Tolerance Scaling",1.e-1);
    scale1_ = parlist.sublist("Step").sublist("Trust Region").sublist("Inexact").sublist("Gradient").get("Relative Tolerance",2.0); 
    // Initialize Trust Region Subproblem Solver Object
    etr_  = StringToETrustRegion(parlist.sublist("Step").sublist("Trust Region").get("Subproblem Solver","Dogleg"));  
    useProjectedGrad_ = parlist.sublist("General").get("Projected Gradient Criticality Measure", false);
    max_fval_         = parlist.sublist("Step").sublist("Line Search").get("Function Evaluation Limit", 20);
    alpha_init_       = parlist.sublist("Step").sublist("Line Search").get("Initial Step Size", 1.0);
    trustRegion_ = TrustRegionFactory<Real>(parlist);
    // Secant Object
    esec_ = StringToESecant(parlist.sublist("General").sublist("Secant").get("Type","Limited-Memory BFGS"));
    useSecantPrecond_ = parlist.sublist("General").sublist("Secant").get("Use as Preconditioner", false);
    useSecantHessVec_ = parlist.sublist("General").sublist("Secant").get("Use as Hessian", false);
    secant_ = SecantFactory<Real>(parlist);
    // Changing Objective Functions
    softUp_ = parlist.sublist("General").get("Variable Objective Function",false);
    // Scale for epsilon active sets
    scaleEps_ = parlist.sublist("General").get("Scale for Epsilon Active Sets",1.0);
  }

  TrustRegionStep( Teuchos::RCP<Secant<Real> > &secant, Teuchos::ParameterList &parlist ) 
    : Step<Real>(),
      secant_(secant), trustRegion_(Teuchos::null),
      xnew_(Teuchos::null), xold_(Teuchos::null), gp_(Teuchos::null),
      etr_(TRUSTREGION_DOGLEG), esec_(SECANT_USERDEFINED),
      useSecantHessVec_(false), useSecantPrecond_(false),
      useProjectedGrad_(false),
      TRflag_(0), TR_nfval_(0), TR_ngrad_(0),
      CGflag_(0), CGiter_(0),
      delMax_(1.e4),
      alpha_init_(1.), max_fval_(20),
      scale0_(1.), scale1_(1.),
      softUp_(false), scaleEps_(1.) {
    Teuchos::RCP<StepState<Real> > step_state = Step<Real>::getState();
    // Trust-Region Parameters
    step_state->searchSize = parlist.sublist("Step").sublist("Trust Region").get("Initial Radius", -1.0);
    delMax_                = parlist.sublist("Step").sublist("Trust Region").get("Maximum Radius", 1000.0);
    // Inexactness Information
    useInexact_.clear();
    useInexact_.push_back(parlist.sublist("General").get("Inexact Objective Function", false));
    useInexact_.push_back(parlist.sublist("General").get("Inexact Gradient", false));
    useInexact_.push_back(parlist.sublist("General").get("Inexact Hessian-Times-A-Vector", false));
    // Trust-Region Inexactness Parameters
    scale0_ = parlist.sublist("Step").sublist("Trust Region").sublist("Inexact").sublist("Gradient").get("Tolerance Scaling",1.e-1);
    scale1_ = parlist.sublist("Step").sublist("Trust Region").sublist("Inexact").sublist("Gradient").get("Relative Tolerance",2.0); 
    // Initialize Trust Region Subproblem Solver Object
    etr_  = StringToETrustRegion(parlist.sublist("Step").sublist("Trust Region").get("Subproblem Solver","Dogleg"));  
    useProjectedGrad_ = parlist.sublist("General").get("Projected Gradient Criticality Measure", false);
    max_fval_         = parlist.sublist("Step").sublist("Line Search").get("Function Evaluation Limit", 20);
    alpha_init_       = parlist.sublist("Step").sublist("Line Search").get("Initial Step Size", 1.0);
    trustRegion_ = TrustRegionFactory<Real>(parlist);
    // Secant Object
    useSecantPrecond_ = parlist.sublist("General").sublist("Secant").get("Use as Preconditioner", false);
    useSecantHessVec_ = parlist.sublist("General").sublist("Secant").get("Use as Hessian", false);
    // Changing Objective Functions
    softUp_ = parlist.sublist("General").get("Variable Objective Function",false);
    // Scale for epsilon active sets
    scaleEps_ = parlist.sublist("General").get("Scale for Epsilon Active Sets",1.0);
  }

  void initialize( Vector<Real> &x, const Vector<Real> &s, const Vector<Real> &g, 
                   Objective<Real> &obj, BoundConstraint<Real> &con, 
                   AlgorithmState<Real> &algo_state ) {
    Teuchos::RCP<StepState<Real> > step_state = Step<Real>::getState();

    trustRegion_->initialize(x,s,g);

    Real htol = std::sqrt(ROL_EPSILON);
    Real ftol = 0.1*ROL_OVERFLOW; 

    step_state->descentVec  = s.clone();
    step_state->gradientVec = g.clone();

    if ( con.isActivated() ) {
      con.project(x);
      xnew_ = x.clone();
      xold_ = x.clone();
    }
 
    if ( con.isActivated() || secant_ != Teuchos::null ) {
      gp_ = g.clone();
    }

    // Update approximate gradient and approximate objective function.
    obj.update(x,true,algo_state.iter);    
    updateGradient(x,obj,con,algo_state);
    algo_state.snorm = 1.e10;
    algo_state.value = obj.value(x,ftol); 
    algo_state.nfval++;

    // Evaluate Objective Function at Cauchy Point
    if ( step_state->searchSize <= 0.0 ) {
      Teuchos::RCP<Vector<Real> > Bg = g.clone();
      if ( useSecantHessVec_ ) {
        secant_->applyB(*Bg,(step_state->gradientVec)->dual(),x);
      }
      else {
        obj.hessVec(*Bg,(step_state->gradientVec)->dual(),x,htol);
      }
      Real gBg = Bg->dot(*(step_state->gradientVec));
      Real alpha = 1.0;
      if ( gBg > ROL_EPSILON ) {
        alpha = algo_state.gnorm*algo_state.gnorm/gBg;
      }
      // Evaluate the objective function at the Cauchy point
      Teuchos::RCP<Vector<Real> > cp = s.clone();
      cp->set((step_state->gradientVec)->dual()); 
      cp->scale(-alpha);
      Teuchos::RCP<Vector<Real> > xcp = x.clone();
      xcp->set(x);
      xcp->plus(*cp);
      if ( con.isActivated() ) {
        con.project(*xcp);
      }
      obj.update(*xcp);
      Real fnew = obj.value(*xcp,ftol); // MUST DO SOMETHING HERE WITH FTOL
      algo_state.nfval++;
      // Perform cubic interpolation to determine initial trust region radius
      Real gs = cp->dot((step_state->gradientVec)->dual());
      Real a  = fnew - algo_state.value - gs - 0.5*alpha*alpha*gBg;
      if ( std::abs(a) < ROL_EPSILON ) { 
        // a = 0 implies the objective is quadratic in the negative gradient direction
        step_state->searchSize = std::min(alpha*algo_state.gnorm,delMax_);
      }
      else {
        Real b  = 0.5*alpha*alpha*gBg;
        Real c  = gs;
        if ( b*b-3.0*a*c > ROL_EPSILON ) {
          // There is at least one critical point
          Real t1 = (-b-std::sqrt(b*b-3.0*a*c))/(3.0*a);
          Real t2 = (-b+std::sqrt(b*b-3.0*a*c))/(3.0*a);
          if ( 6.0*a*t1 + 2.0*b > 0.0 ) {
            // t1 is the minimizer
            step_state->searchSize = std::min(t1*alpha*algo_state.gnorm,delMax_);
          }
          else {
            // t2 is the minimizer
            step_state->searchSize = std::min(t2*alpha*algo_state.gnorm,delMax_);
          }
        }
        else {
          step_state->searchSize = std::min(alpha*algo_state.gnorm,delMax_);
        }
      }
    }
  }

  void compute( Vector<Real> &s, const Vector<Real> &x, Objective<Real> &obj, BoundConstraint<Real> &con, 
                AlgorithmState<Real> &algo_state ) {
    Teuchos::RCP<StepState<Real> > step_state = Step<Real>::getState();

    Real eps = 0.0;
    if ( con.isActivated() ) {
      eps = scaleEps_*algo_state.gnorm;
    }
    ProjectedObjective<Real> pObj(obj,con,secant_,useSecantPrecond_,useSecantHessVec_,eps);

    CGflag_ = 0;
    CGiter_ = 0;
    trustRegion_->run(s,algo_state.snorm,step_state->searchSize,CGflag_,CGiter_,
                            x,*(step_state->gradientVec),algo_state.gnorm,pObj);

    if ( con.isActivated() ) {
      xnew_->set(x);
      xnew_->plus(s);
      con.project(*xnew_);
      s.set(*xnew_);
      s.axpy(-1.0,x);
    }
  }

  void update( Vector<Real> &x, const Vector<Real> &s, Objective<Real> &obj, BoundConstraint<Real> &con, 
               AlgorithmState<Real> &algo_state ) {
    Teuchos::RCP<StepState<Real> > state = Step<Real>::getState();

    Real tol = std::sqrt(ROL_EPSILON);

    Real eps = 0.0;
    if ( con.isActivated() ) {
      eps = algo_state.gnorm;
    }
    ProjectedObjective<Real> pObj(obj,con,secant_,useSecantPrecond_,useSecantHessVec_,eps);

    // Store previous step for constraint computations
    if ( con.isActivated() ) {
      xold_->set(x);
    }

    // Update trust-region information;
    // Performs a hard update on the objective function
    TRflag_   = 0;
    TR_nfval_ = 0;
    TR_ngrad_ = 0;
    Real fold = algo_state.value;
    Real fnew = 0.0;
    algo_state.iter++;
    trustRegion_->update(x,fnew,state->searchSize,TR_nfval_,TR_ngrad_,TRflag_,
                         s,algo_state.snorm,fold,*(state->gradientVec),algo_state.iter,pObj);
    algo_state.nfval += TR_nfval_;
    algo_state.ngrad += TR_ngrad_;
    algo_state.value  = fnew;

    // If step is accepted ...
    // Compute new gradient and update secant storage
    if ( TRflag_ == 0 || TRflag_ == 1 ) {  
      // Perform line search (smoothing) to ensure decrease 
      if ( con.isActivated() ) {
        // Compute new gradient
        obj.gradient(*gp_,x,tol); // MUST DO SOMETHING HERE WITH TOL
        algo_state.ngrad++;
        // Compute smoothed step
        Real alpha = 1.0;
        xnew_->set(x);
        xnew_->axpy(-alpha*alpha_init_,gp_->dual());
        con.project(*xnew_);
        // Compute new objective value
        if ( softUp_ ) {
          obj.update(*xnew_);
        }
        else {
          obj.update(*xnew_,true,algo_state.iter);
        }
        Real ftmp = obj.value(*xnew_,tol); // MUST DO SOMETHING HERE WITH TOL
        algo_state.nfval++;
        // Perform smoothing
        int cnt = 0;
        alpha = 1.0/alpha_init_;
        while ( (fnew-ftmp) <= 1.e-4*(fnew-fold) ) { 
          xnew_->set(x);
          xnew_->axpy(-alpha*alpha_init_,gp_->dual());
          con.project(*xnew_);
          if ( softUp_ ) {
            obj.update(*xnew_,false,algo_state.iter);
          }
          else {
            obj.update(*xnew_,true,algo_state.iter);
          }
          ftmp = obj.value(*xnew_,tol); // MUST DO SOMETHING HERE WITH TOL
          algo_state.nfval++;
          if ( cnt >= max_fval_ ) {
            break;
          }
          alpha *= 0.5;
          cnt++;
        }
        // Store objective function and iteration information
        fnew = ftmp;
        x.set(*xnew_);
      }
      else {
        if (softUp_) {
          pObj.update(x,true,algo_state.iter);
        }
      }

      // Store previous gradient for secant update
      if ( secant_ != Teuchos::null ) {
        gp_->set(*(state->gradientVec));
      }

      // Update objective function and approximate model
      updateGradient(x,obj,con,algo_state);

      // Update secant information
      if ( secant_ != Teuchos::null ) {
        if ( con.isActivated() ) { // Compute new constrained step
          xnew_->set(x);
          xnew_->axpy(-1.0,*xold_);
          secant_->update(*(state->gradientVec),*gp_,*xnew_,algo_state.snorm,algo_state.iter+1);
        }
        else {
          secant_->update(*(state->gradientVec),*gp_,s,algo_state.snorm,algo_state.iter+1);
        }
      }

      // Update algorithm state
      (algo_state.iterateVec)->set(x);
    }
    else { // Step was rejected
      if ( softUp_ ) {
        obj.update(x,true,algo_state.iter);
        fnew = pObj.value(x,tol);
        algo_state.nfval++;
        algo_state.value = fnew; 
      }
    }

  }

  std::string printHeader( void ) const  {
    std::stringstream hist;
    hist << "  ";
    hist << std::setw(6)  << std::left << "iter";
    hist << std::setw(15) << std::left << "value";
    hist << std::setw(15) << std::left << "gnorm";
    hist << std::setw(15) << std::left << "snorm";
    hist << std::setw(15) << std::left << "delta";
    hist << std::setw(10) << std::left << "#fval";
    hist << std::setw(10) << std::left << "#grad";
    hist << std::setw(10) << std::left << "tr_flag";
    if ( etr_ == TRUSTREGION_TRUNCATEDCG ) {
      hist << std::setw(10) << std::left << "iterCG";
      hist << std::setw(10) << std::left << "flagCG";
    }
    hist << "\n";
    return hist.str();
  }

  std::string printName( void ) const {
    std::stringstream hist;
    hist << "\n" << ETrustRegionToString(etr_) << " Trust-Region solver";
    if ( useSecantPrecond_ || useSecantHessVec_ ) {
      if ( useSecantPrecond_ && !useSecantHessVec_ ) {
        hist << " with " << ESecantToString(esec_) << " preconditioning\n";
      }
      else if ( !useSecantPrecond_ && useSecantHessVec_ ) {
        hist << " with " << ESecantToString(esec_) << " Hessian approximation\n";
      }
      else {
        hist << " with " << ESecantToString(esec_) << " preconditioning and Hessian approximation\n";
      }
    }
    else {
      hist << "\n";
    }
    return hist.str();
  }

  std::string print( AlgorithmState<Real> & algo_state, bool print_header = false ) const  {
    const Teuchos::RCP<const StepState<Real> >& step_state = Step<Real>::getStepState();

    std::stringstream hist;
    hist << std::scientific << std::setprecision(6);
    if ( algo_state.iter == 0 ) {
      hist << printName();
    }
    if ( print_header ) {
      hist << printHeader();
    }
    if ( algo_state.iter == 0 ) {
      hist << "  ";
      hist << std::setw(6)  << std::left << algo_state.iter;
      hist << std::setw(15) << std::left << algo_state.value;
      hist << std::setw(15) << std::left << algo_state.gnorm;
      hist << std::setw(15) << std::left << " "; 
      hist << std::setw(15) << std::left << step_state->searchSize; 
      hist << "\n";
    }
    else {
      hist << "  "; 
      hist << std::setw(6)  << std::left << algo_state.iter;  
      hist << std::setw(15) << std::left << algo_state.value; 
      hist << std::setw(15) << std::left << algo_state.gnorm; 
      hist << std::setw(15) << std::left << algo_state.snorm; 
      hist << std::setw(15) << std::left << step_state->searchSize; 
      hist << std::setw(10) << std::left << algo_state.nfval;              
      hist << std::setw(10) << std::left << algo_state.ngrad;              
      hist << std::setw(10) << std::left << TRflag_;              
      if ( etr_ == TRUSTREGION_TRUNCATEDCG ) {
        hist << std::setw(10) << std::left << CGiter_;
        hist << std::setw(10) << std::left << CGflag_;
      }
      hist << "\n";
    }
    return hist.str();
  }

}; // class Step

} // namespace ROL

#endif