ROL
ROL_MoreauYosidaPenaltyStep.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_MOREAUYOSIDAPENALTYSTEP_H
45 #define ROL_MOREAUYOSIDAPENALTYSTEP_H
46 
48 #include "ROL_Types.hpp"
50 #include "ROL_CompositeStep.hpp"
51 #include "ROL_FletcherStep.hpp"
52 #include "ROL_Algorithm.hpp"
53 #include "ROL_ParameterList.hpp"
54 
117 namespace ROL {
118 
119 template <class Real>
120 class MoreauYosidaPenaltyStep : public Step<Real> {
121 private:
122  ROL::Ptr<Algorithm<Real> > algo_;
123  ROL::Ptr<Vector<Real> > x_;
124  ROL::Ptr<Vector<Real> > g_;
125  ROL::Ptr<Vector<Real> > l_;
126  ROL::Ptr<BoundConstraint<Real> > bnd_;
127 
129  Real gLnorm_;
130  Real tau_;
131  bool print_;
133 
134  ROL::ParameterList parlist_;
137 
139  std::string stepname_;
140 
141  void updateState(const Vector<Real> &x, const Vector<Real> &l,
142  Objective<Real> &obj,
144  AlgorithmState<Real> &algo_state) {
146  = dynamic_cast<MoreauYosidaPenalty<Real>&>(obj);
147  Real zerotol = std::sqrt(ROL_EPSILON<Real>());
148  ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
149  // Update objective and constraint.
150  myPen.update(x,true,algo_state.iter);
151  con.update(x,true,algo_state.iter);
152  // Compute norm of the gradient of the Lagrangian
153  algo_state.value = myPen.value(x, zerotol);
154  myPen.gradient(*(state->gradientVec), x, zerotol);
155  con.applyAdjointJacobian(*g_,l,x,zerotol);
156  state->gradientVec->plus(*g_);
157  gLnorm_ = (state->gradientVec)->norm();
158  // Compute constraint violation
159  con.value(*(state->constraintVec),x, zerotol);
160  algo_state.cnorm = (state->constraintVec)->norm();
162  algo_state.gnorm = std::max(gLnorm_,compViolation_);
163  // Update state
164  algo_state.nfval++;
165  algo_state.ngrad++;
166  algo_state.ncval++;
167  }
168 
169  void updateState(const Vector<Real> &x,
170  Objective<Real> &obj,
172  AlgorithmState<Real> &algo_state) {
174  = dynamic_cast<MoreauYosidaPenalty<Real>&>(obj);
175  Real zerotol = std::sqrt(ROL_EPSILON<Real>());
176  ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
177  // Update objective and constraint.
178  myPen.update(x,true,algo_state.iter);
179  // Compute norm of the gradient of the Lagrangian
180  algo_state.value = myPen.value(x, zerotol);
181  myPen.gradient(*(state->gradientVec), x, zerotol);
182  gLnorm_ = (state->gradientVec)->norm();
183  // Compute constraint violation
184  algo_state.cnorm = static_cast<Real>(0);
186  algo_state.gnorm = std::max(gLnorm_,compViolation_);
187  // Update state
188  algo_state.nfval++;
189  algo_state.ngrad++;
190  }
191 
192 public:
193 
195  using Step<Real>::compute;
196  using Step<Real>::update;
197 
199 
200  MoreauYosidaPenaltyStep(ROL::ParameterList &parlist)
201  : Step<Real>(), algo_(ROL::nullPtr),
202  x_(ROL::nullPtr), g_(ROL::nullPtr), l_(ROL::nullPtr),
203  tau_(10), print_(false), parlist_(parlist), subproblemIter_(0),
204  hasEquality_(false) {
205  // Parse parameters
206  Real ten(10), oem6(1.e-6), oem8(1.e-8);
207  ROL::ParameterList& steplist = parlist.sublist("Step").sublist("Moreau-Yosida Penalty");
208  Step<Real>::getState()->searchSize = steplist.get("Initial Penalty Parameter",ten);
209  tau_ = steplist.get("Penalty Parameter Growth Factor",ten);
210  updatePenalty_ = steplist.get("Update Penalty",true);
211  print_ = steplist.sublist("Subproblem").get("Print History",false);
212  // Set parameters for step subproblem
213  Real gtol = steplist.sublist("Subproblem").get("Optimality Tolerance",oem8);
214  Real ctol = steplist.sublist("Subproblem").get("Feasibility Tolerance",oem8);
215  Real stol = oem6*std::min(gtol,ctol);
216  int maxit = steplist.sublist("Subproblem").get("Iteration Limit",1000);
217  parlist_.sublist("Status Test").set("Gradient Tolerance", gtol);
218  parlist_.sublist("Status Test").set("Constraint Tolerance", ctol);
219  parlist_.sublist("Status Test").set("Step Tolerance", stol);
220  parlist_.sublist("Status Test").set("Iteration Limit", maxit);
221  // Get step name from parameterlist
222  stepname_ = steplist.sublist("Subproblem").get("Step Type","Composite Step");
224  }
225 
230  AlgorithmState<Real> &algo_state ) {
231  hasEquality_ = true;
232  // Initialize step state
233  ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
234  state->descentVec = x.clone();
235  state->gradientVec = g.clone();
236  state->constraintVec = c.clone();
237  // Initialize additional storage
238  x_ = x.clone();
239  g_ = g.clone();
240  l_ = l.clone();
241  // Project x onto the feasible set
242  if ( bnd.isActivated() ) {
243  bnd.project(x);
244  }
245  // Initialize the algorithm state
246  algo_state.nfval = 0;
247  algo_state.ncval = 0;
248  algo_state.ngrad = 0;
249  updateState(x,l,obj,con,bnd,algo_state);
250  }
251 
256  AlgorithmState<Real> &algo_state ) {
257  // Initialize step state
258  ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
259  state->descentVec = x.clone();
260  state->gradientVec = g.clone();
261  // Initialize additional storage
262  x_ = x.clone();
263  g_ = g.clone();
264  // Project x onto the feasible set
265  if ( bnd.isActivated() ) {
266  bnd.project(x);
267  }
268  // Initialize the algorithm state
269  algo_state.nfval = 0;
270  algo_state.ncval = 0;
271  algo_state.ngrad = 0;
272  updateState(x,obj,bnd,algo_state);
273 
274  bnd_ = ROL::makePtr<BoundConstraint<Real>>();
275  bnd_->deactivate();
276  }
277 
280  void compute( Vector<Real> &s, const Vector<Real> &x, const Vector<Real> &l,
281  Objective<Real> &obj, Constraint<Real> &con,
282  BoundConstraint<Real> &bnd,
283  AlgorithmState<Real> &algo_state ) {
284  //MoreauYosidaPenalty<Real> &myPen
285  // = dynamic_cast<MoreauYosidaPenalty<Real>&>(obj);
286  Real one(1);
287  Ptr<Objective<Real>> penObj;
289  Ptr<Objective<Real>> raw_obj = makePtrFromRef(obj);
290  Ptr<Constraint<Real>> raw_con = makePtrFromRef(con);
291  Ptr<StepState<Real>> state = Step<Real>::getState();
292  penObj = makePtr<AugmentedLagrangian<Real>>(raw_obj,raw_con,l,one,x,*(state->constraintVec),parlist_);
293  }
294  else if (stepType_ == STEP_FLETCHER) {
295  Ptr<Objective<Real>> raw_obj = makePtrFromRef(obj);
296  Ptr<Constraint<Real>> raw_con = makePtrFromRef(con);
297  Ptr<StepState<Real>> state = Step<Real>::getState();
298  penObj = makePtr<Fletcher<Real>>(raw_obj,raw_con,x,*(state->constraintVec),parlist_);
299  }
300  else {
301  penObj = makePtrFromRef(obj);
302  stepname_ = "Composite Step";
304  }
305  algo_ = ROL::makePtr<Algorithm<Real>>(stepname_,parlist_,false);
306  x_->set(x); l_->set(l);
307  algo_->run(*x_,*l_,*penObj,con,print_);
308  s.set(*x_); s.axpy(-one,x);
309  subproblemIter_ = (algo_->getState())->iter;
310  }
311 
316  AlgorithmState<Real> &algo_state ) {
317  Real one(1);
319  = dynamic_cast<MoreauYosidaPenalty<Real>&>(obj);
320  algo_ = ROL::makePtr<Algorithm<Real>>("Trust Region",parlist_,false);
321  x_->set(x);
322  algo_->run(*x_,myPen,*bnd_,print_);
323  s.set(*x_); s.axpy(-one,x);
324  subproblemIter_ = (algo_->getState())->iter;
325  }
326 
327 
333  AlgorithmState<Real> &algo_state ) {
335  = dynamic_cast<MoreauYosidaPenalty<Real>&>(obj);
336  ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
337  state->SPiter = subproblemIter_;
338  state->descentVec->set(s);
339  // Update iterate and Lagrange multiplier
340  x.plus(s);
341  l.set(*l_);
342  // Update objective and constraint
343  algo_state.iter++;
344  con.update(x,true,algo_state.iter);
345  myPen.update(x,true,algo_state.iter);
346  // Update state
347  updateState(x,l,obj,con,bnd,algo_state);
348  // Update multipliers
349  if (updatePenalty_) {
350  state->searchSize *= tau_;
351  }
352  myPen.updateMultipliers(state->searchSize,x);
353  algo_state.nfval += myPen.getNumberFunctionEvaluations() + ((algo_->getState())->nfval);
354  algo_state.ngrad += myPen.getNumberGradientEvaluations() + ((algo_->getState())->ngrad);
355  algo_state.ncval += (algo_->getState())->ncval;
356  algo_state.snorm = s.norm();
357  algo_state.iterateVec->set(x);
358  algo_state.lagmultVec->set(l);
359  }
360 
363  void update( Vector<Real> &x, const Vector<Real> &s,
365  AlgorithmState<Real> &algo_state ) {
367  = dynamic_cast<MoreauYosidaPenalty<Real>&>(obj);
368  ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
369  state->descentVec->set(s);
370  // Update iterate and Lagrange multiplier
371  x.plus(s);
372  // Update objective and constraint
373  algo_state.iter++;
374  myPen.update(x,true,algo_state.iter);
375  // Update state
376  updateState(x,obj,bnd,algo_state);
377  // Update multipliers
378  if (updatePenalty_) {
379  state->searchSize *= tau_;
380  }
381  myPen.updateMultipliers(state->searchSize,x);
382  algo_state.nfval += myPen.getNumberFunctionEvaluations() + ((algo_->getState())->nfval);
383  algo_state.ngrad += myPen.getNumberGradientEvaluations() + ((algo_->getState())->ngrad);
384  algo_state.snorm = s.norm();
385  algo_state.iterateVec->set(x);
386  }
387 
390  std::string printHeader( void ) const {
391  std::stringstream hist;
392  hist << " ";
393  hist << std::setw(6) << std::left << "iter";
394  hist << std::setw(15) << std::left << "fval";
395  if (hasEquality_) {
396  hist << std::setw(15) << std::left << "cnorm";
397  }
398  hist << std::setw(15) << std::left << "gnorm";
399  hist << std::setw(15) << std::left << "ifeas";
400  hist << std::setw(15) << std::left << "snorm";
401  hist << std::setw(10) << std::left << "penalty";
402  hist << std::setw(8) << std::left << "#fval";
403  hist << std::setw(8) << std::left << "#grad";
404  if (hasEquality_) {
405  hist << std::setw(8) << std::left << "#cval";
406  }
407  hist << std::setw(8) << std::left << "subIter";
408  hist << "\n";
409  return hist.str();
410  }
411 
414  std::string printName( void ) const {
415  std::stringstream hist;
416  hist << "\n" << " Moreau-Yosida Penalty solver";
417  hist << "\n";
418  return hist.str();
419  }
420 
423  std::string print( AlgorithmState<Real> &algo_state, bool pHeader = false ) const {
424  std::stringstream hist;
425  hist << std::scientific << std::setprecision(6);
426  if ( algo_state.iter == 0 ) {
427  hist << printName();
428  }
429  if ( pHeader ) {
430  hist << printHeader();
431  }
432  if ( algo_state.iter == 0 ) {
433  hist << " ";
434  hist << std::setw(6) << std::left << algo_state.iter;
435  hist << std::setw(15) << std::left << algo_state.value;
436  if (hasEquality_) {
437  hist << std::setw(15) << std::left << algo_state.cnorm;
438  }
439  hist << std::setw(15) << std::left << gLnorm_;
440  hist << std::setw(15) << std::left << compViolation_;
441  hist << std::setw(15) << std::left << " ";
442  hist << std::scientific << std::setprecision(2);
443  hist << std::setw(10) << std::left << Step<Real>::getStepState()->searchSize;
444  hist << "\n";
445  }
446  else {
447  hist << " ";
448  hist << std::setw(6) << std::left << algo_state.iter;
449  hist << std::setw(15) << std::left << algo_state.value;
450  if (hasEquality_) {
451  hist << std::setw(15) << std::left << algo_state.cnorm;
452  }
453  hist << std::setw(15) << std::left << gLnorm_;
454  hist << std::setw(15) << std::left << compViolation_;
455  hist << std::setw(15) << std::left << algo_state.snorm;
456  hist << std::scientific << std::setprecision(2);
457  hist << std::setw(10) << std::left << Step<Real>::getStepState()->searchSize;
458  hist << std::scientific << std::setprecision(6);
459  hist << std::setw(8) << std::left << algo_state.nfval;
460  hist << std::setw(8) << std::left << algo_state.ngrad;
461  if (hasEquality_) {
462  hist << std::setw(8) << std::left << algo_state.ncval;
463  }
464  hist << std::setw(8) << std::left << subproblemIter_;
465  hist << "\n";
466  }
467  return hist.str();
468  }
469 
470 }; // class MoreauYosidaPenaltyStep
471 
472 } // namespace ROL
473 
474 #endif
Provides the interface to evaluate objective functions.
virtual void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update constraint functions. x is the optimization variable, flag = true if optimization variable is ...
EStep StringToEStep(std::string s)
Definition: ROL_Types.hpp:389
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void plus(const Vector &x)=0
Compute , where .
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:153
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:69
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.
Implements the computation of optimization steps using Moreau-Yosida regularized bound constraints...
void updateMultipliers(Real mu, const ROL::Vector< Real > &x)
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:80
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:143
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:74
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:157
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
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:158
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:209
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:274
Real testComplementarity(const ROL::Vector< Real > &x)
Defines the general constraint operator interface.
virtual void project(Vector< Real > &x)
Project optimization variables onto the bounds.
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).