ROL
ROL_FletcherStep.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_FLETCHERSTEP_H
45#define ROL_FLETCHERSTEP_H
46
47#include "ROL_FletcherBase.hpp"
48#include "ROL_Step.hpp"
51#include "ROL_Types.hpp"
52#include "ROL_ParameterList.hpp"
59namespace ROL {
60
61template <class Real>
62class FletcherStep : public Step<Real> {
63private:
64 ROL::Ptr<Step<Real> > step_;
65 ROL::Ptr<BoundConstraint<Real> > bnd_;
66
67 ROL::ParameterList parlist_;
68
69 ROL::Ptr<Vector<Real> > x_;
70
71 // Lagrange multiplier update
76 // Subproblem information
77 bool print_;
78 std::string subStep_;
79
80 Real delta_;
84
86
87 ROL::Ptr<Vector<Real> > g_;
88
90
91 // For printing output
92 mutable bool isDeltaChanged_;
93 mutable bool isPenaltyChanged_;
94
96
97 mutable int stepHeaderLength_; // For formatting
98
101 Real gnorm = 0.;
102 // Compute norm of projected gradient
103 if (bnd.isActivated()) {
104 x_->set(x);
105 x_->axpy(-1.,g.dual());
106 bnd.project(*x_);
107 x_->axpy(-1.,x);
108 gnorm = x_->norm();
109 }
110 else {
111 gnorm = g.norm();
112 }
113 return gnorm;
114 }
115
116public:
117
118 using Step<Real>::initialize;
119 using Step<Real>::compute;
120 using Step<Real>::update;
121
123
124 FletcherStep(ROL::ParameterList &parlist)
125 : Step<Real>(), bnd_activated_(false), numSuccessSteps_(0),
127 Real zero(0), one(1), two(2), oe8(1.e8), oe1(1.e-1), oem6(1e-6), oem8(1.e-8);
128
129 ROL::ParameterList& sublist = parlist.sublist("Step").sublist("Fletcher");
130 Step<Real>::getState()->searchSize = sublist.get("Penalty Parameter",one);
131 delta_ = sublist.get("Regularization Parameter",zero);
132 deltaMin_ = sublist.get("Min Regularization Parameter",oem8);
133 deltaUpdate_ = sublist.get("Regularization Parameter Decrease Factor", oe1);
134 // penalty parameters
135 penaltyUpdate_ = sublist.get("Penalty Parameter Growth Factor", two);
136 modifyPenalty_ = sublist.get("Modify Penalty Parameter", false);
137 maxPenaltyParam_ = sublist.get("Maximum Penalty Parameter", oe8);
138 minPenaltyParam_ = sublist.get("Minimum Penalty Parameter", oem6);
139
140 subStep_ = sublist.get("Subproblem Solver", "Trust Region");
141
142 parlist_ = parlist;
143 }
144
149 AlgorithmState<Real> &algo_state ) {
150 bnd_ = ROL::makePtr<BoundConstraint<Real>>();
151 bnd_->deactivate();
152 initialize(x,g,l,c,obj,con,*bnd_,algo_state);
153 }
154
159 AlgorithmState<Real> &algo_state ) {
160 // Determine what kind of step
162
163 ROL::ParameterList trlist(parlist_);
164 bool inexactFletcher = trlist.sublist("Step").sublist("Fletcher").get("Inexact Solves", false);
165 if( inexactFletcher ) {
166 trlist.sublist("General").set("Inexact Objective Value", true);
167 trlist.sublist("General").set("Inexact Gradient", true);
168 }
169 if( bnd_activated_ ) {
170 trlist.sublist("Step").sublist("Trust Region").set("Subproblem Model", "Coleman-Li");
171 }
172
173 if ( subStep_ == "Line Search" ) {
174 step_ = makePtr<LineSearchStep<Real>>(trlist);
175 }
176 else {
177 step_ = makePtr<TrustRegionStep<Real>>(trlist);
178 }
179 etr_ = StringToETrustRegion(parlist_.sublist("Step").sublist("Trust Region").get("Subproblem Solver", "Truncated CG"));
180
181 // Initialize class members
182 g_ = g.clone();
183 x_ = x.clone();
184
185 // Rest of initialize
186 FletcherBase<Real>& fletcher = dynamic_cast<FletcherBase<Real>&>(obj);
187
188 tr_algo_state_.iterateVec = x.clone();
189 tr_algo_state_.minIterVec = x.clone();
190 tr_algo_state_.lagmultVec = l.clone();
191
192 step_->initialize(x, g, obj, bnd, tr_algo_state_);
193
194 // Initialize step state
195 ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
196 state->descentVec = x.clone();
197 state->gradientVec = g.clone();
198 state->constraintVec = c.clone();
199 // Initialize the algorithm state
200 algo_state.nfval = 0;
201 algo_state.ncval = 0;
202 algo_state.ngrad = 0;
203
204 algo_state.value = fletcher.getObjectiveValue(x);
205 algo_state.gnorm = computeProjGradientNorm(*(fletcher.getLagrangianGradient(x)),
206 x, bnd);
207 algo_state.aggregateGradientNorm = tr_algo_state_.gnorm;
208
209 state->constraintVec->set(*(fletcher.getConstraintVec(x)));
210 algo_state.cnorm = (state->constraintVec)->norm();
211 // Update evaluation counters
212 algo_state.ncval = fletcher.getNumberConstraintEvaluations();
213 algo_state.nfval = fletcher.getNumberFunctionEvaluations();
214 algo_state.ngrad = fletcher.getNumberGradientEvaluations();
215 }
216
219 void compute( Vector<Real> &s, const Vector<Real> &x, const Vector<Real> &l,
221 AlgorithmState<Real> &algo_state ) {
222 compute(s,x,l,obj,con,*bnd_, algo_state);
223 }
224
227 void compute( Vector<Real> &s, const Vector<Real> &x, const Vector<Real> &l,
229 BoundConstraint<Real> &bnd, AlgorithmState<Real> &algo_state ) {
230 step_->compute( s, x, obj, bnd, tr_algo_state_ );
231 }
232
237 AlgorithmState<Real> &algo_state ) {
238 update(x,l,s,obj,con,*bnd_, algo_state);
239 }
240
246 AlgorithmState<Real> &algo_state ) {
247
248 // This should be in print, but this will not work there
249 isDeltaChanged_ = false;
250 isPenaltyChanged_ = false;
251 bool modified = false;
252
253 FletcherBase<Real> &fletcher = dynamic_cast<FletcherBase<Real>&>(obj);
254 ROL::Ptr<StepState<Real> > fletcherState = Step<Real>::getState();
255 const ROL::Ptr<const StepState<Real> > state = step_->getStepState();
256
257 step_->update(x,s,obj,bnd,tr_algo_state_);
258 numSuccessSteps_ += (state->flag == 0);
259
260 Real gPhiNorm = tr_algo_state_.gnorm;
261 Real cnorm = (fletcherState->constraintVec)->norm();
262 bool too_infeasible = cnorm > static_cast<Real>(100.)*gPhiNorm;
263 bool too_feasible = cnorm < static_cast<Real>(1e-2)*gPhiNorm;
264
265 if( too_infeasible && !modified && modifyPenalty_ && numSuccessSteps_ > 1 ) {
266 Real penaltyParam = Step<Real>::getStepState()->searchSize;
267 if( penaltyParam >= maxPenaltyParam_ ) {
268 // Penalty parameter too large, exit
269 algo_state.flag = true;
270 }
271 penaltyParam *= penaltyUpdate_;
272 penaltyParam = std::min(penaltyParam, maxPenaltyParam_);
273 fletcher.setPenaltyParameter(penaltyParam);
274 Step<Real>::getState()->searchSize = penaltyParam;
275 isPenaltyChanged_ = true;
276 modified = true;
277 }
278
279 if( too_feasible && !modified && modifyPenalty_ && numSuccessSteps_ > 1 ) {
280 Real penaltyParam = Step<Real>::getStepState()->searchSize;
281 if( penaltyParam <= minPenaltyParam_ ) {
282 // Penalty parameter too small, exit (this is unlikely)
283 algo_state.flag = true;
284 }
285 penaltyParam /= penaltyUpdate_;
286 penaltyParam = std::max(penaltyParam, minPenaltyParam_);
287 fletcher.setPenaltyParameter(penaltyParam);
288 Step<Real>::getState()->searchSize = penaltyParam;
289 isPenaltyChanged_ = true;
290 modified = true;
291 }
292
293 if( delta_ > deltaMin_ && !modified ) {
294 Real deltaNext = delta_ * deltaUpdate_;
295 if( gPhiNorm < deltaNext ) {
296 delta_ = deltaNext;
297 fletcher.setDelta(deltaNext);
298 isDeltaChanged_ = true;
299 modified = true;
300 }
301 }
302
303 if( modified ) {
304 // Penalty function has been changed somehow, need to recompute
305 Real tol = static_cast<Real>(1e-12);
306 tr_algo_state_.value = fletcher.value(x, tol);
307 fletcher.gradient(*g_, x, tol);
308
309 tr_algo_state_.nfval++;
310 tr_algo_state_.ngrad++;
311 tr_algo_state_.ncval++;
312 tr_algo_state_.minIter = tr_algo_state_.iter;
313 tr_algo_state_.minValue = tr_algo_state_.value;
314 tr_algo_state_.gnorm = computeProjGradientNorm(*g_, x, bnd);
315 }
316
317 // Update the step and store in state
318 algo_state.iterateVec->set(x);
319 algo_state.iter++;
320
321 fletcherState->descentVec->set(s);
322 fletcherState->gradientVec->set(*(fletcher.getLagrangianGradient(x)));
323 fletcherState->constraintVec->set(*(fletcher.getConstraintVec(x)));
324
325 // Update objective function value
326 algo_state.value = fletcher.getObjectiveValue(x);
327 // Update constraint value
328 algo_state.cnorm = (fletcherState->constraintVec)->norm();
329 // Update the step size
330 algo_state.snorm = tr_algo_state_.snorm;
331 // Compute gradient of the Lagrangian
332 algo_state.gnorm = computeProjGradientNorm(*(fletcherState->gradientVec),
333 x, bnd);
334 // Compute gradient of penalty function
335 algo_state.aggregateGradientNorm = tr_algo_state_.gnorm;
336 // Update evaluation countersgetConstraintVec
337 algo_state.nfval = fletcher.getNumberFunctionEvaluations();
338 algo_state.ngrad = fletcher.getNumberGradientEvaluations();
339 algo_state.ncval = fletcher.getNumberConstraintEvaluations();
340 // Update objective function and constraints
341 // fletcher.update(x,true,algo_state.iter);
342 // Update multipliers
343 algo_state.lagmultVec->set(*(fletcher.getMultiplierVec(x)));
344 }
345
348 std::string printHeader( void ) const {
349 std::stringstream hist;
350 if( subStep_ == "Trust Region" ) {
351 hist << " ";
352 hist << std::setw(6) << std::left << "iter";
353 hist << std::setw(15) << std::left << "merit";
354 hist << std::setw(15) << std::left << "fval";
355 hist << std::setw(15) << std::left << "gpnorm";
356 hist << std::setw(15) << std::left << "gLnorm";
357 hist << std::setw(15) << std::left << "cnorm";
358 hist << std::setw(15) << std::left << "snorm";
359 hist << std::setw(15) << std::left << "tr_radius";
360 hist << std::setw(10) << std::left << "tr_flag";
361 if ( etr_ == TRUSTREGION_TRUNCATEDCG && subStep_ == "Trust Region") {
362 hist << std::setw(10) << std::left << "iterCG";
363 hist << std::setw(10) << std::left << "flagCG";
364 }
365 hist << std::setw(15) << std::left << "penalty";
366 hist << std::setw(15) << std::left << "delta";
367 hist << std::setw(10) << std::left << "#fval";
368 hist << std::setw(10) << std::left << "#grad";
369 hist << std::setw(10) << std::left << "#cval";
370 hist << "\n";
371 }
372 else {
373 std::string stepHeader = step_->printHeader();
374 stepHeaderLength_ = stepHeader.length();
375 hist << stepHeader.substr(0, stepHeaderLength_-1);
376 hist << std::setw(15) << std::left << "fval";
377 hist << std::setw(15) << std::left << "gLnorm";
378 hist << std::setw(15) << std::left << "cnorm";
379 hist << std::setw(15) << std::left << "penalty";
380 hist << std::setw(15) << std::left << "delta";
381 hist << std::setw(10) << std::left << "#cval";
382 hist << "\n";
383 }
384 return hist.str();
385 }
386
389 std::string printName( void ) const {
390 std::stringstream hist;
391 hist << "\n" << " Fletcher solver : " << subStep_;
392 hist << "\n";
393 return hist.str();
394 }
395
398 std::string print( AlgorithmState<Real> &algo_state, bool pHeader = false ) const {
399 std::string stepHist = step_->print( tr_algo_state_, false );
400 stepHist.erase(std::remove(stepHist.end()-3, stepHist.end(),'\n'), stepHist.end());
401 std::string name = step_->printName();
402 size_t pos = stepHist.find(name);
403 if ( pos != std::string::npos ) {
404 stepHist.erase(pos, name.length());
405 }
406
407 std::stringstream hist;
408 hist << std::scientific << std::setprecision(6);
409 if ( algo_state.iter == 0 ) {
410 hist << printName();
411 }
412 if ( pHeader ) {
413 hist << printHeader();
414 }
415
416 std::string penaltyString = getValueString( Step<Real>::getStepState()->searchSize, isPenaltyChanged_ );
417 std::string deltaString = getValueString( delta_, isDeltaChanged_ );
418
419 if( subStep_ == "Trust Region" ) {
420 hist << " ";
421 hist << std::setw(6) << std::left << algo_state.iter;
422 hist << std::setw(15) << std::left << tr_algo_state_.value;
423 hist << std::setw(15) << std::left << algo_state.value;
424 hist << std::setw(15) << std::left << tr_algo_state_.gnorm;
425 hist << std::setw(15) << std::left << algo_state.gnorm;
426 hist << std::setw(15) << std::left << algo_state.cnorm;
427 hist << std::setw(15) << std::left << stepHist.substr(38,15); // snorm
428 hist << std::setw(15) << std::left << stepHist.substr(53,15); // tr_radius
429 hist << std::setw(10) << std::left << (algo_state.iter == 0 ? "" : stepHist.substr(88,10)); // tr_flag
430 if ( etr_ == TRUSTREGION_TRUNCATEDCG && subStep_ == "Trust Region") {
431 hist << std::setw(10) << std::left << (algo_state.iter == 0 ? "" : stepHist.substr(93,10)); // iterCG
432 hist << std::setw(10) << std::left << (algo_state.iter == 0 ? "" : stepHist.substr(103,10)); // flagCG
433 }
434 hist << std::setw(15) << std::left << penaltyString;
435 hist << std::setw(15) << std::left << deltaString;
436 hist << std::setw(10) << std::left << (algo_state.iter == 0 ? "" : stepHist.substr(68,10)); // #fval
437 hist << std::setw(10) << std::left << (algo_state.iter == 0 ? "" : stepHist.substr(78,10)); // #gval
438 hist << std::setw(10) << std::left << algo_state.ncval;
439 hist << "\n";
440 } else {
441 hist << std::setw(stepHeaderLength_-1) << std::left << stepHist;
442 hist << std::setw(15) << std::left << algo_state.value;
443 hist << std::setw(15) << std::left << algo_state.gnorm;
444 hist << std::setw(15) << std::left << algo_state.cnorm;
445 hist << std::setw(15) << std::left << penaltyString;
446 hist << std::setw(15) << std::left << deltaString;
447 hist << std::setw(10) << std::left << algo_state.ncval;
448 hist << "\n";
449 }
450
451 return hist.str();
452 }
453
454 std::string getValueString( const Real value, const bool print ) const {
455 std::stringstream valString;
456 valString << std::scientific << std::setprecision(6);
457 if( print ) {
458 valString << std::setw(15) << std::left << value;
459 } else {
460 valString << std::setw(15) << "";
461 }
462 return valString.str();
463 }
464
470 AlgorithmState<Real> &algo_state ) {}
471
477 AlgorithmState<Real> &algo_state ) {}
478
479}; // class FletcherStep
480
481} // namespace ROL
482
483#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0 zero)()
Contains definitions of custom data types in ROL.
Provides the interface to apply upper and lower bound constraints.
bool isActivated(void) const
Check if bounds are on.
virtual void project(Vector< Real > &x)
Project optimization variables onto the bounds.
Defines the general constraint operator interface.
const Ptr< Vector< Real > > getLagrangianGradient(const Vector< Real > &x)
int getNumberConstraintEvaluations() const
int getNumberGradientEvaluations() const
void setDelta(Real delta)
const Ptr< Vector< Real > > getMultiplierVec(const Vector< Real > &x)
Real getObjectiveValue(const Vector< Real > &x)
int getNumberFunctionEvaluations() const
const Ptr< Vector< Real > > getConstraintVec(const Vector< Real > &x)
void setPenaltyParameter(Real sigma)
Provides the interface to compute Fletcher steps.
void update(Vector< Real > &x, Vector< Real > &l, const Vector< Real > &s, Objective< Real > &obj, Constraint< Real > &con, BoundConstraint< Real > &bnd, AlgorithmState< Real > &algo_state)
Update step, if successful (equality and bound constraints).
std::string printHeader(void) const
Print iterate header.
ROL::Ptr< BoundConstraint< Real > > bnd_
ROL::ParameterList parlist_
void update(Vector< Real > &x, const Vector< Real > &s, Objective< Real > &obj, BoundConstraint< Real > &con, AlgorithmState< Real > &algo_state)
Update step, for bound constraints; here only to satisfy the interface requirements,...
FletcherStep(ROL::ParameterList &parlist)
void compute(Vector< Real > &s, const Vector< Real > &x, const Vector< Real > &l, Objective< Real > &obj, Constraint< Real > &con, AlgorithmState< Real > &algo_state)
Compute step (equality constraint).
void compute(Vector< Real > &s, const Vector< Real > &x, Objective< Real > &obj, BoundConstraint< Real > &con, AlgorithmState< Real > &algo_state)
Compute step for bound constraints; here only to satisfy the interface requirements,...
ROL::Ptr< Step< Real > > step_
void initialize(Vector< Real > &x, const Vector< Real > &g, Vector< Real > &l, const Vector< Real > &c, Objective< Real > &obj, Constraint< Real > &con, AlgorithmState< Real > &algo_state)
Initialize step with equality constraint.
void initialize(Vector< Real > &x, const Vector< Real > &g, Vector< Real > &l, const Vector< Real > &c, Objective< Real > &obj, Constraint< Real > &con, BoundConstraint< Real > &bnd, AlgorithmState< Real > &algo_state)
Initialize step with equality and bound constraints.
std::string getValueString(const Real value, const bool print) const
std::string printName(void) const
Print step name.
Real computeProjGradientNorm(const Vector< Real > &g, const Vector< Real > &x, BoundConstraint< Real > &bnd)
std::string print(AlgorithmState< Real > &algo_state, bool pHeader=false) const
Print iterate status.
ROL::Ptr< Vector< Real > > g_
AlgorithmState< Real > tr_algo_state_
void compute(Vector< Real > &s, const Vector< Real > &x, const Vector< Real > &l, Objective< Real > &obj, Constraint< Real > &con, BoundConstraint< Real > &bnd, AlgorithmState< Real > &algo_state)
Compute step (equality and bound constraints).
void update(Vector< Real > &x, Vector< Real > &l, const Vector< Real > &s, Objective< Real > &obj, Constraint< Real > &con, AlgorithmState< Real > &algo_state)
Update step, if successful (equality constraint).
ROL::Ptr< Vector< Real > > x_
Provides the interface to evaluate objective functions.
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
virtual Real value(const Vector< Real > &x, Real &tol)=0
Compute value.
Provides the interface to compute optimization steps.
Definition ROL_Step.hpp:68
ROL::Ptr< StepState< Real > > getState(void)
Definition ROL_Step.hpp:73
const ROL::Ptr< const StepState< Real > > getStepState(void) const
Get state for step object.
Definition ROL_Step.hpp:211
Defines the linear algebra or vector space interface.
virtual Real norm() const =0
Returns where .
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
ROL::Objective_SerialSimOpt Objective_SimOpt value(const V &u, const V &z, Real &tol) override
ETrustRegion StringToETrustRegion(std::string s)
State for algorithm class. Will be used for restarts.
ROL::Ptr< Vector< Real > > lagmultVec
ROL::Ptr< Vector< Real > > iterateVec