ROL
ROL_CompositeStep.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Rapid Optimization Library (ROL) Package
5// Copyright (2014) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact lead developers:
38// Drew Kouri (dpkouri@sandia.gov) and
39// Denis Ridzal (dridzal@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
44#ifndef ROL_COMPOSITESTEP_H
45#define ROL_COMPOSITESTEP_H
46
47#include "ROL_Types.hpp"
48#include "ROL_Step.hpp"
49#include "ROL_LAPACK.hpp"
50#include "ROL_LinearAlgebra.hpp"
51#include <sstream>
52#include <iomanip>
53
60namespace ROL {
61
62template <class Real>
63class CompositeStep : public Step<Real> {
64private:
65
66 // Vectors used for cloning.
67 ROL::Ptr<Vector<Real> > xvec_;
68 ROL::Ptr<Vector<Real> > gvec_;
69 ROL::Ptr<Vector<Real> > cvec_;
70 ROL::Ptr<Vector<Real> > lvec_;
71
72 // Diagnostic return flags for subalgorithms.
76
77 // Stopping conditions.
80 Real tolCG_;
81 Real tolOSS_;
83
84 // Tolerances and stopping conditions for subalgorithms.
85 Real lmhtol_;
86 Real qntol_;
87 Real pgtol_;
90 Real tntmax_;
91
92 // Trust-region parameters.
93 Real zeta_;
94 Real Delta_;
96 Real eta_;
98
99 Real ared_;
100 Real pred_;
101 Real snorm_;
102 Real nnorm_;
103 Real tnorm_;
104
105 // Output flags.
112
113 // Performance summary.
120
121 template <typename T> int sgn(T val) {
122 return (T(0) < val) - (val < T(0));
123 }
124
125 void printInfoLS(const std::vector<Real> &res) const {
126 if (infoLS_) {
127 std::stringstream hist;
128 hist << std::scientific << std::setprecision(8);
129 hist << "\n Augmented System Solver:\n";
130 hist << " True Residual\n";
131 for (unsigned j=0; j<res.size(); j++) {
132 hist << " " << std::left << std::setw(14) << res[j] << "\n";
133 }
134 hist << "\n";
135 std::cout << hist.str();
136 }
137 }
138
139 Real setTolOSS(const Real intol) const {
140 return tolOSSfixed_ ? tolOSS_ : intol;
141 }
142
143public:
144 using Step<Real>::initialize;
145 using Step<Real>::compute;
146 using Step<Real>::update;
147
148 virtual ~CompositeStep() {}
149
150 CompositeStep( ROL::ParameterList & parlist ) : Step<Real>() {
151 //ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState();
152 flagCG_ = 0;
153 flagAC_ = 0;
154 iterCG_ = 0;
155
156 ROL::ParameterList& steplist = parlist.sublist("Step").sublist("Composite Step");
157
158 //maxiterOSS_ = steplist.sublist("Optimality System Solver").get("Iteration Limit", 50);
159 tolOSS_ = steplist.sublist("Optimality System Solver").get("Nominal Relative Tolerance", 1e-8);
160 tolOSSfixed_ = steplist.sublist("Optimality System Solver").get("Fix Tolerance", true);
161
162 maxiterCG_ = steplist.sublist("Tangential Subproblem Solver").get("Iteration Limit", 20);
163 tolCG_ = steplist.sublist("Tangential Subproblem Solver").get("Relative Tolerance", 1e-2);
164 Delta_ = steplist.get("Initial Radius", 1e2);
165 useConHess_ = steplist.get("Use Constraint Hessian", true);
166
167 int outLvl = steplist.get("Output Level", 0);
168
170 qntol_ = tolOSS_;
171 pgtol_ = tolOSS_;
174 tntmax_ = 2.0;
175
176 zeta_ = 0.8;
177 penalty_ = 1.0;
178 eta_ = 1e-8;
179
180 snorm_ = 0.0;
181 nnorm_ = 0.0;
182 tnorm_ = 0.0;
183
184 infoALL_ = false;
185 if (outLvl > 0) {
186 infoALL_ = true;
187 }
188 infoQN_ = false;
189 infoLM_ = false;
190 infoTS_ = false;
191 infoAC_ = false;
192 infoLS_ = false;
198
199 totalIterCG_ = 0;
200 totalProj_ = 0;
201 totalNegCurv_ = 0;
202 totalRef_ = 0;
203 totalCallLS_ = 0;
204 totalIterLS_ = 0;
205 }
206
211 AlgorithmState<Real> &algo_state ) {
212 //ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState();
213 ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
214 state->descentVec = x.clone();
215 state->gradientVec = g.clone();
216 state->constraintVec = c.clone();
217
218 xvec_ = x.clone();
219 gvec_ = g.clone();
220 lvec_ = l.clone();
221 cvec_ = c.clone();
222
223 ROL::Ptr<Vector<Real> > ajl = gvec_->clone();
224 ROL::Ptr<Vector<Real> > gl = gvec_->clone();
225
226 algo_state.nfval = 0;
227 algo_state.ncval = 0;
228 algo_state.ngrad = 0;
229
230 Real zerotol = std::sqrt(ROL_EPSILON<Real>());
231
232 // Update objective and constraint.
233 obj.update(x,true,algo_state.iter);
234 algo_state.value = obj.value(x, zerotol);
235 algo_state.nfval++;
236 con.update(x,true,algo_state.iter);
237 con.value(*cvec_, x, zerotol);
238 algo_state.cnorm = cvec_->norm();
239 algo_state.ncval++;
240 obj.gradient(*gvec_, x, zerotol);
241
242 // Compute gradient of Lagrangian at new multiplier guess.
243 computeLagrangeMultiplier(l, x, *gvec_, con);
244 con.applyAdjointJacobian(*ajl, l, x, zerotol);
245 gl->set(*gvec_); gl->plus(*ajl);
246 algo_state.ngrad++;
247 algo_state.gnorm = gl->norm();
248 }
249
252 void compute( Vector<Real> &s, const Vector<Real> &x, const Vector<Real> &l,
254 AlgorithmState<Real> &algo_state ) {
255 //ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState();
256 Real zerotol = std::sqrt(ROL_EPSILON<Real>());
257 Real f = 0.0;
258 ROL::Ptr<Vector<Real> > n = xvec_->clone();
259 ROL::Ptr<Vector<Real> > c = cvec_->clone();
260 ROL::Ptr<Vector<Real> > t = xvec_->clone();
261 ROL::Ptr<Vector<Real> > tCP = xvec_->clone();
262 ROL::Ptr<Vector<Real> > g = gvec_->clone();
263 ROL::Ptr<Vector<Real> > gf = gvec_->clone();
264 ROL::Ptr<Vector<Real> > Wg = xvec_->clone();
265 ROL::Ptr<Vector<Real> > ajl = gvec_->clone();
266
267 Real f_new = 0.0;
268 ROL::Ptr<Vector<Real> > l_new = lvec_->clone();
269 ROL::Ptr<Vector<Real> > c_new = cvec_->clone();
270 ROL::Ptr<Vector<Real> > g_new = gvec_->clone();
271 ROL::Ptr<Vector<Real> > gf_new = gvec_->clone();
272
273 // Evaluate objective ... should have been stored.
274 f = obj.value(x, zerotol);
275 algo_state.nfval++;
276 // Compute gradient of objective ... should have been stored.
277 obj.gradient(*gf, x, zerotol);
278 // Evaluate constraint ... should have been stored.
279 con.value(*c, x, zerotol);
280
281 // Compute quasi-normal step.
282 computeQuasinormalStep(*n, *c, x, zeta_*Delta_, con);
283
284 // Compute gradient of Lagrangian ... should have been stored.
285 con.applyAdjointJacobian(*ajl, l, x, zerotol);
286 g->set(*gf);
287 g->plus(*ajl);
288 algo_state.ngrad++;
289
290 // Solve tangential subproblem.
291 solveTangentialSubproblem(*t, *tCP, *Wg, x, *g, *n, l, Delta_, obj, con);
293
294 // Check acceptance of subproblem solutions, adjust merit function penalty parameter, ensure global convergence.
295 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);
296 }
297
302 AlgorithmState<Real> &algo_state ) {
303 //ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
304
305 Real zero(0);
306 Real one(1);
307 Real two(2);
308 Real seven(7);
309 Real half(0.5);
310 Real zp9(0.9);
311 Real zp8(0.8);
312 Real em12(1e-12);
313 Real zerotol = std::sqrt(ROL_EPSILON<Real>());//zero;
314 Real ratio(zero);
315
316 ROL::Ptr<Vector<Real> > g = gvec_->clone();
317 ROL::Ptr<Vector<Real> > ajl = gvec_->clone();
318 ROL::Ptr<Vector<Real> > gl = gvec_->clone();
319 ROL::Ptr<Vector<Real> > c = cvec_->clone();
320
321 // Determine if the step gives sufficient reduction in the merit function,
322 // update the trust-region radius.
323 ratio = ared_/pred_;
324 if ((std::abs(ared_) < em12) && std::abs(pred_) < em12) {
325 ratio = one;
326 }
327 if (ratio >= eta_) {
328 x.plus(s);
329 if (ratio >= zp9) {
330 Delta_ = std::max(seven*snorm_, Delta_);
331 }
332 else if (ratio >= zp8) {
333 Delta_ = std::max(two*snorm_, Delta_);
334 }
335 obj.update(x,true,algo_state.iter);
336 con.update(x,true,algo_state.iter);
337 flagAC_ = 1;
338 }
339 else {
340 Delta_ = half*std::max(nnorm_, tnorm_);
341 obj.update(x,false,algo_state.iter);
342 con.update(x,false,algo_state.iter);
343 flagAC_ = 0;
344 } // if (ratio >= eta)
345
346 Real val = obj.value(x, zerotol);
347 algo_state.nfval++;
348 obj.gradient(*g, x, zerotol);
349 computeLagrangeMultiplier(l, x, *g, con);
350 con.applyAdjointJacobian(*ajl, l, x, zerotol);
351 gl->set(*g); gl->plus(*ajl);
352 algo_state.ngrad++;
353 con.value(*c, x, zerotol);
354
355 ROL::Ptr<StepState<Real> > state = Step<Real>::getState();
356 state->gradientVec->set(*gl);
357 state->constraintVec->set(*c);
358
359 algo_state.value = val;
360 algo_state.gnorm = gl->norm();
361 algo_state.cnorm = c->norm();
362 algo_state.iter++;
363 algo_state.snorm = snorm_;
364
365 // Update algorithm state
366 //(algo_state.iterateVec)->set(x);
367 }
368
374 AlgorithmState<Real> &algo_state ) {}
375
381 AlgorithmState<Real> &algo_state ) {}
382
385 std::string printHeader( void ) const {
386 std::stringstream hist;
387 hist << " ";
388 hist << std::setw(6) << std::left << "iter";
389 hist << std::setw(15) << std::left << "fval";
390 hist << std::setw(15) << std::left << "cnorm";
391 hist << std::setw(15) << std::left << "gLnorm";
392 hist << std::setw(15) << std::left << "snorm";
393 hist << std::setw(10) << std::left << "delta";
394 hist << std::setw(10) << std::left << "nnorm";
395 hist << std::setw(10) << std::left << "tnorm";
396 hist << std::setw(8) << std::left << "#fval";
397 hist << std::setw(8) << std::left << "#grad";
398 hist << std::setw(8) << std::left << "iterCG";
399 hist << std::setw(8) << std::left << "flagCG";
400 hist << std::setw(8) << std::left << "accept";
401 hist << std::setw(8) << std::left << "linsys";
402 hist << "\n";
403 return hist.str();
404 }
405
406 std::string printName( void ) const {
407 std::stringstream hist;
408 hist << "\n" << " Composite-step trust-region solver";
409 hist << "\n";
410 return hist.str();
411 }
412
415 std::string print( AlgorithmState<Real> & algo_state, bool pHeader = false ) const {
416 //const ROL::Ptr<const StepState<Real> >& step_state = Step<Real>::getStepState();
417
418 std::stringstream hist;
419 hist << std::scientific << std::setprecision(6);
420 if ( algo_state.iter == 0 ) {
421 hist << printName();
422 }
423 if ( pHeader ) {
424 hist << printHeader();
425 }
426 if ( algo_state.iter == 0 ) {
427 hist << " ";
428 hist << std::setw(6) << std::left << algo_state.iter;
429 hist << std::setw(15) << std::left << algo_state.value;
430 hist << std::setw(15) << std::left << algo_state.cnorm;
431 hist << std::setw(15) << std::left << algo_state.gnorm;
432 hist << "\n";
433 }
434 else {
435 hist << " ";
436 hist << std::setw(6) << std::left << algo_state.iter;
437 hist << std::setw(15) << std::left << algo_state.value;
438 hist << std::setw(15) << std::left << algo_state.cnorm;
439 hist << std::setw(15) << std::left << algo_state.gnorm;
440 hist << std::setw(15) << std::left << algo_state.snorm;
441 hist << std::scientific << std::setprecision(2);
442 hist << std::setw(10) << std::left << Delta_;
443 hist << std::setw(10) << std::left << nnorm_;
444 hist << std::setw(10) << std::left << tnorm_;
445 hist << std::scientific << std::setprecision(6);
446 hist << std::setw(8) << std::left << algo_state.nfval;
447 hist << std::setw(8) << std::left << algo_state.ngrad;
448 hist << std::setw(8) << std::left << iterCG_;
449 hist << std::setw(8) << std::left << flagCG_;
450 hist << std::setw(8) << std::left << flagAC_;
451 hist << std::left << totalCallLS_ << "/" << totalIterLS_;
452 hist << "\n";
453 }
454 return hist.str();
455 }
456
469
470 Real one(1);
471 Real zerotol = std::sqrt(ROL_EPSILON<Real>());
472 std::vector<Real> augiters;
473
474 if (infoLM_) {
475 std::stringstream hist;
476 hist << "\n Lagrange multiplier step\n";
477 std::cout << hist.str();
478 }
479
480 /* Apply adjoint of constraint Jacobian to current multiplier. */
481 ROL::Ptr<Vector<Real> > ajl = gvec_->clone();
482 con.applyAdjointJacobian(*ajl, l, x, zerotol);
483
484 /* Form right-hand side of the augmented system. */
485 ROL::Ptr<Vector<Real> > b1 = gvec_->clone();
486 ROL::Ptr<Vector<Real> > b2 = cvec_->clone();
487 // b1 is the negative gradient of the Lagrangian
488 b1->set(gf); b1->plus(*ajl); b1->scale(-one);
489 // b2 is zero
490 b2->zero();
491
492 /* Declare left-hand side of augmented system. */
493 ROL::Ptr<Vector<Real> > v1 = xvec_->clone();
494 ROL::Ptr<Vector<Real> > v2 = lvec_->clone();
495
496 /* Compute linear solver tolerance. */
497 Real b1norm = b1->norm();
498 Real tol = setTolOSS(lmhtol_*b1norm);
499
500 /* Solve augmented system. */
501 augiters = con.solveAugmentedSystem(*v1, *v2, *b1, *b2, x, tol);
502 totalCallLS_++;
503 totalIterLS_ = totalIterLS_ + augiters.size();
504 printInfoLS(augiters);
505
506 /* Return updated Lagrange multiplier. */
507 // v2 is the multiplier update
508 l.plus(*v2);
509
510 } // computeLagrangeMultiplier
511
512
535 void computeQuasinormalStep(Vector<Real> &n, const Vector<Real> &c, const Vector<Real> &x, Real delta, Constraint<Real> &con) {
536
537 if (infoQN_) {
538 std::stringstream hist;
539 hist << "\n Quasi-normal step\n";
540 std::cout << hist.str();
541 }
542
543 Real zero(0);
544 Real one(1);
545 Real zerotol = std::sqrt(ROL_EPSILON<Real>()); //zero;
546 std::vector<Real> augiters;
547
548 /* Compute Cauchy step nCP. */
549 ROL::Ptr<Vector<Real> > nCP = xvec_->clone();
550 ROL::Ptr<Vector<Real> > nCPdual = gvec_->clone();
551 ROL::Ptr<Vector<Real> > nN = xvec_->clone();
552 ROL::Ptr<Vector<Real> > ctemp = cvec_->clone();
553 ROL::Ptr<Vector<Real> > dualc0 = lvec_->clone();
554 dualc0->set(c.dual());
555 con.applyAdjointJacobian(*nCPdual, *dualc0, x, zerotol);
556 nCP->set(nCPdual->dual());
557 con.applyJacobian(*ctemp, *nCP, x, zerotol);
558
559 Real normsquare_ctemp = ctemp->dot(*ctemp);
560 if (normsquare_ctemp != zero) {
561 nCP->scale( -(nCP->dot(*nCP))/normsquare_ctemp );
562 }
563
564 /* If the Cauchy step nCP is outside the trust region,
565 return the scaled Cauchy step. */
566 Real norm_nCP = nCP->norm();
567 if (norm_nCP >= delta) {
568 n.set(*nCP);
569 n.scale( delta/norm_nCP );
570 if (infoQN_) {
571 std::stringstream hist;
572 hist << " taking partial Cauchy step\n";
573 std::cout << hist.str();
574 }
575 return;
576 }
577
578 /* Compute 'Newton' step, for example, by solving a problem
579 related to finding the minimum norm solution of min || c(x_k)*s + c ||^2. */
580 // Compute tolerance for linear solver.
581 con.applyJacobian(*ctemp, *nCP, x, zerotol);
582 ctemp->plus(c);
583 Real tol = setTolOSS(qntol_*ctemp->norm());
584 // Form right-hand side.
585 ctemp->scale(-one);
586 nCPdual->set(nCP->dual());
587 nCPdual->scale(-one);
588 // Declare left-hand side of augmented system.
589 ROL::Ptr<Vector<Real> > dn = xvec_->clone();
590 ROL::Ptr<Vector<Real> > y = lvec_->clone();
591 // Solve augmented system.
592 augiters = con.solveAugmentedSystem(*dn, *y, *nCPdual, *ctemp, x, tol);
593 totalCallLS_++;
594 totalIterLS_ = totalIterLS_ + augiters.size();
595 printInfoLS(augiters);
596
597 nN->set(*dn);
598 nN->plus(*nCP);
599
600 /* Either take full or partial Newton step, depending on
601 the trust-region constraint. */
602 Real norm_nN = nN->norm();
603 if (norm_nN <= delta) {
604 // Take full feasibility step.
605 n.set(*nN);
606 if (infoQN_) {
607 std::stringstream hist;
608 hist << " taking full Newton step\n";
609 std::cout << hist.str();
610 }
611 }
612 else {
613 // Take convex combination n = nCP+tau*(nN-nCP),
614 // so that ||n|| = delta. In other words, solve
615 // scalar quadratic equation: ||nCP+tau*(nN-nCP)||^2 = delta^2.
616 Real aa = dn->dot(*dn);
617 Real bb = dn->dot(*nCP);
618 Real cc = norm_nCP*norm_nCP - delta*delta;
619 Real tau = (-bb+sqrt(bb*bb-aa*cc))/aa;
620 n.set(*nCP);
621 n.axpy(tau, *dn);
622 if (infoQN_) {
623 std::stringstream hist;
624 hist << " taking dogleg step\n";
625 std::cout << hist.str();
626 }
627 }
628
629 } // computeQuasinormalStep
630
631
646 const Vector<Real> &x, const Vector<Real> &g, const Vector<Real> &n, const Vector<Real> &l,
647 Real delta, Objective<Real> &obj, Constraint<Real> &con) {
648
649 /* Initialization of the CG step. */
650 bool orthocheck = true; // set to true if want to check orthogonality
651 // of Wr and r, otherwise set to false
652 Real tol_ortho = 0.5; // orthogonality measure; represets a bound on norm( \hat{S}, 2), where
653 // \hat{S} is defined in Heinkenschloss/Ridzal., "A Matrix-Free Trust-Region SQP Method"
654 Real S_max = 1.0; // another orthogonality measure; norm(S) needs to be bounded by
655 // a modest constant; norm(S) is small if the approximation of
656 // the null space projector is good
657 Real zero(0);
658 Real one(1);
659 Real zerotol = std::sqrt(ROL_EPSILON<Real>());
660 std::vector<Real> augiters;
661 iterCG_ = 1;
662 flagCG_ = 0;
663 t.zero();
664 tCP.zero();
665 ROL::Ptr<Vector<Real> > r = gvec_->clone();
666 ROL::Ptr<Vector<Real> > pdesc = xvec_->clone();
667 ROL::Ptr<Vector<Real> > tprev = xvec_->clone();
668 ROL::Ptr<Vector<Real> > Wr = xvec_->clone();
669 ROL::Ptr<Vector<Real> > Hp = gvec_->clone();
670 ROL::Ptr<Vector<Real> > xtemp = xvec_->clone();
671 ROL::Ptr<Vector<Real> > gtemp = gvec_->clone();
672 ROL::Ptr<Vector<Real> > ltemp = lvec_->clone();
673 ROL::Ptr<Vector<Real> > czero = cvec_->clone();
674 czero->zero();
675 r->set(g);
676 obj.hessVec(*gtemp, n, x, zerotol);
677 r->plus(*gtemp);
678 if (useConHess_) {
679 con.applyAdjointHessian(*gtemp, l, n, x, zerotol);
680 r->plus(*gtemp);
681 }
682 Real normg = r->norm();
683 Real normWg = zero;
684 Real pHp = zero;
685 Real rp = zero;
686 Real alpha = zero;
687 Real normp = zero;
688 Real normr = zero;
689 Real normt = zero;
690 std::vector<Real> normWr(maxiterCG_+1, zero);
691
692 std::vector<ROL::Ptr<Vector<Real > > > p; // stores search directions
693 std::vector<ROL::Ptr<Vector<Real > > > Hps; // stores duals of hessvec's applied to p's
694 std::vector<ROL::Ptr<Vector<Real > > > rs; // stores duals of residuals
695 std::vector<ROL::Ptr<Vector<Real > > > Wrs; // stores duals of projected residuals
696
697 Real rptol(1e-12);
698
699 if (infoTS_) {
700 std::stringstream hist;
701 hist << "\n Tangential subproblem\n";
702 hist << std::setw(6) << std::right << "iter" << std::setw(18) << "||Wr||/||Wr0||" << std::setw(15) << "||s||";
703 hist << std::setw(15) << "delta" << std::setw(15) << "||c'(x)s||" << "\n";
704 std::cout << hist.str();
705 }
706
707 if (normg == 0) {
708 if (infoTS_) {
709 std::stringstream hist;
710 hist << " >>> Tangential subproblem: Initial gradient is zero! \n";
711 std::cout << hist.str();
712 }
713 iterCG_ = 0; Wg.zero(); flagCG_ = 0;
714 return;
715 }
716
717 /* Start CG loop. */
718 while (iterCG_ < maxiterCG_) {
719
720 // Store tangential Cauchy point (which is the current iterate in the second iteration).
721 if (iterCG_ == 2) {
722 tCP.set(t);
723 }
724
725 // Compute (inexact) projection W*r.
726 if (iterCG_ == 1) {
727 // Solve augmented system.
728 Real tol = setTolOSS(pgtol_);
729 augiters = con.solveAugmentedSystem(*Wr, *ltemp, *r, *czero, x, tol);
730 totalCallLS_++;
731 totalIterLS_ = totalIterLS_ + augiters.size();
732 printInfoLS(augiters);
733
734 Wg.set(*Wr);
735 normWg = Wg.norm();
736 if (orthocheck) {
737 Wrs.push_back(xvec_->clone());
738 (Wrs[iterCG_-1])->set(*Wr);
739 }
740 // Check if done (small initial projected residual).
741 if (normWg == zero) {
742 flagCG_ = 0;
743 iterCG_--;
744 if (infoTS_) {
745 std::stringstream hist;
746 hist << " Initial projected residual is close to zero! \n";
747 std::cout << hist.str();
748 }
749 return;
750 }
751 // Set first residual to projected gradient.
752 // change r->set(Wg);
753 r->set(Wg.dual());
754 if (orthocheck) {
755 rs.push_back(xvec_->clone());
756 // change (rs[0])->set(*r);
757 (rs[0])->set(r->dual());
758 }
759 }
760 else {
761 // Solve augmented system.
762 Real tol = setTolOSS(projtol_);
763 augiters = con.solveAugmentedSystem(*Wr, *ltemp, *r, *czero, x, tol);
764 totalCallLS_++;
765 totalIterLS_ = totalIterLS_ + augiters.size();
766 printInfoLS(augiters);
767
768 if (orthocheck) {
769 Wrs.push_back(xvec_->clone());
770 (Wrs[iterCG_-1])->set(*Wr);
771 }
772 }
773
774 normWr[iterCG_-1] = Wr->norm();
775
776 if (infoTS_) {
777 ROL::Ptr<Vector<Real> > ct = cvec_->clone();
778 con.applyJacobian(*ct, t, x, zerotol);
779 Real linc = ct->norm();
780 std::stringstream hist;
781 hist << std::scientific << std::setprecision(6);
782 hist << std::setw(6) << std::right << iterCG_-1 << std::setw(18) << normWr[iterCG_-1]/normWg << std::setw(15) << t.norm();
783 hist << std::setw(15) << delta << std::setw(15) << linc << "\n";
784 std::cout << hist.str();
785 }
786
787 // Check if done (small relative residual).
788 if (normWr[iterCG_-1]/normWg < tolCG_) {
789 flagCG_ = 0;
790 iterCG_ = iterCG_-1;
791 if (infoTS_) {
792 std::stringstream hist;
793 hist << " || W(g + H*(n+s)) || <= cgtol*|| W(g + H*n)|| \n";
794 std::cout << hist.str();
795 }
796 return;
797 }
798
799 // Check nonorthogonality, one-norm of (WR*R/diag^2 - I)
800 if (orthocheck) {
801 ROL::LA::Matrix<Real> Wrr(iterCG_,iterCG_); // holds matrix Wrs'*rs
802 ROL::LA::Matrix<Real> T(iterCG_,iterCG_); // holds matrix T=(1/diag)*Wrs'*rs*(1/diag)
803 ROL::LA::Matrix<Real> Tm1(iterCG_,iterCG_); // holds matrix Tm1=T-I
804 for (int i=0; i<iterCG_; i++) {
805 for (int j=0; j<iterCG_; j++) {
806 Wrr(i,j) = (Wrs[i])->dot(*rs[j]);
807 T(i,j) = Wrr(i,j)/(normWr[i]*normWr[j]);
808 Tm1(i,j) = T(i,j);
809 if (i==j) {
810 Tm1(i,j) = Tm1(i,j) - one;
811 }
812 }
813 }
814 if (Tm1.normOne() >= tol_ortho) {
815 ROL::LAPACK<int,Real> lapack;
816 std::vector<int> ipiv(iterCG_);
817 int info;
818 std::vector<Real> work(3*iterCG_);
819 // compute inverse of T
820 lapack.GETRF(iterCG_, iterCG_, T.values(), T.stride(), &ipiv[0], &info);
821 lapack.GETRI(iterCG_, T.values(), T.stride(), &ipiv[0], &work[0], 3*iterCG_, &info);
822 Tm1 = T;
823 for (int i=0; i<iterCG_; i++) {
824 Tm1(i,i) = Tm1(i,i) - one;
825 }
826 if (Tm1.normOne() > S_max) {
827 flagCG_ = 4;
828 if (infoTS_) {
829 std::stringstream hist;
830 hist << " large nonorthogonality in W(R)'*R detected \n";
831 std::cout << hist.str();
832 }
833 return;
834 }
835 }
836 }
837
838 // Full orthogonalization.
839 p.push_back(xvec_->clone());
840 (p[iterCG_-1])->set(*Wr);
841 (p[iterCG_-1])->scale(-one);
842 for (int j=1; j<iterCG_; j++) {
843 Real scal = (p[iterCG_-1])->dot(*(Hps[j-1])) / (p[j-1])->dot(*(Hps[j-1]));
844 ROL::Ptr<Vector<Real> > pj = xvec_->clone();
845 pj->set(*p[j-1]);
846 pj->scale(-scal);
847 (p[iterCG_-1])->plus(*pj);
848 }
849
850 // change Hps.push_back(gvec_->clone());
851 Hps.push_back(xvec_->clone());
852 // change obj.hessVec(*(Hps[iterCG_-1]), *(p[iterCG_-1]), x, zerotol);
853 obj.hessVec(*Hp, *(p[iterCG_-1]), x, zerotol);
854 if (useConHess_) {
855 con.applyAdjointHessian(*gtemp, l, *(p[iterCG_-1]), x, zerotol);
856 // change (Hps[iterCG_-1])->plus(*gtemp);
857 Hp->plus(*gtemp);
858 }
859 // "Preconditioning" step.
860 (Hps[iterCG_-1])->set(Hp->dual());
861
862 pHp = (p[iterCG_-1])->dot(*(Hps[iterCG_-1]));
863 // change rp = (p[iterCG_-1])->dot(*r);
864 rp = (p[iterCG_-1])->dot(*(rs[iterCG_-1]));
865
866 normp = (p[iterCG_-1])->norm();
867 normr = r->norm();
868
869 // Negative curvature stopping condition.
870 if (pHp <= 0) {
871 pdesc->set(*(p[iterCG_-1])); // p is the descent direction
872 if ((std::abs(rp) >= rptol*normp*normr) && (sgn(rp) == 1)) {
873 pdesc->scale(-one); // -p is the descent direction
874 }
875 flagCG_ = 2;
876 Real a = pdesc->dot(*pdesc);
877 Real b = pdesc->dot(t);
878 Real c = t.dot(t) - delta*delta;
879 // Positive root of a*theta^2 + 2*b*theta + c = 0.
880 Real theta = (-b + std::sqrt(b*b - a*c)) / a;
881 xtemp->set(*(p[iterCG_-1]));
882 xtemp->scale(theta);
883 t.plus(*xtemp);
884 // Store as tangential Cauchy point if terminating in first iteration.
885 if (iterCG_ == 1) {
886 tCP.set(t);
887 }
888 if (infoTS_) {
889 std::stringstream hist;
890 hist << " negative curvature detected \n";
891 std::cout << hist.str();
892 }
893 return;
894 }
895
896 // Want to enforce nonzero alpha's.
897 if (std::abs(rp) < rptol*normp*normr) {
898 flagCG_ = 5;
899 if (infoTS_) {
900 std::stringstream hist;
901 hist << " Zero alpha due to inexactness. \n";
902 std::cout << hist.str();
903 }
904 return;
905 }
906
907 alpha = - rp/pHp;
908
909 // Iterate update.
910 tprev->set(t);
911 xtemp->set(*(p[iterCG_-1]));
912 xtemp->scale(alpha);
913 t.plus(*xtemp);
914
915 // Trust-region stopping condition.
916 normt = t.norm();
917 if (normt >= delta) {
918 pdesc->set(*(p[iterCG_-1])); // p is the descent direction
919 if (sgn(rp) == 1) {
920 pdesc->scale(-one); // -p is the descent direction
921 }
922 Real a = pdesc->dot(*pdesc);
923 Real b = pdesc->dot(*tprev);
924 Real c = tprev->dot(*tprev) - delta*delta;
925 // Positive root of a*theta^2 + 2*b*theta + c = 0.
926 Real theta = (-b + std::sqrt(b*b - a*c)) / a;
927 xtemp->set(*(p[iterCG_-1]));
928 xtemp->scale(theta);
929 t.set(*tprev);
930 t.plus(*xtemp);
931 // Store as tangential Cauchy point if terminating in first iteration.
932 if (iterCG_ == 1) {
933 tCP.set(t);
934 }
935 flagCG_ = 3;
936 if (infoTS_) {
937 std::stringstream hist;
938 hist << " trust-region condition active \n";
939 std::cout << hist.str();
940 }
941 return;
942 }
943
944 // Residual update.
945 xtemp->set(*(Hps[iterCG_-1]));
946 xtemp->scale(alpha);
947 // change r->plus(*gtemp);
948 r->plus(xtemp->dual());
949 if (orthocheck) {
950 // change rs.push_back(gvec_->clone());
951 rs.push_back(xvec_->clone());
952 // change (rs[iterCG_])->set(*r);
953 (rs[iterCG_])->set(r->dual());
954 }
955
956 iterCG_++;
957
958 } // while (iterCG_ < maxiterCG_)
959
960 flagCG_ = 1;
961 if (infoTS_) {
962 std::stringstream hist;
963 hist << " maximum number of iterations reached \n";
964 std::cout << hist.str();
965 }
966
967 } // solveTangentialSubproblem
968
969
972 void accept(Vector<Real> &s, Vector<Real> &n, Vector<Real> &t, Real f_new, Vector<Real> &c_new,
973 Vector<Real> &gf_new, Vector<Real> &l_new, Vector<Real> &g_new,
974 const Vector<Real> &x, const Vector<Real> &l, Real f, const Vector<Real> &gf, const Vector<Real> &c,
975 const Vector<Real> &g, Vector<Real> &tCP, Vector<Real> &Wg,
976 Objective<Real> &obj, Constraint<Real> &con, AlgorithmState<Real> &algo_state) {
977
978 Real beta = 1e-8; // predicted reduction parameter
979 Real tol_red_tang = 1e-3; // internal reduction factor for tangtol
980 Real tol_red_all = 1e-1; // internal reduction factor for qntol, lmhtol, pgtol, projtol, tangtol
981 //bool glob_refine = true; // true - if subsolver tolerances are adjusted in this routine, keep adjusted values globally
982 // false - if subsolver tolerances are adjusted in this routine, discard adjusted values
983 Real tol_fdiff = 1e-12; // relative objective function difference for ared computation
984 int ct_max = 10; // maximum number of globalization tries
985 Real mintol = 1e-16; // smallest projection tolerance value
986
987 // Determines max value of |rpred|/pred.
988 Real rpred_over_pred = 0.5*(1-eta_);
989
990 if (infoAC_) {
991 std::stringstream hist;
992 hist << "\n Composite step acceptance\n";
993 std::cout << hist.str();
994 }
995
996 Real zero = 0.0;
997 Real one = 1.0;
998 Real two = 2.0;
999 Real half = one/two;
1000 Real zerotol = std::sqrt(ROL_EPSILON<Real>());
1001 std::vector<Real> augiters;
1002
1003 Real pred = zero;
1004 Real ared = zero;
1005 Real rpred = zero;
1006 Real part_pred = zero;
1007 Real linc_preproj = zero;
1008 Real linc_postproj = zero;
1009 Real tangtol_start = zero;
1010 Real tangtol = tangtol_;
1011 //Real projtol = projtol_;
1012 bool flag = false;
1013 int num_proj = 0;
1014 bool try_tCP = false;
1015 Real fdiff = zero;
1016
1017 ROL::Ptr<Vector<Real> > xtrial = xvec_->clone();
1018 ROL::Ptr<Vector<Real> > Jl = gvec_->clone();
1019 ROL::Ptr<Vector<Real> > gfJl = gvec_->clone();
1020 ROL::Ptr<Vector<Real> > Jnc = cvec_->clone();
1021 ROL::Ptr<Vector<Real> > t_orig = xvec_->clone();
1022 ROL::Ptr<Vector<Real> > t_dual = gvec_->clone();
1023 ROL::Ptr<Vector<Real> > Jt_orig = cvec_->clone();
1024 ROL::Ptr<Vector<Real> > t_m_tCP = xvec_->clone();
1025 ROL::Ptr<Vector<Real> > ltemp = lvec_->clone();
1026 ROL::Ptr<Vector<Real> > xtemp = xvec_->clone();
1027 ROL::Ptr<Vector<Real> > rt = cvec_->clone();
1028 ROL::Ptr<Vector<Real> > Hn = gvec_->clone();
1029 ROL::Ptr<Vector<Real> > Hto = gvec_->clone();
1030 ROL::Ptr<Vector<Real> > cxxvec = gvec_->clone();
1031 ROL::Ptr<Vector<Real> > czero = cvec_->clone();
1032 czero->zero();
1033 Real Jnc_normsquared = zero;
1034 Real c_normsquared = zero;
1035
1036 // Compute and store some quantities for later use. Necessary
1037 // because of the function and constraint updates below.
1038 con.applyAdjointJacobian(*Jl, l, x, zerotol);
1039 con.applyJacobian(*Jnc, n, x, zerotol);
1040 Jnc->plus(c);
1041 Jnc_normsquared = Jnc->dot(*Jnc);
1042 c_normsquared = c.dot(c);
1043
1044 for (int ct=0; ct<ct_max; ct++) {
1045
1046 try_tCP = true;
1047 t_m_tCP->set(t);
1048 t_m_tCP->scale(-one);
1049 t_m_tCP->plus(tCP);
1050 if (t_m_tCP->norm() == zero) {
1051 try_tCP = false;
1052 }
1053
1054 t_orig->set(t);
1055 con.applyJacobian(*Jt_orig, *t_orig, x, zerotol);
1056 linc_preproj = Jt_orig->norm();
1057 pred = one;
1058 rpred = two*rpred_over_pred*pred;
1059 flag = false;
1060 num_proj = 1;
1061 tangtol_start = tangtol;
1062
1063 while (std::abs(rpred)/pred > rpred_over_pred) {
1064 // Compute projected tangential step.
1065 if (flag) {
1066 tangtol = tol_red_tang*tangtol;
1067 num_proj++;
1068 if (tangtol < mintol) {
1069 if (infoAC_) {
1070 std::stringstream hist;
1071 hist << "\n The projection of the tangential step cannot be done with sufficient precision.\n";
1072 hist << " Is the quasi-normal step very small? Continuing with no global convergence guarantees.\n";
1073 std::cout << hist.str();
1074 }
1075 break;
1076 }
1077 }
1078 // Solve augmented system.
1079 Real tol = setTolOSS(tangtol);
1080 // change augiters = con.solveAugmentedSystem(t, *ltemp, *t_orig, *czero, x, tol);
1081 t_dual->set(t_orig->dual());
1082 augiters = con.solveAugmentedSystem(t, *ltemp, *t_dual, *czero, x, tol);
1083 totalCallLS_++;
1084 totalIterLS_ = totalIterLS_ + augiters.size();
1085 printInfoLS(augiters);
1086 totalProj_++;
1087 con.applyJacobian(*rt, t, x, zerotol);
1088 linc_postproj = rt->norm();
1089
1090 // Compute composite step.
1091 s.set(t);
1092 s.plus(n);
1093
1094 // Compute some quantities before updating the objective and the constraint.
1095 obj.hessVec(*Hn, n, x, zerotol);
1096 if (useConHess_) {
1097 con.applyAdjointHessian(*cxxvec, l, n, x, zerotol);
1098 Hn->plus(*cxxvec);
1099 }
1100 obj.hessVec(*Hto, *t_orig, x, zerotol);
1101 if (useConHess_) {
1102 con.applyAdjointHessian(*cxxvec, l, *t_orig, x, zerotol);
1103 Hto->plus(*cxxvec);
1104 }
1105
1106 // Compute objective, constraint, etc. values at the trial point.
1107 xtrial->set(x);
1108 xtrial->plus(s);
1109 obj.update(*xtrial,false,algo_state.iter);
1110 con.update(*xtrial,false,algo_state.iter);
1111 f_new = obj.value(*xtrial, zerotol);
1112 obj.gradient(gf_new, *xtrial, zerotol);
1113 con.value(c_new, *xtrial, zerotol);
1114 l_new.set(l);
1115 computeLagrangeMultiplier(l_new, *xtrial, gf_new, con);
1116
1117 // Penalty parameter update.
1118 part_pred = - Wg.dot(*t_orig);
1119 gfJl->set(gf);
1120 gfJl->plus(*Jl);
1121 // change part_pred -= gfJl->dot(n);
1122 part_pred -= n.dot(gfJl->dual());
1123 // change part_pred -= half*Hn->dot(n);
1124 part_pred -= half*n.dot(Hn->dual());
1125 // change part_pred -= half*Hto->dot(*t_orig);
1126 part_pred -= half*t_orig->dot(Hto->dual());
1127 ltemp->set(l_new);
1128 ltemp->axpy(-one, l);
1129 // change part_pred -= Jnc->dot(*ltemp);
1130 part_pred -= Jnc->dot(ltemp->dual());
1131
1132 if ( part_pred < -half*penalty_*(c_normsquared-Jnc_normsquared) ) {
1133 penalty_ = ( -two * part_pred / (c_normsquared-Jnc_normsquared) ) + beta;
1134 }
1135
1136 pred = part_pred + penalty_*(c_normsquared-Jnc_normsquared);
1137
1138 // Computation of rpred.
1139 // change rpred = - ltemp->dot(*rt) - penalty_ * rt->dot(*rt) - two * penalty_ * rt->dot(*Jnc);
1140 rpred = - rt->dot(ltemp->dual()) - penalty_ * rt->dot(*rt) - two * penalty_ * rt->dot(*Jnc);
1141 // change ROL::Ptr<Vector<Real> > lrt = lvec_->clone();
1142 //lrt->set(*rt);
1143 //rpred = - ltemp->dot(*rt) - penalty_ * std::pow(rt->norm(), 2) - two * penalty_ * lrt->dot(*Jnc);
1144 flag = 1;
1145
1146 } // while (std::abs(rpred)/pred > rpred_over_pred)
1147
1148 tangtol = tangtol_start;
1149
1150 // Check if the solution of the tangential subproblem is
1151 // disproportionally large compared to total trial step.
1152 xtemp->set(n);
1153 xtemp->plus(t);
1154 if ( t_orig->norm()/xtemp->norm() < tntmax_ ) {
1155 break;
1156 }
1157 else {
1158 t_m_tCP->set(*t_orig);
1159 t_m_tCP->scale(-one);
1160 t_m_tCP->plus(tCP);
1161 if ((t_m_tCP->norm() > 0) && try_tCP) {
1162 if (infoAC_) {
1163 std::stringstream hist;
1164 hist << " ---> now trying tangential Cauchy point\n";
1165 std::cout << hist.str();
1166 }
1167 t.set(tCP);
1168 }
1169 else {
1170 if (infoAC_) {
1171 std::stringstream hist;
1172 hist << " ---> recomputing quasi-normal step and re-solving tangential subproblem\n";
1173 std::cout << hist.str();
1174 }
1175 totalRef_++;
1176 // Reset global quantities.
1177 obj.update(x);
1178 con.update(x);
1179 /*lmhtol = tol_red_all*lmhtol;
1180 qntol = tol_red_all*qntol;
1181 pgtol = tol_red_all*pgtol;
1182 projtol = tol_red_all*projtol;
1183 tangtol = tol_red_all*tangtol;
1184 if (glob_refine) {
1185 lmhtol_ = lmhtol;
1186 qntol_ = qntol;
1187 pgtol_ = pgtol;
1188 projtol_ = projtol;
1189 tangtol_ = tangtol;
1190 }*/
1191 if (!tolOSSfixed_) {
1192 lmhtol_ *= tol_red_all;
1193 qntol_ *= tol_red_all;
1194 pgtol_ *= tol_red_all;
1195 projtol_ *= tol_red_all;
1196 tangtol_ *= tol_red_all;
1197 }
1198 // Recompute the quasi-normal step.
1199 computeQuasinormalStep(n, c, x, zeta_*Delta_, con);
1200 // Solve tangential subproblem.
1201 solveTangentialSubproblem(t, tCP, Wg, x, g, n, l, Delta_, obj, con);
1203 if (flagCG_ == 1) {
1204 totalNegCurv_++;
1205 }
1206 }
1207 } // else w.r.t. ( t_orig->norm()/xtemp->norm() < tntmax )
1208
1209 } // for (int ct=0; ct<ct_max; ct++)
1210
1211 // Compute actual reduction;
1212 fdiff = f - f_new;
1213 // Heuristic 1: If fdiff is very small compared to f, set it to 0,
1214 // in order to prevent machine precision issues.
1215 Real em24(1e-24);
1216 Real em14(1e-14);
1217 if (std::abs(fdiff / (f+em24)) < tol_fdiff) {
1218 fdiff = em14;
1219 }
1220 // change ared = fdiff + (l.dot(c) - l_new.dot(c_new)) + penalty_*(c.dot(c) - c_new.dot(c_new));
1221 // change ared = fdiff + (l.dot(c) - l_new.dot(c_new)) + penalty_*(std::pow(c.norm(),2) - std::pow(c_new.norm(),2));
1222 ared = fdiff + (c.dot(l.dual()) - c_new.dot(l_new.dual())) + penalty_*(c.dot(c) - c_new.dot(c_new));
1223
1224 // Store actual and predicted reduction.
1225 ared_ = ared;
1226 pred_ = pred;
1227
1228 // Store step and vector norms.
1229 snorm_ = s.norm();
1230 nnorm_ = n.norm();
1231 tnorm_ = t.norm();
1232
1233 // Print diagnostics.
1234 if (infoAC_) {
1235 std::stringstream hist;
1236 hist << "\n Trial step info ...\n";
1237 hist << " n_norm = " << nnorm_ << "\n";
1238 hist << " t_norm = " << tnorm_ << "\n";
1239 hist << " s_norm = " << snorm_ << "\n";
1240 hist << " xtrial_norm = " << xtrial->norm() << "\n";
1241 hist << " f_old = " << f << "\n";
1242 hist << " f_trial = " << f_new << "\n";
1243 hist << " f_old-f_trial = " << f-f_new << "\n";
1244 hist << " ||c_old|| = " << c.norm() << "\n";
1245 hist << " ||c_trial|| = " << c_new.norm() << "\n";
1246 hist << " ||Jac*t_preproj|| = " << linc_preproj << "\n";
1247 hist << " ||Jac*t_postproj|| = " << linc_postproj << "\n";
1248 hist << " ||t_tilde||/||t|| = " << t_orig->norm() / t.norm() << "\n";
1249 hist << " ||t_tilde||/||n+t|| = " << t_orig->norm() / snorm_ << "\n";
1250 hist << " # projections = " << num_proj << "\n";
1251 hist << " penalty param = " << penalty_ << "\n";
1252 hist << " ared = " << ared_ << "\n";
1253 hist << " pred = " << pred_ << "\n";
1254 hist << " ared/pred = " << ared_/pred_ << "\n";
1255 std::cout << hist.str();
1256 }
1257
1258 } // accept
1259
1260}; // class CompositeStep
1261
1262} // namespace ROL
1263
1264#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:68
ROL::Ptr< StepState< Real > > getState(void)
Definition ROL_Step.hpp:73
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.