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
22namespace ROL {
23
24template <class Real>
25class Objective_SimOpt : public Objective<Real> {
26public:
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
#define ROL_NUM_CHECKDERIV_STEPS
Number of steps for derivative checks.
Definition ROL_Types.hpp:40
Provides the interface to evaluate simulation-based objective functions.
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)
void update(const Vector< Real > &x, bool flag=true, int iter=-1)
Update objective function.
virtual void hessVec_22(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
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 > &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_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)
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)
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 > > 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)
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 > > 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)
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)
virtual void update(const Vector< Real > &u, const Vector< Real > &z, UpdateType type, int iter=-1)
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)
Real value(const Vector< Real > &x, Real &tol)
Compute value.
virtual void gradient_2(Vector< Real > &g, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
Compute gradient with respect to second component.
virtual Real value(const Vector< Real > &u, const Vector< Real > &z, Real &tol)=0
Compute value.
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 > > 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)
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)
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.
virtual void hessVec_21(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
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)
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)
void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update objective function.
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)
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)
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)
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 > > 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)
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)
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_12(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &u, const Vector< Real > &z, Real &tol)
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)
Provides the interface to evaluate objective functions.
Defines the linear algebra or vector space interface for simulation-based optimization.
ROL::Ptr< const Vector< Real > > get_2() const
ROL::Ptr< const Vector< Real > > get_1() const
void set_1(const Vector< Real > &vec)
void set_2(const Vector< Real > &vec)
Defines the linear algebra or vector space interface.
virtual Real apply(const Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
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,...
virtual void zero()
Set to zero vector.
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual int dimension() const
Return dimension of the vector space.
virtual ROL::Ptr< Vector > basis(const int i) const
Return i-th basis vector.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
virtual Real dot(const Vector &x) const =0
Compute where .
const double weights[4][5]