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