ROL
ROL_MoreauYosidaPenalty.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_MOREAUYOSIDAPENALTY_H
11 #define ROL_MOREAUYOSIDAPENALTY_H
12 
13 #include "ROL_Objective.hpp"
14 #include "ROL_BoundConstraint.hpp"
15 #include "ROL_Vector.hpp"
16 #include "ROL_Types.hpp"
17 #include "ROL_Ptr.hpp"
18 #include <iostream>
19 
28 namespace ROL {
29 
30 template <class Real>
31 class MoreauYosidaPenalty : public Objective<Real> {
32 private:
33  const ROL::Ptr<Objective<Real> > obj_;
34  const ROL::Ptr<BoundConstraint<Real> > bnd_;
35 
36  ROL::Ptr<Vector<Real> > g_;
37  ROL::Ptr<Vector<Real> > l_;
38  ROL::Ptr<Vector<Real> > u_;
39  ROL::Ptr<Vector<Real> > l1_;
40  ROL::Ptr<Vector<Real> > u1_;
41  ROL::Ptr<Vector<Real> > dl1_;
42  ROL::Ptr<Vector<Real> > du1_;
43  ROL::Ptr<Vector<Real> > xlam_;
44  ROL::Ptr<Vector<Real> > v_;
45  ROL::Ptr<Vector<Real> > dv_;
46  ROL::Ptr<Vector<Real> > dv2_;
47  ROL::Ptr<Vector<Real> > lam_;
48  ROL::Ptr<Vector<Real> > tmp_;
49 
50  Real mu_;
51  Real fval_;
53  int nfval_;
54  int ngval_;
57 
58  void computePenalty(const Vector<Real> &x) {
59  if ( bnd_->isActivated() ) {
60  Real one = 1.0;
61  if ( !isConEvaluated_ ) {
62  xlam_->set(x);
63  xlam_->axpy(one/mu_,*lam_);
64 
65  if ( bnd_->isFeasible(*xlam_) ) {
66  l1_->zero(); dl1_->zero();
67  u1_->zero(); du1_->zero();
68  }
69  else {
70  // Compute lower penalty component
71  l1_->set(*l_);
72  bnd_->pruneLowerInactive(*l1_,*xlam_);
73  tmp_->set(*xlam_);
74  bnd_->pruneLowerInactive(*tmp_,*xlam_);
75  l1_->axpy(-one,*tmp_);
76 
77  // Compute upper penalty component
78  u1_->set(*xlam_);
79  bnd_->pruneUpperInactive(*u1_,*xlam_);
80  tmp_->set(*u_);
81  bnd_->pruneUpperInactive(*tmp_,*xlam_);
82  u1_->axpy(-one,*tmp_);
83 
84  // Compute derivative of lower penalty component
85  dl1_->set(l1_->dual());
86  bnd_->pruneLowerInactive(*dl1_,*xlam_);
87 
88  // Compute derivative of upper penalty component
89  du1_->set(u1_->dual());
90  bnd_->pruneUpperInactive(*du1_,*xlam_);
91  }
92 
93  isConEvaluated_ = true;
94  }
95  }
96  }
97 
99  const ROL::Ptr<ROL::BoundConstraint<Real> > &bnd) {
100  g_ = x.dual().clone();
101  l_ = x.clone();
102  l1_ = x.clone();
103  dl1_ = x.dual().clone();
104  u_ = x.clone();
105  u1_ = x.clone();
106  du1_ = x.dual().clone();
107  xlam_ = x.clone();
108  v_ = x.clone();
109  dv_ = x.dual().clone();
110  dv2_ = x.dual().clone();
111  lam_ = x.clone();
112  tmp_ = x.clone();
113 
114  l_->set(*bnd_->getLowerBound());
115  u_->set(*bnd_->getUpperBound());
116 
117  lam_->zero();
118  //lam_->set(*u_);
119  //lam_->plus(*l_);
120  //lam_->scale(0.5);
121  }
122 
123 public:
125 
126  MoreauYosidaPenalty(const ROL::Ptr<Objective<Real> > &obj,
127  const ROL::Ptr<BoundConstraint<Real> > &bnd,
128  const ROL::Vector<Real> &x,
129  const Real mu = 1e1,
130  const bool updateMultiplier = true,
131  const bool updatePenalty = true)
132  : obj_(obj), bnd_(bnd), mu_(mu),
133  fval_(0), isConEvaluated_(false), nfval_(0), ngval_(0),
134  updateMultiplier_(updateMultiplier), updatePenalty_(updatePenalty) {
135  initialize(x,bnd);
136  }
137 
138  MoreauYosidaPenalty(const ROL::Ptr<Objective<Real> > &obj,
139  const ROL::Ptr<BoundConstraint<Real> > &bnd,
140  const ROL::Vector<Real> &x,
141  ROL::ParameterList &parlist)
142  : obj_(obj), bnd_(bnd),
143  fval_(0), isConEvaluated_(false), nfval_(0), ngval_(0) {
144  initialize(x,bnd);
145  ROL::ParameterList &list = parlist.sublist("Step").sublist("Moreau-Yosida Penalty");
146  updateMultiplier_ = list.get("Update Multiplier",true);
147  updatePenalty_ = list.get("Update Penalty",true);
148  mu_ = list.get("Initial Penalty Parameter",1e1);
149  }
150 
151  MoreauYosidaPenalty(const ROL::Ptr<Objective<Real> > &obj,
152  const ROL::Ptr<BoundConstraint<Real> > &bnd,
153  const ROL::Vector<Real> &x,
154  const ROL::Vector<Real> &lam,
155  ROL::ParameterList &parlist)
156  : obj_(obj), bnd_(bnd),
157  fval_(0), isConEvaluated_(false), nfval_(0), ngval_(0) {
158  initialize(x,bnd);
159  lam_->set(lam);
160  ROL::ParameterList &list = parlist.sublist("Step").sublist("Moreau-Yosida Penalty");
161  updateMultiplier_ = list.get("Update Multiplier",true);
162  updatePenalty_ = list.get("Update Penalty",true);
163  mu_ = list.get("Initial Penalty Parameter",1e1);
164  }
165 
166  void updateMultipliers(Real mu, const ROL::Vector<Real> &x) {
167  if ( bnd_->isActivated() ) {
168  if ( updateMultiplier_ ) {
169  const Real one(1);
170  computePenalty(x);
171  lam_->set(*u1_);
172  lam_->axpy(-one,*l1_);
173  lam_->scale(mu_);
174  }
175  if ( updatePenalty_ ) {
176  mu_ = mu;
177  }
178  }
179  nfval_ = 0; ngval_ = 0;
180  isConEvaluated_ = false;
181  }
182 
183  void reset(const Real mu) {
184  lam_->zero();
185  mu_ = mu;
186  nfval_ = 0; ngval_ = 0;
187  isConEvaluated_ = false;
188  }
189 
191  Real val(0);
192  if (bnd_->isActivated()) {
193  computePenalty(x);
194 
195  tmp_->set(x);
196  tmp_->axpy(static_cast<Real>(-1), *l_);
197  Real lower = mu_*std::abs(tmp_->dot(*l1_));
198 
199  tmp_->set(x);
200  tmp_->axpy(static_cast<Real>(-1), *u_);
201  Real upper = mu_*std::abs(tmp_->dot(*u1_));
202 
203  tmp_->set(x);
204  bnd_->project(*tmp_);
205  tmp_->axpy(static_cast<Real>(-1), x);
206  Real xnorm = tmp_->norm();
207 
208  val = std::max(xnorm,std::max(lower,upper));
209  }
210  return val;
211  }
212 
213  Real getObjectiveValue(void) const {
214  return fval_;
215  }
216 
217  ROL::Ptr<Vector<Real> > getGradient(void) const {
218  return g_;
219  }
220 
222  return nfval_;
223  }
224 
226  return ngval_;
227  }
228 
236  void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
237  obj_->update(x,flag,iter);
238  isConEvaluated_ = false;
239  }
240 
247  Real value( const Vector<Real> &x, Real &tol ) {
248  Real half = 0.5;
249  // Compute objective function value
250  fval_ = obj_->value(x,tol);
251  nfval_++;
252  // Add value of the Moreau-Yosida penalty
253  Real fval = fval_;
254  if ( bnd_->isActivated() ) {
255  computePenalty(x);
256  fval += half*mu_*(l1_->dot(*l1_) + u1_->dot(*u1_));
257  }
258  return fval;
259  }
260 
268  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
269  // Compute gradient of objective function
270  obj_->gradient(*g_,x,tol);
271  ngval_++;
272  g.set(*g_);
273  // Add gradient of the Moreau-Yosida penalty
274  if ( bnd_->isActivated() ) {
275  computePenalty(x);
276  g.axpy(-mu_,*dl1_);
277  g.axpy(mu_,*du1_);
278  }
279  }
280 
289  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
290  // Apply objective Hessian to a vector
291  obj_->hessVec(hv,v,x,tol);
292  // Add Hessian of the Moreau-Yosida penalty
293  if ( bnd_->isActivated() ) {
294  Real one = 1.0;
295  computePenalty(x);
296 
297  v_->set(v);
298  bnd_->pruneLowerActive(*v_,*xlam_);
299  v_->scale(-one);
300  v_->plus(v);
301  dv_->set(v_->dual());
302  dv2_->set(*dv_);
303  bnd_->pruneLowerActive(*dv_,*xlam_);
304  dv_->scale(-one);
305  dv_->plus(*dv2_);
306  hv.axpy(mu_,*dv_);
307 
308  v_->set(v);
309  bnd_->pruneUpperActive(*v_,*xlam_);
310  v_->scale(-one);
311  v_->plus(v);
312  dv_->set(v_->dual());
313  dv2_->set(*dv_);
314  bnd_->pruneUpperActive(*dv_,*xlam_);
315  dv_->scale(-one);
316  dv_->plus(*dv2_);
317  hv.axpy(mu_,*dv_);
318  }
319  }
320 
321 // Definitions for parametrized (stochastic) objective functions
322 public:
323  void setParameter(const std::vector<Real> &param) {
325  obj_->setParameter(param);
326  }
327 }; // class MoreauYosidaPenalty
328 
329 } // namespace ROL
330 
331 #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
ROL::Ptr< Vector< Real > > g_
ROL::Ptr< Vector< Real > > l1_
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
ROL::Ptr< Vector< Real > > du1_
void setParameter(const std::vector< Real > &param)
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:119
Contains definitions of custom data types in ROL.
Real value(const Vector< Real > &x, Real &tol)
Compute value.
ROL::Ptr< Vector< Real > > xlam_
void updateMultipliers(Real mu, const ROL::Vector< Real > &x)
ROL::Ptr< Vector< Real > > dv_
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
void computePenalty(const Vector< Real > &x)
ROL::Ptr< Vector< Real > > dl1_
ROL::Ptr< Vector< Real > > tmp_
const ROL::Ptr< Objective< Real > > obj_
Provides the interface to evaluate the Moreau-Yosida penalty function.
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
MoreauYosidaPenalty(const ROL::Ptr< Objective< Real > > &obj, const ROL::Ptr< BoundConstraint< Real > > &bnd, const ROL::Vector< Real > &x, const ROL::Vector< Real > &lam, ROL::ParameterList &parlist)
Provides the interface to apply upper and lower bound constraints.
ROL::Ptr< Vector< Real > > l_
ROL::Ptr< Vector< Real > > u_
virtual void setParameter(const std::vector< Real > &param)
ROL::Ptr< Vector< Real > > v_
MoreauYosidaPenalty(const ROL::Ptr< Objective< Real > > &obj, const ROL::Ptr< BoundConstraint< Real > > &bnd, const ROL::Vector< Real > &x, const Real mu=1e1, const bool updateMultiplier=true, const bool updatePenalty=true)
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
const ROL::Ptr< BoundConstraint< Real > > bnd_
ROL::Ptr< Vector< Real > > lam_
ROL::Ptr< Vector< Real > > dv2_
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
void initialize(const ROL::Vector< Real > &x, const ROL::Ptr< ROL::BoundConstraint< Real > > &bnd)
ROL::Ptr< Vector< Real > > u1_
Real testComplementarity(const ROL::Vector< Real > &x)
MoreauYosidaPenalty(const ROL::Ptr< Objective< Real > > &obj, const ROL::Ptr< BoundConstraint< Real > > &bnd, const ROL::Vector< Real > &x, ROL::ParameterList &parlist)
ROL::Ptr< Vector< Real > > getGradient(void) const