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 
49 namespace ROL {
50 namespace TypeB {
51 
52 template<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 
130 template<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 
165 template<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
183  obj.update(x,UpdateType::Trial);
184  Real fval = obj.value(x,outTol);
185  return fval;
186 }
187 
188 template<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 
215 template<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 
319 template<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 
329 template<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 
418 template<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 
527 template<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 
796 template<typename Real>
797 void 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 
841 template<typename Real>
842 void 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 
848 template<typename Real>
849 void 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
Provides the interface to evaluate objective functions.
virtual void scale(const Real alpha)=0
Compute where .
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
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 void plus(const Vector &x)=0
Compute , where .
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:153
virtual Real value(const Vector< Real > &x, Real &tol)=0
Compute value.
TrustRegionSPGAlgorithm(ParameterList &list, const Ptr< Secant< Real >> &secant=nullPtr)
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
ESecant StringToESecant(std::string s)
Definition: ROL_Types.hpp:541
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)
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:80
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)
void writeName(std::ostream &os) const override
Print step name.
virtual void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update objective function.
virtual Real dot(const Vector &x) const =0
Compute where .
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &s, Real &tol) override
Apply Hessian approximation to vector.
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Real optimalityCriterion(const Vector< Real > &x, const Vector< Real > &g, Vector< Real > &primal, std::ostream &outStream=std::cout) const
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
std::string NumberToString(T Number)
Definition: ROL_Types.hpp:81
Provides the interface to evaluate trust-region model functions.
void initialize(Vector< Real > &x, const Vector< Real > &g, Real ftol, Objective< Real > &obj, BoundConstraint< Real > &bnd, std::ostream &outStream=std::cout)
void writeHeader(std::ostream &os) const override
Print iterate header.
ESecantMode
Definition: ROL_Secant.hpp:57
Provides interface for and implements limited-memory secant operators.
Definition: ROL_Secant.hpp:79
Provides an interface to check status of optimization algorithms.
std::string ETRFlagToString(ETRFlag trf)
virtual void writeExitStatus(std::ostream &os) const
Provides the interface to apply upper and lower bound constraints.
Real computeGradient(const Vector< Real > &x, Vector< Real > &g, Vector< Real > &pwa, Real del, Objective< Real > &obj, std::ostream &outStream=std::cout) const
void writeOutput(std::ostream &os, bool write_header=false) const override
Print iterate status.
void initialize(const Vector< Real > &x, const Vector< Real > &g)
Real dgpstep(Vector< Real > &s, const Vector< Real > &w, const Vector< Real > &x, const Real alpha, std::ostream &outStream=std::cout) const
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
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...
virtual Real norm() const =0
Returns where .
Real computeValue(Real inTol, Real &outTol, Real pRed, Real &fold, int iter, const Vector< Real > &x, const Vector< Real > &xold, Objective< Real > &obj)