ROL
ROL_FletcherObjectiveE_Def.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_FLETCHEROBJECTIVEEDEF_H
46#define ROL_FLETCHEROBJECTIVEEDEF_H
47
48namespace ROL {
49
50template<typename Real>
52 const ROL::Ptr<Constraint<Real>> &con,
53 const Vector<Real> &xprim,
54 const Vector<Real> &xdual,
55 const Vector<Real> &cprim,
56 const Vector<Real> &cdual,
57 ROL::ParameterList &parlist)
58 : FletcherObjectiveBase<Real>(obj, con, xprim, xdual, cprim, cdual, parlist) {
59 Tv_ = xdual.clone();
60 w_ = xdual.clone();
61 wdual_ = xprim.clone();
62 v_ = cdual.clone();
63 wg_ = xdual.clone();
64 vg_ = cdual.clone();
65
66 xzeros_ = xdual.clone(); xzeros_->zero();
67 czeros_ = cprim.clone(); czeros_->zero();
68}
69
70template<typename Real>
71Real FletcherObjectiveE<Real>::value( const Vector<Real> &x, Real &tol ) {
72 Real val(0);
73 int key(0);
74 bool isComputed = fPhi_->get(val,key);
75 if( isComputed && multSolverError_*cnorm_ <= tol) {
76 tol = multSolverError_*cnorm_;
77 }
78 else {
79 // Reset tolerances
80 Real origTol = tol;
81 Real tol2 = origTol;
82 // Compute penalty function value
83 Real fval = FletcherObjectiveBase<Real>::objValue(x, tol2); tol2 = origTol;
84 multSolverError_ = origTol / (static_cast<Real>(2) * std::max(static_cast<Real>(1), cnorm_));
85 FletcherObjectiveBase<Real>::computeMultipliers(*cdual_, *gLdual_, x, *xdual_, *cprim_, multSolverError_);
86 tol = multSolverError_*cnorm_;
87 //val = fval - cprim_->dot(cdual_->dual());
88 val = fval - cprim_->apply(*cdual_);
89 if( quadPenaltyParameter_ > static_cast<Real>(0) )
90 val += static_cast<Real>(0.5)*quadPenaltyParameter_*(cprim_->dot(*cprim_));
91 // Store new penalty function value
92 fPhi_->set(val,key);
93 }
94 return val;
95}
96
97template<typename Real>
99 int key(0);
100 bool isComputed = gPhi_->get(g,key);
101 if( isComputed && gradSolveError_ <= tol) {
102 tol = gradSolveError_;
103 }
104 else {
105 // Reset tolerances
106 Real origTol = tol;
107 Real tol2 = origTol;
108 // Compute penalty function gradient
109 gradSolveError_ = origTol / static_cast<Real>(2);
110 FletcherObjectiveBase<Real>::computeMultipliers(*cdual_, *gLdual_, x, *xdual_, *cprim_, gradSolveError_);
111 gL_->set(gLdual_->dual());
112 bool refine = isComputed;
113 // gPhi = sum y_i H_i w + sigma w + sum v_i H_i gL - H w + gL
114 solveAugmentedSystem( *wdual_, *vg_, *xzeros_, *cprim_, x, gradSolveError_, refine );
115 gradSolveError_ += multSolverError_;
116 tol = gradSolveError_;
117 wg_->set(wdual_->dual());
118 con_->applyAdjointHessian( g, *cdual_, *wdual_, x, tol2 ); tol2 = origTol;
119 g.axpy( sigma_, *wg_ );
120 obj_->hessVec( *Tv_, *wdual_, x, tol2 ); tol2 = origTol;
121 g.axpy( static_cast<Real>(-1), *Tv_ );
122 con_->applyAdjointHessian( *Tv_, *vg_, *gLdual_, x, tol2 ); tol2 = origTol;
123 g.plus( *Tv_ );
124 g.plus( *gL_ );
125 if( quadPenaltyParameter_ > static_cast<Real>(0) ) {
126 con_->applyAdjointJacobian( *Tv_, *cprim_, x, tol2 ); tol2 = origTol;
127 g.axpy( quadPenaltyParameter_, *Tv_ );
128 }
129 gPhi_->set(g,key);
130 }
131}
132
133template<typename Real>
135 // Reset tolerances
136 Real origTol = tol;
137 Real tol2 = origTol;
138 int key(0);
139 bool isComputed = y_->get(*cdual_,key);
140 if( !isComputed || !useInexact_) {
141 // hessVec tol is always set to ~1e-8. So if inexact linear system solves are used, then
142 // the multipliers will always get re-evaluated to high precision, which defeats the purpose
143 // of computing the objective/gradient approximately.
144 // Thus if inexact linear system solves are used, we will not update the multiplier estimates
145 // to high precision.
146 FletcherObjectiveBase<Real>::computeMultipliers(*cdual_, *gLdual_, x, *xdual_, *cprim_, tol);
147 }
148
149 obj_->hessVec( hv, v, x, tol2 ); tol2 = origTol;
150 con_->applyAdjointHessian( *Tv_, *cdual_, v, x, tol2 ); tol2 = origTol;
151 hv.axpy(static_cast<Real>(-1), *Tv_ );
152
153 tol2 = tol;
154 solveAugmentedSystem( *w_, *v_, hv, *czeros_, x, tol2 ); tol2 = origTol;
155 hv.scale( static_cast<Real>(-1) );
156 hv.plus( *w_ );
157
158 Tv_->set(v.dual());
159 tol2 = tol;
160 solveAugmentedSystem( *w_, *v_, *Tv_, *czeros_, x, tol2 ); tol2 = origTol;
161 hv.axpy(static_cast<Real>(-2)*sigma_, *w_);
162
163 wdual_->set(w_->dual());
164
165 obj_->hessVec( *Tv_, *wdual_, x, tol2 ); tol2 = origTol;
166 hv.plus( *Tv_ );
167 con_->applyAdjointHessian( *Tv_, *cdual_, *wdual_, x, tol2 ); tol2 = origTol;
168 hv.axpy( static_cast<Real>(-1), *Tv_ );
169
170 hv.axpy( static_cast<Real>(2)*sigma_, v );
171
172 if( quadPenaltyParameter_ > static_cast<Real>(0) ) {
173 con_->applyJacobian( *b2_, v, x, tol2 ); tol2 = origTol;
174 con_->applyAdjointJacobian( *Tv_, *b2_, x, tol2 ); tol2 = origTol;
175 hv.axpy( quadPenaltyParameter_, *Tv_ );
176 con_->applyAdjointHessian( *Tv_, *cprim_, v, x, tol2); tol2 = origTol;
177 hv.axpy( -quadPenaltyParameter_, *Tv_ );
178 }
179}
180
181template<typename Real>
183 Vector<Real> &v2,
184 const Vector<Real> &b1,
185 const Vector<Real> &b2,
186 const Vector<Real> &x,
187 Real &tol,
188 bool refine) {
189 // Ignore tol for now
190 ROL::Ptr<LinearOperator<Real>>
191 K = ROL::makePtr<AugSystem>(con_, makePtrFromRef(x), delta_);
192 ROL::Ptr<LinearOperator<Real>>
193 P = ROL::makePtr<AugSystemPrecond>(con_, makePtrFromRef(x), makePtrFromRef(b1));
194
195 vv_->zero();
196 bb_->zero();
197 if( refine ) {
198 // TODO: Make sure this tol is actually ok...
199 Real origTol = tol;
200 w1_->set(v1);
201 w2_->set(v2);
202 K->apply(*bb_, *ww_, tol); tol = origTol;
203 bb_->scale(static_cast<Real>(-1));
204 }
205 b1_->plus(b1);
206 b2_->plus(b2);
207
208 // If inexact, change tolerance
209 if( useInexact_ ) krylov_->resetAbsoluteTolerance(tol);
210
211 //con_->solveAugmentedSystem(*v1_,*v2_,*b1_,*b2_,x,tol);
212 flagKrylov_ = 0;
213 tol = krylov_->run(*vv_,*K,*bb_,*P,iterKrylov_,flagKrylov_);
214
215 if( refine ) {
216 v1.plus(*v1_);
217 v2.plus(*v2_);
218 } else {
219 v1.set(*v1_);
220 v2.set(*v2_);
221 }
222}
223
224} // namespace ROL
225
226#endif
const Ptr< Obj > obj_
Defines the general constraint operator interface.
Real objValue(const Vector< Real > &x, Real &tol)
void computeMultipliers(Vector< Real > &y, Vector< Real > &gL, const Vector< Real > &x, Vector< Real > &g, Vector< Real > &c, Real tol)
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) override
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol) override
Apply Hessian approximation to vector.
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol) override
Compute gradient.
Real value(const Vector< Real > &x, Real &tol) override
Compute value.
FletcherObjectiveE(const ROL::Ptr< Objective< Real > > &obj, const ROL::Ptr< Constraint< Real > > &con, const Vector< Real > &xprim, const Vector< Real > &xdual, const Vector< Real > &cprim, const Vector< Real > &cdual, ROL::ParameterList &parlist)
Provides the interface to evaluate objective functions.
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 .