ROL
ROL_NewtonKrylovStep.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 
10 #ifndef ROL_NEWTONKRYLOVSTEP_H
11 #define ROL_NEWTONKRYLOVSTEP_H
12 
13 #include "ROL_Types.hpp"
14 #include "ROL_Step.hpp"
15 
16 #include "ROL_Secant.hpp"
17 #include "ROL_KrylovFactory.hpp"
18 #include "ROL_LinearOperator.hpp"
19 
20 #include <sstream>
21 #include <iomanip>
22 
29 namespace ROL {
30 
31 template <class Real>
32 class NewtonKrylovStep : public Step<Real> {
33 private:
34 
35  ROL::Ptr<Secant<Real> > secant_;
36  ROL::Ptr<Krylov<Real> > krylov_;
37 
40 
41  ROL::Ptr<Vector<Real> > gp_;
42 
45  int verbosity_;
46  const bool computeObj_;
47 
49 
50  std::string krylovName_;
51  std::string secantName_;
52 
53 
54  class HessianNK : public LinearOperator<Real> {
55  private:
56  const ROL::Ptr<Objective<Real> > obj_;
57  const ROL::Ptr<Vector<Real> > x_;
58  public:
59  HessianNK(const ROL::Ptr<Objective<Real> > &obj,
60  const ROL::Ptr<Vector<Real> > &x) : obj_(obj), x_(x) {}
61  void apply(Vector<Real> &Hv, const Vector<Real> &v, Real &tol) const {
62  obj_->hessVec(Hv,v,*x_,tol);
63  }
64  };
65 
66  class PrecondNK : public LinearOperator<Real> {
67  private:
68  const ROL::Ptr<Objective<Real> > obj_;
69  const ROL::Ptr<Vector<Real> > x_;
70  public:
71  PrecondNK(const ROL::Ptr<Objective<Real> > &obj,
72  const ROL::Ptr<Vector<Real> > &x) : obj_(obj), x_(x) {}
73  void apply(Vector<Real> &Hv, const Vector<Real> &v, Real &tol) const {
74  Hv.set(v.dual());
75  }
76  void applyInverse(Vector<Real> &Hv, const Vector<Real> &v, Real &tol) const {
77  obj_->precond(Hv,v,*x_,tol);
78  }
79  };
80 
81 public:
82 
84  using Step<Real>::compute;
85  using Step<Real>::update;
86 
94  NewtonKrylovStep( ROL::ParameterList &parlist, const bool computeObj = true )
95  : Step<Real>(), secant_(ROL::nullPtr), krylov_(ROL::nullPtr),
96  gp_(ROL::nullPtr), iterKrylov_(0), flagKrylov_(0),
97  verbosity_(0), computeObj_(computeObj), useSecantPrecond_(false) {
98  // Parse ParameterList
99  ROL::ParameterList& Glist = parlist.sublist("General");
100  useSecantPrecond_ = Glist.sublist("Secant").get("Use as Preconditioner", false);
101  verbosity_ = Glist.get("Print Verbosity",0);
102  // Initialize Krylov object
103  krylovName_ = Glist.sublist("Krylov").get("Type","Conjugate Gradients");
105  krylov_ = KrylovFactory<Real>(parlist);
106  // Initialize secant object
107  secantName_ = Glist.sublist("Secant").get("Type","Limited-Memory BFGS");
108  esec_ = StringToESecant(secantName_);
109  if ( useSecantPrecond_ ) {
110  secant_ = SecantFactory<Real>(parlist);
111  }
112  }
113 
124  NewtonKrylovStep(ROL::ParameterList &parlist,
125  const ROL::Ptr<Krylov<Real> > &krylov,
126  const ROL::Ptr<Secant<Real> > &secant,
127  const bool computeObj = true)
128  : Step<Real>(), secant_(secant), krylov_(krylov),
130  gp_(ROL::nullPtr), iterKrylov_(0), flagKrylov_(0),
131  verbosity_(0), computeObj_(computeObj), useSecantPrecond_(false) {
132  // Parse ParameterList
133  ROL::ParameterList& Glist = parlist.sublist("General");
134  useSecantPrecond_ = Glist.sublist("Secant").get("Use as Preconditioner", false);
135  verbosity_ = Glist.get("Print Verbosity",0);
136  // Initialize secant object
137  if ( useSecantPrecond_ ) {
138  if(secant_ == ROL::nullPtr ) {
139  secantName_ = Glist.sublist("Secant").get("Type","Limited-Memory BFGS");
141  secant_ = SecantFactory<Real>(parlist);
142  }
143  else {
144  secantName_ = Glist.sublist("Secant").get("User Defined Secant Name",
145  "Unspecified User Defined Secant Method");
146  }
147  }
148  // Initialize Krylov object
149  if ( krylov_ == ROL::nullPtr ) {
150  krylovName_ = Glist.sublist("Krylov").get("Type","Conjugate Gradients");
152  krylov_ = KrylovFactory<Real>(parlist);
153  }
154  else {
155  krylovName_ = Glist.sublist("Krylov").get("User Defined Krylov Name",
156  "Unspecified User Defined Krylov Method");
157  }
158  }
159 
160  void initialize( Vector<Real> &x, const Vector<Real> &s, const Vector<Real> &g,
162  AlgorithmState<Real> &algo_state ) {
163  Step<Real>::initialize(x,s,g,obj,bnd,algo_state);
164  if ( useSecantPrecond_ ) {
165  gp_ = g.clone();
166  }
167  }
168 
169  void compute( Vector<Real> &s, const Vector<Real> &x,
171  AlgorithmState<Real> &algo_state ) {
172  Real one(1);
173  ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState();
174 
175  // Build Hessian and Preconditioner object
176  ROL::Ptr<Objective<Real> > obj_ptr = ROL::makePtrFromRef(obj);
177  ROL::Ptr<LinearOperator<Real> > hessian
178  = ROL::makePtr<HessianNK>(obj_ptr,algo_state.iterateVec);
179  ROL::Ptr<LinearOperator<Real> > precond;
180  if ( useSecantPrecond_ ) {
181  precond = secant_;
182  }
183  else {
184  precond = ROL::makePtr<PrecondNK>(obj_ptr,algo_state.iterateVec);
185  }
186 
187  // Run Krylov method
188  flagKrylov_ = 0;
189  krylov_->run(s,*hessian,*(step_state->gradientVec),*precond,iterKrylov_,flagKrylov_);
190 
191  // Check Krylov flags
192  if ( flagKrylov_ == 2 && iterKrylov_ <= 1 ) {
193  s.set((step_state->gradientVec)->dual());
194  }
195  s.scale(-one);
196  }
197 
198  void update( Vector<Real> &x, const Vector<Real> &s,
200  AlgorithmState<Real> &algo_state ) {
201  Real tol = std::sqrt(ROL_EPSILON<Real>());
202  ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState();
203  step_state->SPiter = iterKrylov_;
204  step_state->SPflag = flagKrylov_;
205 
206  // Update iterate
207  algo_state.iter++;
208  x.plus(s);
209  (step_state->descentVec)->set(s);
210  algo_state.snorm = s.norm();
211 
212  // Compute new gradient
213  if ( useSecantPrecond_ ) {
214  gp_->set(*(step_state->gradientVec));
215  }
216  obj.update(x,true,algo_state.iter);
217  if ( computeObj_ ) {
218  algo_state.value = obj.value(x,tol);
219  algo_state.nfval++;
220  }
221  obj.gradient(*(step_state->gradientVec),x,tol);
222  algo_state.ngrad++;
223 
224  // Update Secant Information
225  if ( useSecantPrecond_ ) {
226  secant_->updateStorage(x,*(step_state->gradientVec),*gp_,s,algo_state.snorm,algo_state.iter+1);
227  }
228 
229  // Update algorithm state
230  (algo_state.iterateVec)->set(x);
231  algo_state.gnorm = step_state->gradientVec->norm();
232  }
233 
234  std::string printHeader( void ) const {
235  std::stringstream hist;
236 
237  if( verbosity_>0 ) {
238  hist << std::string(109,'-') << "\n";
240  hist << " status output definitions\n\n";
241  hist << " iter - Number of iterates (steps taken) \n";
242  hist << " value - Objective function value \n";
243  hist << " gnorm - Norm of the gradient\n";
244  hist << " snorm - Norm of the step (update to optimization vector)\n";
245  hist << " #fval - Cumulative number of times the objective function was evaluated\n";
246  hist << " #grad - Number of times the gradient was computed\n";
247  hist << " iterCG - Number of Krylov iterations used to compute search direction\n";
248  hist << " flagCG - Krylov solver flag" << "\n";
249  hist << std::string(109,'-') << "\n";
250  }
251 
252  hist << " ";
253  hist << std::setw(6) << std::left << "iter";
254  hist << std::setw(15) << std::left << "value";
255  hist << std::setw(15) << std::left << "gnorm";
256  hist << std::setw(15) << std::left << "snorm";
257  hist << std::setw(10) << std::left << "#fval";
258  hist << std::setw(10) << std::left << "#grad";
259  hist << std::setw(10) << std::left << "iterCG";
260  hist << std::setw(10) << std::left << "flagCG";
261  hist << "\n";
262  return hist.str();
263  }
264  std::string printName( void ) const {
265  std::stringstream hist;
266  hist << "\n" << EDescentToString(DESCENT_NEWTONKRYLOV);
267  hist << " using " << krylovName_;
268  if ( useSecantPrecond_ ) {
269  hist << " with " << ESecantToString(esec_) << " preconditioning";
270  }
271  hist << "\n";
272  return hist.str();
273  }
274  std::string print( AlgorithmState<Real> &algo_state, bool print_header = false ) const {
275  std::stringstream hist;
276  hist << std::scientific << std::setprecision(6);
277  if ( algo_state.iter == 0 ) {
278  hist << printName();
279  }
280  if ( print_header ) {
281  hist << printHeader();
282  }
283  if ( algo_state.iter == 0 ) {
284  hist << " ";
285  hist << std::setw(6) << std::left << algo_state.iter;
286  hist << std::setw(15) << std::left << algo_state.value;
287  hist << std::setw(15) << std::left << algo_state.gnorm;
288  hist << "\n";
289  }
290  else {
291  hist << " ";
292  hist << std::setw(6) << std::left << algo_state.iter;
293  hist << std::setw(15) << std::left << algo_state.value;
294  hist << std::setw(15) << std::left << algo_state.gnorm;
295  hist << std::setw(15) << std::left << algo_state.snorm;
296  hist << std::setw(10) << std::left << algo_state.nfval;
297  hist << std::setw(10) << std::left << algo_state.ngrad;
298  hist << std::setw(10) << std::left << iterKrylov_;
299  hist << std::setw(10) << std::left << flagKrylov_;
300  hist << "\n";
301  }
302  return hist.str();
303  }
304 }; // class NewtonKrylovStep
305 
306 } // namespace ROL
307 
308 #endif
Provides the interface to evaluate objective functions.
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: ROL_Vector.hpp:192
virtual void scale(const Real alpha)=0
Compute where .
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
PrecondNK(const ROL::Ptr< Objective< Real > > &obj, const ROL::Ptr< Vector< Real > > &x)
virtual void plus(const Vector &x)=0
Compute , where .
NewtonKrylovStep(ROL::ParameterList &parlist, const ROL::Ptr< Krylov< Real > > &krylov, const ROL::Ptr< Secant< Real > > &secant, const bool computeObj=true)
Constructor.
void update(Vector< Real > &x, const Vector< Real > &s, Objective< Real > &obj, BoundConstraint< Real > &bnd, AlgorithmState< Real > &algo_state)
Update step, if successful.
virtual Real value(const Vector< Real > &x, Real &tol)=0
Compute value.
Provides the interface to compute optimization steps.
Definition: ROL_Step.hpp:34
ROL::Ptr< Secant< Real > > secant_
Secant object (used for quasi-Newton)
int verbosity_
Verbosity level.
Contains definitions of custom data types in ROL.
std::string printName(void) const
Print step name.
ROL::Ptr< Krylov< Real > > krylov_
Krylov solver object (used for inexact Newton)
ESecant StringToESecant(std::string s)
Definition: ROL_Types.hpp:509
std::string EDescentToString(EDescent tr)
Definition: ROL_Types.hpp:386
int flagKrylov_
Termination flag for Krylov method (used for inexact Newton)
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
virtual void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update objective function.
const ROL::Ptr< Vector< Real > > x_
EKrylov
Enumeration of Krylov methods.
EKrylov StringToEKrylov(std::string s)
State for algorithm class. Will be used for restarts.
Definition: ROL_Types.hpp:109
void initialize(Vector< Real > &x, const Vector< Real > &s, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &bnd, AlgorithmState< Real > &algo_state)
Initialize step with bound constraint.
void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
ESecant
Enumeration of secant update algorithms.
Definition: ROL_Types.hpp:452
ROL::Ptr< StepState< Real > > getState(void)
Definition: ROL_Step.hpp:39
HessianNK(const ROL::Ptr< Objective< Real > > &obj, const ROL::Ptr< Vector< Real > > &x)
NewtonKrylovStep(ROL::ParameterList &parlist, const bool computeObj=true)
Constructor.
void compute(Vector< Real > &s, const Vector< Real > &x, Objective< Real > &obj, BoundConstraint< Real > &bnd, AlgorithmState< Real > &algo_state)
Compute step.
Provides interface for and implements limited-memory secant operators.
Definition: ROL_Secant.hpp:45
ROL::Ptr< Vector< Real > > iterateVec
Definition: ROL_Types.hpp:123
std::string printHeader(void) const
Print iterate header.
void applyInverse(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply inverse of linear operator.
ROL::Ptr< Vector< Real > > gp_
Provides definitions for Krylov solvers.
Definition: ROL_Krylov.hpp:24
Provides the interface to apply a linear operator.
const ROL::Ptr< Vector< Real > > x_
Provides the interface to apply upper and lower bound constraints.
Provides the interface to compute optimization steps with projected inexact Newton&#39;s method using lin...
std::string print(AlgorithmState< Real > &algo_state, bool print_header=false) const
Print iterate status.
int iterKrylov_
Number of Krylov iterations (used for inexact Newton)
virtual void initialize(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &con, AlgorithmState< Real > &algo_state)
Initialize step with bound constraint.
Definition: ROL_Step.hpp:54
const ROL::Ptr< Objective< Real > > obj_
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:175
virtual Real norm() const =0
Returns where .
bool useSecantPrecond_
Whether or not a secant approximation is used for preconditioning inexact Newton. ...
std::string ESecantToString(ESecant tr)
Definition: ROL_Types.hpp:461
const ROL::Ptr< Objective< Real > > obj_