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