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
47namespace ROL {
48namespace TypeB {
49
50template<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 // 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
95template<typename Real>
97 const Vector<Real> &g,
98 Objective<Real> &obj,
100 std::ostream &outStream) {
101 const Real one(1);
102 // Initialize data
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
126template<typename Real>
128 const Vector<Real> &g,
129 Objective<Real> &obj,
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);
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
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
290template<typename Real>
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
301template<typename Real>
303 const Vector<Real> &g,
304 Objective<Real> &obj,
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
313template<typename Real>
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
335template<typename Real>
336Real 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,
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
424template<typename Real>
426 const Vector<Real> &v,
427 const Vector<Real> &x,
428 const Vector<Real> &g,
429 Real eps,
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
456template<typename Real>
458 const Vector<Real> &v,
459 const Vector<Real> &x,
460 const Vector<Real> &g,
461 Real eps,
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
487template<typename Real>
488void KelleySachsAlgorithm<Real>::writeHeader( std::ostream& os ) const {
489 std::stringstream hist;
490 if (verbosity_ > 1) {
491 hist << std::string(114,'-') << std::endl;
492 hist << " Kelley-Sachs trust-region method status output definitions" << std::endl << std::endl;
493 hist << " iter - Number of iterates (steps taken)" << std::endl;
494 hist << " value - Objective function value" << std::endl;
495 hist << " gnorm - Norm of the gradient" << std::endl;
496 hist << " snorm - Norm of the step (update to optimization vector)" << std::endl;
497 hist << " delta - Trust-Region radius" << std::endl;
498 hist << " #fval - Number of times the objective function was evaluated" << std::endl;
499 hist << " #grad - Number of times the gradient was computed" << std::endl;
500 hist << " #hess - Number of times the Hessian was applied" << std::endl;
501 hist << std::endl;
502 hist << " tr_flag - Trust-Region flag" << std::endl;
503 for( int flag = TRUtils::SUCCESS; flag != TRUtils::UNDEFINED; ++flag ) {
504 hist << " " << NumberToString(flag) << " - "
505 << TRUtils::ETRFlagToString(static_cast<TRUtils::ETRFlag>(flag)) << std::endl;
506 }
507 hist << std::endl;
508 hist << " iterCG - Number of Truncated CG iterations" << std::endl << std::endl;
509 hist << " flagGC - Trust-Region Truncated CG flag" << std::endl;
510 for( int flag = CG_FLAG_SUCCESS; flag != CG_FLAG_UNDEFINED; ++flag ) {
511 hist << " " << NumberToString(flag) << " - "
512 << ECGFlagToString(static_cast<ECGFlag>(flag)) << std::endl;
513 }
514 hist << std::string(114,'-') << std::endl;
515 }
516 hist << " ";
517 hist << std::setw(6) << std::left << "iter";
518 hist << std::setw(15) << std::left << "value";
519 hist << std::setw(15) << std::left << "gnorm";
520 hist << std::setw(15) << std::left << "snorm";
521 hist << std::setw(15) << std::left << "delta";
522 hist << std::setw(10) << std::left << "#fval";
523 hist << std::setw(10) << std::left << "#grad";
524 hist << std::setw(10) << std::left << "#hess";
525 hist << std::setw(10) << std::left << "tr_flag";
526 hist << std::setw(10) << std::left << "iterCG";
527 hist << std::setw(10) << std::left << "flagCG";
528 hist << std::endl;
529 os << hist.str();
530}
531
532template<typename Real>
533void KelleySachsAlgorithm<Real>::writeName( std::ostream& os ) const {
534 std::stringstream hist;
535 hist << std::endl << "Kelley-Sachs Trust-Region Method (Type B, Bound Constraints)" << std::endl;
536 os << hist.str();
537}
538
539template<typename Real>
540void KelleySachsAlgorithm<Real>::writeOutput( std::ostream& os, bool write_header ) const {
541 std::stringstream hist;
542 hist << std::scientific << std::setprecision(6);
543 if ( state_->iter == 0 ) writeName(os);
544 if ( write_header ) writeHeader(os);
545 if ( state_->iter == 0 ) {
546 hist << " ";
547 hist << std::setw(6) << std::left << state_->iter;
548 hist << std::setw(15) << std::left << state_->value;
549 hist << std::setw(15) << std::left << state_->gnorm;
550 hist << std::setw(15) << std::left << "---";
551 hist << std::setw(15) << std::left << state_->searchSize;
552 hist << std::setw(10) << std::left << state_->nfval;
553 hist << std::setw(10) << std::left << state_->ngrad;
554 hist << std::setw(10) << std::left << nhess_;
555 hist << std::setw(10) << std::left << "---";
556 hist << std::setw(10) << std::left << "---";
557 hist << std::setw(10) << std::left << "---";
558 hist << std::endl;
559 }
560 else {
561 hist << " ";
562 hist << std::setw(6) << std::left << state_->iter;
563 hist << std::setw(15) << std::left << state_->value;
564 hist << std::setw(15) << std::left << state_->gnorm;
565 hist << std::setw(15) << std::left << state_->snorm;
566 hist << std::setw(15) << std::left << state_->searchSize;
567 hist << std::setw(10) << std::left << state_->nfval;
568 hist << std::setw(10) << std::left << state_->ngrad;
569 hist << std::setw(10) << std::left << nhess_;
570 hist << std::setw(10) << std::left << TRflag_;
571 hist << std::setw(10) << std::left << SPiter_;
572 hist << std::setw(10) << std::left << SPflag_;
573 hist << std::endl;
574 }
575 os << hist.str();
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//}
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0 zero)()
Provides the interface to apply upper and lower bound constraints.
void pruneInactive(Vector< Real > &v, const Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the -inactive set.
void pruneActive(Vector< Real > &v, const Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the -active set.
virtual void project(Vector< Real > &x)
Project optimization variables onto the bounds.
Defines the general constraint operator interface.
Provides the interface to evaluate objective functions.
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
virtual Real value(const Vector< Real > &x, Real &tol)=0
Compute value.
virtual void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update objective function.
const Ptr< PolyhedralProjection< Real > > & getPolyhedralProjection()
Get the polyhedral projection object. This is a null pointer if no linear constraints and/or bounds a...
Provides interface for and implements limited-memory secant operators.
Provides an interface to check status of optimization algorithms.
Provides the interface to evaluate trust-region model functions.
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &s, Real &tol) override
Apply Hessian approximation to vector.
virtual void precond(Vector< Real > &Pv, const Vector< Real > &v, const Vector< Real > &s, Real &tol) override
Apply preconditioner to vector.
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.
void initialize(const Vector< Real > &x, const Vector< Real > &g)
virtual void writeExitStatus(std::ostream &os) const
Real trqsol(const Real xtx, const Real ptp, const Real ptx, const Real del) const
void run(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &bnd, std::ostream &outStream=std::cout) override
Run algorithm on bound constrained problems (Type-B). This general interface supports the use of dual...
void writeOutput(std::ostream &os, bool write_header=false) const override
Print iterate status.
KelleySachsAlgorithm(ParameterList &list, const Ptr< Secant< Real > > &secant=nullPtr)
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
void initialize(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &bnd, std::ostream &outStream=std::cout)
void writeName(std::ostream &os) const override
Print step name.
void writeHeader(std::ostream &os) const override
Print iterate header.
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
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
Defines the linear algebra or vector space interface.
virtual Real apply(const Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
virtual Real norm() const =0
Returns where .
virtual void set(const Vector &x)
Set where .
virtual void scale(const Real alpha)=0
Compute where .
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis,...
virtual void plus(const Vector &x)=0
Compute , where .
virtual void zero()
Set to zero vector.
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
std::string ETRFlagToString(ETRFlag trf)
std::string NumberToString(T Number)
Definition ROL_Types.hpp:81
ESecant StringToESecant(std::string s)
@ CG_FLAG_UNDEFINED
@ CG_FLAG_SUCCESS
std::string ECGFlagToString(ECGFlag cgf)