ROL
ROL_TypeB_LinMoreAlgorithm_Def.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_TYPEB_LINMOREALGORITHM_DEF_HPP
45#define ROL_TYPEB_LINMOREALGORITHM_DEF_HPP
46
47namespace ROL {
48namespace TypeB {
49
50template<typename Real>
52 const Ptr<Secant<Real>> &secant) {
53 // Set status test
54 status_->reset();
55 status_->add(makePtr<StatusTest<Real>>(list));
56
57 ParameterList &trlist = list.sublist("Step").sublist("Trust Region");
58 // Trust-Region Parameters
59 state_->searchSize = trlist.get("Initial Radius", -1.0);
60 delMax_ = trlist.get("Maximum Radius", ROL_INF<Real>());
61 eta0_ = trlist.get("Step Acceptance Threshold", 0.05);
62 eta1_ = trlist.get("Radius Shrinking Threshold", 0.05);
63 eta2_ = trlist.get("Radius Growing Threshold", 0.9);
64 gamma0_ = trlist.get("Radius Shrinking Rate (Negative rho)", 0.0625);
65 gamma1_ = trlist.get("Radius Shrinking Rate (Positive rho)", 0.25);
66 gamma2_ = trlist.get("Radius Growing Rate", 2.5);
67 TRsafe_ = trlist.get("Safeguard Size", 100.0);
68 eps_ = TRsafe_*ROL_EPSILON<Real>();
69 interpRad_ = trlist.get("Use Radius Interpolation", false);
70 // Nonmonotone Parameters
71 storageNM_ = trlist.get("Nonmonotone Storage Size", 0);
72 useNM_ = (storageNM_ <= 0 ? false : true);
73 // Krylov Parameters
74 maxit_ = list.sublist("General").sublist("Krylov").get("Iteration Limit", 20);
75 tol1_ = list.sublist("General").sublist("Krylov").get("Absolute Tolerance", 1e-4);
76 tol2_ = list.sublist("General").sublist("Krylov").get("Relative Tolerance", 1e-2);
77 // Algorithm-Specific Parameters
78 ROL::ParameterList &lmlist = trlist.sublist("Lin-More");
79 minit_ = lmlist.get("Maximum Number of Minor Iterations", 10);
80 mu0_ = lmlist.get("Sufficient Decrease Parameter", 1e-2);
81 spexp_ = lmlist.get("Relative Tolerance Exponent", 1.0);
82 spexp_ = std::max(static_cast<Real>(1),std::min(spexp_,static_cast<Real>(2)));
83 redlim_ = lmlist.sublist("Cauchy Point").get("Maximum Number of Reduction Steps", 10);
84 explim_ = lmlist.sublist("Cauchy Point").get("Maximum Number of Expansion Steps", 10);
85 alpha_ = lmlist.sublist("Cauchy Point").get("Initial Step Size", 1.0);
86 normAlpha_ = lmlist.sublist("Cauchy Point").get("Normalize Initial Step Size", false);
87 interpf_ = lmlist.sublist("Cauchy Point").get("Reduction Rate", 0.1);
88 extrapf_ = lmlist.sublist("Cauchy Point").get("Expansion Rate", 10.0);
89 qtol_ = lmlist.sublist("Cauchy Point").get("Decrease Tolerance", 1e-8);
90 interpfPS_ = lmlist.sublist("Projected Search").get("Backtracking Rate", 0.5);
91 pslim_ = lmlist.sublist("Projected Search").get("Maximum Number of Steps", 20);
92 // Output Parameters
93 verbosity_ = list.sublist("General").get("Output Level",0);
94 writeHeader_ = verbosity_ > 2;
95 // Secant Information
96 useSecantPrecond_ = list.sublist("General").sublist("Secant").get("Use as Preconditioner", false);
97 useSecantHessVec_ = list.sublist("General").sublist("Secant").get("Use as Hessian", false);
99 if (useSecantPrecond_ && !useSecantHessVec_) mode = SECANTMODE_INVERSE;
100 else if (useSecantHessVec_ && !useSecantPrecond_) mode = SECANTMODE_FORWARD;
101 // Initialize trust region model
102 model_ = makePtr<TrustRegionModel_U<Real>>(list,secant,mode);
103 if (secant == nullPtr) {
104 esec_ = StringToESecant(list.sublist("General").sublist("Secant").get("Type","Limited-Memory BFGS"));
105 }
106}
107
108template<typename Real>
110 const Vector<Real> &g,
111 Objective<Real> &obj,
113 std::ostream &outStream) {
114 const Real one(1);
115 hasEcon_ = true;
116 if (proj_ == nullPtr) {
117 proj_ = makePtr<PolyhedralProjection<Real>>(makePtrFromRef(bnd));
118 hasEcon_ = false;
119 }
120 // Initialize data
122 nhess_ = 0;
123 // Update approximate gradient and approximate objective function.
124 Real ftol = static_cast<Real>(0.1)*ROL_OVERFLOW<Real>();
125 proj_->project(x,outStream); state_->nproj++;
126 state_->iterateVec->set(x);
127 obj.update(x,UpdateType::Initial,state_->iter);
128 state_->value = obj.value(x,ftol);
129 state_->nfval++;
130 obj.gradient(*state_->gradientVec,x,ftol);
131 state_->ngrad++;
132 state_->stepVec->set(x);
133 state_->stepVec->axpy(-one,state_->gradientVec->dual());
134 proj_->project(*state_->stepVec,outStream); state_->nproj++;
135 state_->stepVec->axpy(-one,x);
136 state_->gnorm = state_->stepVec->norm();
137 state_->snorm = ROL_INF<Real>();
138 // Normalize initial CP step length
139 if (normAlpha_) {
140 alpha_ /= state_->gradientVec->norm();
141 }
142 // Compute initial trust region radius if desired.
143 if ( state_->searchSize <= static_cast<Real>(0) ) {
144 state_->searchSize = state_->gradientVec->norm();
145 }
146 // Initialize null space projection
147 if (hasEcon_) {
148 rcon_ = makePtr<ReducedLinearConstraint<Real>>(proj_->getLinearConstraint(),
149 makePtrFromRef(bnd),
150 makePtrFromRef(x));
151 ns_ = makePtr<NullSpaceOperator<Real>>(rcon_,x,
152 *proj_->getResidual());
153 }
154}
155
156template<typename Real>
158 const Vector<Real> &g,
159 Objective<Real> &obj,
161 std::ostream &outStream ) {
162 const Real zero(0);
163 Real tol0 = std::sqrt(ROL_EPSILON<Real>());
164 Real gfnorm(0), gfnormf(0), tol(0), stol(0), snorm(0);
165 Real ftrial(0), pRed(0), rho(1), q(0), delta(0);
166 int flagCG(0), iterCG(0), maxit(0);
167 // Initialize trust-region data
168 initialize(x,g,obj,bnd,outStream);
169 Ptr<Vector<Real>> s = x.clone();
170 Ptr<Vector<Real>> gmod = g.clone(), gfree = g.clone();
171 Ptr<Vector<Real>> pwa1 = x.clone(), pwa2 = x.clone(), pwa3 = x.clone();
172 Ptr<Vector<Real>> dwa1 = g.clone(), dwa2 = g.clone(), dwa3 = g.clone();
173 // Initialize nonmonotone data
174 Real rhoNM(0), sigmac(0), sigmar(0);
175 Real fr(state_->value), fc(state_->value), fmin(state_->value);
176 TRUtils::ETRFlag TRflagNM;
177 int L(0);
178
179 // Output
180 if (verbosity_ > 0) writeOutput(outStream,true);
181
182 while (status_->check(*state_)) {
183 // Build trust-region model
184 model_->setData(obj,*state_->iterateVec,*state_->gradientVec);
185
186 /**** SOLVE TRUST-REGION SUBPROBLEM ****/
187 // Compute Cauchy point (TRON notation: x = x[1])
188 snorm = dcauchy(*state_->stepVec,alpha_,q,*state_->iterateVec,
189 state_->gradientVec->dual(),state_->searchSize,
190 *model_,*dwa1,*dwa2,outStream); // Solve 1D optimization problem for alpha
191 x.plus(*state_->stepVec); // Set x = x[0] + alpha*g
192 state_->snorm = snorm;
193 delta = state_->searchSize - snorm;
194 pRed = -q;
195
196 // Model gradient at s = x[1] - x[0]
197 gmod->set(*dwa1); // hessVec from Cauchy point computation
198 gmod->plus(*state_->gradientVec);
199 gfree->set(*gmod);
200 //bnd.pruneActive(*gfree,x,zero);
201 pwa1->set(gfree->dual());
202 bnd.pruneActive(*pwa1,x,zero);
203 gfree->set(pwa1->dual());
204 if (hasEcon_) {
205 applyFreePrecond(*pwa1,*gfree,x,*model_,bnd,tol0,*dwa1,*pwa2);
206 gfnorm = pwa1->norm();
207 }
208 else {
209 gfnorm = gfree->norm();
210 }
211 SPiter_ = 0; SPflag_ = 0;
212 if (verbosity_ > 1) {
213 outStream << " Norm of free gradient components: " << gfnorm << std::endl;
214 outStream << std::endl;
215 }
216
217 // Trust-region subproblem solve loop
218 maxit = maxit_;
219 for (int i = 0; i < minit_; ++i) {
220 // Run Truncated CG
221 flagCG = 0; iterCG = 0;
222 gfnormf = zero;
223 tol = std::min(tol1_,tol2_*std::pow(gfnorm,spexp_));
224 stol = tol; //zero;
225 if (gfnorm > zero && delta > zero) {
226 snorm = dtrpcg(*s,flagCG,iterCG,*gfree,x,
227 delta,*model_,bnd,tol,stol,maxit,
228 *pwa1,*dwa1,*pwa2,*dwa2,*pwa3,*dwa3,outStream);
229 maxit -= iterCG;
230 SPiter_ += iterCG;
231 if (verbosity_ > 1) {
232 outStream << " Computation of CG step" << std::endl;
233 outStream << " Current face (i): " << i << std::endl;
234 outStream << " CG step length: " << snorm << std::endl;
235 outStream << " Number of CG iterations: " << iterCG << std::endl;
236 outStream << " CG flag: " << flagCG << std::endl;
237 outStream << " Total number of iterations: " << SPiter_ << std::endl;
238 outStream << std::endl;
239 }
240
241 // Projected search
242 snorm = dprsrch(x,*s,q,gmod->dual(),*model_,bnd,*pwa1,*dwa1,outStream);
243 pRed += -q;
244
245 // Model gradient at s = (x[i+1]-x[i]) - (x[i]-x[0])
246 state_->stepVec->plus(*s);
247 state_->snorm = state_->stepVec->norm();
248 delta = state_->searchSize - state_->snorm;
249 gmod->plus(*dwa1); // gmod = H(x[i+1]-x[i]) + H(x[i]-x[0]) + g
250 gfree->set(*gmod);
251 //bnd.pruneActive(*gfree,x,zero);
252 pwa1->set(gfree->dual());
253 bnd.pruneActive(*pwa1,x,zero);
254 gfree->set(pwa1->dual());
255 if (hasEcon_) {
256 applyFreePrecond(*pwa1,*gfree,x,*model_,bnd,tol0,*dwa1,*pwa2);
257 gfnormf = pwa1->norm();
258 }
259 else {
260 gfnormf = gfree->norm();
261 }
262 if (verbosity_ > 1) {
263 outStream << " Norm of free gradient components: " << gfnormf << std::endl;
264 outStream << std::endl;
265 }
266 }
267
268 // Termination check
269 if (gfnormf <= tol) {
270 SPflag_ = 0;
271 break;
272 }
273 else if (SPiter_ >= maxit_) {
274 SPflag_ = 1;
275 break;
276 }
277 else if (flagCG == 2) {
278 SPflag_ = 2;
279 break;
280 }
281 else if (delta <= zero) {
282 //else if (flagCG == 3 || delta <= zero) {
283 SPflag_ = 3;
284 break;
285 }
286
287 // Update free gradient norm
288 gfnorm = gfnormf;
289 }
290
291 // Compute trial objective value
293 ftrial = obj.value(x,tol0);
294 state_->nfval++;
295
296 // Compute ratio of acutal and predicted reduction
297 TRflag_ = TRUtils::SUCCESS;
298 TRUtils::analyzeRatio<Real>(rho,TRflag_,state_->value,ftrial,pRed,eps_,outStream,verbosity_>1);
299 if (useNM_) {
300 TRUtils::analyzeRatio<Real>(rhoNM,TRflagNM,fr,ftrial,pRed+sigmar,eps_,outStream,verbosity_>1);
301 TRflag_ = (rho < rhoNM ? TRflagNM : TRflag_);
302 rho = (rho < rhoNM ? rhoNM : rho );
303 }
304
305 // Update algorithm state
306 state_->iter++;
307 // Accept/reject step and update trust region radius
308 if ((rho < eta0_ && TRflag_ == TRUtils::SUCCESS) || (TRflag_ >= 2)) { // Step Rejected
309 x.set(*state_->iterateVec);
310 obj.update(x,UpdateType::Revert,state_->iter);
311 if (interpRad_ && (rho < zero && TRflag_ != TRUtils::TRNAN)) {
312 // Negative reduction, interpolate to find new trust-region radius
313 state_->searchSize = TRUtils::interpolateRadius<Real>(*state_->gradientVec,*state_->stepVec,
314 state_->snorm,pRed,state_->value,ftrial,state_->searchSize,gamma0_,gamma1_,eta2_,
315 outStream,verbosity_>1);
316 }
317 else { // Shrink trust-region radius
318 state_->searchSize = gamma1_*std::min(state_->snorm,state_->searchSize);
319 }
320 }
321 else if ((rho >= eta0_ && TRflag_ != TRUtils::NPOSPREDNEG)
322 || (TRflag_ == TRUtils::POSPREDNEG)) { // Step Accepted
323 state_->value = ftrial;
324 obj.update(x,UpdateType::Accept,state_->iter);
325 if (useNM_) {
326 sigmac += pRed; sigmar += pRed;
327 if (ftrial < fmin) { fmin = ftrial; fc = fmin; sigmac = zero; L = 0; }
328 else {
329 L++;
330 if (ftrial > fc) { fc = ftrial; sigmac = zero; }
331 if (L == storageNM_) { fr = fc; sigmar = sigmac; }
332 }
333 }
334 // Increase trust-region radius
335 if (rho >= eta2_) state_->searchSize = std::min(gamma2_*state_->searchSize, delMax_);
336 // Compute gradient at new iterate
337 dwa1->set(*state_->gradientVec);
338 obj.gradient(*state_->gradientVec,x,tol0);
339 state_->ngrad++;
340 state_->gnorm = TypeB::Algorithm<Real>::optimalityCriterion(x,*state_->gradientVec,*pwa1,outStream);
341 state_->iterateVec->set(x);
342 // Update secant information in trust-region model
343 model_->update(x,*state_->stepVec,*dwa1,*state_->gradientVec,
344 state_->snorm,state_->iter);
345 }
346
347 // Update Output
348 if (verbosity_ > 0) writeOutput(outStream,writeHeader_);
349 }
350 if (verbosity_ > 0) TypeB::Algorithm<Real>::writeExitStatus(outStream);
351}
352
353template<typename Real>
355 const Vector<Real> &x, const Real alpha,
356 std::ostream &outStream) const {
357 s.set(x); s.axpy(alpha,w);
358 proj_->project(s,outStream); state_->nproj++;
359 s.axpy(static_cast<Real>(-1),x);
360 return s.norm();
361}
362
363template<typename Real>
365 Real &alpha,
366 Real &q,
367 const Vector<Real> &x,
368 const Vector<Real> &g,
369 const Real del,
371 Vector<Real> &dwa, Vector<Real> &dwa1,
372 std::ostream &outStream) {
373 const Real half(0.5);
374 // const Real zero(0); // Unused
375 Real tol = std::sqrt(ROL_EPSILON<Real>());
376 bool interp = false;
377 Real gs(0), snorm(0);
378 // Compute s = P(x[0] - alpha g[0])
379 snorm = dgpstep(s,g,x,-alpha,outStream);
380 if (snorm > del) {
381 interp = true;
382 }
383 else {
384 model.hessVec(dwa,s,x,tol); nhess_++;
385 gs = s.dot(g);
386 //q = half * s.dot(dwa.dual()) + gs;
387 q = half * s.apply(dwa) + gs;
388 interp = (q > mu0_*gs);
389 }
390 // Either increase or decrease alpha to find approximate Cauchy point
391 int cnt = 0;
392 if (interp) {
393 bool search = true;
394 while (search) {
395 alpha *= interpf_;
396 snorm = dgpstep(s,g,x,-alpha,outStream);
397 if (snorm <= del) {
398 model.hessVec(dwa,s,x,tol); nhess_++;
399 gs = s.dot(g);
400 //q = half * s.dot(dwa.dual()) + gs;
401 q = half * s.apply(dwa) + gs;
402 search = (q > mu0_*gs) && (cnt < redlim_);
403 }
404 cnt++;
405 }
406 }
407 else {
408 bool search = true;
409 Real alphas = alpha;
410 Real qs = q;
411 dwa1.set(dwa);
412 while (search) {
413 alpha *= extrapf_;
414 snorm = dgpstep(s,g,x,-alpha,outStream);
415 if (snorm <= del && cnt < explim_) {
416 model.hessVec(dwa,s,x,tol); nhess_++;
417 gs = s.dot(g);
418 //q = half * s.dot(dwa.dual()) + gs;
419 q = half * s.apply(dwa) + gs;
420 if (q <= mu0_*gs && std::abs(q-qs) > qtol_*std::abs(qs)) {
421 dwa1.set(dwa);
422 search = true;
423 alphas = alpha;
424 qs = q;
425 }
426 else {
427 q = qs;
428 dwa.set(dwa1);
429 search = false;
430 }
431 }
432 else {
433 search = false;
434 }
435 cnt++;
436 }
437 alpha = alphas;
438 snorm = dgpstep(s,g,x,-alpha,outStream);
439 }
440 if (verbosity_ > 1) {
441 outStream << " Cauchy point" << std::endl;
442 outStream << " Step length (alpha): " << alpha << std::endl;
443 outStream << " Step length (alpha*g): " << snorm << std::endl;
444 outStream << " Model decrease (pRed): " << -q << std::endl;
445 if (!interp) {
446 outStream << " Number of extrapolation steps: " << cnt << std::endl;
447 }
448 }
449 return snorm;
450}
451
452template<typename Real>
454 Real &q, const Vector<Real> &g,
457 Vector<Real> &pwa, Vector<Real> &dwa,
458 std::ostream &outStream) {
459 const Real half(0.5);
460 Real tol = std::sqrt(ROL_EPSILON<Real>());
461 Real beta(1), snorm(0), gs(0);
462 int nsteps = 0;
463 // Reduce beta until sufficient decrease is satisfied
464 bool search = true;
465 while (search) {
466 nsteps++;
467 snorm = dgpstep(pwa,s,x,beta,outStream);
468 model.hessVec(dwa,pwa,x,tol); nhess_++;
469 gs = pwa.dot(g);
470 //q = half * pwa.dot(dwa.dual()) + gs;
471 q = half * pwa.apply(dwa) + gs;
472 if (q <= mu0_*gs || nsteps > pslim_) {
473 search = false;
474 }
475 else {
476 beta *= interpfPS_;
477 }
478 }
479 s.set(pwa);
480 x.plus(s);
481 if (verbosity_ > 1) {
482 outStream << std::endl;
483 outStream << " Projected search" << std::endl;
484 outStream << " Step length (beta): " << beta << std::endl;
485 outStream << " Step length (beta*s): " << snorm << std::endl;
486 outStream << " Model decrease (pRed): " << -q << std::endl;
487 outStream << " Number of steps: " << nsteps << std::endl;
488 }
489 return snorm;
490}
491
492template<typename Real>
494 const Real ptp,
495 const Real ptx,
496 const Real del) const {
497 const Real zero(0);
498 Real dsq = del*del;
499 Real rad = ptx*ptx + ptp*(dsq-xtx);
500 rad = std::sqrt(std::max(rad,zero));
501 Real sigma(0);
502 if (ptx > zero) {
503 sigma = (dsq-xtx)/(ptx+rad);
504 }
505 else if (rad > zero) {
506 sigma = (rad-ptx)/ptp;
507 }
508 else {
509 sigma = zero;
510 }
511 return sigma;
512}
513
514template<typename Real>
515Real LinMoreAlgorithm<Real>::dtrpcg(Vector<Real> &w, int &iflag, int &iter,
516 const Vector<Real> &g, const Vector<Real> &x,
517 const Real del, TrustRegionModel_U<Real> &model,
519 const Real tol, const Real stol, const int itermax,
522 std::ostream &outStream) const {
523 // p = step (primal)
524 // q = hessian applied to step p (dual)
525 // t = gradient (dual)
526 // r = preconditioned gradient (primal)
527 Real tol0 = std::sqrt(ROL_EPSILON<Real>());
528 const Real zero(0), one(1), two(2);
529 Real rho(0), kappa(0), beta(0), sigma(0), alpha(0);
530 Real rtr(0), tnorm(0), sMs(0), pMp(0), sMp(0);
531 iter = 0; iflag = 0;
532 // Initialize step
533 w.zero();
534 // Compute residual
535 t.set(g); t.scale(-one);
536 // Preconditioned residual
537 applyFreePrecond(r,t,x,model,bnd,tol0,dwa,pwa);
538 //rho = r.dot(t.dual());
539 rho = r.apply(t);
540 // Initialize direction
541 p.set(r);
542 pMp = (!hasEcon_ ? rho : p.dot(p)); // If no equality constraint, used preconditioned norm
543 // Iterate CG
544 for (iter = 0; iter < itermax; ++iter) {
545 // Apply Hessian to direction dir
546 applyFreeHessian(q,p,x,model,bnd,tol0,pwa);
547 // Compute sigma such that ||s+sigma*dir|| = del
548 //kappa = p.dot(q.dual());
549 kappa = p.apply(q);
550 alpha = (kappa>zero) ? rho/kappa : zero;
551 sigma = dtrqsol(sMs,pMp,sMp,del);
552 // Check for negative curvature or if iterate exceeds trust region
553 if (kappa <= zero || alpha >= sigma) {
554 w.axpy(sigma,p);
555 sMs = del*del;
556 iflag = (kappa<=zero) ? 2 : 3;
557 break;
558 }
559 // Update iterate and residuals
560 w.axpy(alpha,p);
561 t.axpy(-alpha,q);
562 applyFreePrecond(r,t,x,model,bnd,tol0,dwa,pwa);
563 // Exit if residual tolerance is met
564 //rtr = r.dot(t.dual());
565 rtr = r.apply(t);
566 tnorm = t.norm();
567 if (rtr <= stol*stol || tnorm <= tol) {
568 sMs = sMs + two*alpha*sMp + alpha*alpha*pMp;
569 iflag = 0;
570 break;
571 }
572 // Compute p = r + beta * p
573 beta = rtr/rho;
574 p.scale(beta); p.plus(r);
575 rho = rtr;
576 // Update dot products
577 // sMs = <s, inv(M)s> or <s, s> if equality constraint present
578 // sMp = <s, inv(M)p> or <s, p> if equality constraint present
579 // pMp = <p, inv(M)p> or <p, p> if equality constraint present
580 sMs = sMs + two*alpha*sMp + alpha*alpha*pMp;
581 sMp = beta*(sMp + alpha*pMp);
582 pMp = (!hasEcon_ ? rho : p.dot(p)) + beta*beta*pMp;
583 }
584 // Check iteration count
585 if (iter == itermax) {
586 iflag = 1;
587 }
588 if (iflag != 1) {
589 iter++;
590 }
591 return std::sqrt(sMs); // w.norm();
592}
593
594template<typename Real>
596 const Vector<Real> &v,
597 const Vector<Real> &x,
600 Real &tol,
601 Vector<Real> &pwa) const {
602 const Real zero(0);
603 pwa.set(v);
604 bnd.pruneActive(pwa,x,zero);
605 model.hessVec(hv,pwa,x,tol); nhess_++;
606 pwa.set(hv.dual());
607 bnd.pruneActive(pwa,x,zero);
608 hv.set(pwa.dual());
609 //pwa.set(v);
610 //bnd.pruneActive(pwa,x,zero);
611 //model.hessVec(hv,pwa,x,tol); nhess_++;
612 //bnd.pruneActive(hv,x,zero);
613}
614
615template<typename Real>
617 const Vector<Real> &v,
618 const Vector<Real> &x,
621 Real &tol,
622 Vector<Real> &dwa,
623 Vector<Real> &pwa) const {
624 if (!hasEcon_) {
625 const Real zero(0);
626 pwa.set(v.dual());
627 bnd.pruneActive(pwa,x,zero);
628 dwa.set(pwa.dual());
629 model.precond(hv,dwa,x,tol);
630 bnd.pruneActive(hv,x,zero);
631 //dwa.set(v);
632 //bnd.pruneActive(dwa,x,zero);
633 //model.precond(hv,dwa,x,tol);
634 //bnd.pruneActive(hv,x,zero);
635 }
636 else {
637 // Perform null space projection
638 rcon_->setX(makePtrFromRef(x));
639 ns_->update(x);
640 pwa.set(v.dual());
641 ns_->apply(hv,pwa,tol);
642 }
643}
644
645template<typename Real>
646void LinMoreAlgorithm<Real>::writeHeader( std::ostream& os ) const {
647 std::stringstream hist;
648 if (verbosity_ > 1) {
649 hist << std::string(114,'-') << std::endl;
650 hist << " Lin-More trust-region method status output definitions" << std::endl << std::endl;
651 hist << " iter - Number of iterates (steps taken)" << std::endl;
652 hist << " value - Objective function value" << std::endl;
653 hist << " gnorm - Norm of the gradient" << std::endl;
654 hist << " snorm - Norm of the step (update to optimization vector)" << std::endl;
655 hist << " delta - Trust-Region radius" << std::endl;
656 hist << " #fval - Number of times the objective function was evaluated" << std::endl;
657 hist << " #grad - Number of times the gradient was computed" << std::endl;
658 hist << " #hess - Number of times the Hessian was applied" << std::endl;
659 hist << " #proj - Number of times the projection was applied" << std::endl;
660 hist << std::endl;
661 hist << " tr_flag - Trust-Region flag" << std::endl;
662 for( int flag = TRUtils::SUCCESS; flag != TRUtils::UNDEFINED; ++flag ) {
663 hist << " " << NumberToString(flag) << " - "
664 << TRUtils::ETRFlagToString(static_cast<TRUtils::ETRFlag>(flag)) << std::endl;
665 }
666 hist << std::endl;
667 if (minit_ > 0) {
668 hist << " iterCG - Number of Truncated CG iterations" << std::endl << std::endl;
669 hist << " flagGC - Trust-Region Truncated CG flag" << std::endl;
670 for( int flag = CG_FLAG_SUCCESS; flag != CG_FLAG_UNDEFINED; ++flag ) {
671 hist << " " << NumberToString(flag) << " - "
672 << ECGFlagToString(static_cast<ECGFlag>(flag)) << std::endl;
673 }
674 }
675 hist << std::string(114,'-') << std::endl;
676 }
677 hist << " ";
678 hist << std::setw(6) << std::left << "iter";
679 hist << std::setw(15) << std::left << "value";
680 hist << std::setw(15) << std::left << "gnorm";
681 hist << std::setw(15) << std::left << "snorm";
682 hist << std::setw(15) << std::left << "delta";
683 hist << std::setw(10) << std::left << "#fval";
684 hist << std::setw(10) << std::left << "#grad";
685 hist << std::setw(10) << std::left << "#hess";
686 hist << std::setw(10) << std::left << "#proj";
687 hist << std::setw(10) << std::left << "tr_flag";
688 if (minit_ > 0) {
689 hist << std::setw(10) << std::left << "iterCG";
690 hist << std::setw(10) << std::left << "flagCG";
691 }
692 hist << std::endl;
693 os << hist.str();
694}
695
696template<typename Real>
697void LinMoreAlgorithm<Real>::writeName( std::ostream& os ) const {
698 std::stringstream hist;
699 hist << std::endl << "Lin-More Trust-Region Method (Type B, Bound Constraints)" << std::endl;
700 os << hist.str();
701}
702
703template<typename Real>
704void LinMoreAlgorithm<Real>::writeOutput( std::ostream& os, bool write_header ) const {
705 std::stringstream hist;
706 hist << std::scientific << std::setprecision(6);
707 if ( state_->iter == 0 ) writeName(os);
708 if ( write_header ) writeHeader(os);
709 if ( state_->iter == 0 ) {
710 hist << " ";
711 hist << std::setw(6) << std::left << state_->iter;
712 hist << std::setw(15) << std::left << state_->value;
713 hist << std::setw(15) << std::left << state_->gnorm;
714 hist << std::setw(15) << std::left << "---";
715 hist << std::setw(15) << std::left << state_->searchSize;
716 hist << std::setw(10) << std::left << state_->nfval;
717 hist << std::setw(10) << std::left << state_->ngrad;
718 hist << std::setw(10) << std::left << nhess_;
719 hist << std::setw(10) << std::left << state_->nproj;
720 hist << std::setw(10) << std::left << "---";
721 if (minit_ > 0) {
722 hist << std::setw(10) << std::left << "---";
723 hist << std::setw(10) << std::left << "---";
724 }
725 hist << std::endl;
726 }
727 else {
728 hist << " ";
729 hist << std::setw(6) << std::left << state_->iter;
730 hist << std::setw(15) << std::left << state_->value;
731 hist << std::setw(15) << std::left << state_->gnorm;
732 hist << std::setw(15) << std::left << state_->snorm;
733 hist << std::setw(15) << std::left << state_->searchSize;
734 hist << std::setw(10) << std::left << state_->nfval;
735 hist << std::setw(10) << std::left << state_->ngrad;
736 hist << std::setw(10) << std::left << nhess_;
737 hist << std::setw(10) << std::left << state_->nproj;
738 hist << std::setw(10) << std::left << TRflag_;
739 if (minit_ > 0) {
740 hist << std::setw(10) << std::left << SPiter_;
741 hist << std::setw(10) << std::left << SPflag_;
742 }
743 hist << std::endl;
744 }
745 os << hist.str();
746}
747
748} // namespace TypeB
749} // namespace ROL
750
751#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0 zero)()
Provides the interface to apply upper and lower bound constraints.
void pruneActive(Vector< Real > &v, const Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the -active set.
Provides the interface to evaluate objective functions.
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
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 interface for and implements limited-memory secant operators.
Provides an interface to check status of optimization algorithms.
Provides the interface to evaluate trust-region model functions.
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &s, Real &tol) override
Apply Hessian approximation to vector.
virtual void precond(Vector< Real > &Pv, const Vector< Real > &v, const Vector< Real > &s, Real &tol) override
Apply preconditioner to vector.
void initialize(const Vector< Real > &x, const Vector< Real > &g)
Real optimalityCriterion(const Vector< Real > &x, const Vector< Real > &g, Vector< Real > &primal, std::ostream &outStream=std::cout) const
virtual void writeExitStatus(std::ostream &os) const
void initialize(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &bnd, std::ostream &outStream=std::cout)
void applyFreeHessian(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, TrustRegionModel_U< Real > &model, BoundConstraint< Real > &bnd, Real &tol, Vector< Real > &pwa) const
Real dcauchy(Vector< Real > &s, Real &alpha, Real &q, const Vector< Real > &x, const Vector< Real > &g, const Real del, TrustRegionModel_U< Real > &model, Vector< Real > &dwa, Vector< Real > &dwa1, std::ostream &outStream=std::cout)
void writeOutput(std::ostream &os, bool write_header=false) const override
Print iterate status.
void run(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &bnd, std::ostream &outStream=std::cout) override
Run algorithm on bound constrained problems (Type-B). This general interface supports the use of dual...
Real dtrqsol(const Real xtx, const Real ptp, const Real ptx, const Real del) const
Real dtrpcg(Vector< Real > &w, int &iflag, int &iter, const Vector< Real > &g, const Vector< Real > &x, const Real del, TrustRegionModel_U< Real > &model, BoundConstraint< Real > &bnd, const Real tol, const Real stol, const int itermax, Vector< Real > &p, Vector< Real > &q, Vector< Real > &r, Vector< Real > &t, Vector< Real > &pwa, Vector< Real > &dwa, std::ostream &outStream=std::cout) const
void applyFreePrecond(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, TrustRegionModel_U< Real > &model, BoundConstraint< Real > &bnd, Real &tol, Vector< Real > &dwa, Vector< Real > &pwa) const
Real dprsrch(Vector< Real > &x, Vector< Real > &s, Real &q, const Vector< Real > &g, TrustRegionModel_U< Real > &model, BoundConstraint< Real > &bnd, Vector< Real > &pwa, Vector< Real > &dwa, std::ostream &outStream=std::cout)
void writeName(std::ostream &os) const override
Print step name.
Real dgpstep(Vector< Real > &s, const Vector< Real > &w, const Vector< Real > &x, const Real alpha, std::ostream &outStream=std::cout) const
LinMoreAlgorithm(ParameterList &list, const Ptr< Secant< Real > > &secant=nullPtr)
void writeHeader(std::ostream &os) const override
Print iterate header.
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 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 .
std::string ETRFlagToString(ETRFlag trf)
std::string NumberToString(T Number)
Definition ROL_Types.hpp:81
ESecant StringToESecant(std::string s)
ESecantMode
@ SECANTMODE_FORWARD
@ SECANTMODE_INVERSE
@ SECANTMODE_BOTH
@ CG_FLAG_UNDEFINED
@ CG_FLAG_SUCCESS
std::string ECGFlagToString(ECGFlag cgf)