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