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 
13 namespace ROL {
14 namespace TypeB {
15 
16 template<typename Real>
17 LinMoreAlgorithm<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 
92 template<typename Real>
93 void 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 
138 template<typename Real>
139 Real 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 
161 template<typename Real>
162 void 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 
246 template<typename Real>
247 void 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 
447 template<typename Real>
448 Real 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 
457 template<typename Real>
458 Real 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 
546 template<typename Real>
547 Real 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 
586 template<typename Real>
587 Real 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 
608 template<typename Real>
609 Real 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 
688 template<typename Real>
689 void 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 
709 template<typename Real>
710 void 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 
739 template<typename Real>
740 void 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 
790 template<typename Real>
791 void 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 
797 template<typename Real>
798 void 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
ROL::ECGFlagToString
std::string ECGFlagToString(ECGFlag cgf)
Definition: ROL_Types.hpp:801
ROL::Objective
Provides the interface to evaluate objective functions.
Definition: ROL_Objective.hpp:44
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::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::scale
virtual void scale(const Real alpha)=0
Compute where .
ROL::Vector::clone
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
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::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::Vector::plus
virtual void plus(const Vector &x)=0
Compute , where .
ROL::Vector::axpy
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:119
ROL::Objective::value
virtual Real value(const Vector< Real > &x, Real &tol)=0
Compute value.
ROL::TypeB::LinMoreAlgorithm::LinMoreAlgorithm
LinMoreAlgorithm(ParameterList &list, const Ptr< Secant< Real >> &secant=nullPtr)
Definition: ROL_TypeB_LinMoreAlgorithm_Def.hpp:17
ROL::TypeB::Algorithm::writeExitStatus
virtual void writeExitStatus(std::ostream &os) const
Definition: ROL_TypeB_Algorithm_Def.hpp:270
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::StringToESecant
ESecant StringToESecant(std::string s)
Definition: ROL_Types.hpp:513
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::Vector::zero
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:133
ROL::Vector
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
ROL::Objective::update
virtual void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update objective function.
Definition: ROL_Objective.hpp:62
ROL::Vector::dot
virtual Real dot(const Vector &x) const =0
Compute where .
zero
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Definition: ROL_Objective_SerialSimOpt.hpp:111
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::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::TRUtils::ETRFlag
ETRFlag
Enumation of flags used by trust-region solvers.
Definition: ROL_TrustRegionUtilities.hpp:28
ROL::Objective::gradient
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Definition: ROL_ObjectiveDef.hpp:40
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::NumberToString
std::string NumberToString(T Number)
Definition: ROL_Types.hpp:47
ROL::TrustRegionModel_U
Provides the interface to evaluate trust-region model functions.
Definition: ROL_TrustRegionModel_U.hpp:32
ROL::ESecantMode
ESecantMode
Definition: ROL_Secant.hpp:23
ROL::TypeB::LinMoreAlgorithm::writeHeader
void writeHeader(std::ostream &os) const override
Print iterate header.
Definition: ROL_TypeB_LinMoreAlgorithm_Def.hpp:740
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::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::Secant
Provides interface for and implements limited-memory secant operators.
Definition: ROL_Secant.hpp:45
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::StatusTest
Provides an interface to check status of optimization algorithms.
Definition: ROL_StatusTest.hpp:24
ROL::TRUtils::ETRFlagToString
std::string ETRFlagToString(ETRFlag trf)
Definition: ROL_TrustRegionUtilities.hpp:38
ROL::TypeB::LinMoreAlgorithm::writeName
void writeName(std::ostream &os) const override
Print step name.
Definition: ROL_TypeB_LinMoreAlgorithm_Def.hpp:791
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::BoundConstraint
Provides the interface to apply upper and lower bound constraints.
Definition: ROL_BoundConstraint.hpp:39
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::initialize
void initialize(const Vector< Real > &x, const Vector< Real > &g)
Definition: ROL_TypeB_Algorithm_Def.hpp:29
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::Vector::set
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:175
ROL::Vector::norm
virtual Real norm() const =0
Returns where .
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::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::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
Generated by   doxygen 1.8.5