ROL
ROL_TypeB_LinMoreAlgorithm_Def.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Rapid Optimization Library (ROL) Package
5 // Copyright (2014) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact lead developers:
38 // Drew Kouri (dpkouri@sandia.gov) and
39 // Denis Ridzal (dridzal@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
44 #ifndef ROL_TYPEB_LINMOREALGORITHM_DEF_HPP
45 #define ROL_TYPEB_LINMOREALGORITHM_DEF_HPP
46 
47 namespace ROL {
48 namespace TypeB {
49 
50 template<typename Real>
52  const Ptr<Secant<Real>> &secant) {
53  // Set status test
54  status_->reset();
55  status_->add(makePtr<StatusTest<Real>>(list));
56 
57  ParameterList &trlist = list.sublist("Step").sublist("Trust Region");
58  // Trust-Region Parameters
59  state_->searchSize = trlist.get("Initial Radius", -1.0);
60  delMax_ = trlist.get("Maximum Radius", ROL_INF<Real>());
61  eta0_ = trlist.get("Step Acceptance Threshold", 0.05);
62  eta1_ = trlist.get("Radius Shrinking Threshold", 0.05);
63  eta2_ = trlist.get("Radius Growing Threshold", 0.9);
64  gamma0_ = trlist.get("Radius Shrinking Rate (Negative rho)", 0.0625);
65  gamma1_ = trlist.get("Radius Shrinking Rate (Positive rho)", 0.25);
66  gamma2_ = trlist.get("Radius Growing Rate", 2.5);
67  TRsafe_ = trlist.get("Safeguard Size", 100.0);
68  eps_ = TRsafe_*ROL_EPSILON<Real>();
69  interpRad_ = trlist.get("Use Radius Interpolation", false);
70  // Krylov Parameters
71  maxit_ = list.sublist("General").sublist("Krylov").get("Iteration Limit", 20);
72  tol1_ = list.sublist("General").sublist("Krylov").get("Absolute Tolerance", 1e-4);
73  tol2_ = list.sublist("General").sublist("Krylov").get("Relative Tolerance", 1e-2);
74  // Algorithm-Specific Parameters
75  ROL::ParameterList &lmlist = trlist.sublist("Lin-More");
76  minit_ = lmlist.get("Maximum Number of Minor Iterations", 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  interpfPS_ = lmlist.sublist("Projected Search").get("Backtracking Rate", 0.5);
88  pslim_ = lmlist.sublist("Projected Search").get("Maximum Number of Steps", 20);
89  // Output Parameters
90  verbosity_ = list.sublist("General").get("Output Level",0);
91  writeHeader_ = verbosity_ > 2;
92  // Secant Information
93  useSecantPrecond_ = list.sublist("General").sublist("Secant").get("Use as Preconditioner", false);
94  useSecantHessVec_ = list.sublist("General").sublist("Secant").get("Use as Hessian", false);
96  if (useSecantPrecond_ && !useSecantHessVec_) mode = SECANTMODE_INVERSE;
97  else if (useSecantHessVec_ && !useSecantPrecond_) mode = SECANTMODE_FORWARD;
98  // Initialize trust region model
99  model_ = makePtr<TrustRegionModel_U<Real>>(list,secant,mode);
100  if (secant == nullPtr) {
101  esec_ = StringToESecant(list.sublist("General").sublist("Secant").get("Type","Limited-Memory BFGS"));
102  }
103 }
104 
105 template<typename Real>
107  const Vector<Real> &g,
108  Objective<Real> &obj,
110  std::ostream &outStream) {
111  const Real one(1);
112  hasEcon_ = true;
113  if (proj_ == nullPtr) {
114  proj_ = makePtr<PolyhedralProjection<Real>>(makePtrFromRef(bnd));
115  hasEcon_ = false;
116  }
117  // Initialize data
119  nhess_ = 0;
120  // Update approximate gradient and approximate objective function.
121  Real ftol = static_cast<Real>(0.1)*ROL_OVERFLOW<Real>();
122  proj_->project(x,outStream); state_->nproj++;
123  state_->iterateVec->set(x);
124  obj.update(x,UpdateType::Initial,state_->iter);
125  state_->value = obj.value(x,ftol);
126  state_->nfval++;
127  obj.gradient(*state_->gradientVec,x,ftol);
128  state_->ngrad++;
129  state_->stepVec->set(x);
130  state_->stepVec->axpy(-one,state_->gradientVec->dual());
131  proj_->project(*state_->stepVec,outStream); state_->nproj++;
132  state_->stepVec->axpy(-one,x);
133  state_->gnorm = state_->stepVec->norm();
134  state_->snorm = ROL_INF<Real>();
135  // Normalize initial CP step length
136  if (normAlpha_) {
137  alpha_ /= state_->gradientVec->norm();
138  }
139  // Compute initial trust region radius if desired.
140  if ( state_->searchSize <= static_cast<Real>(0) ) {
141  state_->searchSize = state_->gradientVec->norm();
142  }
143  // Initialize null space projection
144  if (hasEcon_) {
145  rcon_ = makePtr<ReducedLinearConstraint<Real>>(proj_->getLinearConstraint(),
146  makePtrFromRef(bnd),
147  makePtrFromRef(x));
148  ns_ = makePtr<NullSpaceOperator<Real>>(rcon_,x,
149  *proj_->getResidual());
150  }
151 }
152 
153 template<typename Real>
155  const Vector<Real> &g,
156  Objective<Real> &obj,
158  std::ostream &outStream ) {
159  const Real zero(0);
160  Real tol0 = std::sqrt(ROL_EPSILON<Real>());
161  Real gfnorm(0), gfnormf(0), tol(0), stol(0), snorm(0);
162  Real ftrial(0), pRed(0), rho(1), q(0);
163  int flagCG(0), iterCG(0), maxit(0);
164  // Initialize trust-region data
165  initialize(x,g,obj,bnd,outStream);
166  Ptr<Vector<Real>> s = x.clone();
167  Ptr<Vector<Real>> gmod = g.clone(), gfree = g.clone();
168  Ptr<Vector<Real>> pwa1 = x.clone(), pwa2 = x.clone(), pwa3 = x.clone();
169  Ptr<Vector<Real>> dwa1 = g.clone(), dwa2 = g.clone(), dwa3 = g.clone();
170 
171  // Output
172  if (verbosity_ > 0) writeOutput(outStream,true);
173 
174  while (status_->check(*state_)) {
175  // Build trust-region model
176  model_->setData(obj,*state_->iterateVec,*state_->gradientVec);
177 
178  /**** SOLVE TRUST-REGION SUBPROBLEM ****/
179  // Compute Cauchy point (TRON notation: x = x[1])
180  snorm = dcauchy(*state_->stepVec,alpha_,q,*state_->iterateVec,
181  state_->gradientVec->dual(),state_->searchSize,
182  *model_,*dwa1,*dwa2,outStream); // Solve 1D optimization problem for alpha
183  x.plus(*state_->stepVec); // Set x = x[0] + alpha*g
184  pRed = -q;
185 
186  // Model gradient at s = x[1] - x[0]
187  gmod->set(*dwa1); // hessVec from Cauchy point computation
188  gmod->plus(*state_->gradientVec);
189  gfree->set(*gmod);
190  //bnd.pruneActive(*gfree,x,zero);
191  pwa1->set(gfree->dual());
192  bnd.pruneActive(*pwa1,x,zero);
193  gfree->set(pwa1->dual());
194  if (hasEcon_) {
195  applyFreePrecond(*pwa1,*gfree,x,*model_,bnd,tol0,*dwa1,*pwa2);
196  gfnorm = pwa1->norm();
197  }
198  else {
199  gfnorm = gfree->norm();
200  }
201  SPiter_ = 0; SPflag_ = 0;
202  if (verbosity_ > 1) {
203  outStream << " Norm of free gradient components: " << gfnorm << std::endl;
204  outStream << std::endl;
205  }
206 
207  // Trust-region subproblem solve loop
208  maxit = maxit_;
209  for (int i = 0; i < minit_; ++i) {
210  // Run Truncated CG
211  flagCG = 0; iterCG = 0;
212  gfnormf = zero;
213  tol = std::min(tol1_,tol2_*std::pow(gfnorm,spexp_));
214  stol = tol; //zero;
215  if (gfnorm > zero) {
216  snorm = dtrpcg(*s,flagCG,iterCG,*gfree,x,
217  state_->searchSize,*model_,bnd,tol,stol,maxit,
218  *pwa1,*dwa1,*pwa2,*dwa2,*pwa3,*dwa3,outStream);
219  maxit -= iterCG;
220  SPiter_ += iterCG;
221  if (verbosity_ > 1) {
222  outStream << " Computation of CG step" << std::endl;
223  outStream << " Current face (i): " << i << std::endl;
224  outStream << " CG step length: " << snorm << std::endl;
225  outStream << " Number of CG iterations: " << iterCG << std::endl;
226  outStream << " CG flag: " << flagCG << std::endl;
227  outStream << " Total number of iterations: " << SPiter_ << std::endl;
228  outStream << std::endl;
229  }
230 
231  // Projected search
232  snorm = dprsrch(x,*s,q,gmod->dual(),*model_,bnd,*pwa1,*dwa1,outStream);
233  pRed += -q;
234 
235  // Model gradient at s = (x[i+1]-x[i]) - (x[i]-x[0])
236  state_->stepVec->plus(*s);
237  gmod->plus(*dwa1); // gmod = H(x[i+1]-x[i]) + H(x[i]-x[0]) + g
238  gfree->set(*gmod);
239  //bnd.pruneActive(*gfree,x,zero);
240  pwa1->set(gfree->dual());
241  bnd.pruneActive(*pwa1,x,zero);
242  gfree->set(pwa1->dual());
243  if (hasEcon_) {
244  applyFreePrecond(*pwa1,*gfree,x,*model_,bnd,tol0,*dwa1,*pwa2);
245  gfnormf = pwa1->norm();
246  }
247  else {
248  gfnormf = gfree->norm();
249  }
250  if (verbosity_ > 1) {
251  outStream << " Norm of free gradient components: " << gfnormf << std::endl;
252  outStream << std::endl;
253  }
254  }
255 
256  // Termination check
257  if (gfnormf <= tol) {
258  SPflag_ = 0;
259  break;
260  }
261  else if (SPiter_ >= maxit_) {
262  SPflag_ = 1;
263  break;
264  }
265  else if (flagCG == 2) {
266  SPflag_ = 2;
267  break;
268  }
269  else if (flagCG == 3) {
270  SPflag_ = 3;
271  break;
272  }
273 
274  // Update free gradient norm
275  gfnorm = gfnormf;
276  }
277  state_->snorm = state_->stepVec->norm();
278 
279  // Compute trial objective value
280  obj.update(x,UpdateType::Trial);
281  ftrial = obj.value(x,tol0);
282  state_->nfval++;
283 
284  // Compute ratio of acutal and predicted reduction
285  TRflag_ = TRUtils::SUCCESS;
286  TRUtils::analyzeRatio<Real>(rho,TRflag_,state_->value,ftrial,pRed,eps_,outStream,verbosity_>1);
287 
288  // Update algorithm state
289  state_->iter++;
290  // Accept/reject step and update trust region radius
291  if ((rho < eta0_ && TRflag_ == TRUtils::SUCCESS) || (TRflag_ >= 2)) { // Step Rejected
292  x.set(*state_->iterateVec);
293  obj.update(x,UpdateType::Revert,state_->iter);
294  if (interpRad_ && (rho < zero && TRflag_ != TRUtils::TRNAN)) {
295  // Negative reduction, interpolate to find new trust-region radius
296  state_->searchSize = TRUtils::interpolateRadius<Real>(*state_->gradientVec,*state_->stepVec,
297  state_->snorm,pRed,state_->value,ftrial,state_->searchSize,gamma0_,gamma1_,eta2_,
298  outStream,verbosity_>1);
299  }
300  else { // Shrink trust-region radius
301  state_->searchSize = gamma1_*std::min(state_->snorm,state_->searchSize);
302  }
303  }
304  else if ((rho >= eta0_ && TRflag_ != TRUtils::NPOSPREDNEG)
305  || (TRflag_ == TRUtils::POSPREDNEG)) { // Step Accepted
306  state_->value = ftrial;
307  obj.update(x,UpdateType::Accept,state_->iter);
308  // Increase trust-region radius
309  if (rho >= eta2_) state_->searchSize = std::min(gamma2_*state_->searchSize, delMax_);
310  // Compute gradient at new iterate
311  dwa1->set(*state_->gradientVec);
312  obj.gradient(*state_->gradientVec,x,tol0);
313  state_->ngrad++;
314  state_->gnorm = TypeB::Algorithm<Real>::optimalityCriterion(x,*state_->gradientVec,*pwa1,outStream);
315  state_->iterateVec->set(x);
316  // Update secant information in trust-region model
317  model_->update(x,*state_->stepVec,*dwa1,*state_->gradientVec,
318  state_->snorm,state_->iter);
319  }
320 
321  // Update Output
322  if (verbosity_ > 0) writeOutput(outStream,writeHeader_);
323  }
324  if (verbosity_ > 0) TypeB::Algorithm<Real>::writeExitStatus(outStream);
325 }
326 
327 template<typename Real>
329  const Vector<Real> &x, const Real alpha,
330  std::ostream &outStream) const {
331  s.set(x); s.axpy(alpha,w);
332  proj_->project(s,outStream); state_->nproj++;
333  s.axpy(static_cast<Real>(-1),x);
334  return s.norm();
335 }
336 
337 template<typename Real>
339  Real &alpha,
340  Real &q,
341  const Vector<Real> &x,
342  const Vector<Real> &g,
343  const Real del,
345  Vector<Real> &dwa, Vector<Real> &dwa1,
346  std::ostream &outStream) {
347  const Real half(0.5);
348  // const Real zero(0); // Unused
349  Real tol = std::sqrt(ROL_EPSILON<Real>());
350  bool interp = false;
351  Real gs(0), snorm(0);
352  // Compute s = P(x[0] - alpha g[0])
353  snorm = dgpstep(s,g,x,-alpha,outStream);
354  if (snorm > del) {
355  interp = true;
356  }
357  else {
358  model.hessVec(dwa,s,x,tol); nhess_++;
359  gs = s.dot(g);
360  //q = half * s.dot(dwa.dual()) + gs;
361  q = half * s.apply(dwa) + gs;
362  interp = (q > mu0_*gs);
363  }
364  // Either increase or decrease alpha to find approximate Cauchy point
365  int cnt = 0;
366  if (interp) {
367  bool search = true;
368  while (search) {
369  alpha *= interpf_;
370  snorm = dgpstep(s,g,x,-alpha,outStream);
371  if (snorm <= del) {
372  model.hessVec(dwa,s,x,tol); nhess_++;
373  gs = s.dot(g);
374  //q = half * s.dot(dwa.dual()) + gs;
375  q = half * s.apply(dwa) + gs;
376  search = (q > mu0_*gs) && (cnt < redlim_);
377  }
378  cnt++;
379  }
380  }
381  else {
382  bool search = true;
383  Real alphas = alpha;
384  Real qs = q;
385  dwa1.set(dwa);
386  while (search) {
387  alpha *= extrapf_;
388  snorm = dgpstep(s,g,x,-alpha,outStream);
389  if (snorm <= del && cnt < explim_) {
390  model.hessVec(dwa,s,x,tol); nhess_++;
391  gs = s.dot(g);
392  //q = half * s.dot(dwa.dual()) + gs;
393  q = half * s.apply(dwa) + gs;
394  if (q <= mu0_*gs && std::abs(q-qs) > qtol_*std::abs(qs)) {
395  dwa1.set(dwa);
396  search = true;
397  alphas = alpha;
398  qs = q;
399  }
400  else {
401  q = qs;
402  dwa.set(dwa1);
403  search = false;
404  }
405  }
406  else {
407  search = false;
408  }
409  cnt++;
410  }
411  alpha = alphas;
412  snorm = dgpstep(s,g,x,-alpha,outStream);
413  }
414  if (verbosity_ > 1) {
415  outStream << " Cauchy point" << std::endl;
416  outStream << " Step length (alpha): " << alpha << std::endl;
417  outStream << " Step length (alpha*g): " << snorm << std::endl;
418  outStream << " Model decrease (pRed): " << -q << std::endl;
419  if (!interp) {
420  outStream << " Number of extrapolation steps: " << cnt << std::endl;
421  }
422  }
423  return snorm;
424 }
425 
426 template<typename Real>
428  Real &q, const Vector<Real> &g,
431  Vector<Real> &pwa, Vector<Real> &dwa,
432  std::ostream &outStream) {
433  const Real half(0.5);
434  Real tol = std::sqrt(ROL_EPSILON<Real>());
435  Real beta(1), snorm(0), gs(0);
436  int nsteps = 0;
437  // Reduce beta until sufficient decrease is satisfied
438  bool search = true;
439  while (search) {
440  nsteps++;
441  snorm = dgpstep(pwa,s,x,beta,outStream);
442  model.hessVec(dwa,pwa,x,tol); nhess_++;
443  gs = pwa.dot(g);
444  //q = half * pwa.dot(dwa.dual()) + gs;
445  q = half * pwa.apply(dwa) + gs;
446  if (q <= mu0_*gs || nsteps > pslim_) {
447  search = false;
448  }
449  else {
450  beta *= interpfPS_;
451  }
452  }
453  s.set(pwa);
454  x.plus(s);
455  if (verbosity_ > 1) {
456  outStream << std::endl;
457  outStream << " Projected search" << std::endl;
458  outStream << " Step length (beta): " << beta << std::endl;
459  outStream << " Step length (beta*s): " << snorm << std::endl;
460  outStream << " Model decrease (pRed): " << -q << std::endl;
461  outStream << " Number of steps: " << nsteps << std::endl;
462  }
463  return snorm;
464 }
465 
466 template<typename Real>
468  const Real ptp,
469  const Real ptx,
470  const Real del) const {
471  const Real zero(0);
472  Real dsq = del*del;
473  Real rad = ptx*ptx + ptp*(dsq-xtx);
474  rad = std::sqrt(std::max(rad,zero));
475  Real sigma(0);
476  if (ptx > zero) {
477  sigma = (dsq-xtx)/(ptx+rad);
478  }
479  else if (rad > zero) {
480  sigma = (rad-ptx)/ptp;
481  }
482  else {
483  sigma = zero;
484  }
485  return sigma;
486 }
487 
488 template<typename Real>
489 Real LinMoreAlgorithm<Real>::dtrpcg(Vector<Real> &w, int &iflag, int &iter,
490  const Vector<Real> &g, const Vector<Real> &x,
491  const Real del, TrustRegionModel_U<Real> &model,
493  const Real tol, const Real stol, const int itermax,
495  Vector<Real> &t, Vector<Real> &pwa, Vector<Real> &dwa,
496  std::ostream &outStream) const {
497  // p = step (primal)
498  // q = hessian applied to step p (dual)
499  // t = gradient (dual)
500  // r = preconditioned gradient (primal)
501  Real tol0 = std::sqrt(ROL_EPSILON<Real>());
502  const Real zero(0), one(1), two(2);
503  Real rho(0), kappa(0), beta(0), sigma(0), alpha(0);
504  Real rtr(0), tnorm(0), sMs(0), pMp(0), sMp(0);
505  iter = 0; iflag = 0;
506  // Initialize step
507  w.zero();
508  // Compute residual
509  t.set(g); t.scale(-one);
510  // Preconditioned residual
511  applyFreePrecond(r,t,x,model,bnd,tol0,dwa,pwa);
512  //rho = r.dot(t.dual());
513  rho = r.apply(t);
514  // Initialize direction
515  p.set(r);
516  pMp = (!hasEcon_ ? rho : p.dot(p)); // If no equality constraint, used preconditioned norm
517  // Iterate CG
518  for (iter = 0; iter < itermax; ++iter) {
519  // Apply Hessian to direction dir
520  applyFreeHessian(q,p,x,model,bnd,tol0,pwa);
521  // Compute sigma such that ||s+sigma*dir|| = del
522  //kappa = p.dot(q.dual());
523  kappa = p.apply(q);
524  alpha = (kappa>zero) ? rho/kappa : zero;
525  sigma = dtrqsol(sMs,pMp,sMp,del);
526  // Check for negative curvature or if iterate exceeds trust region
527  if (kappa <= zero || alpha >= sigma) {
528  w.axpy(sigma,p);
529  sMs = del*del;
530  iflag = (kappa<=zero) ? 2 : 3;
531  break;
532  }
533  // Update iterate and residuals
534  w.axpy(alpha,p);
535  t.axpy(-alpha,q);
536  applyFreePrecond(r,t,x,model,bnd,tol0,dwa,pwa);
537  // Exit if residual tolerance is met
538  //rtr = r.dot(t.dual());
539  rtr = r.apply(t);
540  tnorm = t.norm();
541  if (rtr <= stol*stol || tnorm <= tol) {
542  sMs = sMs + two*alpha*sMp + alpha*alpha*pMp;
543  iflag = 0;
544  break;
545  }
546  // Compute p = r + beta * p
547  beta = rtr/rho;
548  p.scale(beta); p.plus(r);
549  rho = rtr;
550  // Update dot products
551  // sMs = <s, inv(M)s> or <s, s> if equality constraint present
552  // sMp = <s, inv(M)p> or <s, p> if equality constraint present
553  // pMp = <p, inv(M)p> or <p, p> if equality constraint present
554  sMs = sMs + two*alpha*sMp + alpha*alpha*pMp;
555  sMp = beta*(sMp + alpha*pMp);
556  pMp = (!hasEcon_ ? rho : p.dot(p)) + beta*beta*pMp;
557  }
558  // Check iteration count
559  if (iter == itermax) {
560  iflag = 1;
561  }
562  if (iflag != 1) {
563  iter++;
564  }
565  return std::sqrt(sMs); // w.norm();
566 }
567 
568 template<typename Real>
570  const Vector<Real> &v,
571  const Vector<Real> &x,
574  Real &tol,
575  Vector<Real> &pwa) const {
576  const Real zero(0);
577  pwa.set(v);
578  bnd.pruneActive(pwa,x,zero);
579  model.hessVec(hv,pwa,x,tol); nhess_++;
580  pwa.set(hv.dual());
581  bnd.pruneActive(pwa,x,zero);
582  hv.set(pwa.dual());
583  //pwa.set(v);
584  //bnd.pruneActive(pwa,x,zero);
585  //model.hessVec(hv,pwa,x,tol); nhess_++;
586  //bnd.pruneActive(hv,x,zero);
587 }
588 
589 template<typename Real>
591  const Vector<Real> &v,
592  const Vector<Real> &x,
595  Real &tol,
596  Vector<Real> &dwa,
597  Vector<Real> &pwa) const {
598  if (!hasEcon_) {
599  const Real zero(0);
600  pwa.set(v.dual());
601  bnd.pruneActive(pwa,x,zero);
602  dwa.set(pwa.dual());
603  model.precond(hv,dwa,x,tol);
604  bnd.pruneActive(hv,x,zero);
605  //dwa.set(v);
606  //bnd.pruneActive(dwa,x,zero);
607  //model.precond(hv,dwa,x,tol);
608  //bnd.pruneActive(hv,x,zero);
609  }
610  else {
611  // Perform null space projection
612  rcon_->setX(makePtrFromRef(x));
613  ns_->update(x);
614  pwa.set(v.dual());
615  ns_->apply(hv,pwa,tol);
616  }
617 }
618 
619 template<typename Real>
620 void LinMoreAlgorithm<Real>::writeHeader( std::ostream& os ) const {
621  std::stringstream hist;
622  if (verbosity_ > 1) {
623  hist << std::string(114,'-') << std::endl;
624  hist << " Lin-More trust-region method status output definitions" << std::endl << std::endl;
625  hist << " iter - Number of iterates (steps taken)" << std::endl;
626  hist << " value - Objective function value" << std::endl;
627  hist << " gnorm - Norm of the gradient" << std::endl;
628  hist << " snorm - Norm of the step (update to optimization vector)" << std::endl;
629  hist << " delta - Trust-Region radius" << std::endl;
630  hist << " #fval - Number of times the objective function was evaluated" << std::endl;
631  hist << " #grad - Number of times the gradient was computed" << std::endl;
632  hist << " #hess - Number of times the Hessian was applied" << std::endl;
633  hist << " #proj - Number of times the projection was applied" << std::endl;
634  hist << std::endl;
635  hist << " tr_flag - Trust-Region flag" << std::endl;
636  for( int flag = TRUtils::SUCCESS; flag != TRUtils::UNDEFINED; ++flag ) {
637  hist << " " << NumberToString(flag) << " - "
638  << TRUtils::ETRFlagToString(static_cast<TRUtils::ETRFlag>(flag)) << std::endl;
639  }
640  hist << std::endl;
641  if (minit_ > 0) {
642  hist << " iterCG - Number of Truncated CG iterations" << std::endl << std::endl;
643  hist << " flagGC - Trust-Region Truncated CG flag" << std::endl;
644  for( int flag = CG_FLAG_SUCCESS; flag != CG_FLAG_UNDEFINED; ++flag ) {
645  hist << " " << NumberToString(flag) << " - "
646  << ECGFlagToString(static_cast<ECGFlag>(flag)) << std::endl;
647  }
648  }
649  hist << std::string(114,'-') << std::endl;
650  }
651  hist << " ";
652  hist << std::setw(6) << std::left << "iter";
653  hist << std::setw(15) << std::left << "value";
654  hist << std::setw(15) << std::left << "gnorm";
655  hist << std::setw(15) << std::left << "snorm";
656  hist << std::setw(15) << std::left << "delta";
657  hist << std::setw(10) << std::left << "#fval";
658  hist << std::setw(10) << std::left << "#grad";
659  hist << std::setw(10) << std::left << "#hess";
660  hist << std::setw(10) << std::left << "#proj";
661  hist << std::setw(10) << std::left << "tr_flag";
662  if (minit_ > 0) {
663  hist << std::setw(10) << std::left << "iterCG";
664  hist << std::setw(10) << std::left << "flagCG";
665  }
666  hist << std::endl;
667  os << hist.str();
668 }
669 
670 template<typename Real>
671 void LinMoreAlgorithm<Real>::writeName( std::ostream& os ) const {
672  std::stringstream hist;
673  hist << std::endl << "Lin-More Trust-Region Method (Type B, Bound Constraints)" << std::endl;
674  os << hist.str();
675 }
676 
677 template<typename Real>
678 void LinMoreAlgorithm<Real>::writeOutput( std::ostream& os, bool write_header ) const {
679  std::stringstream hist;
680  hist << std::scientific << std::setprecision(6);
681  if ( state_->iter == 0 ) writeName(os);
682  if ( write_header ) writeHeader(os);
683  if ( state_->iter == 0 ) {
684  hist << " ";
685  hist << std::setw(6) << std::left << state_->iter;
686  hist << std::setw(15) << std::left << state_->value;
687  hist << std::setw(15) << std::left << state_->gnorm;
688  hist << std::setw(15) << std::left << "---";
689  hist << std::setw(15) << std::left << state_->searchSize;
690  hist << std::setw(10) << std::left << state_->nfval;
691  hist << std::setw(10) << std::left << state_->ngrad;
692  hist << std::setw(10) << std::left << nhess_;
693  hist << std::setw(10) << std::left << state_->nproj;
694  hist << std::setw(10) << std::left << "---";
695  if (minit_ > 0) {
696  hist << std::setw(10) << std::left << "---";
697  hist << std::setw(10) << std::left << "---";
698  }
699  hist << std::endl;
700  }
701  else {
702  hist << " ";
703  hist << std::setw(6) << std::left << state_->iter;
704  hist << std::setw(15) << std::left << state_->value;
705  hist << std::setw(15) << std::left << state_->gnorm;
706  hist << std::setw(15) << std::left << state_->snorm;
707  hist << std::setw(15) << std::left << state_->searchSize;
708  hist << std::setw(10) << std::left << state_->nfval;
709  hist << std::setw(10) << std::left << state_->ngrad;
710  hist << std::setw(10) << std::left << nhess_;
711  hist << std::setw(10) << std::left << state_->nproj;
712  hist << std::setw(10) << std::left << TRflag_;
713  if (minit_ > 0) {
714  hist << std::setw(10) << std::left << SPiter_;
715  hist << std::setw(10) << std::left << SPflag_;
716  }
717  hist << std::endl;
718  }
719  os << hist.str();
720 }
721 
722 } // namespace TypeB
723 } // namespace ROL
724 
725 #endif
std::string ECGFlagToString(ECGFlag cgf)
Definition: ROL_Types.hpp:829
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
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)
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.
LinMoreAlgorithm(ParameterList &list, const Ptr< Secant< Real >> &secant=nullPtr)
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
ESecant StringToESecant(std::string s)
Definition: ROL_Types.hpp:541
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:167
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:80
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
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)
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 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
ESecantMode
Definition: ROL_Secant.hpp:57
void writeHeader(std::ostream &os) const override
Print iterate header.
Provides interface for and implements limited-memory secant operators.
Definition: ROL_Secant.hpp:79
void initialize(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &bnd, std::ostream &outStream=std::cout)
Provides an interface to check status of optimization algorithms.
std::string ETRFlagToString(ETRFlag trf)
void writeName(std::ostream &os) const override
Print step name.
virtual void writeExitStatus(std::ostream &os) const
void pruneActive(Vector< Real > &v, const Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the -active set.
Provides the interface to apply upper and lower bound constraints.
Real dtrqsol(const Real xtx, const Real ptp, const Real ptx, const Real del) const
Real dgpstep(Vector< Real > &s, const Vector< Real > &w, const Vector< Real > &x, const Real alpha, std::ostream &outStream=std::cout) const
void applyFreePrecond(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, TrustRegionModel_U< Real > &model, BoundConstraint< Real > &bnd, Real &tol, Vector< Real > &dwa, Vector< Real > &pwa) const
void initialize(const Vector< Real > &x, const Vector< Real > &g)
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
virtual Real norm() const =0
Returns where .
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.
virtual void precond(Vector< Real > &Pv, const Vector< Real > &v, const Vector< Real > &s, Real &tol) override
Apply preconditioner to vector.