ROL
ROL_Objective_SimOpt.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_SIMOPT_H
11 #define ROL_OBJECTIVE_SIMOPT_H
12 
13 #include "ROL_Objective.hpp"
14 #include "ROL_Vector_SimOpt.hpp"
15 
22 namespace ROL {
23 
24 template <class Real>
25 class Objective_SimOpt : public Objective<Real> {
26 public:
27 
34  virtual void update( const Vector<Real> &u, const Vector<Real> &z, bool flag = true, int iter = -1 ) {}
35 
36  void update( const Vector<Real> &x, bool flag = true, int iter = -1 ) {
37  const ROL::Vector_SimOpt<Real> &xs = dynamic_cast<const ROL::Vector_SimOpt<Real>&>(
38  dynamic_cast<const ROL::Vector<Real>&>(x));
39  this->update(*(xs.get_1()),*(xs.get_2()),flag,iter);
40  }
41 
42  virtual void update( const Vector<Real> &u, const Vector<Real> &z, UpdateType type, int iter = -1 ) {}
43 
44  void update( const Vector<Real> &x, UpdateType type, int iter = -1 ) {
45  const ROL::Vector_SimOpt<Real> &xs = dynamic_cast<const ROL::Vector_SimOpt<Real>&>(
46  dynamic_cast<const ROL::Vector<Real>&>(x));
47  this->update(*(xs.get_1()),*(xs.get_2()),type,iter);
48  }
49 
50 
53  virtual Real value( const Vector<Real> &u, const Vector<Real> &z, Real &tol ) = 0;
54 
55  Real value( const Vector<Real> &x, Real &tol ) {
56  const ROL::Vector_SimOpt<Real> &xs = dynamic_cast<const ROL::Vector_SimOpt<Real>&>(
57  dynamic_cast<const ROL::Vector<Real>&>(x));
58  return this->value(*(xs.get_1()),*(xs.get_2()),tol);
59  }
60 
61 
64  virtual void gradient_1( Vector<Real> &g, const Vector<Real> &u, const Vector<Real> &z, Real &tol ) {
65  Real ftol = std::sqrt(ROL_EPSILON<Real>());
66  Real h = 0.0;
67  this->update(u,z,UpdateType::Temp);
68  Real v = this->value(u,z,ftol);
69  Real deriv = 0.0;
70  ROL::Ptr<Vector<Real> > unew = u.clone();
71  g.zero();
72  for (int i = 0; i < g.dimension(); i++) {
73  h = u.dot(*u.basis(i))*tol;
74  unew->set(u);
75  unew->axpy(h,*(u.basis(i)));
76  this->update(*unew,z,UpdateType::Temp);
77  deriv = (this->value(*unew,z,ftol) - v)/h;
78  g.axpy(deriv,*(g.basis(i)));
79  }
80  this->update(u,z,UpdateType::Temp);
81  }
84  virtual void gradient_2( Vector<Real> &g, const Vector<Real> &u, const Vector<Real> &z, Real &tol ) {
85  Real ftol = std::sqrt(ROL_EPSILON<Real>());
86  Real h = 0.0;
87  this->update(u,z,UpdateType::Temp);
88  Real v = this->value(u,z,ftol);
89  Real deriv = 0.0;
90  ROL::Ptr<Vector<Real> > znew = z.clone();
91  g.zero();
92  for (int i = 0; i < g.dimension(); i++) {
93  h = z.dot(*z.basis(i))*tol;
94  znew->set(z);
95  znew->axpy(h,*(z.basis(i)));
96  this->update(u,*znew,UpdateType::Temp);
97  deriv = (this->value(u,*znew,ftol) - v)/h;
98  g.axpy(deriv,*(g.basis(i)));
99  }
100  this->update(u,z,UpdateType::Temp);
101  }
102 
103  void gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
105  dynamic_cast<ROL::Vector<Real>&>(g));
106  const ROL::Vector_SimOpt<Real> &xs = dynamic_cast<const ROL::Vector_SimOpt<Real>&>(
107  dynamic_cast<const ROL::Vector<Real>&>(x));
108  ROL::Ptr<Vector<Real> > g1 = gs.get_1()->clone();
109  ROL::Ptr<Vector<Real> > g2 = gs.get_2()->clone();
110  this->gradient_1(*g1,*(xs.get_1()),*(xs.get_2()),tol);
111  this->gradient_2(*g2,*(xs.get_1()),*(xs.get_2()),tol);
112  gs.set_1(*g1);
113  gs.set_2(*g2);
114  }
115 
116 
119  virtual void hessVec_11( Vector<Real> &hv, const Vector<Real> &v,
120  const Vector<Real> &u, const Vector<Real> &z, Real &tol ) {
121  Real gtol = std::sqrt(ROL_EPSILON<Real>());
122  // Compute step length
123  Real h = tol;
124  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
125  h = std::max(1.0,u.norm()/v.norm())*tol;
126  }
127  // Evaluate gradient of first component at (u+hv,z)
128  ROL::Ptr<Vector<Real> > unew = u.clone();
129  unew->set(u);
130  unew->axpy(h,v);
131  this->update(*unew,z,UpdateType::Temp);
132  hv.zero();
133  this->gradient_1(hv,*unew,z,gtol);
134  // Evaluate gradient of first component at (u,z)
135  ROL::Ptr<Vector<Real> > g = hv.clone();
136  this->update(u,z,UpdateType::Temp);
137  this->gradient_1(*g,u,z,gtol);
138  // Compute Newton quotient
139  hv.axpy(-1.0,*g);
140  hv.scale(1.0/h);
141  }
142 
143  virtual void hessVec_12( Vector<Real> &hv, const Vector<Real> &v,
144  const Vector<Real> &u, const Vector<Real> &z, Real &tol ) {
145  Real gtol = std::sqrt(ROL_EPSILON<Real>());
146  // Compute step length
147  Real h = tol;
148  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
149  h = std::max(1.0,u.norm()/v.norm())*tol;
150  }
151  // Evaluate gradient of first component at (u,z+hv)
152  ROL::Ptr<Vector<Real> > znew = z.clone();
153  znew->set(z);
154  znew->axpy(h,v);
155  this->update(u,*znew,UpdateType::Temp);
156  hv.zero();
157  this->gradient_1(hv,u,*znew,gtol);
158  // Evaluate gradient of first component at (u,z)
159  ROL::Ptr<Vector<Real> > g = hv.clone();
160  this->update(u,z,UpdateType::Temp);
161  this->gradient_1(*g,u,z,gtol);
162  // Compute Newton quotient
163  hv.axpy(-1.0,*g);
164  hv.scale(1.0/h);
165  }
166 
167  virtual void hessVec_21( Vector<Real> &hv, const Vector<Real> &v,
168  const Vector<Real> &u, const Vector<Real> &z, Real &tol ) {
169  Real gtol = std::sqrt(ROL_EPSILON<Real>());
170  // Compute step length
171  Real h = tol;
172  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
173  h = std::max(1.0,u.norm()/v.norm())*tol;
174  }
175  // Evaluate gradient of first component at (u+hv,z)
176  ROL::Ptr<Vector<Real> > unew = u.clone();
177  unew->set(u);
178  unew->axpy(h,v);
179  this->update(*unew,z,UpdateType::Temp);
180  hv.zero();
181  this->gradient_2(hv,*unew,z,gtol);
182  // Evaluate gradient of first component at (u,z)
183  ROL::Ptr<Vector<Real> > g = hv.clone();
184  this->update(u,z,UpdateType::Temp);
185  this->gradient_2(*g,u,z,gtol);
186  // Compute Newton quotient
187  hv.axpy(-1.0,*g);
188  hv.scale(1.0/h);
189  }
190 
191  virtual void hessVec_22( Vector<Real> &hv, const Vector<Real> &v,
192  const Vector<Real> &u, const Vector<Real> &z, Real &tol ) {
193  Real gtol = std::sqrt(ROL_EPSILON<Real>());
194  // Compute step length
195  Real h = tol;
196  if (v.norm() > std::sqrt(ROL_EPSILON<Real>())) {
197  h = std::max(1.0,u.norm()/v.norm())*tol;
198  }
199  // Evaluate gradient of first component at (u,z+hv)
200  ROL::Ptr<Vector<Real> > znew = z.clone();
201  znew->set(z);
202  znew->axpy(h,v);
203  this->update(u,*znew,UpdateType::Temp);
204  hv.zero();
205  this->gradient_2(hv,u,*znew,gtol);
206  // Evaluate gradient of first component at (u,z)
207  ROL::Ptr<Vector<Real> > g = hv.clone();
208  this->update(u,z,UpdateType::Temp);
209  this->gradient_2(*g,u,z,gtol);
210  // Compute Newton quotient
211  hv.axpy(-1.0,*g);
212  hv.scale(1.0/h);
213  }
214 
215  void hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
217  dynamic_cast<ROL::Vector<Real>&>(hv));
218  const ROL::Vector_SimOpt<Real> &vs = dynamic_cast<const ROL::Vector_SimOpt<Real>&>(
219  dynamic_cast<const ROL::Vector<Real>&>(v));
220  const ROL::Vector_SimOpt<Real> &xs = dynamic_cast<const ROL::Vector_SimOpt<Real>&>(
221  dynamic_cast<const ROL::Vector<Real>&>(x));
222  ROL::Ptr<Vector<Real> > h11 = (hvs.get_1())->clone();
223  this->hessVec_11(*h11,*(vs.get_1()),*(xs.get_1()),*(xs.get_2()),tol);
224  ROL::Ptr<Vector<Real> > h12 = (hvs.get_1())->clone();
225  this->hessVec_12(*h12,*(vs.get_2()),*(xs.get_1()),*(xs.get_2()),tol);
226  ROL::Ptr<Vector<Real> > h21 = (hvs.get_2())->clone();
227  this->hessVec_21(*h21,*(vs.get_1()),*(xs.get_1()),*(xs.get_2()),tol);
228  ROL::Ptr<Vector<Real> > h22 = (hvs.get_2())->clone();
229  this->hessVec_22(*h22,*(vs.get_2()),*(xs.get_1()),*(xs.get_2()),tol);
230  h11->plus(*h12);
231  hvs.set_1(*h11);
232  h22->plus(*h21);
233  hvs.set_2(*h22);
234  }
235 
236  std::vector<std::vector<Real> > checkGradient_1( const Vector<Real> &u,
237  const Vector<Real> &z,
238  const Vector<Real> &d,
239  const bool printToStream = true,
240  std::ostream & outStream = std::cout,
241  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
242  const int order = 1 ) {
243  return checkGradient_1(u, z, u.dual(), d, printToStream, outStream, numSteps, order);
244  }
245 
246  std::vector<std::vector<Real> > checkGradient_1( const Vector<Real> &u,
247  const Vector<Real> &z,
248  const Vector<Real> &g,
249  const Vector<Real> &d,
250  const bool printToStream,
251  std::ostream & outStream,
252  const int numSteps,
253  const int order ) {
254  std::vector<Real> steps(numSteps);
255  for(int i=0;i<numSteps;++i) {
256  steps[i] = pow(10,-i);
257  }
258 
259  return checkGradient_1(u,z,g,d,steps,printToStream,outStream,order);
260  } // checkGradient_1
261 
262  std::vector<std::vector<Real> > checkGradient_1( const Vector<Real> &u,
263  const Vector<Real> &z,
264  const Vector<Real> &g,
265  const Vector<Real> &d,
266  const std::vector<Real> &steps,
267  const bool printToStream,
268  std::ostream & outStream,
269  const int order ) {
270  ROL_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
271  "Error: finite difference order must be 1,2,3, or 4" );
272 
275 
276  Real tol = std::sqrt(ROL_EPSILON<Real>());
277 
278  int numSteps = steps.size();
279  int numVals = 4;
280  std::vector<Real> tmp(numVals);
281  std::vector<std::vector<Real> > gCheck(numSteps, tmp);
282 
283  // Save the format state of the original outStream.
284  ROL::nullstream oldFormatState;
285  oldFormatState.copyfmt(outStream);
286 
287  // Evaluate objective value at x.
288  this->update(u,z,UpdateType::Temp);
289  Real val = this->value(u,z,tol);
290 
291  // Compute gradient at x.
292  ROL::Ptr<Vector<Real> > gtmp = g.clone();
293  this->gradient_1(*gtmp, u, z, tol);
294  //Real dtg = d.dot(gtmp->dual());
295  Real dtg = d.apply(*gtmp);
296 
297  // Temporary vectors.
298  ROL::Ptr<Vector<Real> > unew = u.clone();
299 
300  for (int i=0; i<numSteps; i++) {
301 
302  Real eta = steps[i];
303 
304  unew->set(u);
305 
306  // Compute gradient, finite-difference gradient, and absolute error.
307  gCheck[i][0] = eta;
308  gCheck[i][1] = dtg;
309 
310  gCheck[i][2] = weights[order-1][0] * val;
311 
312  for(int j=0; j<order; ++j) {
313  // Evaluate at x <- x+eta*c_i*d.
314  unew->axpy(eta*shifts[order-1][j], d);
315 
316  // Only evaluate at shifts where the weight is nonzero
317  if( weights[order-1][j+1] != 0 ) {
318  this->update(*unew,z,UpdateType::Temp);
319  gCheck[i][2] += weights[order-1][j+1] * this->value(*unew,z,tol);
320  }
321  }
322  gCheck[i][2] /= eta;
323 
324  gCheck[i][3] = std::abs(gCheck[i][2] - gCheck[i][1]);
325 
326  if (printToStream) {
327  if (i==0) {
328  outStream << std::right
329  << std::setw(20) << "Step size"
330  << std::setw(20) << "grad'*dir"
331  << std::setw(20) << "FD approx"
332  << std::setw(20) << "abs error"
333  << "\n"
334  << std::setw(20) << "---------"
335  << std::setw(20) << "---------"
336  << std::setw(20) << "---------"
337  << std::setw(20) << "---------"
338  << "\n";
339  }
340  outStream << std::scientific << std::setprecision(11) << std::right
341  << std::setw(20) << gCheck[i][0]
342  << std::setw(20) << gCheck[i][1]
343  << std::setw(20) << gCheck[i][2]
344  << std::setw(20) << gCheck[i][3]
345  << "\n";
346  }
347 
348  }
349 
350  // Reset format state of outStream.
351  outStream.copyfmt(oldFormatState);
352 
353  return gCheck;
354  } // checkGradient_1
355 
356 
357  std::vector<std::vector<Real> > checkGradient_2( const Vector<Real> &u,
358  const Vector<Real> &z,
359  const Vector<Real> &d,
360  const bool printToStream = true,
361  std::ostream & outStream = std::cout,
362  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
363  const int order = 1 ) {
364  return checkGradient_2(u, z, z.dual(), d, printToStream, outStream, numSteps, order);
365  }
366 
367  std::vector<std::vector<Real> > checkGradient_2( const Vector<Real> &u,
368  const Vector<Real> &z,
369  const Vector<Real> &g,
370  const Vector<Real> &d,
371  const bool printToStream,
372  std::ostream & outStream,
373  const int numSteps,
374  const int order ) {
375  std::vector<Real> steps(numSteps);
376  for(int i=0;i<numSteps;++i) {
377  steps[i] = pow(10,-i);
378  }
379 
380  return checkGradient_2(u,z,g,d,steps,printToStream,outStream,order);
381  } // checkGradient_2
382 
383  std::vector<std::vector<Real> > checkGradient_2( const Vector<Real> &u,
384  const Vector<Real> &z,
385  const Vector<Real> &g,
386  const Vector<Real> &d,
387  const std::vector<Real> &steps,
388  const bool printToStream,
389  std::ostream & outStream,
390  const int order ) {
391  ROL_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
392  "Error: finite difference order must be 1,2,3, or 4" );
393 
396 
397  Real tol = std::sqrt(ROL_EPSILON<Real>());
398 
399  int numSteps = steps.size();
400  int numVals = 4;
401  std::vector<Real> tmp(numVals);
402  std::vector<std::vector<Real> > gCheck(numSteps, tmp);
403 
404  // Save the format state of the original outStream.
405  ROL::nullstream oldFormatState;
406  oldFormatState.copyfmt(outStream);
407 
408  // Evaluate objective value at x.
409  this->update(u,z,UpdateType::Temp);
410  Real val = this->value(u,z,tol);
411 
412  // Compute gradient at x.
413  ROL::Ptr<Vector<Real> > gtmp = g.clone();
414  this->gradient_2(*gtmp, u, z, tol);
415  //Real dtg = d.dot(gtmp->dual());
416  Real dtg = d.apply(*gtmp);
417 
418  // Temporary vectors.
419  ROL::Ptr<Vector<Real> > znew = z.clone();
420 
421  for (int i=0; i<numSteps; i++) {
422 
423  Real eta = steps[i];
424 
425  znew->set(z);
426 
427  // Compute gradient, finite-difference gradient, and absolute error.
428  gCheck[i][0] = eta;
429  gCheck[i][1] = dtg;
430 
431  gCheck[i][2] = weights[order-1][0] * val;
432 
433  for(int j=0; j<order; ++j) {
434  // Evaluate at x <- x+eta*c_i*d.
435  znew->axpy(eta*shifts[order-1][j], d);
436 
437  // Only evaluate at shifts where the weight is nonzero
438  if( weights[order-1][j+1] != 0 ) {
439  this->update(u,*znew,UpdateType::Temp);
440  gCheck[i][2] += weights[order-1][j+1] * this->value(u,*znew,tol);
441  }
442  }
443  gCheck[i][2] /= eta;
444 
445  gCheck[i][3] = std::abs(gCheck[i][2] - gCheck[i][1]);
446 
447  if (printToStream) {
448  if (i==0) {
449  outStream << std::right
450  << std::setw(20) << "Step size"
451  << std::setw(20) << "grad'*dir"
452  << std::setw(20) << "FD approx"
453  << std::setw(20) << "abs error"
454  << "\n"
455  << std::setw(20) << "---------"
456  << std::setw(20) << "---------"
457  << std::setw(20) << "---------"
458  << std::setw(20) << "---------"
459  << "\n";
460  }
461  outStream << std::scientific << std::setprecision(11) << std::right
462  << std::setw(20) << gCheck[i][0]
463  << std::setw(20) << gCheck[i][1]
464  << std::setw(20) << gCheck[i][2]
465  << std::setw(20) << gCheck[i][3]
466  << "\n";
467  }
468 
469  }
470 
471  // Reset format state of outStream.
472  outStream.copyfmt(oldFormatState);
473 
474  return gCheck;
475  } // checkGradient_2
476 
477 
478  std::vector<std::vector<Real> > checkHessVec_11( const Vector<Real> &u,
479  const Vector<Real> &z,
480  const Vector<Real> &v,
481  const bool printToStream = true,
482  std::ostream & outStream = std::cout,
483  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
484  const int order = 1 ) {
485 
486  return checkHessVec_11(u, z, u.dual(), v, printToStream, outStream, numSteps, order);
487 
488  }
489 
490  std::vector<std::vector<Real> > checkHessVec_11( const Vector<Real> &u,
491  const Vector<Real> &z,
492  const Vector<Real> &v,
493  const std::vector<Real> &steps,
494  const bool printToStream = true,
495  std::ostream & outStream = std::cout,
496  const int order = 1 ) {
497 
498  return checkHessVec_11(u, z, u.dual(), v, steps, printToStream, outStream, order);
499  }
500 
501 
502  std::vector<std::vector<Real> > checkHessVec_11( const Vector<Real> &u,
503  const Vector<Real> &z,
504  const Vector<Real> &hv,
505  const Vector<Real> &v,
506  const bool printToStream,
507  std::ostream & outStream,
508  const int numSteps,
509  const int order ) {
510  std::vector<Real> steps(numSteps);
511  for(int i=0;i<numSteps;++i) {
512  steps[i] = pow(10,-i);
513  }
514 
515  return checkHessVec_11(u,z,hv,v,steps,printToStream,outStream,order);
516  } // checkHessVec_11
517 
518 
519  std::vector<std::vector<Real> > checkHessVec_11( const Vector<Real> &u,
520  const Vector<Real> &z,
521  const Vector<Real> &hv,
522  const Vector<Real> &v,
523  const std::vector<Real> &steps,
524  const bool printToStream,
525  std::ostream & outStream,
526  const int order ) {
527 
528  ROL_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
529  "Error: finite difference order must be 1,2,3, or 4" );
530 
533 
534 
535  Real tol = std::sqrt(ROL_EPSILON<Real>());
536 
537  int numSteps = steps.size();
538  int numVals = 4;
539  std::vector<Real> tmp(numVals);
540  std::vector<std::vector<Real> > hvCheck(numSteps, tmp);
541 
542  // Save the format state of the original outStream.
543  ROL::nullstream oldFormatState;
544  oldFormatState.copyfmt(outStream);
545 
546  // Compute gradient at x.
547  ROL::Ptr<Vector<Real> > g = hv.clone();
548  this->update(u,z,UpdateType::Temp);
549  this->gradient_1(*g, u, z, tol);
550 
551  // Compute (Hessian at x) times (vector v).
552  ROL::Ptr<Vector<Real> > Hv = hv.clone();
553  this->hessVec_11(*Hv, v, u, z, tol);
554  Real normHv = Hv->norm();
555 
556  // Temporary vectors.
557  ROL::Ptr<Vector<Real> > gdif = hv.clone();
558  ROL::Ptr<Vector<Real> > gnew = hv.clone();
559  ROL::Ptr<Vector<Real> > unew = u.clone();
560 
561  for (int i=0; i<numSteps; i++) {
562 
563  Real eta = steps[i];
564 
565  // Evaluate objective value at x+eta*d.
566  unew->set(u);
567 
568  gdif->set(*g);
569  gdif->scale(weights[order-1][0]);
570 
571  for(int j=0; j<order; ++j) {
572 
573  // Evaluate at x <- x+eta*c_i*d.
574  unew->axpy(eta*shifts[order-1][j], v);
575 
576  // Only evaluate at shifts where the weight is nonzero
577  if( weights[order-1][j+1] != 0 ) {
578  this->update(*unew,z,UpdateType::Temp);
579  this->gradient_1(*gnew, *unew, z, tol);
580  gdif->axpy(weights[order-1][j+1],*gnew);
581  }
582 
583  }
584 
585  gdif->scale(1.0/eta);
586 
587  // Compute norms of hessvec, finite-difference hessvec, and error.
588  hvCheck[i][0] = eta;
589  hvCheck[i][1] = normHv;
590  hvCheck[i][2] = gdif->norm();
591  gdif->axpy(-1.0, *Hv);
592  hvCheck[i][3] = gdif->norm();
593 
594  if (printToStream) {
595  if (i==0) {
596  outStream << std::right
597  << std::setw(20) << "Step size"
598  << std::setw(20) << "norm(Hess*vec)"
599  << std::setw(20) << "norm(FD approx)"
600  << std::setw(20) << "norm(abs error)"
601  << "\n"
602  << std::setw(20) << "---------"
603  << std::setw(20) << "--------------"
604  << std::setw(20) << "---------------"
605  << std::setw(20) << "---------------"
606  << "\n";
607  }
608  outStream << std::scientific << std::setprecision(11) << std::right
609  << std::setw(20) << hvCheck[i][0]
610  << std::setw(20) << hvCheck[i][1]
611  << std::setw(20) << hvCheck[i][2]
612  << std::setw(20) << hvCheck[i][3]
613  << "\n";
614  }
615 
616  }
617 
618  // Reset format state of outStream.
619  outStream.copyfmt(oldFormatState);
620 
621  return hvCheck;
622  } // checkHessVec_11
623 
624 
625  std::vector<std::vector<Real> > checkHessVec_12( const Vector<Real> &u,
626  const Vector<Real> &z,
627  const Vector<Real> &v,
628  const bool printToStream = true,
629  std::ostream & outStream = std::cout,
630  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
631  const int order = 1 ) {
632  return checkHessVec_12(u, z, u.dual(), v, printToStream, outStream, numSteps, order);
633  }
634 
635  std::vector<std::vector<Real> > checkHessVec_12( const Vector<Real> &u,
636  const Vector<Real> &z,
637  const Vector<Real> &v,
638  const std::vector<Real> &steps,
639  const bool printToStream = true,
640  std::ostream & outStream = std::cout,
641  const int order = 1 ) {
642  return checkHessVec_12(u, z, u.dual(), v, steps, printToStream, outStream, order);
643  }
644 
645 
646  std::vector<std::vector<Real> > checkHessVec_12( const Vector<Real> &u,
647  const Vector<Real> &z,
648  const Vector<Real> &hv,
649  const Vector<Real> &v,
650  const bool printToStream,
651  std::ostream & outStream,
652  const int numSteps,
653  const int order ) {
654  std::vector<Real> steps(numSteps);
655  for(int i=0;i<numSteps;++i) {
656  steps[i] = pow(10,-i);
657  }
658 
659  return checkHessVec_12(u,z,hv,v,steps,printToStream,outStream,order);
660  } // checkHessVec_12
661 
662 
663  std::vector<std::vector<Real> > checkHessVec_12( const Vector<Real> &u,
664  const Vector<Real> &z,
665  const Vector<Real> &hv,
666  const Vector<Real> &v,
667  const std::vector<Real> &steps,
668  const bool printToStream,
669  std::ostream & outStream,
670  const int order ) {
671 
672  ROL_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
673  "Error: finite difference order must be 1,2,3, or 4" );
674 
677 
678 
679  Real tol = std::sqrt(ROL_EPSILON<Real>());
680 
681  int numSteps = steps.size();
682  int numVals = 4;
683  std::vector<Real> tmp(numVals);
684  std::vector<std::vector<Real> > hvCheck(numSteps, tmp);
685 
686  // Save the format state of the original outStream.
687  ROL::nullstream oldFormatState;
688  oldFormatState.copyfmt(outStream);
689 
690  // Compute gradient at x.
691  ROL::Ptr<Vector<Real> > g = hv.clone();
692  this->update(u,z,UpdateType::Temp);
693  this->gradient_1(*g, u, z, tol);
694 
695  // Compute (Hessian at x) times (vector v).
696  ROL::Ptr<Vector<Real> > Hv = hv.clone();
697  this->hessVec_12(*Hv, v, u, z, tol);
698  Real normHv = Hv->norm();
699 
700  // Temporary vectors.
701  ROL::Ptr<Vector<Real> > gdif = hv.clone();
702  ROL::Ptr<Vector<Real> > gnew = hv.clone();
703  ROL::Ptr<Vector<Real> > znew = z.clone();
704 
705  for (int i=0; i<numSteps; i++) {
706 
707  Real eta = steps[i];
708 
709  // Evaluate objective value at x+eta*d.
710  znew->set(z);
711 
712  gdif->set(*g);
713  gdif->scale(weights[order-1][0]);
714 
715  for(int j=0; j<order; ++j) {
716 
717  // Evaluate at x <- x+eta*c_i*d.
718  znew->axpy(eta*shifts[order-1][j], v);
719 
720  // Only evaluate at shifts where the weight is nonzero
721  if( weights[order-1][j+1] != 0 ) {
722  this->update(u,*znew,UpdateType::Temp);
723  this->gradient_1(*gnew, u, *znew, tol);
724  gdif->axpy(weights[order-1][j+1],*gnew);
725  }
726 
727  }
728 
729  gdif->scale(1.0/eta);
730 
731  // Compute norms of hessvec, finite-difference hessvec, and error.
732  hvCheck[i][0] = eta;
733  hvCheck[i][1] = normHv;
734  hvCheck[i][2] = gdif->norm();
735  gdif->axpy(-1.0, *Hv);
736  hvCheck[i][3] = gdif->norm();
737 
738  if (printToStream) {
739  if (i==0) {
740  outStream << std::right
741  << std::setw(20) << "Step size"
742  << std::setw(20) << "norm(Hess*vec)"
743  << std::setw(20) << "norm(FD approx)"
744  << std::setw(20) << "norm(abs error)"
745  << "\n"
746  << std::setw(20) << "---------"
747  << std::setw(20) << "--------------"
748  << std::setw(20) << "---------------"
749  << std::setw(20) << "---------------"
750  << "\n";
751  }
752  outStream << std::scientific << std::setprecision(11) << std::right
753  << std::setw(20) << hvCheck[i][0]
754  << std::setw(20) << hvCheck[i][1]
755  << std::setw(20) << hvCheck[i][2]
756  << std::setw(20) << hvCheck[i][3]
757  << "\n";
758  }
759 
760  }
761 
762  // Reset format state of outStream.
763  outStream.copyfmt(oldFormatState);
764 
765  return hvCheck;
766  } // checkHessVec_12
767 
768 
769  std::vector<std::vector<Real> > checkHessVec_21( const Vector<Real> &u,
770  const Vector<Real> &z,
771  const Vector<Real> &v,
772  const bool printToStream = true,
773  std::ostream & outStream = std::cout,
774  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
775  const int order = 1 ) {
776 
777  return checkHessVec_21(u, z, z.dual(), v, printToStream, outStream, numSteps, order);
778 
779  }
780 
781  std::vector<std::vector<Real> > checkHessVec_21( const Vector<Real> &u,
782  const Vector<Real> &z,
783  const Vector<Real> &v,
784  const std::vector<Real> &steps,
785  const bool printToStream = true,
786  std::ostream & outStream = std::cout,
787  const int order = 1 ) {
788 
789  return checkHessVec_21(u, z, z.dual(), v, steps, printToStream, outStream, order);
790  }
791 
792 
793  std::vector<std::vector<Real> > checkHessVec_21( const Vector<Real> &u,
794  const Vector<Real> &z,
795  const Vector<Real> &hv,
796  const Vector<Real> &v,
797  const bool printToStream,
798  std::ostream & outStream,
799  const int numSteps,
800  const int order ) {
801  std::vector<Real> steps(numSteps);
802  for(int i=0;i<numSteps;++i) {
803  steps[i] = pow(10,-i);
804  }
805 
806  return checkHessVec_21(u,z,hv,v,steps,printToStream,outStream,order);
807  } // checkHessVec_21
808 
809 
810  std::vector<std::vector<Real> > checkHessVec_21( const Vector<Real> &u,
811  const Vector<Real> &z,
812  const Vector<Real> &hv,
813  const Vector<Real> &v,
814  const std::vector<Real> &steps,
815  const bool printToStream,
816  std::ostream & outStream,
817  const int order ) {
818 
819  ROL_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
820  "Error: finite difference order must be 1,2,3, or 4" );
821 
824 
825 
826  Real tol = std::sqrt(ROL_EPSILON<Real>());
827 
828  int numSteps = steps.size();
829  int numVals = 4;
830  std::vector<Real> tmp(numVals);
831  std::vector<std::vector<Real> > hvCheck(numSteps, tmp);
832 
833  // Save the format state of the original outStream.
834  ROL::nullstream oldFormatState;
835  oldFormatState.copyfmt(outStream);
836 
837  // Compute gradient at x.
838  ROL::Ptr<Vector<Real> > g = hv.clone();
839  this->update(u,z,UpdateType::Temp);
840  this->gradient_2(*g, u, z, tol);
841 
842  // Compute (Hessian at x) times (vector v).
843  ROL::Ptr<Vector<Real> > Hv = hv.clone();
844  this->hessVec_21(*Hv, v, u, z, tol);
845  Real normHv = Hv->norm();
846 
847  // Temporary vectors.
848  ROL::Ptr<Vector<Real> > gdif = hv.clone();
849  ROL::Ptr<Vector<Real> > gnew = hv.clone();
850  ROL::Ptr<Vector<Real> > unew = u.clone();
851 
852  for (int i=0; i<numSteps; i++) {
853 
854  Real eta = steps[i];
855 
856  // Evaluate objective value at x+eta*d.
857  unew->set(u);
858 
859  gdif->set(*g);
860  gdif->scale(weights[order-1][0]);
861 
862  for(int j=0; j<order; ++j) {
863 
864  // Evaluate at x <- x+eta*c_i*d.
865  unew->axpy(eta*shifts[order-1][j], v);
866 
867  // Only evaluate at shifts where the weight is nonzero
868  if( weights[order-1][j+1] != 0 ) {
869  this->update(*unew,z,UpdateType::Temp);
870  this->gradient_2(*gnew, *unew, z, tol);
871  gdif->axpy(weights[order-1][j+1],*gnew);
872  }
873 
874  }
875 
876  gdif->scale(1.0/eta);
877 
878  // Compute norms of hessvec, finite-difference hessvec, and error.
879  hvCheck[i][0] = eta;
880  hvCheck[i][1] = normHv;
881  hvCheck[i][2] = gdif->norm();
882  gdif->axpy(-1.0, *Hv);
883  hvCheck[i][3] = gdif->norm();
884 
885  if (printToStream) {
886  if (i==0) {
887  outStream << std::right
888  << std::setw(20) << "Step size"
889  << std::setw(20) << "norm(Hess*vec)"
890  << std::setw(20) << "norm(FD approx)"
891  << std::setw(20) << "norm(abs error)"
892  << "\n"
893  << std::setw(20) << "---------"
894  << std::setw(20) << "--------------"
895  << std::setw(20) << "---------------"
896  << std::setw(20) << "---------------"
897  << "\n";
898  }
899  outStream << std::scientific << std::setprecision(11) << std::right
900  << std::setw(20) << hvCheck[i][0]
901  << std::setw(20) << hvCheck[i][1]
902  << std::setw(20) << hvCheck[i][2]
903  << std::setw(20) << hvCheck[i][3]
904  << "\n";
905  }
906 
907  }
908 
909  // Reset format state of outStream.
910  outStream.copyfmt(oldFormatState);
911 
912  return hvCheck;
913  } // checkHessVec_21
914 
915 
916  std::vector<std::vector<Real> > checkHessVec_22( const Vector<Real> &u,
917  const Vector<Real> &z,
918  const Vector<Real> &v,
919  const bool printToStream = true,
920  std::ostream & outStream = std::cout,
921  const int numSteps = ROL_NUM_CHECKDERIV_STEPS,
922  const int order = 1 ) {
923 
924  return checkHessVec_22(u, z, z.dual(), v, printToStream, outStream, numSteps, order);
925 
926  }
927 
928  std::vector<std::vector<Real> > checkHessVec_22( const Vector<Real> &u,
929  const Vector<Real> &z,
930  const Vector<Real> &v,
931  const std::vector<Real> &steps,
932  const bool printToStream = true,
933  std::ostream & outStream = std::cout,
934  const int order = 1 ) {
935 
936  return checkHessVec_22(u, z, z.dual(), v, steps, printToStream, outStream, order);
937  }
938 
939 
940  std::vector<std::vector<Real> > checkHessVec_22( const Vector<Real> &u,
941  const Vector<Real> &z,
942  const Vector<Real> &hv,
943  const Vector<Real> &v,
944  const bool printToStream,
945  std::ostream & outStream,
946  const int numSteps,
947  const int order ) {
948  std::vector<Real> steps(numSteps);
949  for(int i=0;i<numSteps;++i) {
950  steps[i] = pow(10,-i);
951  }
952 
953  return checkHessVec_22(u,z,hv,v,steps,printToStream,outStream,order);
954  } // checkHessVec_22
955 
956 
957  std::vector<std::vector<Real> > checkHessVec_22( const Vector<Real> &u,
958  const Vector<Real> &z,
959  const Vector<Real> &hv,
960  const Vector<Real> &v,
961  const std::vector<Real> &steps,
962  const bool printToStream,
963  std::ostream & outStream,
964  const int order ) {
965 
966  ROL_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
967  "Error: finite difference order must be 1,2,3, or 4" );
968 
971 
972 
973  Real tol = std::sqrt(ROL_EPSILON<Real>());
974 
975  int numSteps = steps.size();
976  int numVals = 4;
977  std::vector<Real> tmp(numVals);
978  std::vector<std::vector<Real> > hvCheck(numSteps, tmp);
979 
980  // Save the format state of the original outStream.
981  ROL::nullstream oldFormatState;
982  oldFormatState.copyfmt(outStream);
983 
984  // Compute gradient at x.
985  ROL::Ptr<Vector<Real> > g = hv.clone();
986  this->update(u,z,UpdateType::Temp);
987  this->gradient_2(*g, u, z, tol);
988 
989  // Compute (Hessian at x) times (vector v).
990  ROL::Ptr<Vector<Real> > Hv = hv.clone();
991  this->hessVec_22(*Hv, v, u, z, tol);
992  Real normHv = Hv->norm();
993 
994  // Temporary vectors.
995  ROL::Ptr<Vector<Real> > gdif = hv.clone();
996  ROL::Ptr<Vector<Real> > gnew = hv.clone();
997  ROL::Ptr<Vector<Real> > znew = z.clone();
998 
999  for (int i=0; i<numSteps; i++) {
1000 
1001  Real eta = steps[i];
1002 
1003  // Evaluate objective value at x+eta*d.
1004  znew->set(z);
1005 
1006  gdif->set(*g);
1007  gdif->scale(weights[order-1][0]);
1008 
1009  for(int j=0; j<order; ++j) {
1010 
1011  // Evaluate at x <- x+eta*c_i*d.
1012  znew->axpy(eta*shifts[order-1][j], v);
1013 
1014  // Only evaluate at shifts where the weight is nonzero
1015  if( weights[order-1][j+1] != 0 ) {
1016  this->update(u,*znew,UpdateType::Temp);
1017  this->gradient_2(*gnew, u, *znew, tol);
1018  gdif->axpy(weights[order-1][j+1],*gnew);
1019  }
1020 
1021  }
1022 
1023  gdif->scale(1.0/eta);
1024 
1025  // Compute norms of hessvec, finite-difference hessvec, and error.
1026  hvCheck[i][0] = eta;
1027  hvCheck[i][1] = normHv;
1028  hvCheck[i][2] = gdif->norm();
1029  gdif->axpy(-1.0, *Hv);
1030  hvCheck[i][3] = gdif->norm();
1031 
1032  if (printToStream) {
1033  if (i==0) {
1034  outStream << std::right
1035  << std::setw(20) << "Step size"
1036  << std::setw(20) << "norm(Hess*vec)"
1037  << std::setw(20) << "norm(FD approx)"
1038  << std::setw(20) << "norm(abs error)"
1039  << "\n"
1040  << std::setw(20) << "---------"
1041  << std::setw(20) << "--------------"
1042  << std::setw(20) << "---------------"
1043  << std::setw(20) << "---------------"
1044  << "\n";
1045  }
1046  outStream << std::scientific << std::setprecision(11) << std::right
1047  << std::setw(20) << hvCheck[i][0]
1048  << std::setw(20) << hvCheck[i][1]
1049  << std::setw(20) << hvCheck[i][2]
1050  << std::setw(20) << hvCheck[i][3]
1051  << "\n";
1052  }
1053 
1054  }
1055 
1056  // Reset format state of outStream.
1057  outStream.copyfmt(oldFormatState);
1058 
1059  return hvCheck;
1060  } // checkHessVec_22
1061 
1062 }; // class Objective_SimOpt
1063 
1064 } // namespace ROL
1065 
1066 #endif
Provides the interface to evaluate objective functions.
Provides the interface to evaluate simulation-based objective functions.
std::vector< std::vector< Real > > checkGradient_2(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &g, const Vector< Real > &d, const std::vector< Real > &steps, const bool printToStream, std::ostream &outStream, const int order)
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
virtual void scale(const Real alpha)=0
Compute where .
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
ROL::Ptr< const Vector< Real > > get_2() const
std::vector< std::vector< Real > > checkGradient_1(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &d, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
virtual int dimension() const
Return dimension of the vector space.
Definition: ROL_Vector.hpp:162
virtual Real apply(const Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
Definition: ROL_Vector.hpp:204
virtual ROL::Ptr< Vector > basis(const int i) const
Return i-th basis vector.
Definition: ROL_Vector.hpp:148
const double weights[4][5]
Definition: ROL_Types.hpp:838
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:119
Defines the linear algebra or vector space interface for simulation-based optimization.
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
std::vector< std::vector< Real > > checkHessVec_21(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &v, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
std::vector< std::vector< Real > > checkHessVec_12(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &v, const std::vector< Real > &steps, const bool printToStream=true, std::ostream &outStream=std::cout, const int order=1)
void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update objective function.
virtual void hessVec_12(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
virtual void update(const Vector< Real > &u, const Vector< Real > &z, UpdateType type, int iter=-1)
void set_1(const Vector< Real > &vec)
std::vector< std::vector< Real > > checkGradient_2(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &d, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
virtual void gradient_2(Vector< Real > &g, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Compute gradient with respect to second component.
std::vector< std::vector< Real > > checkHessVec_21(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &v, const std::vector< Real > &steps, const bool printToStream=true, std::ostream &outStream=std::cout, const int order=1)
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:133
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
std::vector< std::vector< Real > > checkHessVec_22(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &v, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
std::vector< std::vector< Real > > checkHessVec_12(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &hv, const Vector< Real > &v, const bool printToStream, std::ostream &outStream, const int numSteps, const int order)
virtual Real dot(const Vector &x) const =0
Compute where .
std::vector< std::vector< Real > > checkHessVec_12(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &hv, const Vector< Real > &v, const std::vector< Real > &steps, const bool printToStream, std::ostream &outStream, const int order)
std::vector< std::vector< Real > > checkHessVec_11(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &v, const std::vector< Real > &steps, const bool printToStream=true, std::ostream &outStream=std::cout, const int order=1)
basic_nullstream< char, std::char_traits< char >> nullstream
Definition: ROL_Stream.hpp:36
virtual void hessVec_21(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
virtual Real value(const Vector< Real > &u, const Vector< Real > &z, Real &tol)=0
Compute value.
std::vector< std::vector< Real > > checkHessVec_11(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &v, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
std::vector< std::vector< Real > > checkHessVec_12(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &v, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
std::vector< std::vector< Real > > checkHessVec_22(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &hv, const Vector< Real > &v, const bool printToStream, std::ostream &outStream, const int numSteps, const int order)
std::vector< std::vector< Real > > checkHessVec_22(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &v, const std::vector< Real > &steps, const bool printToStream=true, std::ostream &outStream=std::cout, const int order=1)
std::vector< std::vector< Real > > checkHessVec_11(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &hv, const Vector< Real > &v, const std::vector< Real > &steps, const bool printToStream, std::ostream &outStream, const int order)
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
virtual void gradient_1(Vector< Real > &g, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Compute gradient with respect to first component.
virtual void hessVec_22(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
#define ROL_NUM_CHECKDERIV_STEPS
Number of steps for derivative checks.
Definition: ROL_Types.hpp:40
std::vector< std::vector< Real > > checkGradient_1(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &g, const Vector< Real > &d, const std::vector< Real > &steps, const bool printToStream, std::ostream &outStream, const int order)
virtual void update(const Vector< Real > &u, const Vector< Real > &z, bool flag=true, int iter=-1)
Update objective function. u is an iterate, z is an iterate, flag = true if the iterate has changed...
std::vector< std::vector< Real > > checkHessVec_21(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &hv, const Vector< Real > &v, const std::vector< Real > &steps, const bool printToStream, std::ostream &outStream, const int order)
virtual Real norm() const =0
Returns where .
Real value(const Vector< Real > &x, Real &tol)
Compute value.
std::vector< std::vector< Real > > checkHessVec_22(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &hv, const Vector< Real > &v, const std::vector< Real > &steps, const bool printToStream, std::ostream &outStream, const int order)
std::vector< std::vector< Real > > checkGradient_2(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &g, const Vector< Real > &d, const bool printToStream, std::ostream &outStream, const int numSteps, const int order)
virtual void hessVec_11(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Apply Hessian approximation to vector.
std::vector< std::vector< Real > > checkHessVec_21(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &hv, const Vector< Real > &v, const bool printToStream, std::ostream &outStream, const int numSteps, const int order)
std::vector< std::vector< Real > > checkGradient_1(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &g, const Vector< Real > &d, const bool printToStream, std::ostream &outStream, const int numSteps, const int order)
std::vector< std::vector< Real > > checkHessVec_11(const Vector< Real > &u, const Vector< Real > &z, const Vector< Real > &hv, const Vector< Real > &v, const bool printToStream, std::ostream &outStream, const int numSteps, const int order)
void set_2(const Vector< Real > &vec)
ROL::Ptr< const Vector< Real > > get_1() const