ROL
ROL_MoreauYosidaObjective.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_MOREAUYOSIDAOBJECTIVE_H
11 #define ROL_MOREAUYOSIDAOBJECTIVE_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 "ROL_ScalarController.hpp"
19 #include "ROL_ParameterList.hpp"
20 #include <iostream>
21 
30 namespace ROL {
31 
32 template <class Real>
33 class MoreauYosidaObjective : public Objective<Real> {
34 private:
35  const Ptr<Objective<Real>> obj_;
36  const Ptr<BoundConstraint<Real>> bnd_;
37 
38  Ptr<Vector<Real>> l_;
39  Ptr<Vector<Real>> u_;
40  Ptr<Vector<Real>> l1_;
41  Ptr<Vector<Real>> u1_;
42  Ptr<Vector<Real>> dl1_;
43  Ptr<Vector<Real>> du1_;
44  Ptr<Vector<Real>> xlam_;
45  Ptr<Vector<Real>> v_;
46  Ptr<Vector<Real>> dv_;
47  Ptr<Vector<Real>> dv2_;
48  Ptr<Vector<Real>> lam_;
49  Ptr<Vector<Real>> tmp_;
50 
51  Ptr<ScalarController<Real,int>> fval_;
52  Ptr<VectorController<Real,int>> gradient_;
53 
54  Real mu_;
56  int nfval_;
57  int ngrad_;
60 
61  void computePenalty(const Vector<Real> &x) {
62  if ( bnd_->isActivated() ) {
63  Real one = 1.0;
64  if ( !isPenEvaluated_ ) {
65  xlam_->set(x);
66  xlam_->axpy(one/mu_,*lam_);
67 
68  if ( bnd_->isFeasible(*xlam_) ) {
69  l1_->zero(); dl1_->zero();
70  u1_->zero(); du1_->zero();
71  }
72  else {
73  // Compute lower penalty component
74  l1_->set(*l_);
75  bnd_->pruneLowerInactive(*l1_,*xlam_);
76  tmp_->set(*xlam_);
77  bnd_->pruneLowerInactive(*tmp_,*xlam_);
78  l1_->axpy(-one,*tmp_);
79 
80  // Compute upper penalty component
81  u1_->set(*xlam_);
82  bnd_->pruneUpperInactive(*u1_,*xlam_);
83  tmp_->set(*u_);
84  bnd_->pruneUpperInactive(*tmp_,*xlam_);
85  u1_->axpy(-one,*tmp_);
86 
87  // Compute derivative of lower penalty component
88  dl1_->set(l1_->dual());
89  bnd_->pruneLowerInactive(*dl1_,*xlam_);
90 
91  // Compute derivative of upper penalty component
92  du1_->set(u1_->dual());
93  bnd_->pruneUpperInactive(*du1_,*xlam_);
94  }
95 
96  isPenEvaluated_ = true;
97  }
98  }
99  }
100 
101  void initialize(const Vector<Real> &x,
102  const Vector<Real> &g) {
103  fval_ = makePtr<ScalarController<Real,int>>();
104  gradient_ = makePtr<VectorController<Real,int>>();
105 
106  l_ = x.clone();
107  l1_ = x.clone();
108  dl1_ = g.clone();
109  u_ = x.clone();
110  u1_ = x.clone();
111  du1_ = g.clone();
112  xlam_ = x.clone();
113  v_ = x.clone();
114  dv_ = g.clone();
115  dv2_ = g.clone();
116  lam_ = x.clone();
117  tmp_ = x.clone();
118 
119  l_->set(*bnd_->getLowerBound());
120  u_->set(*bnd_->getUpperBound());
121 
122  lam_->zero();
123  //lam_->set(*u_);
124  //lam_->plus(*l_);
125  //lam_->scale(0.5);
126  }
127 
128 public:
130  const Ptr<BoundConstraint<Real>> &bnd,
131  const Vector<Real> &x,
132  const Vector<Real> &g,
133  const Real mu = 1e1,
134  const bool updateMultiplier = true,
135  const bool updatePenalty = true)
136  : obj_(obj), bnd_(bnd), mu_(mu),
137  isPenEvaluated_(false), nfval_(0), ngrad_(0),
139  initialize(x,g);
140  }
141 
143  const Ptr<BoundConstraint<Real>> &bnd,
144  const Vector<Real> &x,
145  const Vector<Real> &g,
146  const Vector<Real> &lam,
147  const Real mu = 1e1,
148  const bool updateMultiplier = true,
149  const bool updatePenalty = true)
151  lam_->set(lam);
152  }
153 
155  const Ptr<BoundConstraint<Real>> &bnd,
156  const Vector<Real> &x,
157  const Vector<Real> &g,
158  ParameterList &parlist)
159  : obj_(obj), bnd_(bnd),
160  isPenEvaluated_(false), nfval_(0), ngrad_(0) {
161  initialize(x,g);
162  ParameterList &list = parlist.sublist("Step").sublist("Moreau-Yosida Penalty");
163  updateMultiplier_ = list.get("Update Multiplier",true);
164  updatePenalty_ = list.get("Update Penalty",true);
165  mu_ = list.get("Initial Penalty Parameter",1e1);
166  }
167 
169  const Ptr<BoundConstraint<Real>> &bnd,
170  const Vector<Real> &x,
171  const Vector<Real> &g,
172  const Vector<Real> &lam,
173  ParameterList &parlist)
174  : MoreauYosidaObjective(obj,bnd,x,g,parlist) {
175  lam_->set(lam);
176  }
177 
179  computePenalty(x);
180  lam_->set(*u1_);
181  lam_->axpy(static_cast<Real>(-1),*l1_);
182  lam_->scale(mu_);
183  isPenEvaluated_ = false;
184  }
185  void updatePenalty(Real mu) {
186  mu_ = mu;
187  isPenEvaluated_ = false;
188  }
189 
190  void updateMultipliers(Real mu, const Vector<Real> &x) {
191  if ( bnd_->isActivated() ) {
192  if ( updateMultiplier_ ) {
193  updateMultiplier(x);
194  //const Real one(1);
195  //computePenalty(x);
196  //lam_->set(*u1_);
197  //lam_->axpy(-one,*l1_);
198  //lam_->scale(mu_);
199  }
200  if ( updatePenalty_ ) {
201  updatePenalty(mu);
202  //mu_ = mu;
203  }
204  }
205  nfval_ = 0; ngrad_ = 0;
206  //isPenEvaluated_ = false;
207  }
208 
209  void reset(const Real mu) {
210  lam_->zero();
211  mu_ = mu;
212  nfval_ = 0; ngrad_ = 0;
213  isPenEvaluated_ = false;
214  }
215 
217  Real val(0);
218  if (bnd_->isActivated()) {
219  computePenalty(x);
220 
221  tmp_->set(x);
222  tmp_->axpy(static_cast<Real>(-1), *l_);
223  Real lower = mu_*std::abs(tmp_->dot(*l1_));
224 
225  tmp_->set(x);
226  tmp_->axpy(static_cast<Real>(-1), *u_);
227  Real upper = mu_*std::abs(tmp_->dot(*u1_));
228 
229  tmp_->set(x);
230  bnd_->project(*tmp_);
231  tmp_->axpy(static_cast<Real>(-1), x);
232  Real xnorm = tmp_->norm();
233 
234  val = std::max(xnorm,std::max(lower,upper));
235  }
236  return val;
237  }
238 
239  Real getObjectiveValue(const Vector<Real> &x, Real &tol) {
240  int key(0);
241  Real val(0);
242  bool isComputed = fval_->get(val,key);
243  if (!isComputed) {
244  val = obj_->value(x,tol); nfval_++;
245  fval_->set(val,key);
246  }
247  return val;
248  }
249 
250  void getObjectiveGradient(Vector<Real> &g, const Vector<Real> &x, Real &tol) {
251  int key(0);
252  bool isComputed = gradient_->get(g,key);
253  if (!isComputed) {
254  obj_->gradient(g,x,tol); ngrad_++;
255  gradient_->set(g,key);
256  }
257  }
258 
260  return nfval_;
261  }
262 
264  return ngrad_;
265  }
266 
274  void update( const Vector<Real> &x, UpdateType type, int iter = -1 ) {
275  obj_->update(x,type,iter);
276  fval_->objectiveUpdate(type);
277  gradient_->objectiveUpdate(type);
278  // Need to do something smart here
279  isPenEvaluated_ = false;
280  }
281 
288  Real value( const Vector<Real> &x, Real &tol ) {
289  // Compute objective function value
290  Real fval = getObjectiveValue(x,tol);
291  // Add value of the Moreau-Yosida penalty
292  if ( bnd_->isActivated() ) {
293  const Real half(0.5);
294  computePenalty(x);
295  fval += half*mu_*(l1_->dot(*l1_) + u1_->dot(*u1_));
296  }
297  return fval;
298  }
299 
307  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
308  // Compute gradient of objective function
309  getObjectiveGradient(g,x,tol);
310  // Add gradient of the Moreau-Yosida penalty
311  if ( bnd_->isActivated() ) {
312  computePenalty(x);
313  g.axpy(-mu_,*dl1_);
314  g.axpy(mu_,*du1_);
315  }
316  }
317 
326  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
327  // Apply objective Hessian to a vector
328  obj_->hessVec(hv,v,x,tol);
329  // Add Hessian of the Moreau-Yosida penalty
330  if ( bnd_->isActivated() ) {
331  const Real one(1);
332  computePenalty(x);
333 
334  v_->set(v);
335  bnd_->pruneLowerActive(*v_,*xlam_);
336  v_->scale(-one);
337  v_->plus(v);
338  dv_->set(v_->dual());
339  dv2_->set(*dv_);
340  bnd_->pruneLowerActive(*dv_,*xlam_);
341  dv_->scale(-one);
342  dv_->plus(*dv2_);
343  hv.axpy(mu_,*dv_);
344 
345  v_->set(v);
346  bnd_->pruneUpperActive(*v_,*xlam_);
347  v_->scale(-one);
348  v_->plus(v);
349  dv_->set(v_->dual());
350  dv2_->set(*dv_);
351  bnd_->pruneUpperActive(*dv_,*xlam_);
352  dv_->scale(-one);
353  dv_->plus(*dv2_);
354  hv.axpy(mu_,*dv_);
355  }
356  }
357 }; // class MoreauYosidaObjective
358 
359 } // namespace ROL
360 
361 #endif
Provides the interface to evaluate objective functions.
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:119
Real getObjectiveValue(const Vector< Real > &x, Real &tol)
Contains definitions of custom data types in ROL.
Real value(const Vector< Real > &x, Real &tol)
Compute value.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
MoreauYosidaObjective(const Ptr< Objective< Real >> &obj, const Ptr< BoundConstraint< Real >> &bnd, const Vector< Real > &x, const Vector< Real > &g, const Vector< Real > &lam, ParameterList &parlist)
void updateMultiplier(const Vector< Real > &x)
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Ptr< ScalarController< Real, int > > fval_
void updateMultipliers(Real mu, const Vector< Real > &x)
void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update Moreau-Yosida penalty function.
Provides the interface to evaluate the Moreau-Yosida penalty function.
void computePenalty(const Vector< Real > &x)
Provides the interface to apply upper and lower bound constraints.
const Ptr< BoundConstraint< Real > > bnd_
MoreauYosidaObjective(const Ptr< Objective< Real >> &obj, const Ptr< BoundConstraint< Real >> &bnd, const Vector< Real > &x, const Vector< Real > &g, const Real mu=1e1, const bool updateMultiplier=true, const bool updatePenalty=true)
MoreauYosidaObjective(const Ptr< Objective< Real >> &obj, const Ptr< BoundConstraint< Real >> &bnd, const Vector< Real > &x, const Vector< Real > &g, ParameterList &parlist)
void initialize(const Vector< Real > &x, const Vector< Real > &g)
void getObjectiveGradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Real testComplementarity(const Vector< Real > &x)
const Ptr< Objective< Real > > obj_
Ptr< VectorController< Real, int > > gradient_
MoreauYosidaObjective(const Ptr< Objective< Real >> &obj, const Ptr< BoundConstraint< Real >> &bnd, const Vector< Real > &x, const Vector< Real > &g, const Vector< Real > &lam, const Real mu=1e1, const bool updateMultiplier=true, const bool updatePenalty=true)