ROL
ROL_InteriorPointPrimalDualResidual.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
45#ifndef ROL_INTERIORPOINT_PRIMALDUAL_RESIDUAL_H
46#define ROL_INTERIORPOINT_PRIMALDUAL_RESIDUAL_H
47
48#include "ROL_Elementwise_Function.hpp"
49#include "ROL_Constraint.hpp"
51#include "ROL_Objective.hpp"
53
54namespace ROL {
55namespace InteriorPoint {
56
73template<class Real> class PrimalDualSymmetrizer;
74
75
76template<class Real>
77class PrimalDualResidual : public Constraint<Real> {
78
79private:
80 typedef Vector<Real> V;
84
85
86 typedef typename PV::size_type size_type;
87
88 ROL::Ptr<OBJ> obj_; // Objective function
89 ROL::Ptr<CON> eqcon_; // Constraint
90 ROL::Ptr<CON> incon_; // Inequality Constraint
91
92 ROL::Ptr<V> qo_; // Storage for optimization variables
93 ROL::Ptr<V> qs_; // Storage for slack variables
94 ROL::Ptr<V> qe_; // Storage for equality multiplier variables
95 ROL::Ptr<V> qi_; // Storage for inequality multiplier variables
96
97 Real mu_; // Penalty parameter
98
99 ROL::Ptr<LinearOperator<Real> > sym_;
100
101 const static size_type OPT = 0; // Optimization vector
102 const static size_type SLACK = 1; // Slack vector
103 const static size_type EQUAL = 2; // Lagrange multipliers for equality constraint
104 const static size_type INEQ = 3; // Lagrange multipliers for inequality constraint
105
106
107
108public:
109
110 PrimalDualResidual( const ROL::Ptr<OBJ> &obj,
111 const ROL::Ptr<CON> &eqcon,
112 const ROL::Ptr<CON> &incon,
113 const V& x ) :
114 obj_(obj), eqcon_(eqcon), incon_(incon), mu_(1.0) {
115
116 // Allocate storage vectors
117 const PV &xpv = dynamic_cast<const PV&>(x);
118
119 qo_ = xpv.get(OPT)->clone();
120 qs_ = xpv.get(SLACK)->clone();
121 qe_ = xpv.get(EQUAL)->clone();
122 qi_ = xpv.get(INEQ)->clone();
123
124 sym_ = ROL::makePtr<PrimalDualSymmetrizer<Real>>(*qs_);
125
126 }
127
128 void value( V &c, const V &x, Real &tol ) {
129
130
131
132 // Downcast to partitioned vectors
133 PV &cpv = dynamic_cast<PV&>(c);
134 const PV &xpv = dynamic_cast<const PV&>(x);
135
136 ROL::Ptr<const V> xo = xpv.get(OPT);
137 ROL::Ptr<const V> xs = xpv.get(SLACK);
138 ROL::Ptr<const V> xe = xpv.get(EQUAL);
139 ROL::Ptr<const V> xi = xpv.get(INEQ);
140
141 c.zero();
142
143 ROL::Ptr<V> co = cpv.get(OPT);
144 ROL::Ptr<V> cs = cpv.get(SLACK);
145 ROL::Ptr<V> ce = cpv.get(EQUAL);
146 ROL::Ptr<V> ci = cpv.get(INEQ);
147
148 // Optimization components
149 obj_->gradient(*co,*xo,tol);
150
151 // Apply adjoint equality Jacobian at xo to xe and store in qo
152 eqcon_->applyAdjointJacobian(*qo_,*xe,*xo,tol);
153 co->axpy(-1.0,*qo_);
154
155 incon_->applyAdjointJacobian(*qo_,*xi,*xo,tol);
156 co->axpy(-1.0,*qo_);
157
158 // Slack components
159 cs->set(*xs);
160
161 Elementwise::Multiply<Real> mult;
162 cs->applyBinary(mult,*xi);
163
164 Elementwise::Fill<Real> fill(-mu_);
165 qs_->applyUnary(fill);
166
167 cs->plus(*qs_); // Sz-e
168
169 // component
170 eqcon_->value(*ce, *xo, tol); // c_e(x)
171
172 // Inequality component
173 incon_->value(*ci, *xo, tol); // c_i(x)-s
174 ci->axpy(-1.0, *xs);
175
176 sym_->update(*xs);
177 sym_->apply(c,c,tol);
178 sym_->applyInverse(c,c,tol);
179
180 }
181
182 void applyJacobian( V &jv, const V &v, const V &x, Real &tol ) {
183
184
185
186 jv.zero();
187
188 // Downcast to partitioned vectors
189 PV &jvpv = dynamic_cast<PV&>(jv);
190 const PV &vpv = dynamic_cast<const PV&>(v);
191 const PV &xpv = dynamic_cast<const PV&>(x);
192
193 ROL::Ptr<V> jvo = jvpv.get(OPT);
194 ROL::Ptr<V> jvs = jvpv.get(SLACK);
195 ROL::Ptr<V> jve = jvpv.get(EQUAL);
196 ROL::Ptr<V> jvi = jvpv.get(INEQ);
197
198 ROL::Ptr<const V> vo = vpv.get(OPT);
199 ROL::Ptr<const V> vs = vpv.get(SLACK);
200 ROL::Ptr<const V> ve = vpv.get(EQUAL);
201 ROL::Ptr<const V> vi = vpv.get(INEQ);
202
203 ROL::Ptr<const V> xo = xpv.get(OPT);
204 ROL::Ptr<const V> xs = xpv.get(SLACK);
205 ROL::Ptr<const V> xe = xpv.get(EQUAL);
206 ROL::Ptr<const V> xi = xpv.get(INEQ);
207
208 // Optimization components
209 obj_->hessVec(*jvo,*vo,*xo,tol);
210
211 eqcon_->applyAdjointHessian(*qo_,*xe,*vo,*xo,tol);
212
213 jvo->axpy(-1.0,*qo_);
214
215 incon_->applyAdjointHessian(*qo_,*xi,*vo,*xo,tol);
216
217 jvo->axpy(-1.0,*qo_);
218
219 eqcon_->applyAdjointJacobian(*qo_,*ve,*xo,tol);
220
221 jvo->axpy(-1.0,*qo_);
222
223 incon_->applyAdjointJacobian(*qo_,*vi,*xo,tol);
224
225 jvo->axpy(-1.0,*qo_);
226
227
228 // Slack components
229 jvs->set(*vs);
230
231 Elementwise::Multiply<Real> mult;
232
233 jvs->applyBinary(mult,*xi);
234
235 qs_->set(*vi);
236
237 qs_->applyBinary(mult,*xs);
238
239 jvs->plus(*qs_);
240
241 // component
242 eqcon_->applyJacobian(*jve,*vo,*xo,tol);
243
244 // Inequality components
245 incon_->applyJacobian(*jvi,*vo,*xo,tol);
246
247 jvi->axpy(-1.0,*vs);
248
249 sym_->update(*xs);
250 sym_->apply(jv,jv,tol);
251 sym_->applyInverse(jv,jv,tol);
252
253 }
254
255 void updatePenalty( Real mu ) {
256 mu_ = mu;
257 }
258
259}; // class PrimalDualResidual
260
261
262
263// Applying this operator to the left- and right-hand-sides, will
264// symmetrize the Primal-Dual Interior-Point KKT system, yielding
265// equation (19.13) from Nocedal & Wright
266
267template<class Real>
269
270 typedef Vector<Real> V;
272
273 typedef typename PV::size_type size_type;
274
275private:
276 ROL::Ptr<V> s_;
277
278 const static size_type OPT = 0; // Optimization vector
279 const static size_type SLACK = 1; // Slack vector
280 const static size_type EQUAL = 2; // Lagrange multipliers for equality constraint
281 const static size_type INEQ = 3; // Lagrange multipliers for inequality constraint
282
283public:
284
286 s_ = s.clone();
287 s_->set(s);
288 }
289
290 void update( const V& s, bool flag = true, int iter = -1 ) {
291 s_->set(s);
292 }
293
294 void apply( V &Hv, const V &v, Real &tol ) const {
295
296
297
298
299 const PV &vpv = dynamic_cast<const PV&>(v);
300 PV &Hvpv = dynamic_cast<PV&>(Hv);
301
302 ROL::Ptr<const V> vo = vpv.get(OPT);
303 ROL::Ptr<const V> vs = vpv.get(SLACK);
304 ROL::Ptr<const V> ve = vpv.get(EQUAL);
305 ROL::Ptr<const V> vi = vpv.get(INEQ);
306
307 ROL::Ptr<V> Hvo = Hvpv.get(OPT);
308 ROL::Ptr<V> Hvs = Hvpv.get(SLACK);
309 ROL::Ptr<V> Hve = Hvpv.get(EQUAL);
310 ROL::Ptr<V> Hvi = Hvpv.get(INEQ);
311
312 Hvo->set(*vo);
313
314 Hvs->set(*vs);
315 Elementwise::Divide<Real> div;
316 Hvs->applyBinary(div,*s_);
317
318 Hve->set(*ve);
319 Hve->scale(-1.0);
320
321 Hvi->set(*vi);
322 Hvi->scale(-1.0);
323
324 }
325
326 void applyInverse( V &Hv, const V &v, Real &tol ) const {
327
328
329
330
331 const PV &vpv = dynamic_cast<const PV&>(v);
332 PV &Hvpv = dynamic_cast<PV&>(Hv);
333
334 ROL::Ptr<const V> vo = vpv.get(OPT);
335 ROL::Ptr<const V> vs = vpv.get(SLACK);
336 ROL::Ptr<const V> ve = vpv.get(EQUAL);
337 ROL::Ptr<const V> vi = vpv.get(INEQ);
338
339 ROL::Ptr<V> Hvo = Hvpv.get(OPT);
340 ROL::Ptr<V> Hvs = Hvpv.get(SLACK);
341 ROL::Ptr<V> Hve = Hvpv.get(EQUAL);
342 ROL::Ptr<V> Hvi = Hvpv.get(INEQ);
343
344 Hvo->set(*vo);
345
346 Hvs->set(*vs);
347 Elementwise::Multiply<Real> mult;
348 Hvs->applyBinary(mult,*s_);
349
350 Hve->set(*ve);
351 Hve->scale(-1.0);
352
353 Hvi->set(*vi);
354 Hvi->scale(-1.0);
355
356 }
357}; // class PrimalDualSymmetrizer
358
359
360} // namespace InteriorPoint
361
362
363
364
365
366} // namespace ROL
367
368
369#endif // ROL_INTERIORPOINT_PRIMALDUAL_RESIDUAL_H
370
Defines the general constraint operator interface.
PrimalDualResidual(const ROL::Ptr< OBJ > &obj, const ROL::Ptr< CON > &eqcon, const ROL::Ptr< CON > &incon, const V &x)
void applyJacobian(V &jv, const V &v, const V &x, Real &tol)
Apply the constraint Jacobian at , , to vector .
void value(V &c, const V &x, Real &tol)
Evaluate the constraint operator at .
void update(const V &s, bool flag=true, int iter=-1)
Update linear operator.
void applyInverse(V &Hv, const V &v, Real &tol) const
Apply inverse of linear operator.
void apply(V &Hv, const V &v, Real &tol) const
Apply linear operator.
Provides the interface to apply a linear operator.
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
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:167
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.