ROL
ROL_PrimalDualInteriorPointResidual.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_PRIMALDUALINTERIORPOINTRESIDUAL_H
11 #define ROL_PRIMALDUALINTERIORPOINTRESIDUAL_H
12 
13 #include "ROL_BoundConstraint.hpp"
14 #include "ROL_Constraint.hpp"
15 #include "ROL_Objective.hpp"
17 #include "ROL_RandomVector.hpp"
18 
19 #include <iostream>
20 
41 namespace ROL {
42 
43 template<class Real>
44 class PrimalDualInteriorPointResidual : public Constraint<Real> {
45 
46  typedef ROL::ParameterList PL;
47 
48  typedef Vector<Real> V;
53 
54  typedef Elementwise::ValueSet<Real> ValueSet;
55 
56  typedef typename PV::size_type size_type;
57 
58 private:
59 
60  static const size_type OPT = 0;
61  static const size_type EQUAL = 1;
62  static const size_type LOWER = 2;
63  static const size_type UPPER = 3;
64 
65  const ROL::Ptr<OBJ> obj_;
66  const ROL::Ptr<CON> con_;
67  const ROL::Ptr<BND> bnd_;
68 
69  ROL::Ptr<const V> x_; // Optimization vector
70  ROL::Ptr<const V> l_; // constraint multiplier
71  ROL::Ptr<const V> zl_; // Lower bound multiplier
72  ROL::Ptr<const V> zu_; // Upper bound multiplier
73 
74  ROL::Ptr<const V> xl_; // Lower bound
75  ROL::Ptr<const V> xu_; // Upper bound
76 
77  const ROL::Ptr<const V> maskL_;
78  const ROL::Ptr<const V> maskU_;
79 
80  ROL::Ptr<V> scratch_; // Scratch vector the same dimension as x
81 
82  Real mu_;
83 
85 
86  Real one_;
87  Real zero_;
88 
89  int nfval_;
90  int ngrad_;
91  int ncval_;
92 
93  Elementwise::Multiply<Real> mult_;
94 
95  class SafeDivide : public Elementwise::BinaryFunction<Real> {
96  public:
97  Real apply( const Real &x, const Real &y ) const {
98  return y != 0 ? x/y : 0;
99  }
100  };
101 
103 
104  class SetZeros : public Elementwise::BinaryFunction<Real> {
105  public:
106  Real apply( const Real &x, const Real &y ) const {
107  return y==1.0 ? 0 : x;
108  }
109  };
110 
112 
113  // Fill in zeros of x with corresponding values of y
114  class InFill : public Elementwise::BinaryFunction<Real> {
115  public:
116  Real apply( const Real &x, const Real &y ) const {
117  return x == 0 ? y : x;
118  }
119  };
120 
122 
123  // Extract the optimization and lagrange multiplier
124  ROL::Ptr<V> getOptMult( V &vec ) {
125  PV &vec_pv = dynamic_cast<PV&>(vec);
126 
127  return CreatePartitioned(vec_pv.get(OPT),vec_pv.get(EQUAL));
128  }
129 
130  // Extract the optimization and lagrange multiplier
131  ROL::Ptr<const V> getOptMult( const V &vec ) {
132  const PV &vec_pv = dynamic_cast<const PV&>(vec);
133 
134  return CreatePartitioned(vec_pv.get(OPT),vec_pv.get(EQUAL));
135  }
136 
137 
138 
139 
140 public:
141 
142  PrimalDualInteriorPointResidual( const ROL::Ptr<OBJ> &obj,
143  const ROL::Ptr<CON> &con,
144  const ROL::Ptr<BND> &bnd,
145  const V &x,
146  const ROL::Ptr<const V> &maskL,
147  const ROL::Ptr<const V> &maskU,
148  ROL::Ptr<V> &scratch,
149  Real mu, bool symmetrize ) :
150  obj_(obj), con_(con), bnd_(bnd), xl_(bnd->getLowerBound()),
151  xu_(bnd->getUpperBound()), maskL_(maskL), maskU_(maskU), scratch_(scratch),
152  mu_(mu), symmetrize_(symmetrize), one_(1.0), zero_(0.0), nfval_(0),
153  ngrad_(0), ncval_(0) {
154 
155  // Get access to the four components
156  const PV &x_pv = dynamic_cast<const PV&>(x);
157 
158  x_ = x_pv.get(OPT);
159  l_ = x_pv.get(EQUAL);
160  zl_ = x_pv.get(LOWER);
161  zu_ = x_pv.get(UPPER);
162 
163  }
164 
165 
166  void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
167 
168  // Get access to the four components
169  const PV &x_pv = dynamic_cast<const PV&>(x);
170 
171  x_ = x_pv.get(OPT);
172  l_ = x_pv.get(EQUAL);
173  zl_ = x_pv.get(LOWER);
174  zu_ = x_pv.get(UPPER);
175 
176  obj_->update(*x_,flag,iter);
177  con_->update(*x_,flag,iter);
178 
179  }
180 
181 
182  // Evaluate the gradient of the Lagrangian
183  void value( V &c, const V &x, Real &tol ) {
184 
185 
186 
187  Elementwise::Shift<Real> subtract_mu(-mu_);
188  Elementwise::Fill<Real> fill_minus_mu(-mu_);
189 
190  const PV &x_pv = dynamic_cast<const PV&>(x);
191  PV &c_pv = dynamic_cast<PV&>(c);
192 
193  x_ = x_pv.get(OPT);
194  l_ = x_pv.get(EQUAL);
195  zl_ = x_pv.get(LOWER);
196  zu_ = x_pv.get(UPPER);
197 
198  ROL::Ptr<V> cx = c_pv.get(OPT);
199  ROL::Ptr<V> cl = c_pv.get(EQUAL);
200  ROL::Ptr<V> czl = c_pv.get(LOWER);
201  ROL::Ptr<V> czu = c_pv.get(UPPER);
202 
203  /********************************************************************************/
204  /* Optimization Components */
205  /********************************************************************************/
206  obj_->gradient(*cx,*x_,tol);
207  ngrad_++;
208 
209  con_->applyAdjointJacobian(*scratch_,*l_,*x_,tol);
210 
211  cx->plus(*scratch_);
212 
213  cx->axpy(-one_,*zl_);
214  cx->plus(*zu_); // cx = g+J'l-zl+zu
215 
216  /********************************************************************************/
217  /* Constraint Components */
218  /********************************************************************************/
219 
220  con_->value(*cl,*x_,tol);
221  ncval_++;
222 
223  /********************************************************************************/
224  /* Lower Bound Components */
225  /********************************************************************************/
226  if( symmetrize_ ) { // -(x-l)+mu/zl
227 
228  czl->applyUnary(fill_minus_mu);
229  czl->applyBinary(div_,*zl_);
230 
231  scratch_->set(*x_);
232  scratch_->axpy(-1.0,*xl_);
233 
234  czl->plus(*scratch_);
235  czl->scale(-1.0);
236  }
237  else { // czl = zl*(x-l)-mu*e
238  czl->set(*x_); // czl = x
239  czl->axpy(-1.0,*xl_); // czl = x-l
240  czl->applyBinary(mult_,*zl_); // czl = zl*(x-l)
241  czl->applyUnary(subtract_mu); // czl = zl*(x-l)-mu*e
242  }
243 
244  // Zero out elements corresponding to infinite lower bounds
245  czl->applyBinary(mult_,*maskL_);
246 
247  /********************************************************************************/
248  /* Upper Bound Components */
249  /********************************************************************************/
250  if( symmetrize_ ) { // -(u-x)+mu/zu
251 
252  czu->applyUnary(fill_minus_mu);
253  czu->applyBinary(div_,*zu_);
254 
255  scratch_->set(*xu_);
256  scratch_->axpy(-1.0, *x_);
257 
258  czu->plus(*scratch_);
259  czu->scale(-1.0);
260  }
261  else { // zu*(u-x)-mu*e
262  czu->set(*xu_); // czu = u
263  czu->axpy(-1.0,*x_); // czu = u-x
264  czu->applyBinary(mult_,*zu_); // czu = zu*(u-x)
265  czu->applyUnary(subtract_mu); // czu = zu*(u-x)-mu*e
266  }
267 
268  // Zero out elements corresponding to infinite upper bounds
269  czu->applyBinary(mult_,*maskU_);
270 
271  }
272 
273  // Evaluate the action of the Hessian of the Lagrangian on a vector
274  //
275  // [ J11 J12 J13 J14 ] [ vx ] [ jvx ] [ J11*vx + J12*vl + J13*vzl + J14*vzu ]
276  // [ J21 0 0 0 ] [ vl ] = [ jvl ] = [ J21*vx ]
277  // [ J31 0 J33 0 ] [ vzl ] [ jvzl ] [ J31*vx + J33*vzl ]
278  // [ J41 0 0 J44 ] [ vzu ] [ jvzu ] [ J41*vx + J44*vzu ]
279  //
280  void applyJacobian( V &jv, const V &v, const V &x, Real &tol ) {
281 
282 
283 
284  PV &jv_pv = dynamic_cast<PV&>(jv);
285  const PV &v_pv = dynamic_cast<const PV&>(v);
286  const PV &x_pv = dynamic_cast<const PV&>(x);
287 
288  // output vector components
289  ROL::Ptr<V> jvx = jv_pv.get(OPT);
290  ROL::Ptr<V> jvl = jv_pv.get(EQUAL);
291  ROL::Ptr<V> jvzl = jv_pv.get(LOWER);
292  ROL::Ptr<V> jvzu = jv_pv.get(UPPER);
293 
294  // input vector components
295  ROL::Ptr<const V> vx = v_pv.get(OPT);
296  ROL::Ptr<const V> vl = v_pv.get(EQUAL);
297  ROL::Ptr<const V> vzl = v_pv.get(LOWER);
298  ROL::Ptr<const V> vzu = v_pv.get(UPPER);
299 
300  x_ = x_pv.get(OPT);
301  l_ = x_pv.get(EQUAL);
302  zl_ = x_pv.get(LOWER);
303  zu_ = x_pv.get(UPPER);
304 
305 
306  /********************************************************************************/
307  /* Optimization Components */
308  /********************************************************************************/
309 
310  obj_->hessVec(*jvx,*vx,*x_,tol);
311  con_->applyAdjointHessian(*scratch_,*l_,*vx,*x_,tol);
312  jvx->plus(*scratch_);
313  con_->applyAdjointJacobian(*scratch_,*vl,*x_,tol);
314  jvx->plus(*scratch_);
315 
316  // H_13 = -I for l_i > -infty
317  scratch_->set(*vzl);
318  scratch_->applyBinary(mult_,*maskL_);
319  jvx->axpy(-1.0,*scratch_);
320 
321  // H_14 = I for u_i < infty
322  scratch_->set(*vzu);
323  scratch_->applyBinary(mult_,*maskU_);
324  jvx->plus(*scratch_);
325 
326  /********************************************************************************/
327  /* Constraint Components */
328  /********************************************************************************/
329 
330  con_->applyJacobian(*jvl,*vx,*x_,tol);
331 
332  /********************************************************************************/
333  /* Lower Bound Components */
334  /********************************************************************************/
335 
336  if( symmetrize_ ) {
337  // czl = x-l-mu/zl
338  // jvzl = -vx - inv(Zl)*(X-L)*vzl
339 
340  jvzl->set(*x_);
341  jvzl->axpy(-1.0,*xl_);
342  jvzl->applyBinary(mult_,*vzl);
343  jvzl->applyBinary(div_,*zl_);
344 
345  jvzl->plus(*vx);
346  jvzl->scale(-1.0);
347 
348  }
349 
350  else {
351  // czl = zl*(x-l)-mu*e
352  // jvzl = Zl*vx + (X-L)*vzl
353 
354  // H_31 = Zl
355  jvzl->set(*vx);
356  jvzl->applyBinary(mult_,*zl_);
357 
358  // H_33 = X-L
359  scratch_->set(*x_);
360  scratch_->axpy(-1.0,*xl_);
361  scratch_->applyBinary(mult_,*vzl);
362 
363  jvzl->plus(*scratch_);
364 
365  }
366 
367  // jvzl[i] = vzl[i] if l[i] = -inf
368  jvzl->applyBinary(mult_,*maskL_);
369  jvzl->applyBinary(inFill_,*vzl);
370 
371  /********************************************************************************/
372  /* Upper Bound Components */
373  /********************************************************************************/
374 
375  if( symmetrize_ ) {
376  // czu = u-x-mu/zu
377  // jvzu = vx - inv(Zu)*(U-X)*vzu
378 
379  jvzu->set(*xu_);
380  jvzu->axpy(-1.0,*x_);
381  jvzu->applyBinary(mult_,*vzu);
382  jvzu->applyBinary(div_,*zu_);
383  jvzu->scale(-1.0);
384  jvzu->plus(*vx);
385 
386  }
387  else {
388  // czu = zu*(u-x)-mu*e
389  // jvzu = -Zu*vx + (U-X)*vzu
390 
391  // H_41 = -Zu
392  scratch_->set(*zu_);
393  scratch_->applyBinary(mult_,*vx);
394 
395  // H_44 = U-X
396  jvzu->set(*xu_);
397  jvzu->axpy(-1.0,*x_);
398  jvzu->applyBinary(mult_,*vzu);
399 
400  jvzu->axpy(-1.0,*scratch_);
401 
402  }
403 
404  // jvzu[i] = vzu[i] if u[i] = inf
405  jvzu->applyBinary(mult_,*maskU_);
406  jvzu->applyBinary(inFill_,*vzu);
407 
408  }
409 
410  // Call this whenever mu changes
411  void reset( const Real mu ) {
412  mu_ = mu;
413  nfval_ = 0;
414  ngrad_ = 0;
415  ncval_ = 0;
416  }
417 
419  return nfval_;
420  }
421 
423  return ngrad_;
424  }
425 
427  return ncval_;
428  }
429 
430 }; // class PrimalDualInteriorPointResidual
431 
432 } // namespace ROL
433 
434 #endif // ROL_PRIMALDUALINTERIORPOINTRESIDUAL_H
Provides the interface to evaluate objective functions.
void value(V &c, const V &x, Real &tol)
Evaluate the constraint operator at .
ROL::Ptr< const Vector< Real > > get(size_type i) const
Defines the linear algebra of vector space on a generic partitioned vector.
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 ...
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
PrimalDualInteriorPointResidual(const ROL::Ptr< OBJ > &obj, const ROL::Ptr< CON > &con, const ROL::Ptr< BND > &bnd, const V &x, const ROL::Ptr< const V > &maskL, const ROL::Ptr< const V > &maskU, ROL::Ptr< V > &scratch, Real mu, bool symmetrize)
void applyJacobian(V &jv, const V &v, const V &x, Real &tol)
Apply the constraint Jacobian at , , to vector .
Provides the interface to apply upper and lower bound constraints.
std::vector< PV >::size_type size_type
Defines the general constraint operator interface.