ROL
ROL_TypeB_TrustRegionSPGAlgorithm_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_TRUSTREGIONSPGALGORITHM_DEF_HPP
45#define ROL_TYPEB_TRUSTREGIONSPGALGORITHM_DEF_HPP
46
47#include <deque>
48
49namespace ROL {
50namespace TypeB {
51
52template<typename Real>
54 const Ptr<Secant<Real>> &secant) {
55 // Set status test
56 status_->reset();
57 status_->add(makePtr<StatusTest<Real>>(list));
58
59 ParameterList &trlist = list.sublist("Step").sublist("Trust Region");
60 // Trust-Region Parameters
61 state_->searchSize = trlist.get("Initial Radius", -1.0);
62 delMax_ = trlist.get("Maximum Radius", ROL_INF<Real>());
63 eta0_ = trlist.get("Step Acceptance Threshold", 0.05);
64 eta1_ = trlist.get("Radius Shrinking Threshold", 0.05);
65 eta2_ = trlist.get("Radius Growing Threshold", 0.9);
66 gamma0_ = trlist.get("Radius Shrinking Rate (Negative rho)", 0.0625);
67 gamma1_ = trlist.get("Radius Shrinking Rate (Positive rho)", 0.25);
68 gamma2_ = trlist.get("Radius Growing Rate", 2.5);
69 TRsafe_ = trlist.get("Safeguard Size", 100.0);
70 eps_ = TRsafe_*ROL_EPSILON<Real>();
71 interpRad_ = trlist.get("Use Radius Interpolation", false);
72 verbosity_ = trlist.sublist("General").get("Output Level", 0);
73 // Algorithm-Specific Parameters
74 ROL::ParameterList &lmlist = trlist.sublist("SPG");
75 useNM_ = lmlist.get("Use Nonmonotone Trust Region", false);
76 maxNM_ = lmlist.get("Maximum Storage Size", 10);
77 mu0_ = lmlist.get("Sufficient Decrease Parameter", 1e-2);
78 spexp_ = lmlist.get("Relative Tolerance Exponent", 1.0);
79 spexp_ = std::max(static_cast<Real>(1),std::min(spexp_,static_cast<Real>(2)));
80 redlim_ = lmlist.sublist("Cauchy Point").get("Maximum Number of Reduction Steps", 10);
81 explim_ = lmlist.sublist("Cauchy Point").get("Maximum Number of Expansion Steps", 10);
82 alpha_ = lmlist.sublist("Cauchy Point").get("Initial Step Size", 1.0);
83 normAlpha_ = lmlist.sublist("Cauchy Point").get("Normalize Initial Step Size", false);
84 interpf_ = lmlist.sublist("Cauchy Point").get("Reduction Rate", 0.1);
85 extrapf_ = lmlist.sublist("Cauchy Point").get("Expansion Rate", 10.0);
86 qtol_ = lmlist.sublist("Cauchy Point").get("Decrease Tolerance", 1e-8);
87 // Spectral projected gradient parameters
88 lambdaMin_ = lmlist.sublist("Solver").get("Minimum Spectral Step Size", 1e-8);
89 lambdaMax_ = lmlist.sublist("Solver").get("Maximum Spectral Step Size", 1e8);
90 gamma_ = lmlist.sublist("Solver").get("Sufficient Decrease Tolerance", 1e-4);
91 maxSize_ = lmlist.sublist("Solver").get("Maximum Storage Size", 10);
92 maxit_ = lmlist.sublist("Solver").get("Iteration Limit", 25);
93 tol1_ = lmlist.sublist("Solver").get("Absolute Tolerance", 1e-4);
94 tol2_ = lmlist.sublist("Solver").get("Relative Tolerance", 1e-2);
95 useMin_ = lmlist.sublist("Solver").get("Use Smallest Model Iterate", true);
96 useNMSP_ = lmlist.sublist("Solver").get("Use Nonmonotone Search", false);
97 // Inexactness Information
98 ParameterList &glist = list.sublist("General");
99 useInexact_.clear();
100 useInexact_.push_back(glist.get("Inexact Objective Function", false));
101 useInexact_.push_back(glist.get("Inexact Gradient", false));
102 useInexact_.push_back(glist.get("Inexact Hessian-Times-A-Vector", false));
103 // Trust-Region Inexactness Parameters
104 ParameterList &ilist = trlist.sublist("Inexact").sublist("Gradient");
105 scale0_ = ilist.get("Tolerance Scaling", static_cast<Real>(0.1));
106 scale1_ = ilist.get("Relative Tolerance", static_cast<Real>(2));
107 // Inexact Function Evaluation Information
108 ParameterList &vlist = trlist.sublist("Inexact").sublist("Value");
109 scale_ = vlist.get("Tolerance Scaling", static_cast<Real>(1.e-1));
110 omega_ = vlist.get("Exponent", static_cast<Real>(0.9));
111 force_ = vlist.get("Forcing Sequence Initial Value", static_cast<Real>(1.0));
112 updateIter_ = vlist.get("Forcing Sequence Update Frequency", static_cast<int>(10));
113 forceFactor_ = vlist.get("Forcing Sequence Reduction Factor", static_cast<Real>(0.1));
114 // Output Parameters
115 verbosity_ = list.sublist("General").get("Output Level",0);
116 writeHeader_ = verbosity_ > 2;
117 // Secant Information
118 useSecantPrecond_ = list.sublist("General").sublist("Secant").get("Use as Preconditioner", false);
119 useSecantHessVec_ = list.sublist("General").sublist("Secant").get("Use as Hessian", false);
121 if (useSecantPrecond_ && !useSecantHessVec_) mode = SECANTMODE_INVERSE;
122 else if (useSecantHessVec_ && !useSecantPrecond_) mode = SECANTMODE_FORWARD;
123 // Initialize trust region model
124 model_ = makePtr<TrustRegionModel_U<Real>>(list,secant,mode);
125 if (secant == nullPtr) {
126 esec_ = StringToESecant(list.sublist("General").sublist("Secant").get("Type","Limited-Memory BFGS"));
127 }
128}
129
130template<typename Real>
132 const Vector<Real> &g,
133 Real ftol,
134 Objective<Real> &obj,
136 std::ostream &outStream) {
137 //const Real one(1);
138 if (proj_ == nullPtr)
139 proj_ = makePtr<PolyhedralProjection<Real>>(makePtrFromRef(bnd));
140 // Initialize data
142 nhess_ = 0;
143 // Update approximate gradient and approximate objective function.
144 proj_->project(x,outStream); state_->nproj++;
145 state_->iterateVec->set(x);
146 obj.update(x,UpdateType::Initial,state_->iter);
147 state_->value = obj.value(x,ftol);
148 state_->nfval++;
149 //obj.gradient(*state_->gradientVec,x,ftol);
150 state_->gnorm = computeGradient(x,*state_->gradientVec,*state_->stepVec,state_->searchSize,obj,outStream);
151 state_->ngrad++;
152 //state_->stepVec->set(x);
153 //state_->stepVec->axpy(-one,state_->gradientVec->dual());
154 //proj_->project(*state_->stepVec,outStream); state_->nproj++;
155 //state_->stepVec->axpy(-one,x);
156 //state_->gnorm = state_->stepVec->norm();
157 state_->snorm = ROL_INF<Real>();
158 // Normalize initial CP step length
159 if (normAlpha_) alpha_ /= state_->gradientVec->norm();
160 // Compute initial trust region radius if desired.
161 if ( state_->searchSize <= static_cast<Real>(0) )
162 state_->searchSize = state_->gradientVec->norm();
163}
164
165template<typename Real>
167 Real &outTol,
168 Real pRed,
169 Real &fold,
170 int iter,
171 const Vector<Real> &x,
172 const Vector<Real> &xold,
173 Objective<Real> &obj) {
174 outTol = std::sqrt(ROL_EPSILON<Real>());
175 if ( useInexact_[0] ) {
176 if (!(iter%updateIter_) && (iter!=0)) force_ *= forceFactor_;
177 const Real one(1);
178 Real eta = static_cast<Real>(0.999)*std::min(eta1_,one-eta2_);
179 outTol = scale_*std::pow(eta*std::min(pRed,force_),one/omega_);
180 if (inTol > outTol) fold = obj.value(xold,outTol);
181 }
182 // Evaluate objective function at new iterate
184 Real fval = obj.value(x,outTol);
185 return fval;
186}
187
188template<typename Real>
190 Vector<Real> &g,
191 Vector<Real> &pwa,
192 Real del,
193 Objective<Real> &obj,
194 std::ostream &outStream) const {
195 Real gnorm(0);
196 if ( useInexact_[1] ) {
197 const Real one(1);
198 Real gtol1 = scale0_*del;
199 Real gtol0 = gtol1 + one;
200 while ( gtol0 > gtol1 ) {
201 obj.gradient(g,x,gtol1);
202 gnorm = TypeB::Algorithm<Real>::optimalityCriterion(x,g,pwa,outStream);
203 gtol0 = gtol1;
204 gtol1 = scale0_*std::min(gnorm,del);
205 }
206 }
207 else {
208 Real gtol = std::sqrt(ROL_EPSILON<Real>());
209 obj.gradient(g,x,gtol);
210 gnorm = TypeB::Algorithm<Real>::optimalityCriterion(x,g,pwa,outStream);
211 }
212 return gnorm;
213}
214
215template<typename Real>
217 const Vector<Real> &g,
218 Objective<Real> &obj,
220 std::ostream &outStream ) {
221 const Real zero(0), one(1);
222 //Real tol0 = std::sqrt(ROL_EPSILON<Real>());
223 Real inTol = static_cast<Real>(0.1)*ROL_OVERFLOW<Real>(), outTol(inTol);
224 Real ftrial(0), fcheck(0), pRed(0), rho(1), q(0);
225 // Initialize trust-region data
226 std::vector<std::string> output;
227 initialize(x,g,inTol,obj,bnd,outStream);
228 Ptr<Vector<Real>> gmod = g.clone();
229 Ptr<Vector<Real>> pwa1 = x.clone(), pwa2 = x.clone();
230 Ptr<Vector<Real>> pwa3 = x.clone(), pwa4 = x.clone();
231 Ptr<Vector<Real>> pwa5 = x.clone(), pwa6 = x.clone();
232 Ptr<Vector<Real>> pwa7 = x.clone();
233 Ptr<Vector<Real>> dwa1 = g.clone(), dwa2 = g.clone();
234
235 // Output
236 if (verbosity_ > 0) writeOutput(outStream,true);
237
238 std::deque<Real> fqueue;
239 if (useNM_) fqueue.push_back(state_->value);
240 while (status_->check(*state_)) {
241 // Build trust-region model
242 model_->setData(obj,*state_->iterateVec,*state_->gradientVec);
243
244 /**** SOLVE TRUST-REGION SUBPROBLEM ****/
245 // Compute Cauchy point (TRON notation: x = x[1])
246 dcauchy(*state_->stepVec,alpha_,q,*state_->iterateVec,
247 state_->gradientVec->dual(),state_->searchSize,
248 *model_,*dwa1,*dwa2,outStream); // Solve 1D optimization problem for alpha
249 x.plus(*state_->stepVec); // Set x = x[0] + alpha*g
250
251 // Model gradient at s = x[1] - x[0]
252 gmod->set(*dwa1); // hessVec from Cauchy point computation
253 gmod->plus(*state_->gradientVec);
254
255 // Apply SPG starting from the Cauchy point
256 dpsg(x,q,*gmod,*state_->iterateVec,state_->searchSize,*model_,
257 *pwa1,*pwa2,*pwa3,*pwa4,*pwa5,*pwa6,*pwa7,*dwa1,outStream);
258 pRed = -q;
259 state_->stepVec->set(x); state_->stepVec->axpy(-one,*state_->iterateVec);
260 state_->snorm = state_->stepVec->norm();
261
262 // Compute trial objective value
263 ftrial = computeValue(inTol,outTol,pRed,state_->value,state_->iter,x,*state_->iterateVec,obj);
264 //obj.update(x,UpdateType::Trial);
265 //ftrial = obj.value(x,tol0);
266 state_->nfval++;
267
268 // Compute ratio of acutal and predicted reduction
269 fcheck = useNM_ ? *std::max_element(fqueue.begin(),fqueue.end()) : state_->value;
270 TRflag_ = TRUtils::SUCCESS;
271 TRUtils::analyzeRatio<Real>(rho,TRflag_,fcheck,ftrial,pRed,eps_,outStream,verbosity_>1);
272 //TRUtils::analyzeRatio<Real>(rho,TRflag_,state_->value,ftrial,pRed,eps_,outStream,verbosity_>1);
273
274 // Update algorithm state
275 state_->iter++;
276 // Accept/reject step and update trust region radius
277 if ((rho < eta0_ && TRflag_ == TRUtils::SUCCESS) || (TRflag_ >= 2)) { // Step Rejected
278 x.set(*state_->iterateVec);
279 obj.update(x,UpdateType::Revert,state_->iter);
280 if (interpRad_ && (rho < zero && TRflag_ != TRUtils::TRNAN)) {
281 // Negative reduction, interpolate to find new trust-region radius
282 state_->searchSize = TRUtils::interpolateRadius<Real>(*state_->gradientVec,*state_->stepVec,
283 state_->snorm,pRed,state_->value,ftrial,state_->searchSize,gamma0_,gamma1_,eta2_,
284 outStream,verbosity_>1);
285 }
286 else { // Shrink trust-region radius
287 state_->searchSize = gamma1_*std::min(state_->snorm,state_->searchSize);
288 }
289 }
290 else if ((rho >= eta0_ && TRflag_ != TRUtils::NPOSPREDNEG)
291 || (TRflag_ == TRUtils::POSPREDNEG)) { // Step Accepted
292 state_->value = ftrial;
293 obj.update(x,UpdateType::Accept,state_->iter);
294 inTol = outTol;
295 // Increase trust-region radius
296 if (rho >= eta2_) state_->searchSize = std::min(gamma2_*state_->searchSize, delMax_);
297 // Compute gradient at new iterate
298 dwa1->set(*state_->gradientVec);
299 //obj.gradient(*state_->gradientVec,x,tol0);
300 //state_->gnorm = TypeB::Algorithm<Real>::optimalityCriterion(x,*state_->gradientVec,*pwa1,outStream);
301 state_->gnorm = computeGradient(x,*state_->gradientVec,*pwa1,state_->searchSize,obj,outStream);
302 state_->ngrad++;
303 state_->iterateVec->set(x);
304 // Update secant information in trust-region model
305 model_->update(x,*state_->stepVec,*dwa1,*state_->gradientVec,
306 state_->snorm,state_->iter);
307 if (useNM_) {
308 if (static_cast<int>(fqueue.size())==maxNM_) fqueue.pop_front();
309 fqueue.push_back(state_->value);
310 }
311 }
312
313 // Update Output
314 if (verbosity_ > 0) writeOutput(outStream,writeHeader_);
315 }
316 if (verbosity_ > 0) TypeB::Algorithm<Real>::writeExitStatus(outStream);
317}
318
319template<typename Real>
321 const Vector<Real> &x, const Real alpha,
322 std::ostream &outStream) const {
323 s.set(x); s.axpy(alpha,w);
324 proj_->project(s,outStream); state_->nproj++;
325 s.axpy(static_cast<Real>(-1),x);
326 return s.norm();
327}
328
329template<typename Real>
331 Real &alpha,
332 Real &q,
333 const Vector<Real> &x,
334 const Vector<Real> &g,
335 const Real del,
337 Vector<Real> &dwa, Vector<Real> &dwa1,
338 std::ostream &outStream) {
339 const Real half(0.5);
340 // const Real zero(0); // Unused
341 Real tol = std::sqrt(ROL_EPSILON<Real>());
342 bool interp = false;
343 Real gs(0), snorm(0);
344 // Compute s = P(x[0] - alpha g[0])
345 snorm = dgpstep(s,g,x,-alpha,outStream);
346 if (snorm > del) {
347 interp = true;
348 }
349 else {
350 model.hessVec(dwa,s,x,tol); nhess_++;
351 gs = s.dot(g);
352 //q = half * s.dot(dwa.dual()) + gs;
353 q = half * s.apply(dwa) + gs;
354 interp = (q > mu0_*gs);
355 }
356 // Either increase or decrease alpha to find approximate Cauchy point
357 int cnt = 0;
358 if (interp) {
359 bool search = true;
360 while (search) {
361 alpha *= interpf_;
362 snorm = dgpstep(s,g,x,-alpha,outStream);
363 if (snorm <= del) {
364 model.hessVec(dwa,s,x,tol); nhess_++;
365 gs = s.dot(g);
366 //q = half * s.dot(dwa.dual()) + gs;
367 q = half * s.apply(dwa) + gs;
368 search = (q > mu0_*gs) && (cnt < redlim_);
369 }
370 cnt++;
371 }
372 }
373 else {
374 bool search = true;
375 Real alphas = alpha;
376 Real qs = q;
377 dwa1.set(dwa);
378 while (search) {
379 alpha *= extrapf_;
380 snorm = dgpstep(s,g,x,-alpha,outStream);
381 if (snorm <= del && cnt < explim_) {
382 model.hessVec(dwa,s,x,tol); nhess_++;
383 gs = s.dot(g);
384 //q = half * s.dot(dwa.dual()) + gs;
385 q = half * s.apply(dwa) + gs;
386 if (q <= mu0_*gs && std::abs(q-qs) > qtol_*std::abs(qs)) {
387 dwa1.set(dwa);
388 search = true;
389 alphas = alpha;
390 qs = q;
391 }
392 else {
393 q = qs;
394 dwa.set(dwa1);
395 search = false;
396 }
397 }
398 else {
399 search = false;
400 }
401 cnt++;
402 }
403 alpha = alphas;
404 snorm = dgpstep(s,g,x,-alpha,outStream);
405 }
406 if (verbosity_ > 1) {
407 outStream << " Cauchy point" << std::endl;
408 outStream << " Step length (alpha): " << alpha << std::endl;
409 outStream << " Step length (alpha*g): " << snorm << std::endl;
410 outStream << " Model decrease (pRed): " << -q << std::endl;
411 if (!interp) {
412 outStream << " Number of extrapolation steps: " << cnt << std::endl;
413 }
414 }
415 return snorm;
416}
417
418template<typename Real>
420 Real &q,
421 Vector<Real> &gmod,
422 const Vector<Real> &x,
423 Real del,
425 Vector<Real> &ymin,
426 Vector<Real> &pwa,
427 Vector<Real> &pwa1,
428 Vector<Real> &pwa2,
429 Vector<Real> &pwa3,
430 Vector<Real> &pwa4,
431 Vector<Real> &pwa5,
432 Vector<Real> &dwa,
433 std::ostream &outStream) {
434 // Use SPG to approximately solve TR subproblem:
435 // min 1/2 <H(y-x), (y-x)> + <g, (y-x)> subject to y\in C, ||y|| \le del
436 //
437 // Inpute:
438 // y = Cauchy step
439 // x = Current iterate
440 // g = Current gradient
441 const Real zero(0), half(0.5), one(1), two(2), eps(std::sqrt(ROL_EPSILON<Real>()));
442 Real tol(std::sqrt(ROL_EPSILON<Real>()));
443 Real alpha(1), sHs(0), alphaTmp(1), mmax(0), qmin(0);
444 std::deque<Real> mqueue; mqueue.push_back(q);
445
446 if (useNMSP_ && useMin_) { qmin = q; ymin.set(y); }
447
448 // Compute initial projected gradient norm
449 pwa1.set(gmod.dual());
450 pwa.set(y); pwa.axpy(-one,pwa1);
451 dproj(pwa,x,del,pwa2,pwa3,pwa4,pwa5,outStream);
452 pwa.axpy(-one,y);
453 Real gnorm = pwa.norm();
454 const Real gtol = std::min(tol1_,tol2_*gnorm);
455
456 // Compute initial step
457 Real lambda = std::max(lambdaMin_,std::min(one/gmod.norm(),lambdaMax_));
458 pwa.set(y); pwa.axpy(-lambda,pwa1); // pwa = y - lambda gmod.dual()
459 dproj(pwa,x,del,pwa2,pwa3,pwa4,pwa5,outStream); // pwa = P(y - lambda gmod.dual())
460 pwa.axpy(-one,y); // pwa = P(y - lambda gmod.dual()) - y = step
461 Real gs = gmod.apply(pwa); // gs = <step, model gradient>
462 Real ss = pwa.dot(pwa); // Norm squared of step
463
464 if (verbosity_ > 1)
465 outStream << " Spectral Projected Gradient" << std::endl;
466
467 SPiter_ = 0;
468 while (SPiter_ < maxit_) {
469 SPiter_++;
470
471 // Evaluate model Hessian
472 model.hessVec(dwa,pwa,x,tol); nhess_++; // dwa = H step
473 sHs = dwa.apply(pwa); // sHs = <step, H step>
474
475 // Perform line search
476 if (useNMSP_) { // Nonmonotone
477 mmax = *std::max_element(mqueue.begin(),mqueue.end());
478 alphaTmp = (-(one-gamma_)*gs + std::sqrt(std::pow((one-gamma_)*gs,two)-two*sHs*(q-mmax)))/sHs;
479 }
480 else { // Exact
481 alphaTmp = -gs/sHs;
482 }
483 alpha = (sHs > zero ? std::min(one,std::max(zero,alphaTmp)) : one);
484
485 // Update model quantities
486 q += alpha * (gs + half * alpha * sHs); // Update model value
487 gmod.axpy(alpha,dwa); // Update model gradient
488 y.axpy(alpha,pwa); // New iterate
489
490 // Update nonmonotone line search information
491 if (useNMSP_) {
492 if (static_cast<int>(mqueue.size())==maxSize_) mqueue.pop_front();
493 mqueue.push_back(q);
494 if (useMin_ && q <= qmin) { qmin = q; ymin.set(y); }
495 }
496
497 // Compute projected gradient norm
498 pwa1.set(gmod.dual());
499 pwa.set(y); pwa.axpy(-one,pwa1);
500 dproj(pwa,x,del,pwa2,pwa3,pwa4,pwa5,outStream);
501 pwa.axpy(-one,y);
502 gnorm = pwa.norm();
503
504 if (verbosity_ > 1) {
505 outStream << std::endl;
506 outStream << " Iterate: " << SPiter_ << std::endl;
507 outStream << " Spectral step length (lambda): " << lambda << std::endl;
508 outStream << " Step length (alpha): " << alpha << std::endl;
509 outStream << " Model decrease (pRed): " << -q << std::endl;
510 outStream << " Optimality criterion: " << gnorm << std::endl;
511 outStream << std::endl;
512 }
513 if (gnorm < gtol) break;
514
515 // Compute new spectral step
516 lambda = (sHs<=eps ? lambdaMax_ : std::max(lambdaMin_,std::min(ss/sHs,lambdaMax_)));
517 pwa.set(y); pwa.axpy(-lambda,pwa1);
518 dproj(pwa,x,del,pwa2,pwa3,pwa4,pwa5,outStream);
519 pwa.axpy(-one,y);
520 gs = gmod.apply(pwa);
521 ss = pwa.dot(pwa);
522 }
523 if (useNMSP_ && useMin_) { q = qmin; y.set(ymin); }
524 SPflag_ = (SPiter_==maxit_) ? 1 : 0;
525}
526
527template<typename Real>
529 const Vector<Real> &x0,
530 Real del,
531 Vector<Real> &y0,
532 Vector<Real> &y1,
533 Vector<Real> &yc,
534 Vector<Real> &pwa,
535 std::ostream &outStream) const {
536 // Solve ||P(t*x0 + (1-t)*(x-x0))-x0|| = del using Brent's method
537 const Real zero(0), half(0.5), one(1), two(2), three(3);
538 const Real eps(ROL_EPSILON<Real>()), tol0(1e1*eps), fudge(1.0-1e-2*sqrt(eps));
539 Real f0(0), f1(0), fc(0), t0(0), t1(1), tc(0), d1(1), d2(1), tol(1);
540 Real p(0), q(0), r(0), s(0), m(0);
541 int cnt(state_->nproj);
542 y1.set(x);
543 proj_->project(y1,outStream); state_->nproj++;
544 pwa.set(y1); pwa.axpy(-one,x0);
545 f1 = pwa.norm();
546 if (f1 <= del) {
547 x.set(y1);
548 return;
549 }
550 y0.set(x0);
551 tc = t0; fc = f0; yc.set(y0);
552 d1 = t1-t0; d2 = d1;
553 int code = 0;
554 while (true) {
555 if (std::abs(fc-del) < std::abs(f1-del)) {
556 t0 = t1; t1 = tc; tc = t0;
557 f0 = f1; f1 = fc; fc = f0;
558 y0.set(y1); y1.set(yc); yc.set(y0);
559 }
560 tol = two*eps*std::abs(t1) + half*tol0;
561 m = half*(tc - t1);
562 if (std::abs(m) <= tol) { code = 1; break; }
563 if ((f1 >= fudge*del && f1 <= del)) break;
564 if (std::abs(d1) < tol || std::abs(f0-del) <= std::abs(f1-del)) {
565 d1 = m; d2 = d1;
566 }
567 else {
568 s = (f1-del)/(f0-del);
569 if (t0 == tc) {
570 p = two*m*s;
571 q = one-s;
572 }
573 else {
574 q = (f0-del)/(fc-del);
575 r = (f1-del)/(fc-del);
576 p = s*(two*m*q*(q-r)-(t1-t0)*(r-one));
577 q = (q-one)*(r-one)*(s-one);
578 }
579 if (p > zero) q = -q;
580 else p = -p;
581 s = d1;
582 d1 = d2;
583 if (two*p < three*m*q-std::abs(tol*q) && p < std::abs(half*s*q)) {
584 d2 = p/q;
585 }
586 else {
587 d1 = m; d2 = d1;
588 }
589 }
590 t0 = t1; f0 = f1; y0.set(y1);
591 if (std::abs(d2) > tol) t1 += d2;
592 else if (m > zero) t1 += tol;
593 else t1 -= tol;
594 y1.set(x); y1.scale(t1); y1.axpy(one-t1,x0);
595 proj_->project(y1,outStream); state_->nproj++;
596 pwa.set(y1); pwa.axpy(-one,x0);
597 f1 = pwa.norm();
598 if ((f1 > del && fc > del) || (f1 <= del && fc <= del)) {
599 tc = t0; fc = f0; yc.set(y0);
600 d1 = t1-t0; d2 = d1;
601 }
602 }
603 if (code==1 && f1>del) x.set(yc);
604 else x.set(y1);
605 if (verbosity_ > 1) {
606 outStream << std::endl;
607 outStream << " Trust-Region Subproblem Projection" << std::endl;
608 outStream << " Number of polyhedral projections: " << state_->nproj-cnt << std::endl;
609 if (code == 1 && f1 > del) {
610 outStream << " Transformed Multiplier: " << tc << std::endl;
611 outStream << " Dual Residual: " << fc-del << std::endl;
612 }
613 else {
614 outStream << " Transformed Multiplier: " << t1 << std::endl;
615 outStream << " Dual Residual: " << f1-del << std::endl;
616 }
617 outStream << " Exit Code: " << code << std::endl;
618 outStream << std::endl;
619 }
620}
621
622// BRACKETING AND BRENTS FOR UNTRANSFORMED MULTIPLIER
623//template<typename Real>
624//void TrustRegionSPGAlgorithm<Real>::dproj(Vector<Real> &x,
625// const Vector<Real> &x0,
626// Real del,
627// Vector<Real> &y0,
628// Vector<Real> &y1,
629// Vector<Real> &yc,
630// Vector<Real> &pwa,
631// std::ostream &outStream) const {
632// // Solve ||P(t*x0 + (1-t)*(x-x0))-x0|| = del using Brent's method
633// const Real zero(0), half(0.5), one(1), two(2), three(3);
634// const Real eps(ROL_EPSILON<Real>()), tol0(1e1*eps), fudge(1.0-1e-2*sqrt(eps));
635// Real f0(0), f1(0), fc(0), u0(0), u1(0), uc(0), t0(1), t1(0), tc(0), d1(1), d2(1), tol(1);
636// Real p(0), q(0), r(0), s(0), m(0);
637// int cnt(state_->nproj);
638// y0.set(x);
639// proj_->project(y0,outStream); state_->nproj++;
640// pwa.set(y0); pwa.axpy(-one,x0);
641// f0 = pwa.norm();
642// if (f0 <= del) {
643// x.set(y0);
644// return;
645// }
646//
647// // Bracketing
648// t1 = static_cast<Real>(1e-1);
649// f1 = one+del;
650// while (f1 >= del) {
651// t1 *= static_cast<Real>(5e-2);
652// y1.set(x); y1.scale(t1); y1.axpy(one-t1,x0);
653// proj_->project(y1,outStream); state_->nproj++;
654// pwa.set(y1); pwa.axpy(-one,x0);
655// f1 = pwa.norm();
656// }
657// u1 = (one-t1)/t1;
658//
659// // Brents
660// uc = u0; tc = t0; fc = f0; yc.set(y0);
661// d1 = u1-u0; d2 = d1;
662// int code = 0;
663// while (true) {
664// if (std::abs(fc-del) < std::abs(f1-del)) {
665// u0 = u1; u1 = uc; uc = u0;
666// t0 = t1; t1 = tc; tc = t0;
667// f0 = f1; f1 = fc; fc = f0;
668// y0.set(y1); y1.set(yc); yc.set(y0);
669// }
670// tol = two*eps*abs(u1) + half*tol0;
671// m = half*(uc - u1);
672// if (std::abs(m) <= tol) { code = 1; break; }
673// if ((f1 >= fudge*del && f1 <= del)) break;
674// if (std::abs(d1) < tol || std::abs(f0-del) <= std::abs(f1-del)) {
675// d1 = m; d2 = d1;
676// }
677// else {
678// s = (f1-del)/(f0-del);
679// if (u0 == uc) {
680// p = two*m*s;
681// q = one-s;
682// }
683// else {
684// q = (f0-del)/(fc-del);
685// r = (f1-del)/(fc-del);
686// p = s*(two*m*q*(q-r)-(u1-u0)*(r-one));
687// q = (q-one)*(r-one)*(s-one);
688// }
689// if (p > zero) q = -q;
690// else p = -p;
691// s = d1;
692// d1 = d2;
693// if (two*p < three*m*q-std::abs(tol*q) && p < std::abs(half*s*q)) {
694// d2 = p/q;
695// }
696// else {
697// d1 = m; d2 = d1;
698// }
699// }
700// u0 = u1; t0 = t1; f0 = f1; y0.set(y1);
701// if (std::abs(d2) > tol) u1 += d2;
702// else if (m > zero) u1 += tol;
703// else u1 -= tol;
704// t1 = one/(one+u1);
705// y1.set(x); y1.scale(t1); y1.axpy(one-t1,x0);
706// proj_->project(y1,outStream); state_->nproj++;
707// pwa.set(y1); pwa.axpy(-one,x0);
708// f1 = pwa.norm();
709// if ((f1 > del && fc > del) || (f1 <= del && fc <= del)) {
710// uc = u0; tc = t0; fc = f0; yc.set(y0);
711// d1 = u1-u0; d2 = d1;
712// }
713// }
714// if (code==1 && f1>del) x.set(yc);
715// else x.set(y1);
716// if (verbosity_ > 1) {
717// outStream << std::endl;
718// outStream << " Trust-Region Subproblem Projection" << std::endl;
719// outStream << " Number of polyhedral projections: " << state_->nproj-cnt << std::endl;
720// if (code == 1 && f1 > del) {
721// outStream << " Multiplier: " << uc << std::endl;
722// outStream << " Transformed Multiplier: " << tc << std::endl;
723// outStream << " Dual Residual: " << fc-del << std::endl;
724// }
725// else {
726// outStream << " Multiplier: " << u1 << std::endl;
727// outStream << " Transformed Multiplier: " << t1 << std::endl;
728// outStream << " Dual Residual: " << f1-del << std::endl;
729// }
730// outStream << " Exit Code: " << code << std::endl;
731// outStream << std::endl;
732// }
733//}
734
735// RIDDERS' METHOD FOR TRUST-REGION PROJECTION
736//template<typename Real>
737//void TrustRegionSPGAlgorithm<Real>::dproj(Vector<Real> &x,
738// const Vector<Real> &x0,
739// Real del,
740// Vector<Real> &y,
741// Vector<Real> &y1,
742// Vector<Real> &yc,
743// Vector<Real> &p,
744// std::ostream &outStream) const {
745// // Solve ||P(t*x0 + (1-t)*(x-x0))-x0|| = del using Ridder's method
746// const Real half(0.5), one(1), tol(1e1*ROL_EPSILON<Real>());
747// const Real fudge(1.0-1e-2*std::sqrt(ROL_EPSILON<Real>()));
748// Real e0(0), e1(0), e2(0), e(0), a0(0), a1(0.5), a2(1), a(0);
749// int cnt(state_->nproj);
750// y.set(x);
751// proj_->project(y,outStream); state_->nproj++;
752// p.set(y); p.axpy(-one,x0);
753// e2 = p.norm();
754// if (e2 <= del) {
755// x.set(y);
756// return;
757// }
758// bool code = 1;
759// while (a2-a0 > tol) {
760// a1 = half*(a0+a2);
761// y.set(x); y.scale(a1); y.axpy(one-a1,x0);
762// proj_->project(y,outStream); state_->nproj++;
763// p.set(y); p.axpy(-one,x0);
764// e1 = p.norm();
765// if (e1 >= fudge*del && e1 <= del) break;
766// a = a1-(a1-a0)*(e1-del)/std::sqrt((e1-del)*(e1-del)-(e0-del)*(e2-del));
767// y.set(x); y.scale(a); y.axpy(one-a,x0);
768// proj_->project(y,outStream); state_->nproj++;
769// p.set(y); p.axpy(-one,x0);
770// e = p.norm();
771// if (e < fudge*del) {
772// if (e1 < fudge*del) { e0 = (a < a1 ? e1 : e); a0 = (a < a1 ? a1 : a); }
773// else { e0 = e; a0 = a; e2 = e1; a2 = a1; };
774// }
775// else if (e > del) {
776// if (e1 < fudge*del) { e0 = e1; a0 = a1; e2 = e; a2 = a; }
777// else { e2 = (a < a1 ? e : e1); a2 = (a < a1 ? a : a1); }
778// }
779// else {
780// code = 0;
781// break; // Exit if fudge*del <= snorm <= del
782// }
783// }
784// x.set(y);
785// if (verbosity_ > 1) {
786// outStream << std::endl;
787// outStream << " Trust-Region Subproblem Projection" << std::endl;
788// outStream << " Number of polyhedral projections: " << state_->nproj-cnt << std::endl;
789// outStream << " Transformed Multiplier: " << a1 << std::endl;
790// outStream << " Dual Residual: " << e1-del << std::endl;
791// outStream << " Exit Code: " << code << std::endl;
792// outStream << std::endl;
793// }
794//}
795
796template<typename Real>
797void TrustRegionSPGAlgorithm<Real>::writeHeader( std::ostream& os ) const {
798 std::stringstream hist;
799 if (verbosity_ > 1) {
800 hist << std::string(114,'-') << std::endl;
801 hist << " SPG trust-region method status output definitions" << std::endl << std::endl;
802 hist << " iter - Number of iterates (steps taken)" << std::endl;
803 hist << " value - Objective function value" << std::endl;
804 hist << " gnorm - Norm of the gradient" << std::endl;
805 hist << " snorm - Norm of the step (update to optimization vector)" << std::endl;
806 hist << " delta - Trust-Region radius" << std::endl;
807 hist << " #fval - Number of times the objective function was evaluated" << std::endl;
808 hist << " #grad - Number of times the gradient was computed" << std::endl;
809 hist << " #hess - Number of times the Hessian was applied" << std::endl;
810 hist << " #proj - Number of times the projection was computed" << std::endl;
811 hist << std::endl;
812 hist << " tr_flag - Trust-Region flag" << std::endl;
813 for( int flag = TRUtils::SUCCESS; flag != TRUtils::UNDEFINED; ++flag ) {
814 hist << " " << NumberToString(flag) << " - "
815 << TRUtils::ETRFlagToString(static_cast<TRUtils::ETRFlag>(flag)) << std::endl;
816 }
817 hist << std::endl;
818 hist << " iterSPG - Number of Spectral Projected Gradient iterations" << std::endl << std::endl;
819 hist << " flagSPG - Trust-Region Truncated CG flag" << std::endl;
820 hist << " 0 - Converged" << std::endl;
821 hist << " 1 - Iteration Limit Exceeded" << std::endl;
822 hist << std::string(114,'-') << std::endl;
823 }
824 hist << " ";
825 hist << std::setw(6) << std::left << "iter";
826 hist << std::setw(15) << std::left << "value";
827 hist << std::setw(15) << std::left << "gnorm";
828 hist << std::setw(15) << std::left << "snorm";
829 hist << std::setw(15) << std::left << "delta";
830 hist << std::setw(10) << std::left << "#fval";
831 hist << std::setw(10) << std::left << "#grad";
832 hist << std::setw(10) << std::left << "#hess";
833 hist << std::setw(10) << std::left << "#proj";
834 hist << std::setw(10) << std::left << "tr_flag";
835 hist << std::setw(10) << std::left << "iterSPG";
836 hist << std::setw(10) << std::left << "flagSPG";
837 hist << std::endl;
838 os << hist.str();
839}
840
841template<typename Real>
842void TrustRegionSPGAlgorithm<Real>::writeName( std::ostream& os ) const {
843 std::stringstream hist;
844 hist << std::endl << "SPG Trust-Region Method (Type B, Bound Constraints)" << std::endl;
845 os << hist.str();
846}
847
848template<typename Real>
849void TrustRegionSPGAlgorithm<Real>::writeOutput( std::ostream& os, bool write_header ) const {
850 std::stringstream hist;
851 hist << std::scientific << std::setprecision(6);
852 if ( state_->iter == 0 ) writeName(os);
853 if ( write_header ) writeHeader(os);
854 if ( state_->iter == 0 ) {
855 hist << " ";
856 hist << std::setw(6) << std::left << state_->iter;
857 hist << std::setw(15) << std::left << state_->value;
858 hist << std::setw(15) << std::left << state_->gnorm;
859 hist << std::setw(15) << std::left << "---";
860 hist << std::setw(15) << std::left << state_->searchSize;
861 hist << std::setw(10) << std::left << state_->nfval;
862 hist << std::setw(10) << std::left << state_->ngrad;
863 hist << std::setw(10) << std::left << nhess_;
864 hist << std::setw(10) << std::left << state_->nproj;
865 hist << std::setw(10) << std::left << "---";
866 hist << std::setw(10) << std::left << "---";
867 hist << std::setw(10) << std::left << "---";
868 hist << std::endl;
869 }
870 else {
871 hist << " ";
872 hist << std::setw(6) << std::left << state_->iter;
873 hist << std::setw(15) << std::left << state_->value;
874 hist << std::setw(15) << std::left << state_->gnorm;
875 hist << std::setw(15) << std::left << state_->snorm;
876 hist << std::setw(15) << std::left << state_->searchSize;
877 hist << std::setw(10) << std::left << state_->nfval;
878 hist << std::setw(10) << std::left << state_->ngrad;
879 hist << std::setw(10) << std::left << nhess_;
880 hist << std::setw(10) << std::left << state_->nproj;
881 hist << std::setw(10) << std::left << TRflag_;
882 hist << std::setw(10) << std::left << SPiter_;
883 hist << std::setw(10) << std::left << SPflag_;
884 hist << std::endl;
885 }
886 os << hist.str();
887}
888
889} // namespace TypeB
890} // namespace ROL
891
892#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Provides the interface to apply upper and lower bound constraints.
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.
Definition: ROL_Secant.hpp:79
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.
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
Real computeValue(Real inTol, Real &outTol, Real pRed, Real &fold, int iter, const Vector< Real > &x, const Vector< Real > &xold, Objective< Real > &obj)
void initialize(Vector< Real > &x, const Vector< Real > &g, Real ftol, Objective< Real > &obj, BoundConstraint< Real > &bnd, std::ostream &outStream=std::cout)
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 dpsg(Vector< Real > &y, Real &q, Vector< Real > &gmod, const Vector< Real > &x, Real del, TrustRegionModel_U< Real > &model, Vector< Real > &ymin, Vector< Real > &pwa, Vector< Real > &pwa1, Vector< Real > &pwa2, Vector< Real > &pwa3, Vector< Real > &pwa4, Vector< Real > &pwa5, Vector< Real > &dwa, std::ostream &outStream=std::cout)
Real computeGradient(const Vector< Real > &x, Vector< Real > &g, Vector< Real > &pwa, Real del, Objective< Real > &obj, std::ostream &outStream=std::cout) const
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...
void writeOutput(std::ostream &os, bool write_header=false) const override
Print iterate status.
TrustRegionSPGAlgorithm(ParameterList &list, const Ptr< Secant< Real > > &secant=nullPtr)
void writeName(std::ostream &os) const override
Print step name.
void writeHeader(std::ostream &os) const override
Print iterate header.
void dproj(Vector< Real > &x, const Vector< Real > &x0, Real del, Vector< Real > &y0, Vector< Real > &y1, Vector< Real > &yc, Vector< Real > &pwa, std::ostream &outStream=std::cout) const
Real dgpstep(Vector< Real > &s, const Vector< Real > &w, const Vector< Real > &x, const Real alpha, std::ostream &outStream=std::cout) const
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
virtual Real apply(const Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
Definition: ROL_Vector.hpp:238
virtual Real norm() const =0
Returns where .
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
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,...
Definition: ROL_Vector.hpp:226
virtual void plus(const Vector &x)=0
Compute , where .
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 .
Definition: ROL_Vector.hpp:153
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)
Definition: ROL_Types.hpp:541
ESecantMode
Definition: ROL_Secant.hpp:57
@ SECANTMODE_FORWARD
Definition: ROL_Secant.hpp:58
@ SECANTMODE_INVERSE
Definition: ROL_Secant.hpp:59
@ SECANTMODE_BOTH
Definition: ROL_Secant.hpp:60