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