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 
58 namespace ROL {
59 
60 template <class Real>
61 class Fletcher : public FletcherBase<Real> {
62 private:
63  // Required for Fletcher penalty function definition
66 
68 
69  // Evaluation counters
73 
74  using FletcherBase<Real>::fPhi_; // value of penalty function
75  using FletcherBase<Real>::gPhi_; // gradient of penalty function
76 
77  using FletcherBase<Real>::y_; // multiplier estimate
78 
79  using FletcherBase<Real>::fval_; // value of objective function
80  using FletcherBase<Real>::g_; // gradient of objective value
81  using FletcherBase<Real>::c_; // constraint value
82  using FletcherBase<Real>::scaledc_; // penaltyParameter_ * c_
83  using FletcherBase<Real>::gL_; // gradient of Lagrangian (g - A*y)
84 
91 
92  using FletcherBase<Real>::delta_; // regularization parameter
93 
95 
96  // Temporaries
97  Ptr<Vector<Real> > Tv_; // Temporary for matvecs
98  Ptr<Vector<Real> > w_; // first component of augmented system solve solution
99  Ptr<Vector<Real> > v_; // second component of augmented system solve solution
100 
101  Ptr<Vector<Real> > xzeros_; // zero vector
102  Ptr<Vector<Real> > czeros_; // zero vector
103 
105 
106  using FletcherBase<Real>::multSolverError_; // Error from augmented system solve in value()
107  using FletcherBase<Real>::gradSolveError_; // Error from augmented system solve in gradient()
108 
109  // For Augmented system solves
119 
120  class AugSystem : public LinearOperator<Real> {
121  private:
122  const Ptr<Constraint<Real> > con_;
123  const Ptr<const Vector<Real> > x_;
124  const Real delta_;
125  public:
126  AugSystem(const Ptr<Constraint<Real> > &con,
127  const Ptr<const Vector<Real> > &x,
128  const Real delta) : con_(con), x_(x), delta_(delta) {}
129 
130  void apply(Vector<Real> &Hv, const Vector<Real> &v, Real &tol) const {
131  PartitionedVector<Real> &Hvp = dynamic_cast<PartitionedVector<Real>&>(Hv);
132  const PartitionedVector<Real> &vp = dynamic_cast<const PartitionedVector<Real>&>(v);
133 
134  con_->applyAdjointJacobian(*(Hvp.get(0)), *(vp.get(1)), *x_, tol);
135  Hvp.get(0)->plus(*(vp.get(0)));
136 
137  con_->applyJacobian(*(Hvp.get(1)), *(vp.get(0)), *x_, tol);
138  Hvp.get(1)->axpy(-delta_*delta_, *(vp.get(1)));
139  }
140  };
141 
142  class AugSystemPrecond : public LinearOperator<Real> {
143  private:
144  const Ptr<Constraint<Real> > con_;
145  const Ptr<const Vector<Real> > x_;
146  public:
148  const Ptr<const Vector<Real> > x) : con_(con), x_(x) {}
149 
150  void apply(Vector<Real> &Hv, const Vector<Real> &v, Real &tol) const {
151  Hv.set(v.dual());
152  }
153  void applyInverse(Vector<Real> &Hv, const Vector<Real> &v, Real &tol) const {
154  Real zero(0);
155  PartitionedVector<Real> &Hvp = dynamic_cast<PartitionedVector<Real>&>(Hv);
156  const PartitionedVector<Real> &vp = dynamic_cast<const PartitionedVector<Real>&>(v);
157 
158  Hvp.set(0, *(vp.get(0)));
159  // Second x should be dual, but unused?
160  con_->applyPreconditioner(*(Hvp.get(1)),*(vp.get(1)),*x_,*x_, zero);
161  }
162  };
163 
164 public:
165  Fletcher(const ROL::Ptr<Objective<Real> > &obj,
166  const ROL::Ptr<Constraint<Real> > &con,
167  const Vector<Real> &optVec,
168  const Vector<Real> &conVec,
169  ROL::ParameterList &parlist)
170  : FletcherBase<Real>(obj, con) {
171 
172  gPhi_ = optVec.dual().clone();
173  y_ = conVec.dual().clone();
174  g_ = optVec.dual().clone();
175  gL_ = optVec.dual().clone();
176  c_ = conVec.clone();
177  scaledc_ = conVec.clone();
178 
179  Tv_ = optVec.dual().clone();
180  w_ = optVec.dual().clone();
181  v_ = conVec.dual().clone();
182 
183  xzeros_ = optVec.dual().clone();
184  xzeros_->zero();
185  czeros_ = conVec.clone();
186  czeros_->zero();
187 
188  v1_ = optVec.dual().clone();
189  v2_ = conVec.dual().clone();
190  vv_ = makePtr<PartitionedVector<Real>>(std::vector<Ptr<Vector<Real>> >({v1_, v2_}));
191 
192  b1_ = optVec.dual().clone();
193  b2_ = conVec.clone();
194  bb_ = makePtr<PartitionedVector<Real>>(std::vector<Ptr<Vector<Real>> >({b1_, b2_}));
195 
196  ROL::ParameterList& sublist = parlist.sublist("Step").sublist("Fletcher");
197  HessianApprox_ = sublist.get("Level of Hessian Approximation", 0);
198  penaltyParameter_ = sublist.get("Penalty Parameter", 1.0);
199 
200  delta_ = sublist.get("Regularization Parameter", 0.0);
201 
202  useInexact_ = sublist.get("Inexact Solves", false);
203 
204  ROL::ParameterList krylovList;
205  Real atol = static_cast<Real>(1e-12);
206  Real rtol = static_cast<Real>(1e-2);
207  krylovList.sublist("General").sublist("Krylov").set("Type", "GMRES");
208  krylovList.sublist("General").sublist("Krylov").set("Absolute Tolerance", atol);
209  krylovList.sublist("General").sublist("Krylov").set("Relative Tolerance", rtol);
210  krylovList.sublist("General").sublist("Krylov").set("Iteration Limit", 200);
211  krylov_ = KrylovFactory<Real>(krylovList);
212  }
213 
214  void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
215  obj_->update(x,flag,iter);
216  con_->update(x,flag,iter);
217  isValueComputed_ = (flag ? false : isValueComputed_);
218  isGradientComputed_ = (flag ? false : isGradientComputed_);
219  isMultiplierComputed_ = (flag ? false : isMultiplierComputed_);
220  isObjValueComputed_ = (flag ? false : isObjValueComputed_);
221  isObjGradComputed_ = (flag ? false : isObjGradComputed_);
222  isConValueComputed_ = (flag ? false : isConValueComputed_);
223  }
224 
225  Real value( const Vector<Real> &x, Real &tol ) {
226  if( isValueComputed_ )
227  return fPhi_;
228 
229  // Reset tolerances
230  Real origTol = tol;
231  Real tol2 = origTol;
232 
233  FletcherBase<Real>::objValue(x, tol2); tol2 = origTol;
234  multSolverError_ = origTol / static_cast<Real>(2);
236  tol = multSolverError_;
237 
238  fPhi_ = fval_ - c_->dot(y_->dual());
239  isValueComputed_ = true;
240 
241  return fPhi_;
242  }
243 
244  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
245  if( isGradientComputed_ && gradSolveError_ <= tol) {
246  tol = gradSolveError_;
247  g.set(*gPhi_);
248  return;
249  }
250 
251  // Reset tolerances
252  Real origTol = tol;
253  Real tol2 = origTol;
254 
255  gradSolveError_ = origTol / static_cast<Real>(2);
257 
258  // gPhi = sum y_i H_i w + sigma w + sum v_i H_i gL - H w + gL
261  tol = gradSolveError_;
262 
263  con_->applyAdjointHessian( *gPhi_, *y_, *w_, x, tol2 ); tol2 = origTol;
264  gPhi_->axpy( penaltyParameter_, *w_ );
265 
266  obj_->hessVec( *Tv_, *w_, x, tol2 ); tol2 = origTol;
267  gPhi_->axpy( static_cast<Real>(-1), *Tv_ );
268 
269  con_->applyAdjointHessian( *Tv_, *v_, *gL_, x, tol2 ); tol2 = origTol;
270  gPhi_->plus( *Tv_ );
271 
272  gPhi_->plus( *gL_ );
273 
274  g.set(*gPhi_);
275 
276  isGradientComputed_ = true;
277  }
278 
279  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
280  // Reset tolerances
281  Real origTol = tol;
282  Real tol2 = origTol;
283 
284  computeMultipliers(x, tol);
285 
286  obj_->hessVec( hv, v, x, tol2 ); tol2 = origTol;
287  con_->applyAdjointHessian( *Tv_, *y_, v, x, tol2 ); tol2 = origTol;
288  hv.axpy(static_cast<Real>(-1), *Tv_ );
289 
290  tol2 = tol;
291  solveAugmentedSystem( *w_, *v_, hv, *czeros_, x, tol2 ); tol2 = origTol;
292  hv.scale( static_cast<Real>(-1) );
293  hv.plus( *w_ );
294 
295  Tv_->set(v);
296  tol2 = tol;
297  solveAugmentedSystem( *w_, *v_, *Tv_, *czeros_, x, tol2 ); tol2 = origTol;
298  hv.axpy(static_cast<Real>(-2)*penaltyParameter_, *w_);
299 
300  obj_->hessVec( *Tv_, *w_, x, tol2 ); tol2 = origTol;
301  hv.plus( *Tv_ );
302  con_->applyAdjointHessian( *Tv_, *y_, *w_, x, tol2 ); tol2 = origTol;
303  hv.axpy( static_cast<Real>(-1), *Tv_ );
304 
305  hv.axpy(static_cast<Real>(2)*penaltyParameter_, v);
306  }
307 
309  Vector<Real> &v2,
310  const Vector<Real> &b1,
311  const Vector<Real> &b2,
312  const Vector<Real> &x,
313  Real &tol) {
314  // Ignore tol for now
315  ROL::Ptr<LinearOperator<Real> > K
316  = ROL::makePtr<AugSystem>(con_, makePtrFromRef(x), delta_);
317  ROL::Ptr<LinearOperator<Real> > P
318  = ROL::makePtr<AugSystemPrecond>(con_, makePtrFromRef(x));
319 
320  v1_->set(v1);
321  v2_->set(v2);
322 
323  b1_->set(b1);
324  b2_->set(b2);
325 
326  // If inexact, change tolerance
327  if( useInexact_ ) {
328  krylov_->resetAbsoluteTolerance(tol);
329  }
330 
331  flagKrylov_ = 0;
332  tol = krylov_->run(*vv_,*K,*bb_,*P,iterKrylov_,flagKrylov_);
333  v1.set(*v1_);
334  v2.set(*v2_);
335  }
336 
337  void computeMultipliers(const Vector<Real>& x, const Real tol) {
338  if( isMultiplierComputed_ && multSolverError_ <= tol) {
339  return;
340  }
341  Real tol2 = tol;
342  FletcherBase<Real>::objGrad(x, tol2); tol2 = tol;
344 
345  multSolverError_ = tol;
347 
348  isMultiplierComputed_ = true;
349  }
350 
351 }; // class Fletcher
352 
353 } // namespace ROL
354 
355 #endif
Provides the interface to evaluate objective functions.
Real value(const Vector< Real > &x, Real &tol)
Compute value.
Ptr< Vector< Real > > c_
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: ROL_Vector.hpp:226
Ptr< Vector< Real > > gPhi_
virtual void scale(const Real alpha)=0
Compute where .
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
const Ptr< Constraint< Real > > con_
ROL::Ptr< const Vector< Real > > get(size_type i) const
virtual void plus(const Vector &x)=0
Compute , where .
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:153
Defines the linear algebra of vector space on a generic partitioned vector.
void objValue(const Vector< Real > &x, Real &tol)
const Ptr< Objective< Real > > obj_
void applyInverse(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply inverse of linear operator.
Contains definitions of custom data types in ROL.
Fletcher(const ROL::Ptr< Objective< Real > > &obj, const ROL::Ptr< Constraint< Real > > &con, const Vector< Real > &optVec, const Vector< Real > &conVec, ROL::ParameterList &parlist)
const Ptr< const Vector< Real > > x_
const Ptr< const Vector< Real > > x_
AugSystemPrecond(const Ptr< Constraint< Real > > con, const Ptr< const Vector< Real > > x)
Ptr< Vector< Real > > czeros_
Ptr< Vector< Real > > Tv_
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:80
Ptr< Krylov< Real > > krylov_
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Ptr< Vector< Real > > gL_
Ptr< Vector< Real > > v1_
Ptr< Vector< Real > > b1_
AugSystem(const Ptr< Constraint< Real > > &con, const Ptr< const Vector< Real > > &x, const Real delta)
Ptr< Vector< Real > > w_
Ptr< Vector< Real > > scaledc_
void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
Ptr< Vector< Real > > g_
void conValue(const Vector< Real > &x, Real &tol)
void set(const V &x)
Set where .
void computeMultipliers(const Vector< Real > &x, const Real tol)
Ptr< Vector< Real > > y_
Ptr< Vector< Real > > xzeros_
Ptr< PartitionedVector< Real > > vv_
Provides the interface to apply a linear operator.
const Ptr< Constraint< Real > > con_
void solveAugmentedSystem(Vector< Real > &v1, Vector< Real > &v2, const Vector< Real > &b1, const Vector< Real > &b2, const Vector< Real > &x, Real &tol)
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
Ptr< Vector< Real > > v2_
Ptr< Vector< Real > > b2_
Ptr< PartitionedVector< Real > > bb_
const Ptr< Constraint< Real > > con_
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
void objGrad(const Vector< Real > &x, Real &tol)
Defines the general constraint operator interface.
void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
Ptr< Vector< Real > > v_