ROL
ROL_CompositeStep.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_COMPOSITESTEP_H
11#define ROL_COMPOSITESTEP_H
12
13#include "ROL_Types.hpp"
14#include "ROL_Step.hpp"
15#include "ROL_LAPACK.hpp"
16#include "ROL_LinearAlgebra.hpp"
17#include <sstream>
18#include <iomanip>
19
26namespace ROL {
27
28template <class Real>
29class CompositeStep : public Step<Real> {
30private:
31
32 // Vectors used for cloning.
33 ROL::Ptr<Vector<Real> > xvec_;
34 ROL::Ptr<Vector<Real> > gvec_;
35 ROL::Ptr<Vector<Real> > cvec_;
36 ROL::Ptr<Vector<Real> > lvec_;
37
38 // Diagnostic return flags for subalgorithms.
42
43 // Stopping conditions.
46 Real tolCG_;
47 Real tolOSS_;
49
50 // Tolerances and stopping conditions for subalgorithms.
51 Real lmhtol_;
52 Real qntol_;
53 Real pgtol_;
56 Real tntmax_;
57
58 // Trust-region parameters.
59 Real zeta_;
60 Real Delta_;
62 Real eta_;
64
65 Real ared_;
66 Real pred_;
67 Real snorm_;
68 Real nnorm_;
69 Real tnorm_;
70
71 // Output flags.
72 bool infoQN_;
73 bool infoLM_;
74 bool infoTS_;
75 bool infoAC_;
76 bool infoLS_;
78
79 // Performance summary.
86
87 template <typename T> int sgn(T val) {
88 return (T(0) < val) - (val < T(0));
89 }
90
91 void printInfoLS(const std::vector<Real> &res) const {
92 if (infoLS_) {
93 std::stringstream hist;
94 hist << std::scientific << std::setprecision(8);
95 hist << "\n Augmented System Solver:\n";
96 hist << " True Residual\n";
97 for (unsigned j=0; j<res.size(); j++) {
98 hist << " " << std::left << std::setw(14) << res[j] << "\n";
99 }
100 hist << "\n";
101 std::cout << hist.str();
102 }
103 }
104
105 Real setTolOSS(const Real intol) const {
106 return tolOSSfixed_ ? tolOSS_ : intol;
107 }
108
109public:
110 using Step<Real>::initialize;
111 using Step<Real>::compute;
112 using Step<Real>::update;
113
114 virtual ~CompositeStep() {}
115
116 CompositeStep( ROL::ParameterList & parlist ) : Step<Real>() {
117 //ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState();
118 flagCG_ = 0;
119 flagAC_ = 0;
120 iterCG_ = 0;
121
122 ROL::ParameterList& steplist = parlist.sublist("Step").sublist("Composite Step");
123
124 //maxiterOSS_ = steplist.sublist("Optimality System Solver").get("Iteration Limit", 50);
125 tolOSS_ = steplist.sublist("Optimality System Solver").get("Nominal Relative Tolerance", 1e-8);
126 tolOSSfixed_ = steplist.sublist("Optimality System Solver").get("Fix Tolerance", true);
127
128 maxiterCG_ = steplist.sublist("Tangential Subproblem Solver").get("Iteration Limit", 20);
129 tolCG_ = steplist.sublist("Tangential Subproblem Solver").get("Relative Tolerance", 1e-2);
130 Delta_ = steplist.get("Initial Radius", 1e2);
131 useConHess_ = steplist.get("Use Constraint Hessian", true);
132
133 int outLvl = steplist.get("Output Level", 0);
134
136 qntol_ = tolOSS_;
137 pgtol_ = tolOSS_;
140 tntmax_ = 2.0;
141
142 zeta_ = 0.8;
143 penalty_ = 1.0;
144 eta_ = 1e-8;
145
146 snorm_ = 0.0;
147 nnorm_ = 0.0;
148 tnorm_ = 0.0;
149
150 infoALL_ = false;
151 if (outLvl > 0) {
152 infoALL_ = true;
153 }
154 infoQN_ = false;
155 infoLM_ = false;
156 infoTS_ = false;
157 infoAC_ = false;
158 infoLS_ = false;
164
165 totalIterCG_ = 0;
166 totalProj_ = 0;
167 totalNegCurv_ = 0;
168 totalRef_ = 0;
169 totalCallLS_ = 0;
170 totalIterLS_ = 0;
171 }
172
177 AlgorithmState<Real> &algo_state ) {
178 //ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState();
179 ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
180 state->descentVec = x.clone();
181 state->gradientVec = g.clone();
182 state->constraintVec = c.clone();
183
184 xvec_ = x.clone();
185 gvec_ = g.clone();
186 lvec_ = l.clone();
187 cvec_ = c.clone();
188
189 ROL::Ptr<Vector<Real> > ajl = gvec_->clone();
190 ROL::Ptr<Vector<Real> > gl = gvec_->clone();
191
192 algo_state.nfval = 0;
193 algo_state.ncval = 0;
194 algo_state.ngrad = 0;
195
196 Real zerotol = std::sqrt(ROL_EPSILON<Real>());
197
198 // Update objective and constraint.
199 obj.update(x,true,algo_state.iter);
200 algo_state.value = obj.value(x, zerotol);
201 algo_state.nfval++;
202 con.update(x,true,algo_state.iter);
203 con.value(*cvec_, x, zerotol);
204 algo_state.cnorm = cvec_->norm();
205 algo_state.ncval++;
206 obj.gradient(*gvec_, x, zerotol);
207
208 // Compute gradient of Lagrangian at new multiplier guess.
209 computeLagrangeMultiplier(l, x, *gvec_, con);
210 con.applyAdjointJacobian(*ajl, l, x, zerotol);
211 gl->set(*gvec_); gl->plus(*ajl);
212 algo_state.ngrad++;
213 algo_state.gnorm = gl->norm();
214 }
215
218 void compute( Vector<Real> &s, const Vector<Real> &x, const Vector<Real> &l,
220 AlgorithmState<Real> &algo_state ) {
221 //ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState();
222 Real zerotol = std::sqrt(ROL_EPSILON<Real>());
223 Real f = 0.0;
224 ROL::Ptr<Vector<Real> > n = xvec_->clone();
225 ROL::Ptr<Vector<Real> > c = cvec_->clone();
226 ROL::Ptr<Vector<Real> > t = xvec_->clone();
227 ROL::Ptr<Vector<Real> > tCP = xvec_->clone();
228 ROL::Ptr<Vector<Real> > g = gvec_->clone();
229 ROL::Ptr<Vector<Real> > gf = gvec_->clone();
230 ROL::Ptr<Vector<Real> > Wg = xvec_->clone();
231 ROL::Ptr<Vector<Real> > ajl = gvec_->clone();
232
233 Real f_new = 0.0;
234 ROL::Ptr<Vector<Real> > l_new = lvec_->clone();
235 ROL::Ptr<Vector<Real> > c_new = cvec_->clone();
236 ROL::Ptr<Vector<Real> > g_new = gvec_->clone();
237 ROL::Ptr<Vector<Real> > gf_new = gvec_->clone();
238
239 // Evaluate objective ... should have been stored.
240 f = obj.value(x, zerotol);
241 algo_state.nfval++;
242 // Compute gradient of objective ... should have been stored.
243 obj.gradient(*gf, x, zerotol);
244 // Evaluate constraint ... should have been stored.
245 con.value(*c, x, zerotol);
246
247 // Compute quasi-normal step.
248 computeQuasinormalStep(*n, *c, x, zeta_*Delta_, con);
249
250 // Compute gradient of Lagrangian ... should have been stored.
251 con.applyAdjointJacobian(*ajl, l, x, zerotol);
252 g->set(*gf);
253 g->plus(*ajl);
254 algo_state.ngrad++;
255
256 // Solve tangential subproblem.
257 solveTangentialSubproblem(*t, *tCP, *Wg, x, *g, *n, l, Delta_, obj, con);
259
260 // Check acceptance of subproblem solutions, adjust merit function penalty parameter, ensure global convergence.
261 accept(s, *n, *t, f_new, *c_new, *gf_new, *l_new, *g_new, x, l, f, *gf, *c, *g, *tCP, *Wg, obj, con, algo_state);
262 }
263
268 AlgorithmState<Real> &algo_state ) {
269 //ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
270
271 Real zero(0);
272 Real one(1);
273 Real two(2);
274 Real seven(7);
275 Real half(0.5);
276 Real zp9(0.9);
277 Real zp8(0.8);
278 Real em12(1e-12);
279 Real zerotol = std::sqrt(ROL_EPSILON<Real>());//zero;
280 Real ratio(zero);
281
282 ROL::Ptr<Vector<Real> > g = gvec_->clone();
283 ROL::Ptr<Vector<Real> > ajl = gvec_->clone();
284 ROL::Ptr<Vector<Real> > gl = gvec_->clone();
285 ROL::Ptr<Vector<Real> > c = cvec_->clone();
286
287 // Determine if the step gives sufficient reduction in the merit function,
288 // update the trust-region radius.
289 ratio = ared_/pred_;
290 if ((std::abs(ared_) < em12) && std::abs(pred_) < em12) {
291 ratio = one;
292 }
293 if (ratio >= eta_) {
294 x.plus(s);
295 if (ratio >= zp9) {
296 Delta_ = std::max(seven*snorm_, Delta_);
297 }
298 else if (ratio >= zp8) {
299 Delta_ = std::max(two*snorm_, Delta_);
300 }
301 obj.update(x,true,algo_state.iter);
302 con.update(x,true,algo_state.iter);
303 flagAC_ = 1;
304 }
305 else {
306 Delta_ = half*std::max(nnorm_, tnorm_);
307 obj.update(x,false,algo_state.iter);
308 con.update(x,false,algo_state.iter);
309 flagAC_ = 0;
310 } // if (ratio >= eta)
311
312 Real val = obj.value(x, zerotol);
313 algo_state.nfval++;
314 obj.gradient(*g, x, zerotol);
315 computeLagrangeMultiplier(l, x, *g, con);
316 con.applyAdjointJacobian(*ajl, l, x, zerotol);
317 gl->set(*g); gl->plus(*ajl);
318 algo_state.ngrad++;
319 con.value(*c, x, zerotol);
320
321 ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
322 state->gradientVec->set(*gl);
323 state->constraintVec->set(*c);
324
325 algo_state.value = val;
326 algo_state.gnorm = gl->norm();
327 algo_state.cnorm = c->norm();
328 algo_state.iter++;
329 algo_state.snorm = snorm_;
330
331 // Update algorithm state
332 //(algo_state.iterateVec)->set(x);
333 }
334
340 AlgorithmState<Real> &algo_state ) {}
341
347 AlgorithmState<Real> &algo_state ) {}
348
351 std::string printHeader( void ) const {
352 std::stringstream hist;
353 hist << " ";
354 hist << std::setw(6) << std::left << "iter";
355 hist << std::setw(15) << std::left << "fval";
356 hist << std::setw(15) << std::left << "cnorm";
357 hist << std::setw(15) << std::left << "gLnorm";
358 hist << std::setw(15) << std::left << "snorm";
359 hist << std::setw(10) << std::left << "delta";
360 hist << std::setw(10) << std::left << "nnorm";
361 hist << std::setw(10) << std::left << "tnorm";
362 hist << std::setw(8) << std::left << "#fval";
363 hist << std::setw(8) << std::left << "#grad";
364 hist << std::setw(8) << std::left << "iterCG";
365 hist << std::setw(8) << std::left << "flagCG";
366 hist << std::setw(8) << std::left << "accept";
367 hist << std::setw(8) << std::left << "linsys";
368 hist << "\n";
369 return hist.str();
370 }
371
372 std::string printName( void ) const {
373 std::stringstream hist;
374 hist << "\n" << " Composite-step trust-region solver";
375 hist << "\n";
376 return hist.str();
377 }
378
381 std::string print( AlgorithmState<Real> & algo_state, bool pHeader = false ) const {
382 //const ROL::Ptr<const StepState<Real> >& step_state = Step<Real>::getStepState();
383
384 std::stringstream hist;
385 hist << std::scientific << std::setprecision(6);
386 if ( algo_state.iter == 0 ) {
387 hist << printName();
388 }
389 if ( pHeader ) {
390 hist << printHeader();
391 }
392 if ( algo_state.iter == 0 ) {
393 hist << " ";
394 hist << std::setw(6) << std::left << algo_state.iter;
395 hist << std::setw(15) << std::left << algo_state.value;
396 hist << std::setw(15) << std::left << algo_state.cnorm;
397 hist << std::setw(15) << std::left << algo_state.gnorm;
398 hist << "\n";
399 }
400 else {
401 hist << " ";
402 hist << std::setw(6) << std::left << algo_state.iter;
403 hist << std::setw(15) << std::left << algo_state.value;
404 hist << std::setw(15) << std::left << algo_state.cnorm;
405 hist << std::setw(15) << std::left << algo_state.gnorm;
406 hist << std::setw(15) << std::left << algo_state.snorm;
407 hist << std::scientific << std::setprecision(2);
408 hist << std::setw(10) << std::left << Delta_;
409 hist << std::setw(10) << std::left << nnorm_;
410 hist << std::setw(10) << std::left << tnorm_;
411 hist << std::scientific << std::setprecision(6);
412 hist << std::setw(8) << std::left << algo_state.nfval;
413 hist << std::setw(8) << std::left << algo_state.ngrad;
414 hist << std::setw(8) << std::left << iterCG_;
415 hist << std::setw(8) << std::left << flagCG_;
416 hist << std::setw(8) << std::left << flagAC_;
417 hist << std::left << totalCallLS_ << "/" << totalIterLS_;
418 hist << "\n";
419 }
420 return hist.str();
421 }
422
435
436 Real one(1);
437 Real zerotol = std::sqrt(ROL_EPSILON<Real>());
438 std::vector<Real> augiters;
439
440 if (infoLM_) {
441 std::stringstream hist;
442 hist << "\n Lagrange multiplier step\n";
443 std::cout << hist.str();
444 }
445
446 /* Apply adjoint of constraint Jacobian to current multiplier. */
447 ROL::Ptr<Vector<Real> > ajl = gvec_->clone();
448 con.applyAdjointJacobian(*ajl, l, x, zerotol);
449
450 /* Form right-hand side of the augmented system. */
451 ROL::Ptr<Vector<Real> > b1 = gvec_->clone();
452 ROL::Ptr<Vector<Real> > b2 = cvec_->clone();
453 // b1 is the negative gradient of the Lagrangian
454 b1->set(gf); b1->plus(*ajl); b1->scale(-one);
455 // b2 is zero
456 b2->zero();
457
458 /* Declare left-hand side of augmented system. */
459 ROL::Ptr<Vector<Real> > v1 = xvec_->clone();
460 ROL::Ptr<Vector<Real> > v2 = lvec_->clone();
461
462 /* Compute linear solver tolerance. */
463 Real b1norm = b1->norm();
464 Real tol = setTolOSS(lmhtol_*b1norm);
465
466 /* Solve augmented system. */
467 augiters = con.solveAugmentedSystem(*v1, *v2, *b1, *b2, x, tol);
468 totalCallLS_++;
469 totalIterLS_ = totalIterLS_ + augiters.size();
470 printInfoLS(augiters);
471
472 /* Return updated Lagrange multiplier. */
473 // v2 is the multiplier update
474 l.plus(*v2);
475
476 } // computeLagrangeMultiplier
477
478
501 void computeQuasinormalStep(Vector<Real> &n, const Vector<Real> &c, const Vector<Real> &x, Real delta, Constraint<Real> &con) {
502
503 if (infoQN_) {
504 std::stringstream hist;
505 hist << "\n Quasi-normal step\n";
506 std::cout << hist.str();
507 }
508
509 Real zero(0);
510 Real one(1);
511 Real zerotol = std::sqrt(ROL_EPSILON<Real>()); //zero;
512 std::vector<Real> augiters;
513
514 /* Compute Cauchy step nCP. */
515 ROL::Ptr<Vector<Real> > nCP = xvec_->clone();
516 ROL::Ptr<Vector<Real> > nCPdual = gvec_->clone();
517 ROL::Ptr<Vector<Real> > nN = xvec_->clone();
518 ROL::Ptr<Vector<Real> > ctemp = cvec_->clone();
519 ROL::Ptr<Vector<Real> > dualc0 = lvec_->clone();
520 dualc0->set(c.dual());
521 con.applyAdjointJacobian(*nCPdual, *dualc0, x, zerotol);
522 nCP->set(nCPdual->dual());
523 con.applyJacobian(*ctemp, *nCP, x, zerotol);
524
525 Real normsquare_ctemp = ctemp->dot(*ctemp);
526 if (normsquare_ctemp != zero) {
527 nCP->scale( -(nCP->dot(*nCP))/normsquare_ctemp );
528 }
529
530 /* If the Cauchy step nCP is outside the trust region,
531 return the scaled Cauchy step. */
532 Real norm_nCP = nCP->norm();
533 if (norm_nCP >= delta) {
534 n.set(*nCP);
535 n.scale( delta/norm_nCP );
536 if (infoQN_) {
537 std::stringstream hist;
538 hist << " taking partial Cauchy step\n";
539 std::cout << hist.str();
540 }
541 return;
542 }
543
544 /* Compute 'Newton' step, for example, by solving a problem
545 related to finding the minimum norm solution of min || c(x_k)*s + c ||^2. */
546 // Compute tolerance for linear solver.
547 con.applyJacobian(*ctemp, *nCP, x, zerotol);
548 ctemp->plus(c);
549 Real tol = setTolOSS(qntol_*ctemp->norm());
550 // Form right-hand side.
551 ctemp->scale(-one);
552 nCPdual->set(nCP->dual());
553 nCPdual->scale(-one);
554 // Declare left-hand side of augmented system.
555 ROL::Ptr<Vector<Real> > dn = xvec_->clone();
556 ROL::Ptr<Vector<Real> > y = lvec_->clone();
557 // Solve augmented system.
558 augiters = con.solveAugmentedSystem(*dn, *y, *nCPdual, *ctemp, x, tol);
559 totalCallLS_++;
560 totalIterLS_ = totalIterLS_ + augiters.size();
561 printInfoLS(augiters);
562
563 nN->set(*dn);
564 nN->plus(*nCP);
565
566 /* Either take full or partial Newton step, depending on
567 the trust-region constraint. */
568 Real norm_nN = nN->norm();
569 if (norm_nN <= delta) {
570 // Take full feasibility step.
571 n.set(*nN);
572 if (infoQN_) {
573 std::stringstream hist;
574 hist << " taking full Newton step\n";
575 std::cout << hist.str();
576 }
577 }
578 else {
579 // Take convex combination n = nCP+tau*(nN-nCP),
580 // so that ||n|| = delta. In other words, solve
581 // scalar quadratic equation: ||nCP+tau*(nN-nCP)||^2 = delta^2.
582 Real aa = dn->dot(*dn);
583 Real bb = dn->dot(*nCP);
584 Real cc = norm_nCP*norm_nCP - delta*delta;
585 Real tau = (-bb+sqrt(bb*bb-aa*cc))/aa;
586 n.set(*nCP);
587 n.axpy(tau, *dn);
588 if (infoQN_) {
589 std::stringstream hist;
590 hist << " taking dogleg step\n";
591 std::cout << hist.str();
592 }
593 }
594
595 } // computeQuasinormalStep
596
597
612 const Vector<Real> &x, const Vector<Real> &g, const Vector<Real> &n, const Vector<Real> &l,
613 Real delta, Objective<Real> &obj, Constraint<Real> &con) {
614
615 /* Initialization of the CG step. */
616 bool orthocheck = true; // set to true if want to check orthogonality
617 // of Wr and r, otherwise set to false
618 Real tol_ortho = 0.5; // orthogonality measure; represets a bound on norm( \hat{S}, 2), where
619 // \hat{S} is defined in Heinkenschloss/Ridzal., "A Matrix-Free Trust-Region SQP Method"
620 Real S_max = 1.0; // another orthogonality measure; norm(S) needs to be bounded by
621 // a modest constant; norm(S) is small if the approximation of
622 // the null space projector is good
623 Real zero(0);
624 Real one(1);
625 Real zerotol = std::sqrt(ROL_EPSILON<Real>());
626 std::vector<Real> augiters;
627 iterCG_ = 1;
628 flagCG_ = 0;
629 t.zero();
630 tCP.zero();
631 ROL::Ptr<Vector<Real> > r = gvec_->clone();
632 ROL::Ptr<Vector<Real> > pdesc = xvec_->clone();
633 ROL::Ptr<Vector<Real> > tprev = xvec_->clone();
634 ROL::Ptr<Vector<Real> > Wr = xvec_->clone();
635 ROL::Ptr<Vector<Real> > Hp = gvec_->clone();
636 ROL::Ptr<Vector<Real> > xtemp = xvec_->clone();
637 ROL::Ptr<Vector<Real> > gtemp = gvec_->clone();
638 ROL::Ptr<Vector<Real> > ltemp = lvec_->clone();
639 ROL::Ptr<Vector<Real> > czero = cvec_->clone();
640 czero->zero();
641 r->set(g);
642 obj.hessVec(*gtemp, n, x, zerotol);
643 r->plus(*gtemp);
644 if (useConHess_) {
645 con.applyAdjointHessian(*gtemp, l, n, x, zerotol);
646 r->plus(*gtemp);
647 }
648 Real normg = r->norm();
649 Real normWg = zero;
650 Real pHp = zero;
651 Real rp = zero;
652 Real alpha = zero;
653 Real normp = zero;
654 Real normr = zero;
655 Real normt = zero;
656 std::vector<Real> normWr(maxiterCG_+1, zero);
657
658 std::vector<ROL::Ptr<Vector<Real > > > p; // stores search directions
659 std::vector<ROL::Ptr<Vector<Real > > > Hps; // stores duals of hessvec's applied to p's
660 std::vector<ROL::Ptr<Vector<Real > > > rs; // stores duals of residuals
661 std::vector<ROL::Ptr<Vector<Real > > > Wrs; // stores duals of projected residuals
662
663 Real rptol(1e-12);
664
665 if (infoTS_) {
666 std::stringstream hist;
667 hist << "\n Tangential subproblem\n";
668 hist << std::setw(6) << std::right << "iter" << std::setw(18) << "||Wr||/||Wr0||" << std::setw(15) << "||s||";
669 hist << std::setw(15) << "delta" << std::setw(15) << "||c'(x)s||" << "\n";
670 std::cout << hist.str();
671 }
672
673 if (normg == 0) {
674 if (infoTS_) {
675 std::stringstream hist;
676 hist << " >>> Tangential subproblem: Initial gradient is zero! \n";
677 std::cout << hist.str();
678 }
679 iterCG_ = 0; Wg.zero(); flagCG_ = 0;
680 return;
681 }
682
683 /* Start CG loop. */
684 while (iterCG_ < maxiterCG_) {
685
686 // Store tangential Cauchy point (which is the current iterate in the second iteration).
687 if (iterCG_ == 2) {
688 tCP.set(t);
689 }
690
691 // Compute (inexact) projection W*r.
692 if (iterCG_ == 1) {
693 // Solve augmented system.
694 Real tol = setTolOSS(pgtol_);
695 augiters = con.solveAugmentedSystem(*Wr, *ltemp, *r, *czero, x, tol);
696 totalCallLS_++;
697 totalIterLS_ = totalIterLS_ + augiters.size();
698 printInfoLS(augiters);
699
700 Wg.set(*Wr);
701 normWg = Wg.norm();
702 if (orthocheck) {
703 Wrs.push_back(xvec_->clone());
704 (Wrs[iterCG_-1])->set(*Wr);
705 }
706 // Check if done (small initial projected residual).
707 if (normWg == zero) {
708 flagCG_ = 0;
709 iterCG_--;
710 if (infoTS_) {
711 std::stringstream hist;
712 hist << " Initial projected residual is close to zero! \n";
713 std::cout << hist.str();
714 }
715 return;
716 }
717 // Set first residual to projected gradient.
718 // change r->set(Wg);
719 r->set(Wg.dual());
720 if (orthocheck) {
721 rs.push_back(xvec_->clone());
722 // change (rs[0])->set(*r);
723 (rs[0])->set(r->dual());
724 }
725 }
726 else {
727 // Solve augmented system.
728 Real tol = setTolOSS(projtol_);
729 augiters = con.solveAugmentedSystem(*Wr, *ltemp, *r, *czero, x, tol);
730 totalCallLS_++;
731 totalIterLS_ = totalIterLS_ + augiters.size();
732 printInfoLS(augiters);
733
734 if (orthocheck) {
735 Wrs.push_back(xvec_->clone());
736 (Wrs[iterCG_-1])->set(*Wr);
737 }
738 }
739
740 normWr[iterCG_-1] = Wr->norm();
741
742 if (infoTS_) {
743 ROL::Ptr<Vector<Real> > ct = cvec_->clone();
744 con.applyJacobian(*ct, t, x, zerotol);
745 Real linc = ct->norm();
746 std::stringstream hist;
747 hist << std::scientific << std::setprecision(6);
748 hist << std::setw(6) << std::right << iterCG_-1 << std::setw(18) << normWr[iterCG_-1]/normWg << std::setw(15) << t.norm();
749 hist << std::setw(15) << delta << std::setw(15) << linc << "\n";
750 std::cout << hist.str();
751 }
752
753 // Check if done (small relative residual).
754 if (normWr[iterCG_-1]/normWg < tolCG_) {
755 flagCG_ = 0;
756 iterCG_ = iterCG_-1;
757 if (infoTS_) {
758 std::stringstream hist;
759 hist << " || W(g + H*(n+s)) || <= cgtol*|| W(g + H*n)|| \n";
760 std::cout << hist.str();
761 }
762 return;
763 }
764
765 // Check nonorthogonality, one-norm of (WR*R/diag^2 - I)
766 if (orthocheck) {
767 ROL::LA::Matrix<Real> Wrr(iterCG_,iterCG_); // holds matrix Wrs'*rs
768 ROL::LA::Matrix<Real> T(iterCG_,iterCG_); // holds matrix T=(1/diag)*Wrs'*rs*(1/diag)
769 ROL::LA::Matrix<Real> Tm1(iterCG_,iterCG_); // holds matrix Tm1=T-I
770 for (int i=0; i<iterCG_; i++) {
771 for (int j=0; j<iterCG_; j++) {
772 Wrr(i,j) = (Wrs[i])->dot(*rs[j]);
773 T(i,j) = Wrr(i,j)/(normWr[i]*normWr[j]);
774 Tm1(i,j) = T(i,j);
775 if (i==j) {
776 Tm1(i,j) = Tm1(i,j) - one;
777 }
778 }
779 }
780 if (Tm1.normOne() >= tol_ortho) {
781 ROL::LAPACK<int,Real> lapack;
782 std::vector<int> ipiv(iterCG_);
783 int info;
784 std::vector<Real> work(3*iterCG_);
785 // compute inverse of T
786 lapack.GETRF(iterCG_, iterCG_, T.values(), T.stride(), &ipiv[0], &info);
787 lapack.GETRI(iterCG_, T.values(), T.stride(), &ipiv[0], &work[0], 3*iterCG_, &info);
788 Tm1 = T;
789 for (int i=0; i<iterCG_; i++) {
790 Tm1(i,i) = Tm1(i,i) - one;
791 }
792 if (Tm1.normOne() > S_max) {
793 flagCG_ = 4;
794 if (infoTS_) {
795 std::stringstream hist;
796 hist << " large nonorthogonality in W(R)'*R detected \n";
797 std::cout << hist.str();
798 }
799 return;
800 }
801 }
802 }
803
804 // Full orthogonalization.
805 p.push_back(xvec_->clone());
806 (p[iterCG_-1])->set(*Wr);
807 (p[iterCG_-1])->scale(-one);
808 for (int j=1; j<iterCG_; j++) {
809 Real scal = (p[iterCG_-1])->dot(*(Hps[j-1])) / (p[j-1])->dot(*(Hps[j-1]));
810 ROL::Ptr<Vector<Real> > pj = xvec_->clone();
811 pj->set(*p[j-1]);
812 pj->scale(-scal);
813 (p[iterCG_-1])->plus(*pj);
814 }
815
816 // change Hps.push_back(gvec_->clone());
817 Hps.push_back(xvec_->clone());
818 // change obj.hessVec(*(Hps[iterCG_-1]), *(p[iterCG_-1]), x, zerotol);
819 obj.hessVec(*Hp, *(p[iterCG_-1]), x, zerotol);
820 if (useConHess_) {
821 con.applyAdjointHessian(*gtemp, l, *(p[iterCG_-1]), x, zerotol);
822 // change (Hps[iterCG_-1])->plus(*gtemp);
823 Hp->plus(*gtemp);
824 }
825 // "Preconditioning" step.
826 (Hps[iterCG_-1])->set(Hp->dual());
827
828 pHp = (p[iterCG_-1])->dot(*(Hps[iterCG_-1]));
829 // change rp = (p[iterCG_-1])->dot(*r);
830 rp = (p[iterCG_-1])->dot(*(rs[iterCG_-1]));
831
832 normp = (p[iterCG_-1])->norm();
833 normr = r->norm();
834
835 // Negative curvature stopping condition.
836 if (pHp <= 0) {
837 pdesc->set(*(p[iterCG_-1])); // p is the descent direction
838 if ((std::abs(rp) >= rptol*normp*normr) && (sgn(rp) == 1)) {
839 pdesc->scale(-one); // -p is the descent direction
840 }
841 flagCG_ = 2;
842 Real a = pdesc->dot(*pdesc);
843 Real b = pdesc->dot(t);
844 Real c = t.dot(t) - delta*delta;
845 // Positive root of a*theta^2 + 2*b*theta + c = 0.
846 Real theta = (-b + std::sqrt(b*b - a*c)) / a;
847 xtemp->set(*(p[iterCG_-1]));
848 xtemp->scale(theta);
849 t.plus(*xtemp);
850 // Store as tangential Cauchy point if terminating in first iteration.
851 if (iterCG_ == 1) {
852 tCP.set(t);
853 }
854 if (infoTS_) {
855 std::stringstream hist;
856 hist << " negative curvature detected \n";
857 std::cout << hist.str();
858 }
859 return;
860 }
861
862 // Want to enforce nonzero alpha's.
863 if (std::abs(rp) < rptol*normp*normr) {
864 flagCG_ = 5;
865 if (infoTS_) {
866 std::stringstream hist;
867 hist << " Zero alpha due to inexactness. \n";
868 std::cout << hist.str();
869 }
870 return;
871 }
872
873 alpha = - rp/pHp;
874
875 // Iterate update.
876 tprev->set(t);
877 xtemp->set(*(p[iterCG_-1]));
878 xtemp->scale(alpha);
879 t.plus(*xtemp);
880
881 // Trust-region stopping condition.
882 normt = t.norm();
883 if (normt >= delta) {
884 pdesc->set(*(p[iterCG_-1])); // p is the descent direction
885 if (sgn(rp) == 1) {
886 pdesc->scale(-one); // -p is the descent direction
887 }
888 Real a = pdesc->dot(*pdesc);
889 Real b = pdesc->dot(*tprev);
890 Real c = tprev->dot(*tprev) - delta*delta;
891 // Positive root of a*theta^2 + 2*b*theta + c = 0.
892 Real theta = (-b + std::sqrt(b*b - a*c)) / a;
893 xtemp->set(*(p[iterCG_-1]));
894 xtemp->scale(theta);
895 t.set(*tprev);
896 t.plus(*xtemp);
897 // Store as tangential Cauchy point if terminating in first iteration.
898 if (iterCG_ == 1) {
899 tCP.set(t);
900 }
901 flagCG_ = 3;
902 if (infoTS_) {
903 std::stringstream hist;
904 hist << " trust-region condition active \n";
905 std::cout << hist.str();
906 }
907 return;
908 }
909
910 // Residual update.
911 xtemp->set(*(Hps[iterCG_-1]));
912 xtemp->scale(alpha);
913 // change r->plus(*gtemp);
914 r->plus(xtemp->dual());
915 if (orthocheck) {
916 // change rs.push_back(gvec_->clone());
917 rs.push_back(xvec_->clone());
918 // change (rs[iterCG_])->set(*r);
919 (rs[iterCG_])->set(r->dual());
920 }
921
922 iterCG_++;
923
924 } // while (iterCG_ < maxiterCG_)
925
926 flagCG_ = 1;
927 if (infoTS_) {
928 std::stringstream hist;
929 hist << " maximum number of iterations reached \n";
930 std::cout << hist.str();
931 }
932
933 } // solveTangentialSubproblem
934
935
938 void accept(Vector<Real> &s, Vector<Real> &n, Vector<Real> &t, Real f_new, Vector<Real> &c_new,
939 Vector<Real> &gf_new, Vector<Real> &l_new, Vector<Real> &g_new,
940 const Vector<Real> &x, const Vector<Real> &l, Real f, const Vector<Real> &gf, const Vector<Real> &c,
941 const Vector<Real> &g, Vector<Real> &tCP, Vector<Real> &Wg,
942 Objective<Real> &obj, Constraint<Real> &con, AlgorithmState<Real> &algo_state) {
943
944 Real beta = 1e-8; // predicted reduction parameter
945 Real tol_red_tang = 1e-3; // internal reduction factor for tangtol
946 Real tol_red_all = 1e-1; // internal reduction factor for qntol, lmhtol, pgtol, projtol, tangtol
947 //bool glob_refine = true; // true - if subsolver tolerances are adjusted in this routine, keep adjusted values globally
948 // false - if subsolver tolerances are adjusted in this routine, discard adjusted values
949 Real tol_fdiff = 1e-12; // relative objective function difference for ared computation
950 int ct_max = 10; // maximum number of globalization tries
951 Real mintol = 1e-16; // smallest projection tolerance value
952
953 // Determines max value of |rpred|/pred.
954 Real rpred_over_pred = 0.5*(1-eta_);
955
956 if (infoAC_) {
957 std::stringstream hist;
958 hist << "\n Composite step acceptance\n";
959 std::cout << hist.str();
960 }
961
962 Real zero = 0.0;
963 Real one = 1.0;
964 Real two = 2.0;
965 Real half = one/two;
966 Real zerotol = std::sqrt(ROL_EPSILON<Real>());
967 std::vector<Real> augiters;
968
969 Real pred = zero;
970 Real ared = zero;
971 Real rpred = zero;
972 Real part_pred = zero;
973 Real linc_preproj = zero;
974 Real linc_postproj = zero;
975 Real tangtol_start = zero;
976 Real tangtol = tangtol_;
977 //Real projtol = projtol_;
978 bool flag = false;
979 int num_proj = 0;
980 bool try_tCP = false;
981 Real fdiff = zero;
982
983 ROL::Ptr<Vector<Real> > xtrial = xvec_->clone();
984 ROL::Ptr<Vector<Real> > Jl = gvec_->clone();
985 ROL::Ptr<Vector<Real> > gfJl = gvec_->clone();
986 ROL::Ptr<Vector<Real> > Jnc = cvec_->clone();
987 ROL::Ptr<Vector<Real> > t_orig = xvec_->clone();
988 ROL::Ptr<Vector<Real> > t_dual = gvec_->clone();
989 ROL::Ptr<Vector<Real> > Jt_orig = cvec_->clone();
990 ROL::Ptr<Vector<Real> > t_m_tCP = xvec_->clone();
991 ROL::Ptr<Vector<Real> > ltemp = lvec_->clone();
992 ROL::Ptr<Vector<Real> > xtemp = xvec_->clone();
993 ROL::Ptr<Vector<Real> > rt = cvec_->clone();
994 ROL::Ptr<Vector<Real> > Hn = gvec_->clone();
995 ROL::Ptr<Vector<Real> > Hto = gvec_->clone();
996 ROL::Ptr<Vector<Real> > cxxvec = gvec_->clone();
997 ROL::Ptr<Vector<Real> > czero = cvec_->clone();
998 czero->zero();
999 Real Jnc_normsquared = zero;
1000 Real c_normsquared = zero;
1001
1002 // Compute and store some quantities for later use. Necessary
1003 // because of the function and constraint updates below.
1004 con.applyAdjointJacobian(*Jl, l, x, zerotol);
1005 con.applyJacobian(*Jnc, n, x, zerotol);
1006 Jnc->plus(c);
1007 Jnc_normsquared = Jnc->dot(*Jnc);
1008 c_normsquared = c.dot(c);
1009
1010 for (int ct=0; ct<ct_max; ct++) {
1011
1012 try_tCP = true;
1013 t_m_tCP->set(t);
1014 t_m_tCP->scale(-one);
1015 t_m_tCP->plus(tCP);
1016 if (t_m_tCP->norm() == zero) {
1017 try_tCP = false;
1018 }
1019
1020 t_orig->set(t);
1021 con.applyJacobian(*Jt_orig, *t_orig, x, zerotol);
1022 linc_preproj = Jt_orig->norm();
1023 pred = one;
1024 rpred = two*rpred_over_pred*pred;
1025 flag = false;
1026 num_proj = 1;
1027 tangtol_start = tangtol;
1028
1029 while (std::abs(rpred)/pred > rpred_over_pred) {
1030 // Compute projected tangential step.
1031 if (flag) {
1032 tangtol = tol_red_tang*tangtol;
1033 num_proj++;
1034 if (tangtol < mintol) {
1035 if (infoAC_) {
1036 std::stringstream hist;
1037 hist << "\n The projection of the tangential step cannot be done with sufficient precision.\n";
1038 hist << " Is the quasi-normal step very small? Continuing with no global convergence guarantees.\n";
1039 std::cout << hist.str();
1040 }
1041 break;
1042 }
1043 }
1044 // Solve augmented system.
1045 Real tol = setTolOSS(tangtol);
1046 // change augiters = con.solveAugmentedSystem(t, *ltemp, *t_orig, *czero, x, tol);
1047 t_dual->set(t_orig->dual());
1048 augiters = con.solveAugmentedSystem(t, *ltemp, *t_dual, *czero, x, tol);
1049 totalCallLS_++;
1050 totalIterLS_ = totalIterLS_ + augiters.size();
1051 printInfoLS(augiters);
1052 totalProj_++;
1053 con.applyJacobian(*rt, t, x, zerotol);
1054 linc_postproj = rt->norm();
1055
1056 // Compute composite step.
1057 s.set(t);
1058 s.plus(n);
1059
1060 // Compute some quantities before updating the objective and the constraint.
1061 obj.hessVec(*Hn, n, x, zerotol);
1062 if (useConHess_) {
1063 con.applyAdjointHessian(*cxxvec, l, n, x, zerotol);
1064 Hn->plus(*cxxvec);
1065 }
1066 obj.hessVec(*Hto, *t_orig, x, zerotol);
1067 if (useConHess_) {
1068 con.applyAdjointHessian(*cxxvec, l, *t_orig, x, zerotol);
1069 Hto->plus(*cxxvec);
1070 }
1071
1072 // Compute objective, constraint, etc. values at the trial point.
1073 xtrial->set(x);
1074 xtrial->plus(s);
1075 obj.update(*xtrial,false,algo_state.iter);
1076 con.update(*xtrial,false,algo_state.iter);
1077 f_new = obj.value(*xtrial, zerotol);
1078 obj.gradient(gf_new, *xtrial, zerotol);
1079 con.value(c_new, *xtrial, zerotol);
1080 l_new.set(l);
1081 computeLagrangeMultiplier(l_new, *xtrial, gf_new, con);
1082
1083 // Penalty parameter update.
1084 part_pred = - Wg.dot(*t_orig);
1085 gfJl->set(gf);
1086 gfJl->plus(*Jl);
1087 // change part_pred -= gfJl->dot(n);
1088 part_pred -= n.dot(gfJl->dual());
1089 // change part_pred -= half*Hn->dot(n);
1090 part_pred -= half*n.dot(Hn->dual());
1091 // change part_pred -= half*Hto->dot(*t_orig);
1092 part_pred -= half*t_orig->dot(Hto->dual());
1093 ltemp->set(l_new);
1094 ltemp->axpy(-one, l);
1095 // change part_pred -= Jnc->dot(*ltemp);
1096 part_pred -= Jnc->dot(ltemp->dual());
1097
1098 if ( part_pred < -half*penalty_*(c_normsquared-Jnc_normsquared) ) {
1099 penalty_ = ( -two * part_pred / (c_normsquared-Jnc_normsquared) ) + beta;
1100 }
1101
1102 pred = part_pred + penalty_*(c_normsquared-Jnc_normsquared);
1103
1104 // Computation of rpred.
1105 // change rpred = - ltemp->dot(*rt) - penalty_ * rt->dot(*rt) - two * penalty_ * rt->dot(*Jnc);
1106 rpred = - rt->dot(ltemp->dual()) - penalty_ * rt->dot(*rt) - two * penalty_ * rt->dot(*Jnc);
1107 // change ROL::Ptr<Vector<Real> > lrt = lvec_->clone();
1108 //lrt->set(*rt);
1109 //rpred = - ltemp->dot(*rt) - penalty_ * std::pow(rt->norm(), 2) - two * penalty_ * lrt->dot(*Jnc);
1110 flag = 1;
1111
1112 } // while (std::abs(rpred)/pred > rpred_over_pred)
1113
1114 tangtol = tangtol_start;
1115
1116 // Check if the solution of the tangential subproblem is
1117 // disproportionally large compared to total trial step.
1118 xtemp->set(n);
1119 xtemp->plus(t);
1120 if ( t_orig->norm()/xtemp->norm() < tntmax_ ) {
1121 break;
1122 }
1123 else {
1124 t_m_tCP->set(*t_orig);
1125 t_m_tCP->scale(-one);
1126 t_m_tCP->plus(tCP);
1127 if ((t_m_tCP->norm() > 0) && try_tCP) {
1128 if (infoAC_) {
1129 std::stringstream hist;
1130 hist << " ---> now trying tangential Cauchy point\n";
1131 std::cout << hist.str();
1132 }
1133 t.set(tCP);
1134 }
1135 else {
1136 if (infoAC_) {
1137 std::stringstream hist;
1138 hist << " ---> recomputing quasi-normal step and re-solving tangential subproblem\n";
1139 std::cout << hist.str();
1140 }
1141 totalRef_++;
1142 // Reset global quantities.
1143 obj.update(x);
1144 con.update(x);
1145 /*lmhtol = tol_red_all*lmhtol;
1146 qntol = tol_red_all*qntol;
1147 pgtol = tol_red_all*pgtol;
1148 projtol = tol_red_all*projtol;
1149 tangtol = tol_red_all*tangtol;
1150 if (glob_refine) {
1151 lmhtol_ = lmhtol;
1152 qntol_ = qntol;
1153 pgtol_ = pgtol;
1154 projtol_ = projtol;
1155 tangtol_ = tangtol;
1156 }*/
1157 if (!tolOSSfixed_) {
1158 lmhtol_ *= tol_red_all;
1159 qntol_ *= tol_red_all;
1160 pgtol_ *= tol_red_all;
1161 projtol_ *= tol_red_all;
1162 tangtol_ *= tol_red_all;
1163 }
1164 // Recompute the quasi-normal step.
1165 computeQuasinormalStep(n, c, x, zeta_*Delta_, con);
1166 // Solve tangential subproblem.
1167 solveTangentialSubproblem(t, tCP, Wg, x, g, n, l, Delta_, obj, con);
1169 if (flagCG_ == 1) {
1170 totalNegCurv_++;
1171 }
1172 }
1173 } // else w.r.t. ( t_orig->norm()/xtemp->norm() < tntmax )
1174
1175 } // for (int ct=0; ct<ct_max; ct++)
1176
1177 // Compute actual reduction;
1178 fdiff = f - f_new;
1179 // Heuristic 1: If fdiff is very small compared to f, set it to 0,
1180 // in order to prevent machine precision issues.
1181 Real em24(1e-24);
1182 Real em14(1e-14);
1183 if (std::abs(fdiff / (f+em24)) < tol_fdiff) {
1184 fdiff = em14;
1185 }
1186 // change ared = fdiff + (l.dot(c) - l_new.dot(c_new)) + penalty_*(c.dot(c) - c_new.dot(c_new));
1187 // change ared = fdiff + (l.dot(c) - l_new.dot(c_new)) + penalty_*(std::pow(c.norm(),2) - std::pow(c_new.norm(),2));
1188 ared = fdiff + (c.dot(l.dual()) - c_new.dot(l_new.dual())) + penalty_*(c.dot(c) - c_new.dot(c_new));
1189
1190 // Store actual and predicted reduction.
1191 ared_ = ared;
1192 pred_ = pred;
1193
1194 // Store step and vector norms.
1195 snorm_ = s.norm();
1196 nnorm_ = n.norm();
1197 tnorm_ = t.norm();
1198
1199 // Print diagnostics.
1200 if (infoAC_) {
1201 std::stringstream hist;
1202 hist << "\n Trial step info ...\n";
1203 hist << " n_norm = " << nnorm_ << "\n";
1204 hist << " t_norm = " << tnorm_ << "\n";
1205 hist << " s_norm = " << snorm_ << "\n";
1206 hist << " xtrial_norm = " << xtrial->norm() << "\n";
1207 hist << " f_old = " << f << "\n";
1208 hist << " f_trial = " << f_new << "\n";
1209 hist << " f_old-f_trial = " << f-f_new << "\n";
1210 hist << " ||c_old|| = " << c.norm() << "\n";
1211 hist << " ||c_trial|| = " << c_new.norm() << "\n";
1212 hist << " ||Jac*t_preproj|| = " << linc_preproj << "\n";
1213 hist << " ||Jac*t_postproj|| = " << linc_postproj << "\n";
1214 hist << " ||t_tilde||/||t|| = " << t_orig->norm() / t.norm() << "\n";
1215 hist << " ||t_tilde||/||n+t|| = " << t_orig->norm() / snorm_ << "\n";
1216 hist << " # projections = " << num_proj << "\n";
1217 hist << " penalty param = " << penalty_ << "\n";
1218 hist << " ared = " << ared_ << "\n";
1219 hist << " pred = " << pred_ << "\n";
1220 hist << " ared/pred = " << ared_/pred_ << "\n";
1221 std::cout << hist.str();
1222 }
1223
1224 } // accept
1225
1226}; // class CompositeStep
1227
1228} // namespace ROL
1229
1230#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Contains definitions of custom data types in ROL.
Provides the interface to apply upper and lower bound constraints.
Implements the computation of optimization steps with composite-step trust-region methods.
void initialize(Vector< Real > &x, const Vector< Real > &g, Vector< Real > &l, const Vector< Real > &c, Objective< Real > &obj, Constraint< Real > &con, AlgorithmState< Real > &algo_state)
Initialize step.
void compute(Vector< Real > &s, const Vector< Real > &x, Objective< Real > &obj, BoundConstraint< Real > &con, AlgorithmState< Real > &algo_state)
Compute step for bound constraints; here only to satisfy the interface requirements,...
void update(Vector< Real > &x, const Vector< Real > &s, Objective< Real > &obj, BoundConstraint< Real > &con, AlgorithmState< Real > &algo_state)
Update step, for bound constraints; here only to satisfy the interface requirements,...
ROL::Ptr< Vector< Real > > xvec_
ROL::Ptr< Vector< Real > > lvec_
Real setTolOSS(const Real intol) const
std::string printName(void) const
Print step name.
void computeLagrangeMultiplier(Vector< Real > &l, const Vector< Real > &x, const Vector< Real > &gf, Constraint< Real > &con)
Compute Lagrange multipliers by solving the least-squares problem minimizing the gradient of the Lagr...
ROL::Ptr< Vector< Real > > gvec_
void printInfoLS(const std::vector< Real > &res) const
void compute(Vector< Real > &s, const Vector< Real > &x, const Vector< Real > &l, Objective< Real > &obj, Constraint< Real > &con, AlgorithmState< Real > &algo_state)
Compute step.
ROL::Ptr< Vector< Real > > cvec_
void computeQuasinormalStep(Vector< Real > &n, const Vector< Real > &c, const Vector< Real > &x, Real delta, Constraint< Real > &con)
Compute quasi-normal step by minimizing the norm of the linearized constraint.
void accept(Vector< Real > &s, Vector< Real > &n, Vector< Real > &t, Real f_new, Vector< Real > &c_new, Vector< Real > &gf_new, Vector< Real > &l_new, Vector< Real > &g_new, const Vector< Real > &x, const Vector< Real > &l, Real f, const Vector< Real > &gf, const Vector< Real > &c, const Vector< Real > &g, Vector< Real > &tCP, Vector< Real > &Wg, Objective< Real > &obj, Constraint< Real > &con, AlgorithmState< Real > &algo_state)
Check acceptance of subproblem solutions, adjust merit function penalty parameter,...
void solveTangentialSubproblem(Vector< Real > &t, Vector< Real > &tCP, Vector< Real > &Wg, const Vector< Real > &x, const Vector< Real > &g, const Vector< Real > &n, const Vector< Real > &l, Real delta, Objective< Real > &obj, Constraint< Real > &con)
Solve tangential subproblem.
CompositeStep(ROL::ParameterList &parlist)
std::string printHeader(void) const
Print iterate header.
std::string print(AlgorithmState< Real > &algo_state, bool pHeader=false) const
Print iterate status.
void update(Vector< Real > &x, Vector< Real > &l, const Vector< Real > &s, Objective< Real > &obj, Constraint< Real > &con, AlgorithmState< Real > &algo_state)
Update step, if successful.
Defines the general constraint operator interface.
virtual void value(Vector< Real > &c, const Vector< Real > &x, Real &tol)=0
Evaluate the constraint operator at .
virtual void applyAdjointHessian(Vector< Real > &ahuv, const Vector< Real > &u, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the derivative of the adjoint of the constraint Jacobian at to vector in direction ,...
virtual void applyAdjointJacobian(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the adjoint of the the constraint Jacobian at , , to vector .
virtual void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update constraint function.
virtual void applyJacobian(Vector< Real > &jv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the constraint Jacobian at , , to vector .
virtual std::vector< Real > solveAugmentedSystem(Vector< Real > &v1, Vector< Real > &v2, const Vector< Real > &b1, const Vector< Real > &b2, const Vector< Real > &x, Real &tol)
Approximately solves the augmented system
Provides the interface to evaluate objective functions.
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
virtual Real value(const Vector< Real > &x, Real &tol)=0
Compute value.
virtual void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update objective function.
Provides the interface to compute optimization steps.
Definition ROL_Step.hpp:34
ROL::Ptr< StepState< Real > > getState(void)
Definition ROL_Step.hpp:39
Defines the linear algebra or vector space interface.
virtual Real norm() const =0
Returns where .
virtual void set(const Vector &x)
Set 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 plus(const Vector &x)=0
Compute , where .
virtual void zero()
Set to zero vector.
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
virtual Real dot(const Vector &x) const =0
Compute where .
State for algorithm class. Will be used for restarts.