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),
138  updateMultiplier_(updateMultiplier), updatePenalty_(updatePenalty) {
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)
150  : MoreauYosidaObjective(obj,bnd,x,g,mu,updateMultiplier,updatePenalty) {
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 
178  void updateMultipliers(Real mu, const Vector<Real> &x) {
179  if ( bnd_->isActivated() ) {
180  if ( updateMultiplier_ ) {
181  const Real one(1);
182  computePenalty(x);
183  lam_->set(*u1_);
184  lam_->axpy(-one,*l1_);
185  lam_->scale(mu_);
186  }
187  if ( updatePenalty_ ) {
188  mu_ = mu;
189  }
190  }
191  nfval_ = 0; ngrad_ = 0;
192  isPenEvaluated_ = false;
193  }
194 
195  void reset(const Real mu) {
196  lam_->zero();
197  mu_ = mu;
198  nfval_ = 0; ngrad_ = 0;
199  isPenEvaluated_ = false;
200  }
201 
203  Real val(0);
204  if (bnd_->isActivated()) {
205  computePenalty(x);
206 
207  tmp_->set(x);
208  tmp_->axpy(static_cast<Real>(-1), *l_);
209  Real lower = mu_*std::abs(tmp_->dot(*l1_));
210 
211  tmp_->set(x);
212  tmp_->axpy(static_cast<Real>(-1), *u_);
213  Real upper = mu_*std::abs(tmp_->dot(*u1_));
214 
215  tmp_->set(x);
216  bnd_->project(*tmp_);
217  tmp_->axpy(static_cast<Real>(-1), x);
218  Real xnorm = tmp_->norm();
219 
220  val = std::max(xnorm,std::max(lower,upper));
221  }
222  return val;
223  }
224 
225  Real getObjectiveValue(const Vector<Real> &x, Real &tol) {
226  int key(0);
227  Real val(0);
228  bool isComputed = fval_->get(val,key);
229  if (!isComputed) {
230  val = obj_->value(x,tol); nfval_++;
231  fval_->set(val,key);
232  }
233  return val;
234  }
235 
236  void getObjectiveGradient(Vector<Real> &g, const Vector<Real> &x, Real &tol) {
237  int key(0);
238  bool isComputed = gradient_->get(g,key);
239  if (!isComputed) {
240  obj_->gradient(g,x,tol); ngrad_++;
241  gradient_->set(g,key);
242  }
243  }
244 
246  return nfval_;
247  }
248 
250  return ngrad_;
251  }
252 
260  void update( const Vector<Real> &x, UpdateType type, int iter = -1 ) {
261  obj_->update(x,type,iter);
262  fval_->objectiveUpdate(type);
263  gradient_->objectiveUpdate(type);
264  // Need to do something smart here
265  isPenEvaluated_ = false;
266  }
267 
274  Real value( const Vector<Real> &x, Real &tol ) {
275  // Compute objective function value
276  Real fval = getObjectiveValue(x,tol);
277  // Add value of the Moreau-Yosida penalty
278  if ( bnd_->isActivated() ) {
279  const Real half(0.5);
280  computePenalty(x);
281  fval += half*mu_*(l1_->dot(*l1_) + u1_->dot(*u1_));
282  }
283  return fval;
284  }
285 
293  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
294  // Compute gradient of objective function
295  getObjectiveGradient(g,x,tol);
296  // Add gradient of the Moreau-Yosida penalty
297  if ( bnd_->isActivated() ) {
298  computePenalty(x);
299  g.axpy(-mu_,*dl1_);
300  g.axpy(mu_,*du1_);
301  }
302  }
303 
312  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
313  // Apply objective Hessian to a vector
314  obj_->hessVec(hv,v,x,tol);
315  // Add Hessian of the Moreau-Yosida penalty
316  if ( bnd_->isActivated() ) {
317  const Real one(1);
318  computePenalty(x);
319 
320  v_->set(v);
321  bnd_->pruneLowerActive(*v_,*xlam_);
322  v_->scale(-one);
323  v_->plus(v);
324  dv_->set(v_->dual());
325  dv2_->set(*dv_);
326  bnd_->pruneLowerActive(*dv_,*xlam_);
327  dv_->scale(-one);
328  dv_->plus(*dv2_);
329  hv.axpy(mu_,*dv_);
330 
331  v_->set(v);
332  bnd_->pruneUpperActive(*v_,*xlam_);
333  v_->scale(-one);
334  v_->plus(v);
335  dv_->set(v_->dual());
336  dv2_->set(*dv_);
337  bnd_->pruneUpperActive(*dv_,*xlam_);
338  dv_->scale(-one);
339  dv_->plus(*dv2_);
340  hv.axpy(mu_,*dv_);
341  }
342  }
343 }; // class MoreauYosidaObjective
344 
345 } // namespace ROL
346 
347 #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 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)