ROL
ROL_ObjectiveDef.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#ifndef ROL_OBJECTIVE_DEF_H
45#define ROL_OBJECTIVE_DEF_H
46
51namespace ROL {
52
53template<typename Real>
54Real Objective<Real>::dirDeriv( const Vector<Real> &x, const Vector<Real> &d, Real &tol) {
55 if (dual_ == nullPtr) dual_ = x.dual().clone();
56 gradient(*dual_,x,tol);
57 //return d.dot(dual_->dual());
58 return d.apply(*dual_);
59 //Real dnorm = d.norm(), zero(0);
60 //if ( dnorm == zero ) {
61 // return zero;
62 //}
63 //Real cbrteps = std::cbrt(ROL_EPSILON<Real>()), one(1), v0(0), v1(0);
64 //Real xnorm = x.norm(), h = cbrteps * std::max(xnorm/dnorm,one);
65 //v0 = value(x,tol);
66 //prim_->set(x); prim_->axpy(h, d);
67 //update(*prim_,UpdateType::Temp);
68 //v1 = value(*prim_,tol);
69 //update(x,UpdateType::Revert);
70 //return (v1 - v0) / h;
71}
72
73template<typename Real>
74void Objective<Real>::gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
75 if (prim_ == nullPtr) prim_ = x.clone();
76 if (basis_ == nullPtr) basis_ = x.clone();
77
78 const Real cbrteps = std::cbrt(ROL_EPSILON<Real>()), zero(0), one(1);
79 Real f0 = value(x,tol), h(0), xi(0), gi(0);
80 g.zero();
81 for (int i = 0; i < x.dimension(); i++) {
82 basis_->set(*x.basis(i));
83 xi = x.dot(*basis_);
84 h = cbrteps * std::max(std::abs(xi),one) * (xi < zero ? -one : one);
85 prim_->set(x); prim_->axpy(h,*basis_);
86 h = prim_->dot(*basis_) - xi;
88 gi = (value(*prim_,tol) - f0) / h;
89 g.axpy(gi,*g.basis(i));
90 }
92}
93
94template<typename Real>
95void Objective<Real>::hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
96 const Real zero(0), vnorm = v.norm();
97 // Get Step Length
98 if ( vnorm == zero ) {
99 hv.zero();
100 }
101 else {
102 if (prim_ == nullPtr) prim_ = x.clone();
103 if (dual_ == nullPtr) dual_ = hv.clone();
104
105 //Real h = 2.0/(v.norm()*v.norm())*tol;
106 const Real one(1), h(std::max(one,x.norm()/vnorm)*tol);
107
108 gradient(*dual_,x,tol); // Compute gradient at x
109 prim_->set(x); prim_->axpy(h,v); // Set prim = x + hv
110 update(*prim_,UpdateType::Temp); // Temporarily update objective at x + hv
111 gradient(hv,*prim_,tol); // Compute gradient at x + hv
112 hv.axpy(-one,*dual_); // Compute difference (f'(x+hv)-f'(x))
113 hv.scale(one/h); // Compute Newton quotient (f'(x+hv)-f'(x))/h
114 update(x,UpdateType::Revert); // Reset objective to x
115 }
116}
117
118template<typename Real>
119std::vector<std::vector<Real>> Objective<Real>::checkGradient( const Vector<Real> &x,
120 const Vector<Real> &g,
121 const Vector<Real> &d,
122 const bool printToStream,
123 std::ostream & outStream,
124 const int numSteps,
125 const int order ) {
126
127 const Real ten(10);
128 std::vector<Real> steps(numSteps);
129 for(int i=0;i<numSteps;++i) {
130 steps[i] = pow(ten,static_cast<Real>(-i));
131 }
132
133 return checkGradient(x,g,d,steps,printToStream,outStream,order);
134
135} // checkGradient
136
137template<typename Real>
138std::vector<std::vector<Real>> Objective<Real>::checkGradient( const Vector<Real> &x,
139 const Vector<Real> &g,
140 const Vector<Real> &d,
141 const std::vector<Real> &steps,
142 const bool printToStream,
143 std::ostream & outStream,
144 const int order ) {
145
146 ROL_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
147 "Error: finite difference order must be 1,2,3, or 4" );
148
151
152 Real tol = std::sqrt(ROL_EPSILON<Real>());
153
154 int numSteps = steps.size();
155 int numVals = 4;
156 std::vector<Real> tmp(numVals);
157 std::vector<std::vector<Real>> gCheck(numSteps, tmp);
158
159 // Save the format state of the original outStream.
160 nullstream oldFormatState;
161 oldFormatState.copyfmt(outStream);
162
163 // Evaluate objective value at x.
165 Real val = value(x,tol);
166
167 // Compute gradient at x.
168 Ptr<Vector<Real>> gtmp = g.clone();
169 gradient(*gtmp, x, tol);
170 //Real dtg = d.dot(gtmp->dual());
171 Real dtg = d.apply(*gtmp);
172
173 // Temporary vectors.
174 Ptr<Vector<Real>> xnew = x.clone();
175
176 for (int i=0; i<numSteps; i++) {
177
178 Real eta = steps[i];
179
180 xnew->set(x);
181
182 // Compute gradient, finite-difference gradient, and absolute error.
183 gCheck[i][0] = eta;
184 gCheck[i][1] = dtg;
185
186 gCheck[i][2] = weights[order-1][0] * val;
187
188 for(int j=0; j<order; ++j) {
189 // Evaluate at x <- x+eta*c_i*d.
190 xnew->axpy(eta*shifts[order-1][j], d);
191
192 // Only evaluate at shifts where the weight is nonzero
193 if( weights[order-1][j+1] != 0 ) {
195 gCheck[i][2] += weights[order-1][j+1] * this->value(*xnew,tol);
196 }
197 }
198
199 gCheck[i][2] /= eta;
200
201 gCheck[i][3] = std::abs(gCheck[i][2] - gCheck[i][1]);
202
203 if (printToStream) {
204 if (i==0) {
205 outStream << std::right
206 << std::setw(20) << "Step size"
207 << std::setw(20) << "grad'*dir"
208 << std::setw(20) << "FD approx"
209 << std::setw(20) << "abs error"
210 << "\n"
211 << std::setw(20) << "---------"
212 << std::setw(20) << "---------"
213 << std::setw(20) << "---------"
214 << std::setw(20) << "---------"
215 << "\n";
216 }
217 outStream << std::scientific << std::setprecision(11) << std::right
218 << std::setw(20) << gCheck[i][0]
219 << std::setw(20) << gCheck[i][1]
220 << std::setw(20) << gCheck[i][2]
221 << std::setw(20) << gCheck[i][3]
222 << "\n";
223 }
224
225 }
226
227 // Reset format state of outStream.
228 outStream.copyfmt(oldFormatState);
229
230 return gCheck;
231} // checkGradient
232
233template<typename Real>
234std::vector<std::vector<Real>> Objective<Real>::checkHessVec( const Vector<Real> &x,
235 const Vector<Real> &hv,
236 const Vector<Real> &v,
237 const bool printToStream,
238 std::ostream & outStream,
239 const int numSteps,
240 const int order ) {
241 const Real ten(10);
242 std::vector<Real> steps(numSteps);
243 for(int i=0;i<numSteps;++i) {
244 steps[i] = pow(ten,static_cast<Real>(-i));
245 }
246
247 return checkHessVec(x,hv,v,steps,printToStream,outStream,order);
248} // checkHessVec
249
250
251
252template<typename Real>
253std::vector<std::vector<Real>> Objective<Real>::checkHessVec( const Vector<Real> &x,
254 const Vector<Real> &hv,
255 const Vector<Real> &v,
256 const std::vector<Real> &steps,
257 const bool printToStream,
258 std::ostream & outStream,
259 const int order ) {
260
261 ROL_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
262 "Error: finite difference order must be 1,2,3, or 4" );
263
266
267 const Real one(1);
268 Real tol = std::sqrt(ROL_EPSILON<Real>());
269
270 int numSteps = steps.size();
271 int numVals = 4;
272 std::vector<Real> tmp(numVals);
273 std::vector<std::vector<Real>> hvCheck(numSteps, tmp);
274
275 // Save the format state of the original outStream.
276 nullstream oldFormatState;
277 oldFormatState.copyfmt(outStream);
278
279 // Compute gradient at x.
280 Ptr<Vector<Real>> g = hv.clone();
282 gradient(*g, x, tol);
283
284 // Compute (Hessian at x) times (vector v).
285 Ptr<Vector<Real>> Hv = hv.clone();
286 hessVec(*Hv, v, x, tol);
287 Real normHv = Hv->norm();
288
289 // Temporary vectors.
290 Ptr<Vector<Real>> gdif = hv.clone();
291 Ptr<Vector<Real>> gnew = hv.clone();
292 Ptr<Vector<Real>> xnew = x.clone();
293
294 for (int i=0; i<numSteps; i++) {
295 Real eta = steps[i];
296 // Evaluate objective value at x+eta*d.
297 xnew->set(x);
298 gdif->set(*g);
299 gdif->scale(weights[order-1][0]);
300 for (int j=0; j<order; ++j) {
301 // Evaluate at x <- x+eta*c_i*d.
302 xnew->axpy(eta*shifts[order-1][j], v);
303 // Only evaluate at shifts where the weight is nonzero
304 if ( weights[order-1][j+1] != 0 ) {
306 gradient(*gnew, *xnew, tol);
307 gdif->axpy(weights[order-1][j+1],*gnew);
308 }
309 }
310 gdif->scale(one/eta);
311
312 // Compute norms of hessvec, finite-difference hessvec, and error.
313 hvCheck[i][0] = eta;
314 hvCheck[i][1] = normHv;
315 hvCheck[i][2] = gdif->norm();
316 gdif->axpy(-one, *Hv);
317 hvCheck[i][3] = gdif->norm();
318
319 if (printToStream) {
320 if (i==0) {
321 outStream << std::right
322 << std::setw(20) << "Step size"
323 << std::setw(20) << "norm(Hess*vec)"
324 << std::setw(20) << "norm(FD approx)"
325 << std::setw(20) << "norm(abs error)"
326 << "\n"
327 << std::setw(20) << "---------"
328 << std::setw(20) << "--------------"
329 << std::setw(20) << "---------------"
330 << std::setw(20) << "---------------"
331 << "\n";
332 }
333 outStream << std::scientific << std::setprecision(11) << std::right
334 << std::setw(20) << hvCheck[i][0]
335 << std::setw(20) << hvCheck[i][1]
336 << std::setw(20) << hvCheck[i][2]
337 << std::setw(20) << hvCheck[i][3]
338 << "\n";
339 }
340
341 }
342
343 // Reset format state of outStream.
344 outStream.copyfmt(oldFormatState);
345
346 return hvCheck;
347} // checkHessVec
348
349template<typename Real>
350std::vector<Real> Objective<Real>::checkHessSym( const Vector<Real> &x,
351 const Vector<Real> &hv,
352 const Vector<Real> &v,
353 const Vector<Real> &w,
354 const bool printToStream,
355 std::ostream & outStream ) {
356
357 Real tol = std::sqrt(ROL_EPSILON<Real>());
358
359 // Compute (Hessian at x) times (vector v).
360 Ptr<Vector<Real>> h = hv.clone();
362 hessVec(*h, v, x, tol);
363 //Real wHv = w.dot(h->dual());
364 Real wHv = w.apply(*h);
365
366 hessVec(*h, w, x, tol);
367 //Real vHw = v.dot(h->dual());
368 Real vHw = v.apply(*h);
369
370 std::vector<Real> hsymCheck(3, 0);
371
372 hsymCheck[0] = wHv;
373 hsymCheck[1] = vHw;
374 hsymCheck[2] = std::abs(vHw-wHv);
375
376 // Save the format state of the original outStream.
377 nullstream oldFormatState;
378 oldFormatState.copyfmt(outStream);
379
380 if (printToStream) {
381 outStream << std::right
382 << std::setw(20) << "<w, H(x)v>"
383 << std::setw(20) << "<v, H(x)w>"
384 << std::setw(20) << "abs error"
385 << "\n";
386 outStream << std::scientific << std::setprecision(11) << std::right
387 << std::setw(20) << hsymCheck[0]
388 << std::setw(20) << hsymCheck[1]
389 << std::setw(20) << hsymCheck[2]
390 << "\n";
391 }
392
393 // Reset format state of outStream.
394 outStream.copyfmt(oldFormatState);
395
396 return hsymCheck;
397
398} // checkHessSym
399
400} // namespace ROL
401
402#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
virtual std::vector< Real > checkHessSym(const Vector< Real > &x, const Vector< Real > &v, const Vector< Real > &w, const bool printToStream=true, std::ostream &outStream=std::cout)
Hessian symmetry check.
virtual Real dirDeriv(const Vector< Real > &x, const Vector< Real > &d, Real &tol)
Compute directional derivative.
virtual std::vector< std::vector< Real > > checkGradient(const Vector< Real > &x, const Vector< Real > &d, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Finite-difference gradient check.
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
virtual std::vector< std::vector< Real > > checkHessVec(const Vector< Real > &x, const Vector< Real > &v, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Finite-difference Hessian-applied-to-vector check.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
virtual Real apply(const Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
Definition: ROL_Vector.hpp:238
virtual Real norm() const =0
Returns 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,...
Definition: ROL_Vector.hpp:226
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.
virtual int dimension() const
Return dimension of the vector space.
Definition: ROL_Vector.hpp:196
virtual ROL::Ptr< Vector > basis(const int i) const
Return i-th basis vector.
Definition: ROL_Vector.hpp:182
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:153
virtual Real dot(const Vector &x) const =0
Compute where .
const double weights[4][5]
Definition: ROL_Types.hpp:868
ROL::Objective_SerialSimOpt Objective_SimOpt value(const V &u, const V &z, Real &tol) override
virtual void update(const Vector< Real > &u, const Vector< Real > &z, bool flag=true, int iter=-1) override