ROL
ROL_MoreauYosidaPenaltyStep.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_MOREAUYOSIDAPENALTYSTEP_H
11 #define ROL_MOREAUYOSIDAPENALTYSTEP_H
12 
14 #include "ROL_Types.hpp"
16 #include "ROL_CompositeStep.hpp"
17 #include "ROL_FletcherStep.hpp"
18 #include "ROL_BundleStep.hpp"
19 #include "ROL_TrustRegionStep.hpp"
20 #include "ROL_LineSearchStep.hpp"
21 #include "ROL_Algorithm.hpp"
23 #include "ROL_BundleStatusTest.hpp"
24 #include "ROL_ParameterList.hpp"
25 
88 namespace ROL {
89 
90 template<class Real>
91 class AugmentedLagrangianStep;
92 
93 template <class Real>
94 class MoreauYosidaPenaltyStep : public Step<Real> {
95 private:
96  ROL::Ptr<StatusTest<Real>> status_;
97  ROL::Ptr<Step<Real>> step_;
98  ROL::Ptr<Algorithm<Real>> algo_;
99  ROL::Ptr<Vector<Real>> x_;
100  ROL::Ptr<Vector<Real>> g_;
101  ROL::Ptr<Vector<Real>> l_;
102  ROL::Ptr<BoundConstraint<Real>> bnd_;
103 
105  Real gLnorm_;
106  Real tau_;
107  bool print_;
109 
110  ROL::ParameterList parlist_;
113 
115  std::string stepname_;
116 
117  void updateState(const Vector<Real> &x, const Vector<Real> &l,
118  Objective<Real> &obj,
120  AlgorithmState<Real> &algo_state) {
122  = dynamic_cast<MoreauYosidaPenalty<Real>&>(obj);
123  Real zerotol = std::sqrt(ROL_EPSILON<Real>());
124  ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
125  // Update objective and constraint.
126  myPen.update(x,true,algo_state.iter);
127  con.update(x,true,algo_state.iter);
128  // Compute norm of the gradient of the Lagrangian
129  algo_state.value = myPen.value(x, zerotol);
130  myPen.gradient(*(state->gradientVec), x, zerotol);
131  con.applyAdjointJacobian(*g_,l,x,zerotol);
132  state->gradientVec->plus(*g_);
133  gLnorm_ = (state->gradientVec)->norm();
134  // Compute constraint violation
135  con.value(*(state->constraintVec),x, zerotol);
136  algo_state.cnorm = (state->constraintVec)->norm();
138  algo_state.gnorm = std::max(gLnorm_,compViolation_);
139  // Update state
140  algo_state.nfval++;
141  algo_state.ngrad++;
142  algo_state.ncval++;
143  }
144 
145  void updateState(const Vector<Real> &x,
146  Objective<Real> &obj,
148  AlgorithmState<Real> &algo_state) {
150  = dynamic_cast<MoreauYosidaPenalty<Real>&>(obj);
151  Real zerotol = std::sqrt(ROL_EPSILON<Real>());
152  ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
153  // Update objective and constraint.
154  myPen.update(x,true,algo_state.iter);
155  // Compute norm of the gradient of the Lagrangian
156  algo_state.value = myPen.value(x, zerotol);
157  myPen.gradient(*(state->gradientVec), x, zerotol);
158  gLnorm_ = (state->gradientVec)->norm();
159  // Compute constraint violation
160  algo_state.cnorm = static_cast<Real>(0);
162  algo_state.gnorm = std::max(gLnorm_,compViolation_);
163  // Update state
164  algo_state.nfval++;
165  algo_state.ngrad++;
166  }
167 
168 public:
169 
171  using Step<Real>::compute;
172  using Step<Real>::update;
173 
175 
176  MoreauYosidaPenaltyStep(ROL::ParameterList &parlist)
177  : Step<Real>(), algo_(ROL::nullPtr),
178  x_(ROL::nullPtr), g_(ROL::nullPtr), l_(ROL::nullPtr),
179  tau_(10), print_(false), parlist_(parlist), subproblemIter_(0),
180  hasEquality_(false) {
181  // Parse parameters
182  Real ten(10), oem6(1.e-6), oem8(1.e-8);
183  ROL::ParameterList& steplist = parlist.sublist("Step").sublist("Moreau-Yosida Penalty");
184  Step<Real>::getState()->searchSize = steplist.get("Initial Penalty Parameter",ten);
185  tau_ = steplist.get("Penalty Parameter Growth Factor",ten);
186  updatePenalty_ = steplist.get("Update Penalty",true);
187  print_ = steplist.sublist("Subproblem").get("Print History",false);
188  // Set parameters for step subproblem
189  Real gtol = steplist.sublist("Subproblem").get("Optimality Tolerance",oem8);
190  Real ctol = steplist.sublist("Subproblem").get("Feasibility Tolerance",oem8);
191  Real stol = oem6*std::min(gtol,ctol);
192  int maxit = steplist.sublist("Subproblem").get("Iteration Limit",1000);
193  parlist_.sublist("Status Test").set("Gradient Tolerance", gtol);
194  parlist_.sublist("Status Test").set("Constraint Tolerance", ctol);
195  parlist_.sublist("Status Test").set("Step Tolerance", stol);
196  parlist_.sublist("Status Test").set("Iteration Limit", maxit);
197  // Get step name from parameterlist
198  stepname_ = steplist.sublist("Subproblem").get("Step Type","Composite Step");
200  }
201 
206  AlgorithmState<Real> &algo_state ) {
207  hasEquality_ = true;
208  // Initialize step state
209  ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
210  state->descentVec = x.clone();
211  state->gradientVec = g.clone();
212  state->constraintVec = c.clone();
213  // Initialize additional storage
214  x_ = x.clone();
215  g_ = g.clone();
216  l_ = l.clone();
217  // Project x onto the feasible set
218  if ( bnd.isActivated() ) {
219  bnd.project(x);
220  }
221  // Initialize the algorithm state
222  algo_state.nfval = 0;
223  algo_state.ncval = 0;
224  algo_state.ngrad = 0;
225  updateState(x,l,obj,con,bnd,algo_state);
226  }
227 
232  AlgorithmState<Real> &algo_state ) {
233  // Initialize step state
234  ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
235  state->descentVec = x.clone();
236  state->gradientVec = g.clone();
237  // Initialize additional storage
238  x_ = x.clone();
239  g_ = g.clone();
240  // Project x onto the feasible set
241  if ( bnd.isActivated() ) {
242  bnd.project(x);
243  }
244  // Initialize the algorithm state
245  algo_state.nfval = 0;
246  algo_state.ncval = 0;
247  algo_state.ngrad = 0;
248  updateState(x,obj,bnd,algo_state);
249 
250  bnd_ = ROL::makePtr<BoundConstraint<Real>>();
251  bnd_->deactivate();
252  }
253 
256  void compute( Vector<Real> &s, const Vector<Real> &x, const Vector<Real> &l,
257  Objective<Real> &obj, Constraint<Real> &con,
258  BoundConstraint<Real> &bnd,
259  AlgorithmState<Real> &algo_state ) {
260  //MoreauYosidaPenalty<Real> &myPen
261  // = dynamic_cast<MoreauYosidaPenalty<Real>&>(obj);
262  Real one(1);
263  Ptr<Objective<Real>> penObj;
265  Ptr<Objective<Real>> raw_obj = makePtrFromRef(obj);
266  Ptr<Constraint<Real>> raw_con = makePtrFromRef(con);
267  Ptr<StepState<Real>> state = Step<Real>::getState();
268  penObj = makePtr<AugmentedLagrangian<Real>>(raw_obj,raw_con,l,one,x,*(state->constraintVec),parlist_);
269  step_ = makePtr<AugmentedLagrangianStep<Real>>(parlist_);
270  }
271  else if (stepType_ == STEP_FLETCHER) {
272  Ptr<Objective<Real>> raw_obj = makePtrFromRef(obj);
273  Ptr<Constraint<Real>> raw_con = makePtrFromRef(con);
274  Ptr<StepState<Real>> state = Step<Real>::getState();
275  penObj = makePtr<Fletcher<Real>>(raw_obj,raw_con,x,*(state->constraintVec),parlist_);
276  step_ = makePtr<FletcherStep<Real>>(parlist_);
277  }
278  else {
279  penObj = makePtrFromRef(obj);
280  stepname_ = "Composite Step";
282  step_ = makePtr<CompositeStep<Real>>(parlist_);
283  }
284  status_ = makePtr<ConstraintStatusTest<Real>>(parlist_);
285  algo_ = ROL::makePtr<Algorithm<Real>>(step_,status_,false);
286  x_->set(x); l_->set(l);
287  algo_->run(*x_,*l_,*penObj,con,print_);
288  s.set(*x_); s.axpy(-one,x);
289  subproblemIter_ = (algo_->getState())->iter;
290  }
291 
296  AlgorithmState<Real> &algo_state ) {
297  Real one(1);
299  = dynamic_cast<MoreauYosidaPenalty<Real>&>(obj);
300  if (stepType_ == STEP_BUNDLE) {
301  status_ = makePtr<BundleStatusTest<Real>>(parlist_);
302  step_ = makePtr<BundleStep<Real>>(parlist_);
303  }
304  else if (stepType_ == STEP_LINESEARCH) {
305  status_ = makePtr<StatusTest<Real>>(parlist_);
306  step_ = makePtr<LineSearchStep<Real>>(parlist_);
307  }
308  else {
309  status_ = makePtr<StatusTest<Real>>(parlist_);
310  step_ = makePtr<TrustRegionStep<Real>>(parlist_);
311  }
312  algo_ = ROL::makePtr<Algorithm<Real>>(step_,status_,false);
313  x_->set(x);
314  algo_->run(*x_,myPen,*bnd_,print_);
315  s.set(*x_); s.axpy(-one,x);
316  subproblemIter_ = (algo_->getState())->iter;
317  }
318 
319 
325  AlgorithmState<Real> &algo_state ) {
327  = dynamic_cast<MoreauYosidaPenalty<Real>&>(obj);
328  ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
329  state->SPiter = subproblemIter_;
330  state->descentVec->set(s);
331  // Update iterate and Lagrange multiplier
332  x.plus(s);
333  l.set(*l_);
334  // Update objective and constraint
335  algo_state.iter++;
336  con.update(x,true,algo_state.iter);
337  myPen.update(x,true,algo_state.iter);
338  // Update state
339  updateState(x,l,obj,con,bnd,algo_state);
340  // Update multipliers
341  if (updatePenalty_) {
342  state->searchSize *= tau_;
343  }
344  myPen.updateMultipliers(state->searchSize,x);
345  algo_state.nfval += myPen.getNumberFunctionEvaluations() + ((algo_->getState())->nfval);
346  algo_state.ngrad += myPen.getNumberGradientEvaluations() + ((algo_->getState())->ngrad);
347  algo_state.ncval += (algo_->getState())->ncval;
348  algo_state.snorm = s.norm();
349  algo_state.iterateVec->set(x);
350  algo_state.lagmultVec->set(l);
351  }
352 
355  void update( Vector<Real> &x, const Vector<Real> &s,
357  AlgorithmState<Real> &algo_state ) {
359  = dynamic_cast<MoreauYosidaPenalty<Real>&>(obj);
360  ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
361  state->descentVec->set(s);
362  // Update iterate and Lagrange multiplier
363  x.plus(s);
364  // Update objective and constraint
365  algo_state.iter++;
366  myPen.update(x,true,algo_state.iter);
367  // Update state
368  updateState(x,obj,bnd,algo_state);
369  // Update multipliers
370  if (updatePenalty_) {
371  state->searchSize *= tau_;
372  }
373  myPen.updateMultipliers(state->searchSize,x);
374  algo_state.nfval += myPen.getNumberFunctionEvaluations() + ((algo_->getState())->nfval);
375  algo_state.ngrad += myPen.getNumberGradientEvaluations() + ((algo_->getState())->ngrad);
376  algo_state.snorm = s.norm();
377  algo_state.iterateVec->set(x);
378  }
379 
382  std::string printHeader( void ) const {
383  std::stringstream hist;
384  hist << " ";
385  hist << std::setw(6) << std::left << "iter";
386  hist << std::setw(15) << std::left << "fval";
387  if (hasEquality_) {
388  hist << std::setw(15) << std::left << "cnorm";
389  }
390  hist << std::setw(15) << std::left << "gnorm";
391  hist << std::setw(15) << std::left << "ifeas";
392  hist << std::setw(15) << std::left << "snorm";
393  hist << std::setw(10) << std::left << "penalty";
394  hist << std::setw(8) << std::left << "#fval";
395  hist << std::setw(8) << std::left << "#grad";
396  if (hasEquality_) {
397  hist << std::setw(8) << std::left << "#cval";
398  }
399  hist << std::setw(8) << std::left << "subIter";
400  hist << "\n";
401  return hist.str();
402  }
403 
406  std::string printName( void ) const {
407  std::stringstream hist;
408  hist << "\n" << " Moreau-Yosida Penalty solver";
409  hist << "\n";
410  return hist.str();
411  }
412 
415  std::string print( AlgorithmState<Real> &algo_state, bool pHeader = false ) const {
416  std::stringstream hist;
417  hist << std::scientific << std::setprecision(6);
418  if ( algo_state.iter == 0 ) {
419  hist << printName();
420  }
421  if ( pHeader ) {
422  hist << printHeader();
423  }
424  if ( algo_state.iter == 0 ) {
425  hist << " ";
426  hist << std::setw(6) << std::left << algo_state.iter;
427  hist << std::setw(15) << std::left << algo_state.value;
428  if (hasEquality_) {
429  hist << std::setw(15) << std::left << algo_state.cnorm;
430  }
431  hist << std::setw(15) << std::left << gLnorm_;
432  hist << std::setw(15) << std::left << compViolation_;
433  hist << std::setw(15) << std::left << " ";
434  hist << std::scientific << std::setprecision(2);
435  hist << std::setw(10) << std::left << Step<Real>::getStepState()->searchSize;
436  hist << "\n";
437  }
438  else {
439  hist << " ";
440  hist << std::setw(6) << std::left << algo_state.iter;
441  hist << std::setw(15) << std::left << algo_state.value;
442  if (hasEquality_) {
443  hist << std::setw(15) << std::left << algo_state.cnorm;
444  }
445  hist << std::setw(15) << std::left << gLnorm_;
446  hist << std::setw(15) << std::left << compViolation_;
447  hist << std::setw(15) << std::left << algo_state.snorm;
448  hist << std::scientific << std::setprecision(2);
449  hist << std::setw(10) << std::left << Step<Real>::getStepState()->searchSize;
450  hist << std::scientific << std::setprecision(6);
451  hist << std::setw(8) << std::left << algo_state.nfval;
452  hist << std::setw(8) << std::left << algo_state.ngrad;
453  if (hasEquality_) {
454  hist << std::setw(8) << std::left << algo_state.ncval;
455  }
456  hist << std::setw(8) << std::left << subproblemIter_;
457  hist << "\n";
458  }
459  return hist.str();
460  }
461 
462 }; // class MoreauYosidaPenaltyStep
463 
464 } // namespace ROL
465 
466 #endif
Provides the interface to evaluate objective functions.
EStep StringToEStep(std::string s)
Definition: ROL_Types.hpp:357
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void plus(const Vector &x)=0
Compute , where .
virtual void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update constraint function.
void updateState(const Vector< Real > &x, const Vector< Real > &l, Objective< Real > &obj, Constraint< Real > &con, BoundConstraint< Real > &bnd, AlgorithmState< Real > &algo_state)
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:119
bool isActivated(void) const
Check if bounds are on.
std::string printHeader(void) const
Print iterate header.
Provides the interface to compute optimization steps.
Definition: ROL_Step.hpp:34
MoreauYosidaPenaltyStep(ROL::ParameterList &parlist)
Contains definitions of custom data types in ROL.
Real value(const Vector< Real > &x, Real &tol)
Compute value.
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 constraint.
void updateMultipliers(Real mu, const ROL::Vector< Real > &x)
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
virtual void value(Vector< Real > &c, const Vector< Real > &x, Real &tol)=0
Evaluate the constraint operator at .
ROL::Ptr< Algorithm< Real > > algo_
ROL::Ptr< BoundConstraint< Real > > bnd_
State for algorithm class. Will be used for restarts.
Definition: ROL_Types.hpp:109
ROL::Ptr< StatusTest< Real > > status_
void update(Vector< Real > &x, const Vector< Real > &s, Objective< Real > &obj, BoundConstraint< Real > &bnd, AlgorithmState< Real > &algo_state)
Update step, for bound constraints.
ROL::Ptr< StepState< Real > > getState(void)
Definition: ROL_Step.hpp:39
Provides the interface to evaluate the Moreau-Yosida penalty function.
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).
ROL::Ptr< Vector< Real > > iterateVec
Definition: ROL_Types.hpp:123
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
virtual void project(Vector< Real > &x)
Project optimization variables onto the bounds.
void compute(Vector< Real > &s, const Vector< Real > &x, Objective< Real > &obj, BoundConstraint< Real > &bnd, AlgorithmState< Real > &algo_state)
Compute step for bound constraints.
void initialize(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &bnd, AlgorithmState< Real > &algo_state)
Initialize step without equality constraint.
Provides the interface to apply upper and lower bound constraints.
virtual void applyAdjointJacobian(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the adjoint of the the constraint Jacobian at , , to vector .
std::string printName(void) const
Print step name.
ROL::Ptr< Vector< Real > > lagmultVec
Definition: ROL_Types.hpp:124
void updateState(const Vector< Real > &x, Objective< Real > &obj, BoundConstraint< Real > &bnd, AlgorithmState< Real > &algo_state)
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update Moreau-Yosida penalty function.
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:175
virtual Real norm() const =0
Returns where .
std::string print(AlgorithmState< Real > &algo_state, bool pHeader=false) const
Print iterate status.
EStep
Enumeration of step types.
Definition: ROL_Types.hpp:242
Real testComplementarity(const ROL::Vector< Real > &x)
Defines the general constraint operator interface.
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).