ROL
ROL_TypeB_KelleySachsAlgorithm_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_KELLEYSACHSALGORITHM_DEF_HPP
45 #define ROL_TYPEB_KELLEYSACHSALGORITHM_DEF_HPP
46 
47 namespace ROL {
48 namespace TypeB {
49 
50 template<typename Real>
51 KelleySachsAlgorithm<Real>::KelleySachsAlgorithm(ParameterList &list,
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  // Krylov Parameters
70  maxit_ = list.sublist("General").sublist("Krylov").get("Iteration Limit", 20);
71  tol1_ = list.sublist("General").sublist("Krylov").get("Absolute Tolerance", 1e-4);
72  tol2_ = list.sublist("General").sublist("Krylov").get("Relative Tolerance", 1e-2);
73  // Algorithm-Specific Parameters
74  minit_ = trlist.sublist("Kelley-Sachs").get("Maximum Number of Smoothing Iterations", 20);
75  mu0_ = trlist.sublist("Kelley-Sachs").get("Sufficient Decrease Parameter", 1e-4);
76  mu1_ = trlist.sublist("Kelley-Sachs").get("Post-Smoothing Decrease Parameter", 0.9999);
77  eps0_ = trlist.sublist("Kelley-Sachs").get("Binding Set Tolerance", 1e-3);
78  beta_ = trlist.sublist("Kelley-Sachs").get("Post-Smoothing Backtracking Rate", 1e-2);
79  alpha0_ = trlist.sublist("Kelley-Sachs").get("Initial Post-Smoothing Step Size", 1.0);
80  // Output Parameters
81  verbosity_ = list.sublist("General").get("Output Level",0);
82  writeHeader_ = verbosity_ > 2;
83  // Secant Information
84  useSecantPrecond_ = list.sublist("General").sublist("Secant").get("Use as Preconditioner", false);
85  useSecantHessVec_ = list.sublist("General").sublist("Secant").get("Use as Hessian", false);
86  // Initialize trust region model
87  model_ = makePtr<TrustRegionModel_U<Real>>(list,secant);
88  useSecantPrecond_ = list.sublist("General").sublist("Secant").get("Use as Preconditioner", false);
89  useSecantHessVec_ = list.sublist("General").sublist("Secant").get("Use as Hessian", false);
90  if (secant == nullPtr) {
91  esec_ = StringToESecant(list.sublist("General").sublist("Secant").get("Type","Limited-Memory BFGS"));
92  }
93 }
94 
95 template<typename Real>
96 void KelleySachsAlgorithm<Real>::initialize(Vector<Real> &x,
97  const Vector<Real> &g,
98  Objective<Real> &obj,
99  BoundConstraint<Real> &bnd,
100  std::ostream &outStream) {
101  const Real one(1);
102  // Initialize data
103  TypeB::Algorithm<Real>::initialize(x,g);
104  nhess_ = 0;
105  // Update approximate gradient and approximate objective function.
106  Real ftol = static_cast<Real>(0.1)*ROL_OVERFLOW<Real>();
107  bnd.project(x);
108  state_->iterateVec->set(x);
109  obj.update(x,UpdateType::Initial,state_->iter);
110  state_->value = obj.value(x,ftol);
111  state_->nfval++;
112  obj.gradient(*state_->gradientVec,x,ftol);
113  state_->ngrad++;
114  state_->stepVec->set(x);
115  state_->stepVec->axpy(-one,state_->gradientVec->dual());
116  bnd.project(*state_->stepVec);
117  state_->stepVec->axpy(-one,x);
118  state_->gnorm = state_->stepVec->norm();
119  state_->snorm = ROL_INF<Real>();
120  // Compute initial trust region radius if desired.
121  if ( state_->searchSize <= static_cast<Real>(0) ) {
122  state_->searchSize = state_->gradientVec->norm();
123  }
124 }
125 
126 template<typename Real>
127 void KelleySachsAlgorithm<Real>::run(Vector<Real> &x,
128  const Vector<Real> &g,
129  Objective<Real> &obj,
130  BoundConstraint<Real> &bnd,
131  std::ostream &outStream ) {
132  const Real one(1), xeps(ROL_EPSILON<Real>());
133  Real ftol = std::sqrt(ROL_EPSILON<Real>());
134  Real gfnorm(0), gfnormf(0), tol(0), stol(0), eps(0);
135  Real ftrial(0), fnew(0), pRed(0), rho(1), alpha(1), lambda(1);
136  int cnt(0);
137  // Initialize trust-region data
138  initialize(x,g,obj,bnd,outStream);
139  Ptr<Vector<Real>> s = x.clone(), gfree = g.clone();
140  Ptr<Vector<Real>> pwa1 = x.clone(), pwa2 = x.clone(), pwa3 = x.clone();
141  Ptr<Vector<Real>> dwa1 = g.clone(), dwa2 = g.clone(), dwa3 = g.clone();
142 
143  // Output
144  if (verbosity_ > 0) writeOutput(outStream,true);
145 
146  // Compute initial free gradient
147  gfree->set(*state_->gradientVec);
148  //bnd.pruneActive(*gfree,*state_->gradientVec,x,xeps,eps);
149  //gfnorm = gfree->norm();
150  pwa1->set(gfree->dual());
151  bnd.pruneActive(*pwa1,state_->gradientVec->dual(),x,xeps,eps);
152  gfree->set(pwa1->dual());
153  gfnorm = gfree->norm();
154 
155  // Initial tolerances
156  eps = std::min(eps0_,std::sqrt(state_->gnorm));
157  tol = tol2_*std::sqrt(state_->gnorm);
158  stol = std::min(tol1_,tol*gfnorm);
159 
160  while (status_->check(*state_)) {
161  // Build trust-region model
162  model_->setData(obj,*state_->iterateVec,*state_->gradientVec,ftol);
163 
164  /**** SOLVE TRUST-REGION SUBPROBLEM ****/
165  pRed = trpcg(*s,SPflag_,SPiter_,*gfree,x,*state_->gradientVec,
166  state_->searchSize,*model_,bnd,eps,stol,maxit_,
167  *pwa1,*dwa1,*pwa2,*dwa2,*pwa3,*dwa3,outStream);
168  x.plus(*s);
169  bnd.project(x);
170 
171  // Compute trial objective value
172  obj.update(x,UpdateType::Trial);
173  ftrial = obj.value(x,ftol); state_->nfval++;
174 
175  // Compute ratio of acutal and predicted reduction
176  TRflag_ = TRUtils::SUCCESS;
177  TRUtils::analyzeRatio<Real>(rho,TRflag_,state_->value,ftrial,pRed,eps_,outStream,verbosity_>1);
178 
179  // Check sufficient decrease
180  if ( rho >= eta0_ ) {
181  lambda = std::min(one, state_->searchSize/gfnorm);
182  // Compute Scaled Measure || x - P( x - lam * PI(g) ) ||
183  pwa1->set(*state_->iterateVec);
184  pwa1->axpy(-lambda,gfree->dual());
185  bnd.project(*pwa1);
186  pwa1->axpy(-one,*state_->iterateVec);
187  gfnormf = pwa1->norm();
188  // Sufficient decrease?
189  if (state_->value-ftrial < mu0_*gfnormf*state_->gnorm) {
190  TRflag_ = TRUtils::QMINSUFDEC;
191  }
192  if ( verbosity_ > 1 ) {
193  outStream << " Norm of projected free gradient: " << gfnormf << std::endl;
194  outStream << " Decrease lower bound (constraints): " << mu0_*gfnormf*state_->gnorm << std::endl;
195  outStream << " Trust-region flag (constraints): " << TRflag_ << std::endl;
196  outStream << std::endl;
197  }
198  }
199 
200  // Update algorithm state
201  state_->iter++;
202  // Accept/reject step and update trust region radius
203  if ((rho < eta0_ && TRflag_ == TRUtils::SUCCESS) || (TRflag_ >= 2)) {
204  state_->stepVec->set(x);
205  state_->stepVec->axpy(-one,*state_->iterateVec);
206  state_->snorm = state_->stepVec->norm();
207  x.set(*state_->iterateVec);
208  obj.update(x,UpdateType::Revert,state_->iter);
209  // Decrease trust-region radius
210  state_->searchSize = gamma1_*std::min(state_->snorm,state_->searchSize);
211  }
212  else if ((rho >= eta0_ && TRflag_ != TRUtils::NPOSPREDNEG)
213  || (TRflag_ == TRUtils::POSPREDNEG)) {
214  if (rho >= eta0_ && rho < eta1_) {
215  // Decrease trust-region radius
216  state_->searchSize = gamma1_*std::min(state_->snorm,state_->searchSize);
217  }
218  else if (rho >= eta2_) {
219  // Increase trust-region radius
220  state_->searchSize = std::min(delMax_,gamma2_*state_->searchSize);
221  }
222  // Projected search
223  cnt = 0;
224  alpha = one;
225  obj.gradient(*dwa1,x,ftol); state_->ngrad++;
226  pwa2->set(dwa1->dual());
227  pwa1->set(x); pwa1->axpy(-alpha/alpha0_,*pwa2);
228  bnd.project(*pwa1);
229  obj.update(*pwa1,UpdateType::Trial);
230  fnew = obj.value(*pwa1,ftol); state_->nfval++;
231  while ((fnew-ftrial) >= mu1_*(state_->value-ftrial)) {
232  alpha *= beta_;
233  pwa1->set(x); pwa1->axpy(-alpha/alpha0_,*pwa2);
234  bnd.project(*pwa1);
235  obj.update(*pwa1,UpdateType::Trial);
236  fnew = obj.value(*pwa1,ftol); state_->nfval++;
237  if ( cnt >= minit_ ) break;
238  cnt++;
239  }
240  state_->stepVec->set(*pwa1);
241  state_->stepVec->axpy(-one,*state_->iterateVec);
242  state_->snorm = state_->stepVec->norm();
243  // Store objective function and iteration information
244  if (std::isnan(fnew)) {
245  TRflag_ = TRUtils::TRNAN;
246  rho = -one;
247  x.set(*state_->iterateVec);
248  obj.update(x,UpdateType::Revert,state_->iter);
249  // Decrease trust-region radius
250  state_->searchSize = gamma1_*std::min(state_->snorm,state_->searchSize);
251  }
252  else {
253  // Update objective function information
254  x.set(*pwa1);
255  state_->iterateVec->set(x);
256  state_->value = fnew;
257  dwa1->set(*state_->gradientVec);
258  obj.update(x,UpdateType::Accept,state_->iter);
259  obj.gradient(*state_->gradientVec,x,ftol); state_->ngrad++;
260  // Compute free gradient
261  gfree->set(*state_->gradientVec);
262  //bnd.pruneActive(*gfree,*state_->gradientVec,x,xeps,eps);
263  //gfnorm = gfree->norm();
264  pwa1->set(gfree->dual());
265  bnd.pruneActive(*pwa1,state_->gradientVec->dual(),x,xeps,eps);
266  gfree->set(pwa1->dual());
267  gfnorm = gfree->norm();
268  // Compute criticality measure
269  pwa1->set(x);
270  pwa1->axpy(-one,state_->gradientVec->dual());
271  bnd.project(*pwa1);
272  pwa1->axpy(-one,x);
273  state_->gnorm = pwa1->norm();
274  // Update secant information in trust-region model
275  model_->update(x,*state_->stepVec,*dwa1,*state_->gradientVec,
276  state_->snorm,state_->iter);
277  // Update tolerances
278  eps = std::min(eps0_,std::sqrt(state_->gnorm));
279  tol = tol2_*std::sqrt(state_->gnorm);
280  stol = std::min(tol1_,tol*gfnorm);
281  }
282  }
283 
284  // Update Output
285  if (verbosity_ > 0) writeOutput(outStream,writeHeader_);
286  }
287  if (verbosity_ > 0) TypeB::Algorithm<Real>::writeExitStatus(outStream);
288 }
289 
290 template<typename Real>
291 void KelleySachsAlgorithm<Real>::run( Problem<Real> &problem,
292  std::ostream &outStream ) {
293  if (problem.getPolyhedralProjection() == nullPtr) {
294  TypeB::Algorithm<Real>::run(problem,outStream);
295  }
296  else {
297  throw Exception::NotImplemented(">>> TypeB::KelleySachsAlgorithm::run : This algorithm cannot solve problems with linear equality constraints!");
298  }
299 }
300 
301 template<typename Real>
302 void KelleySachsAlgorithm<Real>::run( Vector<Real> &x,
303  const Vector<Real> &g,
304  Objective<Real> &obj,
305  BoundConstraint<Real> &bnd,
306  Constraint<Real> &linear_econ,
307  Vector<Real> &linear_emul,
308  const Vector<Real> &linear_eres,
309  std::ostream &outStream ) {
310  throw Exception::NotImplemented(">>> TypeB::KelleySachsAlgorithm::run : This algorithm cannot solve problems with linear equality constraints!");
311 }
312 
313 template<typename Real>
314 Real KelleySachsAlgorithm<Real>::trqsol(const Real xtx,
315  const Real ptp,
316  const Real ptx,
317  const Real del) const {
318  const Real zero(0);
319  Real dsq = del*del;
320  Real rad = ptx*ptx + ptp*(dsq-xtx);
321  rad = std::sqrt(std::max(rad,zero));
322  Real sigma(0);
323  if (ptx > zero) {
324  sigma = (dsq-xtx)/(ptx+rad);
325  }
326  else if (rad > zero) {
327  sigma = (rad-ptx)/ptp;
328  }
329  else {
330  sigma = zero;
331  }
332  return sigma;
333 }
334 
335 template<typename Real>
336 Real KelleySachsAlgorithm<Real>::trpcg(Vector<Real> &w, int &iflag, int &iter,
337  const Vector<Real> &g, const Vector<Real> &x,
338  const Vector<Real> &g0,
339  const Real del, TrustRegionModel_U<Real> &model,
340  BoundConstraint<Real> &bnd, const Real eps,
341  const Real tol, const int itermax,
342  Vector<Real> &p, Vector<Real> &q, Vector<Real> &r,
343  Vector<Real> &t, Vector<Real> &pwa, Vector<Real> &dwa,
344  std::ostream &outStream) const {
345  // p = step (primal)
346  // q = hessian applied to step p (dual)
347  // t = gradient (dual)
348  // r = preconditioned gradient (primal)
349  Real ftol = std::sqrt(ROL_EPSILON<Real>());
350  const Real zero(0), half(0.5), one(1), two(2);
351  Real rho(0), kappa(0), beta(0), sigma(0), alpha(0), pRed(0);
352  Real rtr(0), tnorm(0), rnorm0(0), sMs(0), pMp(0), sMp(0);
353  iter = 0; iflag = 0;
354  // Initialize step
355  w.zero();
356  // Compute residual
357  t.set(g); t.scale(-one);
358  // Preconditioned residual
359  applyFreePrecond(r,t,x,g0,eps,model,bnd,ftol,pwa,dwa);
360  //rho = r.dot(t.dual());
361  rho = r.apply(t);
362  rnorm0 = std::sqrt(rho);
363  if ( rnorm0 == zero ) {
364  return zero;
365  }
366  // Initialize direction
367  p.set(r);
368  pMp = rho;
369  // Iterate CG
370  for (iter = 0; iter < itermax; ++iter) {
371  // Apply Hessian to direction dir
372  applyFreeHessian(q,p,x,g0,eps,model,bnd,ftol,pwa,dwa);
373  // Compute sigma such that ||s+sigma*dir|| = del
374  //kappa = p.dot(q.dual());
375  kappa = p.apply(q);
376  alpha = (kappa>zero) ? rho/kappa : zero;
377  sigma = trqsol(sMs,pMp,sMp,del);
378  // Check for negative curvature or if iterate exceeds trust region
379  if (kappa <= zero || alpha >= sigma) {
380  w.axpy(sigma,p);
381  sMs = del*del;
382  iflag = (kappa<=zero) ? 2 : 3;
383  break;
384  }
385  pRed += half*alpha*rho;
386  // Update iterate and residuals
387  w.axpy(alpha,p);
388  t.axpy(-alpha,q);
389  applyFreePrecond(r,t,x,g0,eps,model,bnd,ftol,pwa,dwa);
390  // Exit if residual tolerance is met
391  //rtr = r.dot(t.dual());
392  rtr = r.apply(t);
393  tnorm = t.norm();
394  if (rtr <= tol*tol || tnorm <= tol) {
395  sMs = sMs + two*alpha*sMp + alpha*alpha*pMp;
396  iflag = 0;
397  break;
398  }
399  // Compute p = r + beta * p
400  beta = rtr/rho;
401  p.scale(beta); p.plus(r);
402  rho = rtr;
403  // Update dot products
404  // sMs = <s, inv(M)s> or <s, s> if equality constraint present
405  // sMp = <s, inv(M)p> or <s, p> if equality constraint present
406  // pMp = <p, inv(M)p> or <p, p> if equality constraint present
407  sMs = sMs + two*alpha*sMp + alpha*alpha*pMp;
408  sMp = beta*(sMp + alpha*pMp);
409  pMp = rho + beta*beta*pMp;
410  }
411  // Check iteration count
412  if (iter == itermax) {
413  iflag = 1;
414  }
415  if (iflag > 1) {
416  pRed += sigma*(rho-half*sigma*kappa);
417  }
418  if (iflag != 1) {
419  iter++;
420  }
421  return pRed;
422 }
423 
424 template<typename Real>
425 void KelleySachsAlgorithm<Real>::applyFreeHessian(Vector<Real> &hv,
426  const Vector<Real> &v,
427  const Vector<Real> &x,
428  const Vector<Real> &g,
429  Real eps,
430  TrustRegionModel_U<Real> &model,
431  BoundConstraint<Real> &bnd,
432  Real &tol,
433  Vector<Real> &pwa,
434  Vector<Real> &dwa) const {
435  const Real xeps(ROL_EPSILON<Real>());
436  pwa.set(v);
437  bnd.pruneActive(pwa,g.dual(),x,xeps,eps);
438  model.hessVec(hv,pwa,x,tol); nhess_++;
439  pwa.set(hv.dual());
440  bnd.pruneActive(pwa,g.dual(),x,xeps,eps);
441  hv.set(pwa.dual());
442  pwa.set(v);
443  bnd.pruneInactive(pwa,g.dual(),x,xeps,eps);
444  hv.plus(pwa.dual());
445  //pwa.set(v);
446  //bnd.pruneActive(pwa,g,x,xeps,eps);
447  //model.hessVec(hv,pwa,x,tol); nhess_++;
448  //bnd.pruneActive(hv,g,x,xeps,eps);
449  //pwa.set(v);
450  //bnd.pruneInactive(pwa,g,x,xeps,eps);
451  //dwa.set(pwa.dual());
452  //bnd.pruneInactive(dwa,g,x,xeps,eps);
453  //hv.plus(dwa);
454 }
455 
456 template<typename Real>
457 void KelleySachsAlgorithm<Real>::applyFreePrecond(Vector<Real> &hv,
458  const Vector<Real> &v,
459  const Vector<Real> &x,
460  const Vector<Real> &g,
461  Real eps,
462  TrustRegionModel_U<Real> &model,
463  BoundConstraint<Real> &bnd,
464  Real &tol,
465  Vector<Real> &pwa,
466  Vector<Real> &dwa) const {
467  const Real xeps(ROL_EPSILON<Real>());
468  pwa.set(v.dual());
469  bnd.pruneActive(pwa,g.dual(),x,xeps,eps);
470  dwa.set(pwa.dual());
471  model.precond(hv,dwa,x,tol);
472  bnd.pruneActive(hv,g.dual(),x,xeps,eps);
473  pwa.set(v.dual());
474  bnd.pruneInactive(pwa,g.dual(),x,xeps,eps);
475  hv.plus(pwa);
476  //dwa.set(v);
477  //bnd.pruneActive(dwa,g,x,xeps,eps);
478  //model.precond(hv,dwa,x,tol);
479  //bnd.pruneActive(hv,g,x,xeps,eps);
480  //dwa.set(v);
481  //bnd.pruneInactive(dwa,g,x,xeps,eps);
482  //pwa.set(dwa.dual());
483  //bnd.pruneInactive(pwa,g,x,xeps,eps);
484  //hv.plus(pwa);
485 }
486 
487 template<typename Real>
488 void KelleySachsAlgorithm<Real>::writeHeader( std::ostream& os ) const {
489  std::ios_base::fmtflags osFlags(os.flags());
490  if (verbosity_ > 1) {
491  os << std::string(114,'-') << std::endl;
492  os << " Kelley-Sachs trust-region method status output definitions" << std::endl << std::endl;
493  os << " iter - Number of iterates (steps taken)" << std::endl;
494  os << " value - Objective function value" << std::endl;
495  os << " gnorm - Norm of the gradient" << std::endl;
496  os << " snorm - Norm of the step (update to optimization vector)" << std::endl;
497  os << " delta - Trust-Region radius" << std::endl;
498  os << " #fval - Number of times the objective function was evaluated" << std::endl;
499  os << " #grad - Number of times the gradient was computed" << std::endl;
500  os << " #hess - Number of times the Hessian was applied" << std::endl;
501  os << std::endl;
502  os << " tr_flag - Trust-Region flag" << std::endl;
503  for( int flag = TRUtils::SUCCESS; flag != TRUtils::UNDEFINED; ++flag ) {
504  os << " " << NumberToString(flag) << " - "
505  << TRUtils::ETRFlagToString(static_cast<TRUtils::ETRFlag>(flag)) << std::endl;
506  }
507  os << std::endl;
508  os << " iterCG - Number of Truncated CG iterations" << std::endl << std::endl;
509  os << " flagGC - Trust-Region Truncated CG flag" << std::endl;
510  for( int flag = CG_FLAG_SUCCESS; flag != CG_FLAG_UNDEFINED; ++flag ) {
511  os << " " << NumberToString(flag) << " - "
512  << ECGFlagToString(static_cast<ECGFlag>(flag)) << std::endl;
513  }
514  os << std::string(114,'-') << std::endl;
515  }
516  os << " ";
517  os << std::setw(6) << std::left << "iter";
518  os << std::setw(15) << std::left << "value";
519  os << std::setw(15) << std::left << "gnorm";
520  os << std::setw(15) << std::left << "snorm";
521  os << std::setw(15) << std::left << "delta";
522  os << std::setw(10) << std::left << "#fval";
523  os << std::setw(10) << std::left << "#grad";
524  os << std::setw(10) << std::left << "#hess";
525  os << std::setw(10) << std::left << "tr_flag";
526  os << std::setw(10) << std::left << "iterCG";
527  os << std::setw(10) << std::left << "flagCG";
528  os << std::endl;
529  os.flags(osFlags);
530 }
531 
532 template<typename Real>
533 void KelleySachsAlgorithm<Real>::writeName( std::ostream& os ) const {
534  std::ios_base::fmtflags osFlags(os.flags());
535  os << std::endl << "Kelley-Sachs Trust-Region Method (Type B, Bound Constraints)" << std::endl;
536  os.flags(osFlags);
537 }
538 
539 template<typename Real>
540 void KelleySachsAlgorithm<Real>::writeOutput( std::ostream& os, bool write_header ) const {
541  std::ios_base::fmtflags osFlags(os.flags());
542  os << std::scientific << std::setprecision(6);
543  if ( state_->iter == 0 ) writeName(os);
544  if ( write_header ) writeHeader(os);
545  if ( state_->iter == 0 ) {
546  os << " ";
547  os << std::setw(6) << std::left << state_->iter;
548  os << std::setw(15) << std::left << state_->value;
549  os << std::setw(15) << std::left << state_->gnorm;
550  os << std::setw(15) << std::left << "---";
551  os << std::setw(15) << std::left << state_->searchSize;
552  os << std::setw(10) << std::left << state_->nfval;
553  os << std::setw(10) << std::left << state_->ngrad;
554  os << std::setw(10) << std::left << nhess_;
555  os << std::setw(10) << std::left << "---";
556  os << std::setw(10) << std::left << "---";
557  os << std::setw(10) << std::left << "---";
558  os << std::endl;
559  }
560  else {
561  os << " ";
562  os << std::setw(6) << std::left << state_->iter;
563  os << std::setw(15) << std::left << state_->value;
564  os << std::setw(15) << std::left << state_->gnorm;
565  os << std::setw(15) << std::left << state_->snorm;
566  os << std::setw(15) << std::left << state_->searchSize;
567  os << std::setw(10) << std::left << state_->nfval;
568  os << std::setw(10) << std::left << state_->ngrad;
569  os << std::setw(10) << std::left << nhess_;
570  os << std::setw(10) << std::left << TRflag_;
571  os << std::setw(10) << std::left << SPiter_;
572  os << std::setw(10) << std::left << SPflag_;
573  os << std::endl;
574  }
575  os.flags(osFlags);
576 }
577 
578 } // namespace TypeB
579 } // namespace ROL
580 
581 #endif
582 
583 // ORIGINAL KELLEY-SACHS ALGORITHM
584 //template<typename Real>
585 //std::vector<std::string> KelleySachsAlgorithm<Real>::run(Vector<Real> &x,
586 // const Vector<Real> &g,
587 // Objective<Real> &obj,
588 // BoundConstraint<Real> &bnd,
589 // std::ostream &outStream ) {
590 // const Real one(1);
591 // Real ftol = std::sqrt(ROL_EPSILON<Real>()), tol1(1e-2);
592 // Real gfnorm(0), gfnormf(0), tol(0), stol(0), eps(0);
593 // Real ftrial(0), fnew(0), pRed(0), rho(1), alpha(1), lambda(1);
594 // int rflag(0), cnt(0);
595 // // Initialize trust-region data
596 // initialize(x,g,obj,bnd,outStream);
597 // Ptr<Vector<Real>> s = x.clone(), gfree = g.clone();
598 // Ptr<Vector<Real>> pwa1 = x.clone(), pwa2 = x.clone(), pwa3 = x.clone();
599 // Ptr<Vector<Real>> dwa1 = g.clone(), dwa2 = g.clone(), dwa3 = g.clone();
600 //
601 // // Output
602 // if (verbosity_ > 0) writeOutput(outStream,true);
603 //
604 // // Initial tolerances
605 // rflag = 0;
606 // eps = std::min(eps0_,std::sqrt(state_->gnorm));
607 // tol = tol1*std::min(tol0_,std::sqrt(state_->gnorm));
608 //
609 // // Compute initial free gradient
610 // gfree->set(*state_->gradientVec);
611 // bnd.pruneActive(*gfree,*state_->gradientVec,x,eps);
612 // gfnorm = gfree->norm();
613 // stol = tol*gfnorm;
614 //
615 // while (status_->check(*state_)) {
616 // // Build trust-region model
617 // model_->setData(obj,*state_->iterateVec,*state_->gradientVec);
618 //
619 // /**** SOLVE TRUST-REGION SUBPROBLEM ****/
620 // pRed = trpcg(*s,SPflag_,SPiter_,*gfree,x,*state_->gradientVec,
621 // state_->searchSize,*model_,bnd,eps,stol,maxit_,
622 // *pwa1,*dwa1,*pwa2,*dwa2,*pwa3,*dwa3,outStream);
623 // x.plus(*s);
624 // bnd.project(x);
625 //
626 // // Compute trial objective value
627 // obj.update(x,false);
628 // ftrial = obj.value(x,ftol); state_->nfval++;
629 //
630 // // Compute ratio of acutal and predicted reduction
631 // TRflag_ = TRUtils::SUCCESS;
632 // TRUtils::analyzeRatio<Real>(rho,TRflag_,state_->value,ftrial,pRed,eps_,outStream,verbosity_>1);
633 //
634 // // Check sufficient decrease
635 // if ( rho >= eta0_ ) {
636 // lambda = std::min(one, state_->searchSize/gfnorm);
637 // // Compute Scaled Measure || x - P( x - lam * PI(g) ) ||
638 // pwa1->set(*state_->iterateVec);
639 // pwa1->axpy(-lambda,gfree->dual());
640 // bnd.project(*pwa1);
641 // pwa1->axpy(-one,*state_->iterateVec);
642 // gfnormf = pwa1->norm();
643 // // Sufficient decrease?
644 // if (state_->value-ftrial < mu0_*gfnormf*state_->gnorm) {
645 // TRflag_ = TRUtils::QMINSUFDEC;
646 // }
647 // if ( verbosity_ > 1 ) {
648 // outStream << " Norm of projected free gradient: " << gfnormf << std::endl;
649 // outStream << " Decrease lower bound (constraints): " << mu0_*gfnormf*state_->gnorm << std::endl;
650 // outStream << " Trust-region flag (constraints): " << TRflag_ << std::endl;
651 // outStream << std::endl;
652 // }
653 // }
654 //
655 // // Update algorithm state
656 // state_->iter++;
657 // // Accept/reject step and update trust region radius
658 // if ((rho < eta0_ && TRflag_ == TRUtils::SUCCESS) || (TRflag_ >= 2)) {
659 // state_->stepVec->set(x);
660 // state_->stepVec->axpy(-one,*state_->iterateVec);
661 // state_->snorm = state_->stepVec->norm();
662 // rflag = 1;
663 // x.set(*state_->iterateVec);
664 // obj.update(x,false,state_->iter);
665 // // Decrease trust-region radius
666 // state_->searchSize = gamma1_*state_->searchSize;
667 // }
668 // else if ((rho >= eta0_ && TRflag_ != TRUtils::NPOSPREDNEG)
669 // || (TRflag_ == TRUtils::POSPREDNEG)) {
670 // if (rho >= eta2_ && rflag == 0 && state_->searchSize != delMax_) {
671 // state_->stepVec->set(x);
672 // state_->stepVec->axpy(-one,*state_->iterateVec);
673 // state_->snorm = state_->stepVec->norm();
674 // x.set(*state_->iterateVec);
675 // obj.update(x,false,state_->iter);
676 // // Increase trust-region radius
677 // state_->searchSize = std::min(delMax_,gamma2_*state_->searchSize);
678 // }
679 // else {
680 // if (rho >= eta0_ && rho < eta1_) {
681 // // Decrease trust-region radius
682 // state_->searchSize = gamma1_*state_->searchSize;
683 // }
684 // // Projected search
685 // cnt = 0;
686 // alpha = one;
687 // obj.gradient(*dwa1,x,ftol); state_->ngrad++;
688 // pwa2->set(dwa1->dual());
689 // pwa1->set(x); pwa1->axpy(-alpha/alpha0_,*pwa2);
690 // bnd.project(*pwa1);
691 // obj.update(*pwa1,false);
692 // fnew = obj.value(*pwa1,ftol); state_->nfval++;
693 // while ((fnew-ftrial) >= mu1_*(state_->value-ftrial)) {
694 // alpha *= beta_;
695 // pwa1->set(x); pwa1->axpy(-alpha/alpha0_,*pwa2);
696 // bnd.project(*pwa1);
697 // obj.update(*pwa1,false);
698 // fnew = obj.value(*pwa1,ftol); state_->nfval++;
699 // if ( cnt >= minit_ ) break;
700 // cnt++;
701 // }
702 // state_->stepVec->set(*pwa1);
703 // state_->stepVec->axpy(-one,*state_->iterateVec);
704 // state_->snorm = state_->stepVec->norm();
705 // // Store objective function and iteration information
706 // if (std::isnan(fnew)) {
707 // TRflag_ = TRUtils::TRNAN;
708 // rho = -one;
709 // x.set(*state_->iterateVec);
710 // obj.update(x,false,state_->iter);
711 // // Decrease trust-region radius
712 // state_->searchSize = gamma1_*state_->searchSize;
713 // }
714 // else {
715 // // Update objective function information
716 // x.set(*pwa1);
717 // state_->iterateVec->set(x);
718 // state_->value = fnew;
719 // dwa1->set(*state_->gradientVec);
720 // obj.update(x,true,state_->iter);
721 // obj.gradient(*state_->gradientVec,x,ftol); state_->ngrad++;
722 // // Compute free gradient
723 // gfree->set(*state_->gradientVec);
724 // bnd.pruneActive(*gfree,*state_->gradientVec,x,eps);
725 // gfnorm = gfree->norm();
726 // stol = tol*gfnorm;
727 // // Compute criticality measure
728 // pwa1->set(x);
729 // pwa1->axpy(-one,state_->gradientVec->dual());
730 // bnd.project(*pwa1);
731 // pwa1->axpy(-one,x);
732 // state_->gnorm = pwa1->norm();
733 // // Update secant information in trust-region model
734 // model_->update(x,*state_->stepVec,*dwa1,*state_->gradientVec,
735 // state_->snorm,state_->iter);
736 // // Update tolerances
737 // rflag = 0;
738 // eps = std::min(eps0_,std::sqrt(state_->gnorm));
739 // tol = tol1*std::min(tol0_,std::sqrt(state_->gnorm));
740 // }
741 // }
742 // }
743 //
744 // // Update Output
745 // if (verbosity_ > 0) writeOutput(outStream,writeHeader_);
746 // }
747 // if (verbosity_ > 0) TypeB::Algorithm<Real>::writeExitStatus(outStream);
748 //}
ROL::ECGFlagToString
std::string ECGFlagToString(ECGFlag cgf)
Definition: ROL_Types.hpp:831
ROL::Objective
Provides the interface to evaluate objective functions.
Definition: ROL_Objective.hpp:78
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:226
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:238
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:153
ROL::Objective::value
virtual Real value(const Vector< Real > &x, Real &tol)=0
Compute value.
ROL::TypeB::KelleySachsAlgorithm::writeOutput
void writeOutput(std::ostream &os, const bool write_header=false) const override
Print iterate status.
Definition: ROL_TypeB_KelleySachsAlgorithm_Def.hpp:540
ROL::TypeB::Algorithm::writeExitStatus
virtual void writeExitStatus(std::ostream &os) const
Definition: ROL_TypeB_Algorithm_Def.hpp:304
ROL::StringToESecant
ESecant StringToESecant(std::string s)
Definition: ROL_Types.hpp:543
ROL::Vector::zero
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:167
ROL::Vector
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:80
ROL::Objective::update
virtual void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update objective function.
Definition: ROL_Objective.hpp:96
zero
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Definition: ROL_Objective_SerialSimOpt.hpp:145
ROL::TypeB::KelleySachsAlgorithm::KelleySachsAlgorithm
KelleySachsAlgorithm(ParameterList &list, const Ptr< Secant< Real >> &secant=nullPtr)
Definition: ROL_TypeB_KelleySachsAlgorithm_Def.hpp:51
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:184
ROL::TypeB::KelleySachsAlgorithm::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_KelleySachsAlgorithm_Def.hpp:127
ROL::Objective::gradient
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Definition: ROL_ObjectiveDef.hpp:74
ROL::NumberToString
std::string NumberToString(T Number)
Definition: ROL_Types.hpp:81
ROL::TypeB::KelleySachsAlgorithm::initialize
void initialize(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &bnd, std::ostream &outStream=std::cout)
Definition: ROL_TypeB_KelleySachsAlgorithm_Def.hpp:96
ROL::TrustRegionModel_U
Provides the interface to evaluate trust-region model functions.
Definition: ROL_TrustRegionModel_U.hpp:66
ROL::TypeB::KelleySachsAlgorithm::applyFreePrecond
void applyFreePrecond(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, const Vector< Real > &g, Real eps, TrustRegionModel_U< Real > &model, BoundConstraint< Real > &bnd, Real &tol, Vector< Real > &pwa, Vector< Real > &dwa) const
Definition: ROL_TypeB_KelleySachsAlgorithm_Def.hpp:457
ROL::TypeB::KelleySachsAlgorithm::applyFreeHessian
void applyFreeHessian(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, const Vector< Real > &g, Real eps, TrustRegionModel_U< Real > &model, BoundConstraint< Real > &bnd, Real &tol, Vector< Real > &pwa, Vector< Real > &dwa) const
Definition: ROL_TypeB_KelleySachsAlgorithm_Def.hpp:425
ROL::Secant
Provides interface for and implements limited-memory secant operators.
Definition: ROL_Secant.hpp:79
ROL::TypeB::KelleySachsAlgorithm::trqsol
Real trqsol(const Real xtx, const Real ptp, const Real ptx, const Real del) const
Definition: ROL_TypeB_KelleySachsAlgorithm_Def.hpp:314
ROL::StatusTest
Provides an interface to check status of optimization algorithms.
Definition: ROL_StatusTest.hpp:58
ROL::TRUtils::ETRFlagToString
std::string ETRFlagToString(ETRFlag trf)
Definition: ROL_TrustRegionUtilities.hpp:72
ROL::BoundConstraint::project
virtual void project(Vector< Real > &x)
Project optimization variables onto the bounds.
Definition: ROL_BoundConstraint_Def.hpp:74
ROL::TypeB::KelleySachsAlgorithm::writeHeader
void writeHeader(std::ostream &os) const override
Print iterate header.
Definition: ROL_TypeB_KelleySachsAlgorithm_Def.hpp:488
ROL::BoundConstraint::pruneInactive
void pruneInactive(Vector< Real > &v, const Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the -inactive set.
Definition: ROL_BoundConstraint_Def.hpp:263
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:203
ROL::TypeB::KelleySachsAlgorithm::writeName
void writeName(std::ostream &os) const override
Print step name.
Definition: ROL_TypeB_KelleySachsAlgorithm_Def.hpp:533
ROL::BoundConstraint
Provides the interface to apply upper and lower bound constraints.
Definition: ROL_BoundConstraint.hpp:73
ROL::Problem::getPolyhedralProjection
const Ptr< PolyhedralProjection< Real > > & getPolyhedralProjection()
Get the polyhedral projection object. This is a null pointer if no linear constraints and/or bounds a...
Definition: function/prox/ROL_Problem_Def.hpp:517
ROL::TypeB::Algorithm::initialize
void initialize(const Vector< Real > &x, const Vector< Real > &g)
Definition: ROL_TypeB_Algorithm_Def.hpp:63
ROL::TypeB::Algorithm::run
virtual void run(Problem< Real > &problem, std::ostream &outStream=std::cout)
Run algorithm on bound constrained problems (Type-B). This is the primary Type-B interface.
Definition: ROL_TypeB_Algorithm_Def.hpp:107
ROL::Vector::set
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
ROL::Vector::norm
virtual Real norm() const =0
Returns where .
ROL::TypeB::KelleySachsAlgorithm::trpcg
Real trpcg(Vector< Real > &w, int &iflag, int &iter, const Vector< Real > &g, const Vector< Real > &x, const Vector< Real > &g0, const Real del, TrustRegionModel_U< Real > &model, BoundConstraint< Real > &bnd, Real eps, const Real tol, 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_KelleySachsAlgorithm_Def.hpp:336
ROL::Constraint
Defines the general constraint operator interface.
Definition: ROL_Constraint.hpp:86
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:192
Generated by   doxygen 1.8.5