44 #ifndef ROL_TYPEB_LINMOREALGORITHM_DEF_HPP 45 #define ROL_TYPEB_LINMOREALGORITHM_DEF_HPP 50 template<
typename Real>
57 ParameterList &trlist = list.sublist(
"Step").sublist(
"Trust Region");
59 state_->searchSize = trlist.get(
"Initial Radius", -1.0);
60 delMax_ = trlist.get(
"Maximum Radius", ROL_INF<Real>());
61 eta0_ = trlist.get(
"Step Acceptance Threshold", 0.05);
62 eta1_ = trlist.get(
"Radius Shrinking Threshold", 0.05);
63 eta2_ = trlist.get(
"Radius Growing Threshold", 0.9);
64 gamma0_ = trlist.get(
"Radius Shrinking Rate (Negative rho)", 0.0625);
65 gamma1_ = trlist.get(
"Radius Shrinking Rate (Positive rho)", 0.25);
66 gamma2_ = trlist.get(
"Radius Growing Rate", 2.5);
67 TRsafe_ = trlist.get(
"Safeguard Size", 100.0);
68 eps_ = TRsafe_*ROL_EPSILON<Real>();
69 interpRad_ = trlist.get(
"Use Radius Interpolation",
false);
71 maxit_ = list.sublist(
"General").sublist(
"Krylov").get(
"Iteration Limit", 20);
72 tol1_ = list.sublist(
"General").sublist(
"Krylov").get(
"Absolute Tolerance", 1e-4);
73 tol2_ = list.sublist(
"General").sublist(
"Krylov").get(
"Relative Tolerance", 1e-2);
75 ROL::ParameterList &lmlist = trlist.sublist(
"Lin-More");
76 minit_ = lmlist.get(
"Maximum Number of Minor Iterations", 10);
77 mu0_ = lmlist.get(
"Sufficient Decrease Parameter", 1e-2);
78 spexp_ = lmlist.get(
"Relative Tolerance Exponent", 1.0);
79 spexp_ = std::max(static_cast<Real>(1),std::min(spexp_,static_cast<Real>(2)));
80 redlim_ = lmlist.sublist(
"Cauchy Point").get(
"Maximum Number of Reduction Steps", 10);
81 explim_ = lmlist.sublist(
"Cauchy Point").get(
"Maximum Number of Expansion Steps", 10);
82 alpha_ = lmlist.sublist(
"Cauchy Point").get(
"Initial Step Size", 1.0);
83 normAlpha_ = lmlist.sublist(
"Cauchy Point").get(
"Normalize Initial Step Size",
false);
84 interpf_ = lmlist.sublist(
"Cauchy Point").get(
"Reduction Rate", 0.1);
85 extrapf_ = lmlist.sublist(
"Cauchy Point").get(
"Expansion Rate", 10.0);
86 qtol_ = lmlist.sublist(
"Cauchy Point").get(
"Decrease Tolerance", 1e-8);
87 interpfPS_ = lmlist.sublist(
"Projected Search").get(
"Backtracking Rate", 0.5);
88 pslim_ = lmlist.sublist(
"Projected Search").get(
"Maximum Number of Steps", 20);
90 verbosity_ = list.sublist(
"General").get(
"Output Level",0);
91 writeHeader_ = verbosity_ > 2;
93 useSecantPrecond_ = list.sublist(
"General").sublist(
"Secant").get(
"Use as Preconditioner",
false);
94 useSecantHessVec_ = list.sublist(
"General").sublist(
"Secant").get(
"Use as Hessian",
false);
99 model_ = makePtr<TrustRegionModel_U<Real>>(list,secant,mode);
100 if (secant == nullPtr) {
101 esec_ =
StringToESecant(list.sublist(
"General").sublist(
"Secant").get(
"Type",
"Limited-Memory BFGS"));
105 template<
typename Real>
110 std::ostream &outStream) {
113 if (proj_ == nullPtr) {
114 proj_ = makePtr<PolyhedralProjection<Real>>(makePtrFromRef(bnd));
121 Real ftol =
static_cast<Real
>(0.1)*ROL_OVERFLOW<Real>();
122 proj_->project(x,outStream); state_->nproj++;
123 state_->iterateVec->set(x);
125 state_->value = obj.
value(x,ftol);
127 obj.
gradient(*state_->gradientVec,x,ftol);
129 state_->stepVec->set(x);
130 state_->stepVec->axpy(-one,state_->gradientVec->dual());
131 proj_->project(*state_->stepVec,outStream); state_->nproj++;
132 state_->stepVec->axpy(-one,x);
133 state_->gnorm = state_->stepVec->norm();
134 state_->snorm = ROL_INF<Real>();
137 alpha_ /= state_->gradientVec->norm();
140 if ( state_->searchSize <= static_cast<Real>(0) ) {
141 state_->searchSize = state_->gradientVec->norm();
145 rcon_ = makePtr<ReducedLinearConstraint<Real>>(proj_->getLinearConstraint(),
148 ns_ = makePtr<NullSpaceOperator<Real>>(rcon_,x,
149 *proj_->getResidual());
153 template<
typename Real>
158 std::ostream &outStream ) {
160 Real tol0 = std::sqrt(ROL_EPSILON<Real>());
161 Real gfnorm(0), gfnormf(0), tol(0), stol(0), snorm(0);
162 Real ftrial(0), pRed(0), rho(1), q(0);
163 int flagCG(0), iterCG(0), maxit(0);
165 initialize(x,g,obj,bnd,outStream);
166 Ptr<Vector<Real>> s = x.
clone();
167 Ptr<Vector<Real>> gmod = g.
clone(), gfree = g.
clone();
168 Ptr<Vector<Real>> pwa1 = x.
clone(), pwa2 = x.
clone(), pwa3 = x.
clone();
169 Ptr<Vector<Real>> dwa1 = g.
clone(), dwa2 = g.
clone(), dwa3 = g.
clone();
172 if (verbosity_ > 0) writeOutput(outStream,
true);
174 while (status_->check(*state_)) {
176 model_->setData(obj,*state_->iterateVec,*state_->gradientVec);
180 snorm = dcauchy(*state_->stepVec,alpha_,q,*state_->iterateVec,
181 state_->gradientVec->dual(),state_->searchSize,
182 *model_,*dwa1,*dwa2,outStream);
183 x.
plus(*state_->stepVec);
188 gmod->plus(*state_->gradientVec);
191 pwa1->set(gfree->dual());
193 gfree->set(pwa1->dual());
195 applyFreePrecond(*pwa1,*gfree,x,*model_,bnd,tol0,*dwa1,*pwa2);
196 gfnorm = pwa1->norm();
199 gfnorm = gfree->norm();
201 SPiter_ = 0; SPflag_ = 0;
202 if (verbosity_ > 1) {
203 outStream <<
" Norm of free gradient components: " << gfnorm << std::endl;
204 outStream << std::endl;
209 for (
int i = 0; i < minit_; ++i) {
211 flagCG = 0; iterCG = 0;
213 tol = std::min(tol1_,tol2_*std::pow(gfnorm,spexp_));
216 snorm = dtrpcg(*s,flagCG,iterCG,*gfree,x,
217 state_->searchSize,*model_,bnd,tol,stol,maxit,
218 *pwa1,*dwa1,*pwa2,*dwa2,*pwa3,*dwa3,outStream);
221 if (verbosity_ > 1) {
222 outStream <<
" Computation of CG step" << std::endl;
223 outStream <<
" Current face (i): " << i << std::endl;
224 outStream <<
" CG step length: " << snorm << std::endl;
225 outStream <<
" Number of CG iterations: " << iterCG << std::endl;
226 outStream <<
" CG flag: " << flagCG << std::endl;
227 outStream <<
" Total number of iterations: " << SPiter_ << std::endl;
228 outStream << std::endl;
232 snorm = dprsrch(x,*s,q,gmod->
dual(),*model_,bnd,*pwa1,*dwa1,outStream);
236 state_->stepVec->plus(*s);
240 pwa1->set(gfree->dual());
241 bnd.pruneActive(*pwa1,x,
zero);
242 gfree->set(pwa1->dual());
244 applyFreePrecond(*pwa1,*gfree,x,*model_,bnd,tol0,*dwa1,*pwa2);
245 gfnormf = pwa1->norm();
248 gfnormf = gfree->norm();
250 if (verbosity_ > 1) {
251 outStream <<
" Norm of free gradient components: " << gfnormf << std::endl;
252 outStream << std::endl;
257 if (gfnormf <= tol) {
261 else if (SPiter_ >= maxit_) {
265 else if (flagCG == 2) {
269 else if (flagCG == 3) {
277 state_->snorm = state_->stepVec->norm();
281 ftrial = obj.
value(x,tol0);
286 TRUtils::analyzeRatio<Real>(rho,TRflag_,state_->value,ftrial,pRed,eps_,outStream,verbosity_>1);
292 x.
set(*state_->iterateVec);
296 state_->searchSize = TRUtils::interpolateRadius<Real>(*state_->gradientVec,*state_->stepVec,
297 state_->snorm,pRed,state_->value,ftrial,state_->searchSize,gamma0_,gamma1_,eta2_,
298 outStream,verbosity_>1);
301 state_->searchSize = gamma1_*std::min(state_->snorm,state_->searchSize);
306 state_->value = ftrial;
309 if (rho >= eta2_) state_->searchSize = std::min(gamma2_*state_->searchSize, delMax_);
311 dwa1->set(*state_->gradientVec);
312 obj.
gradient(*state_->gradientVec,x,tol0);
315 state_->iterateVec->set(x);
317 model_->update(x,*state_->stepVec,*dwa1,*state_->gradientVec,
318 state_->snorm,state_->iter);
322 if (verbosity_ > 0) writeOutput(outStream,writeHeader_);
327 template<
typename Real>
330 std::ostream &outStream)
const {
332 proj_->project(s,outStream); state_->nproj++;
333 s.
axpy(static_cast<Real>(-1),x);
337 template<
typename Real>
346 std::ostream &outStream) {
347 const Real half(0.5);
349 Real tol = std::sqrt(ROL_EPSILON<Real>());
351 Real gs(0), snorm(0);
353 snorm = dgpstep(s,g,x,-alpha,outStream);
358 model.
hessVec(dwa,s,x,tol); nhess_++;
361 q = half * s.
apply(dwa) + gs;
362 interp = (q > mu0_*gs);
370 snorm = dgpstep(s,g,x,-alpha,outStream);
372 model.
hessVec(dwa,s,x,tol); nhess_++;
375 q = half * s.
apply(dwa) + gs;
376 search = (q > mu0_*gs) && (cnt < redlim_);
388 snorm = dgpstep(s,g,x,-alpha,outStream);
389 if (snorm <= del && cnt < explim_) {
390 model.
hessVec(dwa,s,x,tol); nhess_++;
393 q = half * s.
apply(dwa) + gs;
394 if (q <= mu0_*gs && std::abs(q-qs) > qtol_*std::abs(qs)) {
412 snorm = dgpstep(s,g,x,-alpha,outStream);
414 if (verbosity_ > 1) {
415 outStream <<
" Cauchy point" << std::endl;
416 outStream <<
" Step length (alpha): " << alpha << std::endl;
417 outStream <<
" Step length (alpha*g): " << snorm << std::endl;
418 outStream <<
" Model decrease (pRed): " << -q << std::endl;
420 outStream <<
" Number of extrapolation steps: " << cnt << std::endl;
426 template<
typename Real>
432 std::ostream &outStream) {
433 const Real half(0.5);
434 Real tol = std::sqrt(ROL_EPSILON<Real>());
435 Real beta(1), snorm(0), gs(0);
441 snorm = dgpstep(pwa,s,x,beta,outStream);
442 model.
hessVec(dwa,pwa,x,tol); nhess_++;
445 q = half * pwa.
apply(dwa) + gs;
446 if (q <= mu0_*gs || nsteps > pslim_) {
455 if (verbosity_ > 1) {
456 outStream << std::endl;
457 outStream <<
" Projected search" << std::endl;
458 outStream <<
" Step length (beta): " << beta << std::endl;
459 outStream <<
" Step length (beta*s): " << snorm << std::endl;
460 outStream <<
" Model decrease (pRed): " << -q << std::endl;
461 outStream <<
" Number of steps: " << nsteps << std::endl;
466 template<
typename Real>
470 const Real del)
const {
473 Real rad = ptx*ptx + ptp*(dsq-xtx);
474 rad = std::sqrt(std::max(rad,
zero));
477 sigma = (dsq-xtx)/(ptx+rad);
479 else if (rad >
zero) {
480 sigma = (rad-ptx)/ptp;
488 template<
typename Real>
493 const Real tol,
const Real stol,
const int itermax,
496 std::ostream &outStream)
const {
501 Real tol0 = std::sqrt(ROL_EPSILON<Real>());
502 const Real
zero(0), one(1), two(2);
503 Real rho(0), kappa(0), beta(0), sigma(0), alpha(0);
504 Real rtr(0), tnorm(0), sMs(0), pMp(0), sMp(0);
511 applyFreePrecond(r,t,x,model,bnd,tol0,dwa,pwa);
516 pMp = (!hasEcon_ ? rho : p.
dot(p));
518 for (iter = 0; iter < itermax; ++iter) {
520 applyFreeHessian(q,p,x,model,bnd,tol0,pwa);
524 alpha = (kappa>
zero) ? rho/kappa :
zero;
525 sigma = dtrqsol(sMs,pMp,sMp,del);
527 if (kappa <= zero || alpha >= sigma) {
530 iflag = (kappa<=
zero) ? 2 : 3;
536 applyFreePrecond(r,t,x,model,bnd,tol0,dwa,pwa);
541 if (rtr <= stol*stol || tnorm <= tol) {
542 sMs = sMs + two*alpha*sMp + alpha*alpha*pMp;
554 sMs = sMs + two*alpha*sMp + alpha*alpha*pMp;
555 sMp = beta*(sMp + alpha*pMp);
556 pMp = (!hasEcon_ ? rho : p.
dot(p)) + beta*beta*pMp;
559 if (iter == itermax) {
565 return std::sqrt(sMs);
568 template<
typename Real>
579 model.
hessVec(hv,pwa,x,tol); nhess_++;
589 template<
typename Real>
612 rcon_->setX(makePtrFromRef(x));
615 ns_->apply(hv,pwa,tol);
619 template<
typename Real>
621 std::stringstream hist;
622 if (verbosity_ > 1) {
623 hist << std::string(114,
'-') << std::endl;
624 hist <<
" Lin-More trust-region method status output definitions" << std::endl << std::endl;
625 hist <<
" iter - Number of iterates (steps taken)" << std::endl;
626 hist <<
" value - Objective function value" << std::endl;
627 hist <<
" gnorm - Norm of the gradient" << std::endl;
628 hist <<
" snorm - Norm of the step (update to optimization vector)" << std::endl;
629 hist <<
" delta - Trust-Region radius" << std::endl;
630 hist <<
" #fval - Number of times the objective function was evaluated" << std::endl;
631 hist <<
" #grad - Number of times the gradient was computed" << std::endl;
632 hist <<
" #hess - Number of times the Hessian was applied" << std::endl;
633 hist <<
" #proj - Number of times the projection was applied" << std::endl;
635 hist <<
" tr_flag - Trust-Region flag" << std::endl;
642 hist <<
" iterCG - Number of Truncated CG iterations" << std::endl << std::endl;
643 hist <<
" flagGC - Trust-Region Truncated CG flag" << std::endl;
649 hist << std::string(114,
'-') << std::endl;
652 hist << std::setw(6) << std::left <<
"iter";
653 hist << std::setw(15) << std::left <<
"value";
654 hist << std::setw(15) << std::left <<
"gnorm";
655 hist << std::setw(15) << std::left <<
"snorm";
656 hist << std::setw(15) << std::left <<
"delta";
657 hist << std::setw(10) << std::left <<
"#fval";
658 hist << std::setw(10) << std::left <<
"#grad";
659 hist << std::setw(10) << std::left <<
"#hess";
660 hist << std::setw(10) << std::left <<
"#proj";
661 hist << std::setw(10) << std::left <<
"tr_flag";
663 hist << std::setw(10) << std::left <<
"iterCG";
664 hist << std::setw(10) << std::left <<
"flagCG";
670 template<
typename Real>
672 std::stringstream hist;
673 hist << std::endl <<
"Lin-More Trust-Region Method (Type B, Bound Constraints)" << std::endl;
677 template<
typename Real>
679 std::stringstream hist;
680 hist << std::scientific << std::setprecision(6);
681 if ( state_->iter == 0 ) writeName(os);
682 if ( write_header ) writeHeader(os);
683 if ( state_->iter == 0 ) {
685 hist << std::setw(6) << std::left << state_->iter;
686 hist << std::setw(15) << std::left << state_->value;
687 hist << std::setw(15) << std::left << state_->gnorm;
688 hist << std::setw(15) << std::left <<
"---";
689 hist << std::setw(15) << std::left << state_->searchSize;
690 hist << std::setw(10) << std::left << state_->nfval;
691 hist << std::setw(10) << std::left << state_->ngrad;
692 hist << std::setw(10) << std::left << nhess_;
693 hist << std::setw(10) << std::left << state_->nproj;
694 hist << std::setw(10) << std::left <<
"---";
696 hist << std::setw(10) << std::left <<
"---";
697 hist << std::setw(10) << std::left <<
"---";
703 hist << std::setw(6) << std::left << state_->iter;
704 hist << std::setw(15) << std::left << state_->value;
705 hist << std::setw(15) << std::left << state_->gnorm;
706 hist << std::setw(15) << std::left << state_->snorm;
707 hist << std::setw(15) << std::left << state_->searchSize;
708 hist << std::setw(10) << std::left << state_->nfval;
709 hist << std::setw(10) << std::left << state_->ngrad;
710 hist << std::setw(10) << std::left << nhess_;
711 hist << std::setw(10) << std::left << state_->nproj;
712 hist << std::setw(10) << std::left << TRflag_;
714 hist << std::setw(10) << std::left << SPiter_;
715 hist << std::setw(10) << std::left << SPflag_;
std::string ECGFlagToString(ECGFlag cgf)
Provides the interface to evaluate objective functions.
virtual void scale(const Real alpha)=0
Compute where .
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual Real apply(const Vector< Real > &x) const
Apply to a dual vector. This is equivalent to the call .
Real dcauchy(Vector< Real > &s, Real &alpha, Real &q, const Vector< Real > &x, const Vector< Real > &g, const Real del, TrustRegionModel_U< Real > &model, Vector< Real > &dwa, Vector< Real > &dwa1, std::ostream &outStream=std::cout)
virtual void plus(const Vector &x)=0
Compute , where .
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
virtual Real value(const Vector< Real > &x, Real &tol)=0
Compute value.
LinMoreAlgorithm(ParameterList &list, const Ptr< Secant< Real >> &secant=nullPtr)
Real dtrpcg(Vector< Real > &w, int &iflag, int &iter, const Vector< Real > &g, const Vector< Real > &x, const Real del, TrustRegionModel_U< Real > &model, BoundConstraint< Real > &bnd, const Real tol, const Real stol, const int itermax, Vector< Real > &p, Vector< Real > &q, Vector< Real > &r, Vector< Real > &t, Vector< Real > &pwa, Vector< Real > &dwa, std::ostream &outStream=std::cout) const
ESecant StringToESecant(std::string s)
virtual void zero()
Set to zero vector.
Defines the linear algebra or vector space interface.
virtual void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update objective function.
virtual Real dot(const Vector &x) const =0
Compute where .
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &s, Real &tol) override
Apply Hessian approximation to vector.
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Real optimalityCriterion(const Vector< Real > &x, const Vector< Real > &g, Vector< Real > &primal, std::ostream &outStream=std::cout) const
Real dprsrch(Vector< Real > &x, Vector< Real > &s, Real &q, const Vector< Real > &g, TrustRegionModel_U< Real > &model, BoundConstraint< Real > &bnd, Vector< Real > &pwa, Vector< Real > &dwa, std::ostream &outStream=std::cout)
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
std::string NumberToString(T Number)
Provides the interface to evaluate trust-region model functions.
void applyFreeHessian(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, TrustRegionModel_U< Real > &model, BoundConstraint< Real > &bnd, Real &tol, Vector< Real > &pwa) const
void writeHeader(std::ostream &os) const override
Print iterate header.
Provides interface for and implements limited-memory secant operators.
void initialize(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &bnd, std::ostream &outStream=std::cout)
Provides an interface to check status of optimization algorithms.
std::string ETRFlagToString(ETRFlag trf)
void writeName(std::ostream &os) const override
Print step name.
virtual void writeExitStatus(std::ostream &os) const
void pruneActive(Vector< Real > &v, const Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the -active set.
Provides the interface to apply upper and lower bound constraints.
Real dtrqsol(const Real xtx, const Real ptp, const Real ptx, const Real del) const
Real dgpstep(Vector< Real > &s, const Vector< Real > &w, const Vector< Real > &x, const Real alpha, std::ostream &outStream=std::cout) const
void applyFreePrecond(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, TrustRegionModel_U< Real > &model, BoundConstraint< Real > &bnd, Real &tol, Vector< Real > &dwa, Vector< Real > &pwa) const
void initialize(const Vector< Real > &x, const Vector< Real > &g)
virtual void set(const Vector &x)
Set where .
virtual Real norm() const =0
Returns where .
void run(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &bnd, std::ostream &outStream=std::cout) override
Run algorithm on bound constrained problems (Type-B). This general interface supports the use of dual...
void writeOutput(std::ostream &os, bool write_header=false) const override
Print iterate status.
virtual void precond(Vector< Real > &Pv, const Vector< Real > &v, const Vector< Real > &s, Real &tol) override
Apply preconditioner to vector.