ROL
ROL_ProjectedNewtonKrylovStep.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Rapid Optimization Library (ROL) Package
5 // Copyright (2014) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact lead developers:
38 // Drew Kouri (dpkouri@sandia.gov) and
39 // Denis Ridzal (dridzal@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
44 #ifndef ROL_PROJECTEDNEWTONKRYLOVSTEP_H
45 #define ROL_PROJECTEDNEWTONKRYLOVSTEP_H
46 
47 #include "ROL_Types.hpp"
48 #include "ROL_Step.hpp"
49 
50 #include "ROL_Secant.hpp"
51 #include "ROL_Krylov.hpp"
52 #include "ROL_LinearOperator.hpp"
53 
54 #include <sstream>
55 #include <iomanip>
56 
64 namespace ROL {
65 
66 template <class Real>
67 class ProjectedNewtonKrylovStep : public Step<Real> {
68 private:
69 
70  Teuchos::RCP<Secant<Real> > secant_;
71  Teuchos::RCP<Krylov<Real> > krylov_;
72 
75 
76  Teuchos::RCP<Vector<Real> > gp_;
77  Teuchos::RCP<Vector<Real> > d_;
78 
81  int verbosity_;
82  const bool computeObj_;
83 
86 
87  class HessianPNK : public LinearOperator<Real> {
88  private:
89  const Teuchos::RCP<Objective<Real> > obj_;
90  const Teuchos::RCP<BoundConstraint<Real> > bnd_;
91  const Teuchos::RCP<Vector<Real> > x_;
92  const Teuchos::RCP<Vector<Real> > g_;
93  Teuchos::RCP<Vector<Real> > v_;
94  Real eps_;
95  public:
96  HessianPNK(const Teuchos::RCP<Objective<Real> > &obj,
97  const Teuchos::RCP<BoundConstraint<Real> > &bnd,
98  const Teuchos::RCP<Vector<Real> > &x,
99  const Teuchos::RCP<Vector<Real> > &g,
100  Real eps = 0 )
101  : obj_(obj), bnd_(bnd), x_(x), g_(g), eps_(eps) {
102  v_ = x_->clone();
103  }
104  void apply(Vector<Real> &Hv, const Vector<Real> &v, Real &tol) const {
105  v_->set(v);
106  bnd_->pruneActive(*v_,*g_,*x_,eps_);
107  obj_->hessVec(Hv,*v_,*x_,tol);
108  bnd_->pruneActive(Hv,*g_,*x_,eps_);
109  v_->set(v);
110  bnd_->pruneInactive(*v_,*g_,*x_,eps_);
111  Hv.plus(v_->dual());
112  }
113  };
114 
115  class PrecondPNK : public LinearOperator<Real> {
116  private:
117  const Teuchos::RCP<Objective<Real> > obj_;
118  const Teuchos::RCP<Secant<Real> > secant_;
119  const Teuchos::RCP<BoundConstraint<Real> > bnd_;
120  const Teuchos::RCP<Vector<Real> > x_;
121  const Teuchos::RCP<Vector<Real> > g_;
122  Teuchos::RCP<Vector<Real> > v_;
123  Real eps_;
124  const bool useSecant_;
125  public:
126  PrecondPNK(const Teuchos::RCP<Objective<Real> > &obj,
127  const Teuchos::RCP<BoundConstraint<Real> > &bnd,
128  const Teuchos::RCP<Vector<Real> > &x,
129  const Teuchos::RCP<Vector<Real> > &g,
130  Real eps = 0 )
131  : obj_(obj), bnd_(bnd), x_(x), g_(g), eps_(eps), useSecant_(false) {
132  v_ = x_->clone();
133  }
134  PrecondPNK(const Teuchos::RCP<Secant<Real> > &secant,
135  const Teuchos::RCP<BoundConstraint<Real> > &bnd,
136  const Teuchos::RCP<Vector<Real> > &x,
137  const Teuchos::RCP<Vector<Real> > &g,
138  Real eps = 0 )
139  : secant_(secant), bnd_(bnd), x_(x), g_(g), eps_(eps), useSecant_(true) {
140  v_ = x_->clone();
141  }
142  void apply(Vector<Real> &Hv, const Vector<Real> &v, Real &tol) const {
143  Hv.set(v.dual());
144  }
145  void applyInverse(Vector<Real> &Hv, const Vector<Real> &v, Real &tol) const {
146  v_->set(v);
147  bnd_->pruneActive(*v_,*g_,*x_,eps_);
148  if ( useSecant_ ) {
149  secant_->applyH(Hv,*v_);
150  }
151  else {
152  obj_->precond(Hv,*v_,*x_,tol);
153  }
154  bnd_->pruneActive(Hv,*g_,*x_,eps_);
155  v_->set(v);
156  bnd_->pruneInactive(*v_,*g_,*x_,eps_);
157  Hv.plus(v_->dual());
158  }
159  };
160 
161 public:
162 
164  using Step<Real>::compute;
165  using Step<Real>::update;
166 
174  ProjectedNewtonKrylovStep( Teuchos::ParameterList &parlist, const bool computeObj = true )
175  : Step<Real>(), secant_(Teuchos::null), krylov_(Teuchos::null),
176  gp_(Teuchos::null), d_(Teuchos::null),
178  computeObj_(computeObj), useSecantPrecond_(false) {
179  // Parse ParameterList
180  Teuchos::ParameterList& Glist = parlist.sublist("General");
181  useSecantPrecond_ = Glist.sublist("Secant").get("Use as Preconditioner", false);
182  useProjectedGrad_ = Glist.get("Projected Gradient Criticality Measure", false);
183  verbosity_ = Glist.get("Print Verbosity",0);
184  // Initialize Krylov object
185  ekv_ = StringToEKrylov(Glist.sublist("Krylov").get("Type","Conjugate Gradients"));
186  krylov_ = KrylovFactory<Real>(parlist);
187  // Initialize secant object
188  esec_ = StringToESecant(Glist.sublist("Secant").get("Type","Limited-Memory BFGS"));
189  if ( useSecantPrecond_ ) {
190  secant_ = SecantFactory<Real>(parlist);
191  }
192  }
193 
204  ProjectedNewtonKrylovStep(Teuchos::ParameterList &parlist,
205  const Teuchos::RCP<Krylov<Real> > &krylov,
206  const Teuchos::RCP<Secant<Real> > &secant,
207  const bool computeObj = true)
208  : Step<Real>(), secant_(secant), krylov_(krylov),
210  gp_(Teuchos::null), d_(Teuchos::null),
212  computeObj_(computeObj), useSecantPrecond_(false) {
213  // Parse ParameterList
214  Teuchos::ParameterList& Glist = parlist.sublist("General");
215  useSecantPrecond_ = Glist.sublist("Secant").get("Use as Preconditioner", false);
216  useProjectedGrad_ = Glist.get("Projected Gradient Criticality Measure", false);
217  verbosity_ = Glist.get("Print Verbosity",0);
218  // Initialize secant object
219  if ( useSecantPrecond_ && secant_ == Teuchos::null ) {
220  esec_ = StringToESecant(Glist.sublist("Secant").get("Type","Limited-Memory BFGS"));
221  secant_ = SecantFactory<Real>(parlist);
222  }
223  // Initialize Krylov object
224  if ( krylov_ == Teuchos::null ) {
225  ekv_ = StringToEKrylov(Glist.sublist("Krylov").get("Type","Conjugate Gradients"));
226  krylov_ = KrylovFactory<Real>(parlist);
227  }
228  }
229 
230  void initialize( Vector<Real> &x, const Vector<Real> &s, const Vector<Real> &g,
232  AlgorithmState<Real> &algo_state ) {
233  Step<Real>::initialize(x,s,g,obj,bnd,algo_state);
234  gp_ = g.clone();
235  d_ = s.clone();
236  }
237 
238  void compute( Vector<Real> &s, const Vector<Real> &x,
240  AlgorithmState<Real> &algo_state ) {
241  Real one(1);
242  Teuchos::RCP<StepState<Real> > step_state = Step<Real>::getState();
243 
244  // Build Hessian and Preconditioner object
245  Teuchos::RCP<Objective<Real> > obj_ptr = Teuchos::rcpFromRef(obj);
246  Teuchos::RCP<BoundConstraint<Real> > bnd_ptr = Teuchos::rcpFromRef(bnd);
247  Teuchos::RCP<LinearOperator<Real> > hessian
248  = Teuchos::rcp(new HessianPNK(obj_ptr,bnd_ptr,algo_state.iterateVec,
249  step_state->gradientVec,algo_state.gnorm));
250  Teuchos::RCP<LinearOperator<Real> > precond;
251  if (useSecantPrecond_) {
252  precond = Teuchos::rcp(new PrecondPNK(secant_,bnd_ptr,
253  algo_state.iterateVec,step_state->gradientVec,algo_state.gnorm));
254  }
255  else {
256  precond = Teuchos::rcp(new PrecondPNK(obj_ptr,bnd_ptr,
257  algo_state.iterateVec,step_state->gradientVec,algo_state.gnorm));
258  }
259 
260  // Run Krylov method
261  flagKrylov_ = 0;
262  krylov_->run(s,*hessian,*(step_state->gradientVec),*precond,iterKrylov_,flagKrylov_);
263 
264  // Check Krylov flags
265  if ( flagKrylov_ == 2 && iterKrylov_ <= 1 ) {
266  s.set((step_state->gradientVec)->dual());
267  }
268  s.scale(-one);
269  }
270 
271  void update( Vector<Real> &x, const Vector<Real> &s,
273  AlgorithmState<Real> &algo_state ) {
274  Real tol = std::sqrt(ROL_EPSILON<Real>()), one(1);
275  Teuchos::RCP<StepState<Real> > step_state = Step<Real>::getState();
276 
277  // Update iterate and store previous step
278  algo_state.iter++;
279  d_->set(x);
280  x.plus(s);
281  bnd.project(x);
282  (step_state->descentVec)->set(x);
283  (step_state->descentVec)->axpy(-one,*d_);
284  algo_state.snorm = s.norm();
285 
286  // Compute new gradient
287  if ( useSecantPrecond_ ) {
288  gp_->set(*(step_state->gradientVec));
289  }
290  obj.update(x,true,algo_state.iter);
291  if ( computeObj_ ) {
292  algo_state.value = obj.value(x,tol);
293  algo_state.nfval++;
294  }
295  obj.gradient(*(step_state->gradientVec),x,tol);
296  algo_state.ngrad++;
297 
298  // Update Secant Information
299  if ( useSecantPrecond_ ) {
300  secant_->updateStorage(x,*(step_state->gradientVec),*gp_,s,algo_state.snorm,algo_state.iter+1);
301  }
302 
303  // Update algorithm state
304  (algo_state.iterateVec)->set(x);
305  if ( useProjectedGrad_ ) {
306  gp_->set(*(step_state->gradientVec));
307  bnd.computeProjectedGradient( *gp_, x );
308  algo_state.gnorm = gp_->norm();
309  }
310  else {
311  d_->set(x);
312  d_->axpy(-one,(step_state->gradientVec)->dual());
313  bnd.project(*d_);
314  d_->axpy(-one,x);
315  algo_state.gnorm = d_->norm();
316  }
317  }
318 
319  std::string printHeader( void ) const {
320  std::stringstream hist;
321 
322  if( verbosity_>0 ) {
323  hist << std::string(109,'-') << "\n";
325  hist << " status output definitions\n\n";
326  hist << " iter - Number of iterates (steps taken) \n";
327  hist << " value - Objective function value \n";
328  hist << " gnorm - Norm of the gradient\n";
329  hist << " snorm - Norm of the step (update to optimization vector)\n";
330  hist << " #fval - Cumulative number of times the objective function was evaluated\n";
331  hist << " #grad - Number of times the gradient was computed\n";
332  hist << " iterCG - Number of Krylov iterations used to compute search direction\n";
333  hist << " flagCG - Krylov solver flag" << "\n";
334  hist << std::string(109,'-') << "\n";
335  }
336 
337  hist << " ";
338  hist << std::setw(6) << std::left << "iter";
339  hist << std::setw(15) << std::left << "value";
340  hist << std::setw(15) << std::left << "gnorm";
341  hist << std::setw(15) << std::left << "snorm";
342  hist << std::setw(10) << std::left << "#fval";
343  hist << std::setw(10) << std::left << "#grad";
344  hist << std::setw(10) << std::left << "iterCG";
345  hist << std::setw(10) << std::left << "flagCG";
346  hist << "\n";
347  return hist.str();
348  }
349  std::string printName( void ) const {
350  std::stringstream hist;
351  hist << "\n" << EDescentToString(DESCENT_NEWTONKRYLOV);
352  hist << " using " << EKrylovToString(ekv_);
353  if ( useSecantPrecond_ ) {
354  hist << " with " << ESecantToString(esec_) << " preconditioning";
355  }
356  hist << "\n";
357  return hist.str();
358  }
359  std::string print( AlgorithmState<Real> &algo_state, bool print_header = false ) const {
360  std::stringstream hist;
361  hist << std::scientific << std::setprecision(6);
362  if ( algo_state.iter == 0 ) {
363  hist << printName();
364  }
365  if ( print_header ) {
366  hist << printHeader();
367  }
368  if ( algo_state.iter == 0 ) {
369  hist << " ";
370  hist << std::setw(6) << std::left << algo_state.iter;
371  hist << std::setw(15) << std::left << algo_state.value;
372  hist << std::setw(15) << std::left << algo_state.gnorm;
373  hist << "\n";
374  }
375  else {
376  hist << " ";
377  hist << std::setw(6) << std::left << algo_state.iter;
378  hist << std::setw(15) << std::left << algo_state.value;
379  hist << std::setw(15) << std::left << algo_state.gnorm;
380  hist << std::setw(15) << std::left << algo_state.snorm;
381  hist << std::setw(10) << std::left << algo_state.nfval;
382  hist << std::setw(10) << std::left << algo_state.ngrad;
383  hist << std::setw(10) << std::left << iterKrylov_;
384  hist << std::setw(10) << std::left << flagKrylov_;
385  hist << "\n";
386  }
387  return hist.str();
388  }
389 }; // class ProjectedNewtonKrylovStep
390 
391 } // namespace ROL
392 
393 #endif
Provides the interface to evaluate objective functions.
bool useSecantPrecond_
Whether or not a secant approximation is used for preconditioning inexact Newton. ...
virtual void scale(const Real alpha)=0
Compute where .
void applyInverse(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply inverse of linear operator.
void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
void update(Vector< Real > &x, const Vector< Real > &s, Objective< Real > &obj, BoundConstraint< Real > &bnd, AlgorithmState< Real > &algo_state)
Update step, if successful.
virtual void plus(const Vector &x)=0
Compute , where .
PrecondPNK(const Teuchos::RCP< Secant< Real > > &secant, const Teuchos::RCP< BoundConstraint< Real > > &bnd, const Teuchos::RCP< Vector< Real > > &x, const Teuchos::RCP< Vector< Real > > &g, Real eps=0)
virtual Real value(const Vector< Real > &x, Real &tol)=0
Compute value.
Provides the interface to compute optimization steps.
Definition: ROL_Step.hpp:69
std::string printName(void) const
Print step name.
Teuchos::RCP< StepState< Real > > getState(void)
Definition: ROL_Step.hpp:74
Contains definitions of custom data types in ROL.
void compute(Vector< Real > &s, const Vector< Real > &x, Objective< Real > &obj, BoundConstraint< Real > &bnd, AlgorithmState< Real > &algo_state)
Compute step.
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
ESecant StringToESecant(std::string s)
Definition: ROL_Types.hpp:477
std::string EDescentToString(EDescent tr)
Definition: ROL_Types.hpp:354
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
Provides the interface to compute optimization steps with projected inexact ProjectedNewton&#39;s method ...
EKrylov
Enumeration of Krylov methods.
Definition: ROL_Types.hpp:493
EKrylov StringToEKrylov(std::string s)
Definition: ROL_Types.hpp:546
State for algorithm class. Will be used for restarts.
Definition: ROL_Types.hpp:91
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
bool useProjectedGrad_
Whether or not to use to the projected gradient criticality measure.
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:213
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:420
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.
int flagKrylov_
Termination flag for Krylov method (used for inexact Newton)
ProjectedNewtonKrylovStep(Teuchos::ParameterList &parlist, const Teuchos::RCP< Krylov< Real > > &krylov, const Teuchos::RCP< Secant< Real > > &secant, const bool computeObj=true)
Constructor.
Provides interface for and implements limited-memory secant operators.
Definition: ROL_Secant.hpp:70
std::string printHeader(void) const
Print iterate header.
std::string print(AlgorithmState< Real > &algo_state, bool print_header=false) const
Print iterate status.
const Teuchos::RCP< BoundConstraint< Real > > bnd_
Teuchos::RCP< Secant< Real > > secant_
Secant object (used for quasi-Newton)
Provides definitions for Krylov solvers.
Definition: ROL_Krylov.hpp:57
Provides the interface to apply a linear operator.
PrecondPNK(const Teuchos::RCP< Objective< Real > > &obj, const Teuchos::RCP< BoundConstraint< Real > > &bnd, const Teuchos::RCP< Vector< Real > > &x, const Teuchos::RCP< Vector< Real > > &g, Real eps=0)
Provides the interface to apply upper and lower bound constraints.
void computeProjectedGradient(Vector< Real > &g, const Vector< Real > &x)
Compute projected gradient.
std::string EKrylovToString(EKrylov tr)
Definition: ROL_Types.hpp:501
int iterKrylov_
Number of Krylov iterations (used for inexact Newton)
Teuchos::RCP< Krylov< Real > > krylov_
Krylov solver object (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:89
ProjectedNewtonKrylovStep(Teuchos::ParameterList &parlist, const bool computeObj=true)
Constructor.
Teuchos::RCP< Vector< Real > > iterateVec
Definition: ROL_Types.hpp:105
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:196
HessianPNK(const Teuchos::RCP< Objective< Real > > &obj, const Teuchos::RCP< BoundConstraint< Real > > &bnd, const Teuchos::RCP< Vector< Real > > &x, const Teuchos::RCP< Vector< Real > > &g, Real eps=0)
virtual Real norm() const =0
Returns where .
virtual void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
const Teuchos::RCP< Objective< Real > > obj_
std::string ESecantToString(ESecant tr)
Definition: ROL_Types.hpp:429
const Teuchos::RCP< BoundConstraint< Real > > bnd_
virtual void project(Vector< Real > &x)
Project optimization variables onto the bounds.