ROL
ROL_Fletcher.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_FLETCHER_H
46#define ROL_FLETCHER_H
47
48#include "ROL_FletcherBase.hpp"
49#include "ROL_Objective.hpp"
50#include "ROL_Constraint.hpp"
51#include "ROL_Vector.hpp"
52#include "ROL_Types.hpp"
53#include "ROL_Ptr.hpp"
54#include "ROL_Krylov.hpp"
56#include <iostream>
57
58namespace ROL {
59
60template <class Real>
61class Fletcher : public FletcherBase<Real> {
62private:
63 // Required for Fletcher penalty function definition
64 using FletcherBase<Real>::obj_;
65 using FletcherBase<Real>::con_;
66
69
70 // Evaluation counters
71 using FletcherBase<Real>::nfval_;
72 using FletcherBase<Real>::ngval_;
73 using FletcherBase<Real>::ncval_;
74
75 using FletcherBase<Real>::fPhi_; // value of penalty function
76 using FletcherBase<Real>::gPhi_; // gradient of penalty function
77
78 using FletcherBase<Real>::y_; // multiplier estimate
79
80 using FletcherBase<Real>::fval_; // value of objective function
81 using FletcherBase<Real>::g_; // gradient of objective value
82 using FletcherBase<Real>::c_; // constraint value
83 using FletcherBase<Real>::scaledc_; // penaltyParameter_ * c_
84 using FletcherBase<Real>::gL_; // gradient of Lagrangian (g - A*y)
85
86 using FletcherBase<Real>::cnorm_; // norm of constraint violation
87
94
95 using FletcherBase<Real>::delta_; // regularization parameter
96
98
99 // Temporaries
100 Ptr<Vector<Real> > Tv_; // Temporary for matvecs
101 Ptr<Vector<Real> > w_; // first component of augmented system solve solution
102 Ptr<Vector<Real> > v_; // second component of augmented system solve solution
103 Ptr<Vector<Real> > wg_; // first component of augmented system solve solution for gradient
104 Ptr<Vector<Real> > vg_; // second component of augmented system solve solution for gradient
105
106 Ptr<Vector<Real> > xzeros_; // zero vector
107 Ptr<Vector<Real> > czeros_; // zero vector
108
109 using FletcherBase<Real>::useInexact_;
110
111 using FletcherBase<Real>::multSolverError_; // Error from augmented system solve in value()
112 using FletcherBase<Real>::gradSolveError_; // Error from augmented system solve in gradient()
113
114 // For Augmented system solves
115 using FletcherBase<Real>::krylov_;
116 using FletcherBase<Real>::iterKrylov_;
117 using FletcherBase<Real>::flagKrylov_;
118 using FletcherBase<Real>::v1_;
119 using FletcherBase<Real>::v2_;
120 using FletcherBase<Real>::vv_;
121 using FletcherBase<Real>::w1_;
122 using FletcherBase<Real>::w2_;
123 using FletcherBase<Real>::ww_;
124 using FletcherBase<Real>::b1_;
125 using FletcherBase<Real>::b2_;
126 using FletcherBase<Real>::bb_;
127
128 class AugSystem : public LinearOperator<Real> {
129 private:
130 const Ptr<Constraint<Real> > con_;
131 const Ptr<const Vector<Real> > x_;
132 const Real delta_;
133 public:
134 AugSystem(const Ptr<Constraint<Real> > &con,
135 const Ptr<const Vector<Real> > &x,
136 const Real delta) : con_(con), x_(x), delta_(delta) {}
137
138 void apply(Vector<Real> &Hv, const Vector<Real> &v, Real &tol) const {
139 PartitionedVector<Real> &Hvp = dynamic_cast<PartitionedVector<Real>&>(Hv);
140 const PartitionedVector<Real> &vp = dynamic_cast<const PartitionedVector<Real>&>(v);
141
142 con_->applyAdjointJacobian(*(Hvp.get(0)), *(vp.get(1)), *x_, tol);
143 Hvp.get(0)->plus(*(vp.get(0)));
144
145 con_->applyJacobian(*(Hvp.get(1)), *(vp.get(0)), *x_, tol);
146 Hvp.get(1)->axpy(-delta_*delta_, *(vp.get(1)));
147 }
148 };
149
150 class AugSystemPrecond : public LinearOperator<Real> {
151 private:
152 const Ptr<Constraint<Real> > con_;
153 const Ptr<const Vector<Real> > x_;
154 public:
156 const Ptr<const Vector<Real> > x) : con_(con), x_(x) {}
157
158 void apply(Vector<Real> &Hv, const Vector<Real> &v, Real &tol) const {
159 Hv.set(v.dual());
160 }
161 void applyInverse(Vector<Real> &Hv, const Vector<Real> &v, Real &tol) const {
162 Real zero(0);
163 PartitionedVector<Real> &Hvp = dynamic_cast<PartitionedVector<Real>&>(Hv);
164 const PartitionedVector<Real> &vp = dynamic_cast<const PartitionedVector<Real>&>(v);
165
166 Hvp.set(0, *(vp.get(0)));
167 // Second x should be dual, but unused?
168 con_->applyPreconditioner(*(Hvp.get(1)),*(vp.get(1)),*x_,*x_, zero);
169 }
170 };
171
172public:
173 Fletcher(const ROL::Ptr<Objective<Real> > &obj,
174 const ROL::Ptr<Constraint<Real> > &con,
175 const Vector<Real> &optVec,
176 const Vector<Real> &conVec,
177 ROL::ParameterList &parlist)
178 : FletcherBase<Real>(obj, con) {
179
180 gPhi_ = optVec.dual().clone();
181 y_ = conVec.dual().clone();
182 g_ = optVec.dual().clone();
183 gL_ = optVec.dual().clone();
184 c_ = conVec.clone();
185 scaledc_ = conVec.clone();
186
187 Tv_ = optVec.dual().clone();
188 w_ = optVec.dual().clone();
189 v_ = conVec.dual().clone();
190 wg_ = optVec.dual().clone();
191 vg_ = conVec.dual().clone();
192
193 xzeros_ = optVec.dual().clone();
194 xzeros_->zero();
195 czeros_ = conVec.clone();
196 czeros_->zero();
197
198 v1_ = optVec.dual().clone();
199 v2_ = conVec.dual().clone();
200 vv_ = makePtr<PartitionedVector<Real>>(std::vector<Ptr<Vector<Real>> >({v1_, v2_}));
201
202 w1_ = optVec.dual().clone();
203 w2_ = conVec.dual().clone();
204 ww_ = makePtr<PartitionedVector<Real>>(std::vector<Ptr<Vector<Real>> >({w1_, w2_}));
205
206 b1_ = optVec.dual().clone();
207 b2_ = conVec.clone();
208 bb_ = makePtr<PartitionedVector<Real>>(std::vector<Ptr<Vector<Real>> >({b1_, b2_}));
209
210 ROL::ParameterList& sublist = parlist.sublist("Step").sublist("Fletcher");
211 HessianApprox_ = sublist.get("Level of Hessian Approximation", 0);
212 penaltyParameter_ = sublist.get("Penalty Parameter", 1.0);
213 quadPenaltyParameter_ = sublist.get("Quadratic Penalty Parameter", 0.0);
214
215 delta_ = sublist.get("Regularization Parameter", 0.0);
216
217 useInexact_ = sublist.get("Inexact Solves", false);
218
219 ROL::ParameterList krylovList;
220 Real atol = static_cast<Real>(1e-12);
221 Real rtol = static_cast<Real>(1e-2);
222 krylovList.sublist("General").sublist("Krylov").set("Type", "GMRES");
223 krylovList.sublist("General").sublist("Krylov").set("Absolute Tolerance", atol);
224 krylovList.sublist("General").sublist("Krylov").set("Relative Tolerance", rtol);
225 krylovList.sublist("General").sublist("Krylov").set("Iteration Limit", 200);
226 krylov_ = KrylovFactory<Real>(krylovList);
227 }
228
229 void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
230 obj_->update(x,flag,iter);
231 con_->update(x,flag,iter);
232 isValueComputed_ = (flag ? false : isValueComputed_);
233 isGradientComputed_ = (flag ? false : isGradientComputed_);
235 isObjValueComputed_ = (flag ? false : isObjValueComputed_);
236 isObjGradComputed_ = (flag ? false : isObjGradComputed_);
237 isConValueComputed_ = (flag ? false : isConValueComputed_);
238 }
239
240 Real value( const Vector<Real> &x, Real &tol ) {
243 return fPhi_;
244 }
245
246 Real zero(0);
247
248 // Reset tolerances
249 Real origTol = tol;
250 Real tol2 = origTol;
251
252 FletcherBase<Real>::objValue(x, tol2); tol2 = origTol;
253 multSolverError_ = origTol / (static_cast<Real>(2) * std::max(static_cast<Real>(1), cnorm_));
256
257 fPhi_ = fval_ - c_->dot(y_->dual());
258
260 fPhi_ = fPhi_ + Real(0.5)*quadPenaltyParameter_*(c_->dot(c_->dual()));
261 }
262
263 isValueComputed_ = true;
264
265 return fPhi_;
266 }
267
268 void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
269 if( isGradientComputed_ && gradSolveError_ <= tol) {
270 tol = gradSolveError_;
271 g.set(*gPhi_);
272 return;
273 }
274
275 Real zero(0);
276
277 // Reset tolerances
278 Real origTol = tol;
279 Real tol2 = origTol;
280
281 gradSolveError_ = origTol / static_cast<Real>(2);
283
284 bool refine = isGradientComputed_;
285
286 // gPhi = sum y_i H_i w + sigma w + sum v_i H_i gL - H w + gL
289 tol = gradSolveError_;
290
291 con_->applyAdjointHessian( *gPhi_, *y_, *wg_, x, tol2 ); tol2 = origTol;
292 gPhi_->axpy( penaltyParameter_, *wg_ );
293
294 obj_->hessVec( *Tv_, *wg_, x, tol2 ); tol2 = origTol;
295 gPhi_->axpy( static_cast<Real>(-1), *Tv_ );
296
297 con_->applyAdjointHessian( *Tv_, *vg_, *gL_, x, tol2 ); tol2 = origTol;
298 gPhi_->plus( *Tv_ );
299
300 gPhi_->plus( *gL_ );
301
303 con_->applyAdjointJacobian( *Tv_, *c_, x, tol2 ); tol2 = origTol;
305 }
306
307 g.set(*gPhi_);
308 isGradientComputed_ = true;
309 }
310
311 void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
312 Real zero(0);
313
314 // Reset tolerances
315 Real origTol = tol;
316 Real tol2 = origTol;
317
319 // hessVec tol is always set to ~1e-8. So if inexact linear system solves are used, then
320 // the multipliers will always get re-evaluated to high precision, which defeats the purpose
321 // of computing the objective/gradient approximately.
322 // Thus if inexact linear system solves are used, we will not update the multiplier estimates
323 // to high precision.
324 computeMultipliers(x, tol);
325 }
326
327 obj_->hessVec( hv, v, x, tol2 ); tol2 = origTol;
328 con_->applyAdjointHessian( *Tv_, *y_, v, x, tol2 ); tol2 = origTol;
329 hv.axpy(static_cast<Real>(-1), *Tv_ );
330
331 tol2 = tol;
332 solveAugmentedSystem( *w_, *v_, hv, *czeros_, x, tol2 ); tol2 = origTol;
333 hv.scale( static_cast<Real>(-1) );
334 hv.plus( *w_ );
335
336 Tv_->set(v);
337 tol2 = tol;
338 solveAugmentedSystem( *w_, *v_, *Tv_, *czeros_, x, tol2 ); tol2 = origTol;
339 hv.axpy(static_cast<Real>(-2)*penaltyParameter_, *w_);
340
341 obj_->hessVec( *Tv_, *w_, x, tol2 ); tol2 = origTol;
342 hv.plus( *Tv_ );
343 con_->applyAdjointHessian( *Tv_, *y_, *w_, x, tol2 ); tol2 = origTol;
344 hv.axpy( static_cast<Real>(-1), *Tv_ );
345
346 hv.axpy( static_cast<Real>(2)*penaltyParameter_, v );
347
349 con_->applyJacobian( *b2_, v, x, tol2 ); tol2 = origTol;
350 con_->applyAdjointJacobian( *Tv_, *b2_, x, tol2 ); tol2 = origTol;
352 con_->applyAdjointHessian( *Tv_, *c_, v, x, tol2); tol2 = origTol;
354 }
355
356 }
357
359 Vector<Real> &v2,
360 const Vector<Real> &b1,
361 const Vector<Real> &b2,
362 const Vector<Real> &x,
363 Real &tol,
364 bool refine = false) {
365 // Ignore tol for now
366 ROL::Ptr<LinearOperator<Real> > K
367 = ROL::makePtr<AugSystem>(con_, makePtrFromRef(x), delta_);
368 ROL::Ptr<LinearOperator<Real> > P
369 = ROL::makePtr<AugSystemPrecond>(con_, makePtrFromRef(x));
370
371 b1_->set(b1);
372 b2_->set(b2);
373
374 if( refine ) {
375 // TODO: Make sure this tol is actually ok...
376 Real origTol = tol;
377 w1_->set(v1);
378 w2_->set(v2);
379 K->apply(*vv_, *ww_, tol); tol = origTol;
380
381 b1_->axpy( static_cast<Real>(-1), *v1_ );
382 b2_->axpy( static_cast<Real>(-1), *v2_ );
383 }
384
385 v1_->zero();
386 v2_->zero();
387
388 // If inexact, change tolerance
389 if( useInexact_ ) {
390 krylov_->resetAbsoluteTolerance(tol);
391 }
392
393 flagKrylov_ = 0;
394 tol = krylov_->run(*vv_,*K,*bb_,*P,iterKrylov_,flagKrylov_);
395
396 if( refine ) {
397 v1.plus(*v1_);
398 v2.plus(*v2_);
399 } else {
400 v1.set(*v1_);
401 v2.set(*v2_);
402 }
403 }
404
405 void computeMultipliers(const Vector<Real>& x, const Real tol) {
407 return;
408 }
409
410 if( !isMultiplierComputed_ ) {
411 Real tol2 = tol;
412 FletcherBase<Real>::objGrad(x, tol2); tol2 = tol;
414 cnorm_ = c_->norm();
415 }
416
417 bool refine = isMultiplierComputed_;
418
419 multSolverError_ = tol;
421
423 }
424
425}; // class Fletcher
426
427} // namespace ROL
428
429#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0 zero)()
Contains definitions of custom data types in ROL.
Defines the general constraint operator interface.
Ptr< Vector< Real > > v1_
Ptr< Vector< Real > > b2_
Ptr< Vector< Real > > c_
Ptr< Vector< Real > > w1_
const Ptr< Constraint< Real > > con_
void conValue(const Vector< Real > &x, Real &tol)
Ptr< Vector< Real > > scaledc_
Ptr< Vector< Real > > v2_
const Ptr< Objective< Real > > obj_
Ptr< Vector< Real > > gPhi_
Ptr< Vector< Real > > w2_
void objGrad(const Vector< Real > &x, Real &tol)
void objValue(const Vector< Real > &x, Real &tol)
Ptr< Vector< Real > > g_
Ptr< Vector< Real > > y_
Ptr< PartitionedVector< Real > > vv_
Ptr< Vector< Real > > gL_
Ptr< PartitionedVector< Real > > ww_
Ptr< PartitionedVector< Real > > bb_
Ptr< Krylov< Real > > krylov_
Ptr< Vector< Real > > b1_
const Ptr< const Vector< Real > > x_
const Ptr< Constraint< Real > > con_
AugSystemPrecond(const Ptr< Constraint< Real > > con, const Ptr< const Vector< Real > > x)
void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
void applyInverse(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply inverse of linear operator.
const Ptr< Constraint< Real > > con_
void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
const Ptr< const Vector< Real > > x_
AugSystem(const Ptr< Constraint< Real > > &con, const Ptr< const Vector< Real > > &x, const Real delta)
Ptr< Vector< Real > > v_
Fletcher(const ROL::Ptr< Objective< Real > > &obj, const ROL::Ptr< Constraint< Real > > &con, const Vector< Real > &optVec, const Vector< Real > &conVec, ROL::ParameterList &parlist)
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
Ptr< Vector< Real > > w_
void solveAugmentedSystem(Vector< Real > &v1, Vector< Real > &v2, const Vector< Real > &b1, const Vector< Real > &b2, const Vector< Real > &x, Real &tol, bool refine=false)
Ptr< Vector< Real > > wg_
Real value(const Vector< Real > &x, Real &tol)
Compute value.
Ptr< Vector< Real > > czeros_
void computeMultipliers(const Vector< Real > &x, const Real tol)
Ptr< Vector< Real > > Tv_
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Ptr< Vector< Real > > vg_
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Ptr< Vector< Real > > xzeros_
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.
ROL::Ptr< const Vector< Real > > get(size_type i) const
void set(const V &x)
Set where .
Defines the linear algebra or vector space interface.
virtual void set(const Vector &x)
Set where .
virtual void scale(const Real alpha)=0
Compute where .
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
virtual void plus(const Vector &x)=0
Compute , where .
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .