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 
17 namespace ROL {
18 
19 template<typename Real>
20 Real 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 
39 template<typename Real>
40 void 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 
60 template<typename Real>
61 void 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 
84 template<typename Real>
85 void 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 
106 template<typename Real>
107 std::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 
125 template<typename Real>
126 std::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 
221 template<typename Real>
222 std::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 
240 template<typename Real>
241 std::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 
337 template<typename Real>
338 std::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 
388 template<typename Real>
389 std::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
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::scale
virtual void scale(const Real alpha)=0
Compute where .
ROL::Vector::clone
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
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::Vector::dimension
virtual int dimension() const
Return dimension of the vector space.
Definition: ROL_Vector.hpp:162
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::basis
virtual ROL::Ptr< Vector > basis(const int i) const
Return i-th basis vector.
Definition: ROL_Vector.hpp:148
ROL::Finite_Difference_Arrays::weights
const double weights[4][5]
Definition: ROL_Types.hpp:838
ROL::Vector::axpy
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:119
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
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::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::Finite_Difference_Arrays::shifts
const int shifts[4][4]
Definition: ROL_Types.hpp:831
ROL::value
ROL::Objective_SimOpt value
ROL::Vector::zero
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:133
ROL::Vector
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
ROL::Vector::dot
virtual Real dot(const Vector &x) const =0
Compute where .
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
zero
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Definition: ROL_Objective_SerialSimOpt.hpp:111
ROL::details::nullstream
basic_nullstream< char, std::char_traits< char >> nullstream
Definition: ROL_Stream.hpp:36
ROL::Objective::gradient
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Definition: ROL_ObjectiveDef.hpp:40
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::norm
virtual Real norm() const =0
Returns where .
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
Generated by   doxygen 1.8.5