ROL
interiorpoint/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 }
256
257 void updatePenalty( Real mu ) {
258 mu_ = mu;
259 }
260
261}; // class PrimalDualResidual
262
263
264
265// Applying this operator to the left- and right-hand-sides, will
266// symmetrize the Primal-Dual Interior-Point KKT system, yielding
267// equation (19.13) from Nocedal & Wright
268
269template<class Real>
270class PrimalDualSymmetrizer : public LinearOperator<Real> {
271
272 typedef Vector<Real> V;
274
275 typedef typename PV::size_type size_type;
276
277private:
278 ROL::Ptr<V> s_;
279
280 const static size_type OPT = 0; // Optimization vector
281 const static size_type SLACK = 1; // Slack vector
282 const static size_type EQUAL = 2; // Lagrange multipliers for equality constraint
283 const static size_type INEQ = 3; // Lagrange multipliers for inequality constraint
284
285public:
286
288 s_ = s.clone();
289 s_->set(s);
290 }
291
292 void update( const V& s, bool flag = true, int iter = -1 ) {
293 s_->set(s);
294 }
295
296 void apply( V &Hv, const V &v, Real &tol ) const {
297
298
299
300
301 const PV &vpv = dynamic_cast<const PV&>(v);
302 PV &Hvpv = dynamic_cast<PV&>(Hv);
303
304 ROL::Ptr<const V> vo = vpv.get(OPT);
305 ROL::Ptr<const V> vs = vpv.get(SLACK);
306 ROL::Ptr<const V> ve = vpv.get(EQUAL);
307 ROL::Ptr<const V> vi = vpv.get(INEQ);
308
309 ROL::Ptr<V> Hvo = Hvpv.get(OPT);
310 ROL::Ptr<V> Hvs = Hvpv.get(SLACK);
311 ROL::Ptr<V> Hve = Hvpv.get(EQUAL);
312 ROL::Ptr<V> Hvi = Hvpv.get(INEQ);
313
314 Hvo->set(*vo);
315
316 Hvs->set(*vs);
317 Elementwise::Divide<Real> div;
318 Hvs->applyBinary(div,*s_);
319
320 Hve->set(*ve);
321 Hve->scale(-1.0);
322
323 Hvi->set(*vi);
324 Hvi->scale(-1.0);
325
326 }
327
328 void applyInverse( V &Hv, const V &v, Real &tol ) const {
329
330
331
332
333 const PV &vpv = dynamic_cast<const PV&>(v);
334 PV &Hvpv = dynamic_cast<PV&>(Hv);
335
336 ROL::Ptr<const V> vo = vpv.get(OPT);
337 ROL::Ptr<const V> vs = vpv.get(SLACK);
338 ROL::Ptr<const V> ve = vpv.get(EQUAL);
339 ROL::Ptr<const V> vi = vpv.get(INEQ);
340
341 ROL::Ptr<V> Hvo = Hvpv.get(OPT);
342 ROL::Ptr<V> Hvs = Hvpv.get(SLACK);
343 ROL::Ptr<V> Hve = Hvpv.get(EQUAL);
344 ROL::Ptr<V> Hvi = Hvpv.get(INEQ);
345
346 Hvo->set(*vo);
347
348 Hvs->set(*vs);
349 Elementwise::Multiply<Real> mult;
350 Hvs->applyBinary(mult,*s_);
351
352 Hve->set(*ve);
353 Hve->scale(-1.0);
354
355 Hvi->set(*vi);
356 Hvi->scale(-1.0);
357
358 }
359}; // class PrimalDualSymmetrizer
360
361
362} // namespace InteriorPoint
363
364
365
366
367
368} // namespace ROL
369
370
371#endif // ROL_INTERIORPOINT_PRIMALDUAL_RESIDUAL_H
372
Defines the general constraint operator interface.
Express the Primal-Dual Interior Point gradient as an equality constraint.
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.