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