44#ifndef ROL_TYPEB_LINMOREALGORITHM_DEF_HPP
45#define ROL_TYPEB_LINMOREALGORITHM_DEF_HPP
50template<
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 storageNM_ = trlist.get(
"Nonmonotone Storage Size", 0);
72 useNM_ = (storageNM_ <= 0 ? false :
true);
74 maxit_ = list.sublist(
"General").sublist(
"Krylov").get(
"Iteration Limit", 20);
75 tol1_ = list.sublist(
"General").sublist(
"Krylov").get(
"Absolute Tolerance", 1e-4);
76 tol2_ = list.sublist(
"General").sublist(
"Krylov").get(
"Relative Tolerance", 1e-2);
78 ROL::ParameterList &lmlist = trlist.sublist(
"Lin-More");
79 minit_ = lmlist.get(
"Maximum Number of Minor Iterations", 10);
80 mu0_ = lmlist.get(
"Sufficient Decrease Parameter", 1e-2);
81 spexp_ = lmlist.get(
"Relative Tolerance Exponent", 1.0);
82 spexp_ = std::max(
static_cast<Real
>(1),std::min(spexp_,
static_cast<Real
>(2)));
83 redlim_ = lmlist.sublist(
"Cauchy Point").get(
"Maximum Number of Reduction Steps", 10);
84 explim_ = lmlist.sublist(
"Cauchy Point").get(
"Maximum Number of Expansion Steps", 10);
85 alpha_ = lmlist.sublist(
"Cauchy Point").get(
"Initial Step Size", 1.0);
86 normAlpha_ = lmlist.sublist(
"Cauchy Point").get(
"Normalize Initial Step Size",
false);
87 interpf_ = lmlist.sublist(
"Cauchy Point").get(
"Reduction Rate", 0.1);
88 extrapf_ = lmlist.sublist(
"Cauchy Point").get(
"Expansion Rate", 10.0);
89 qtol_ = lmlist.sublist(
"Cauchy Point").get(
"Decrease Tolerance", 1e-8);
90 interpfPS_ = lmlist.sublist(
"Projected Search").get(
"Backtracking Rate", 0.5);
91 pslim_ = lmlist.sublist(
"Projected Search").get(
"Maximum Number of Steps", 20);
93 verbosity_ = list.sublist(
"General").get(
"Output Level",0);
94 writeHeader_ = verbosity_ > 2;
96 useSecantPrecond_ = list.sublist(
"General").sublist(
"Secant").get(
"Use as Preconditioner",
false);
97 useSecantHessVec_ = list.sublist(
"General").sublist(
"Secant").get(
"Use as Hessian",
false);
102 model_ = makePtr<TrustRegionModel_U<Real>>(list,secant,mode);
103 if (secant == nullPtr) {
104 esec_ =
StringToESecant(list.sublist(
"General").sublist(
"Secant").get(
"Type",
"Limited-Memory BFGS"));
108template<
typename Real>
113 std::ostream &outStream) {
116 if (proj_ == nullPtr) {
117 proj_ = makePtr<PolyhedralProjection<Real>>(makePtrFromRef(bnd));
124 Real ftol =
static_cast<Real
>(0.1)*ROL_OVERFLOW<Real>();
125 proj_->project(x,outStream); state_->nproj++;
126 state_->iterateVec->set(x);
128 state_->value = obj.
value(x,ftol);
130 obj.
gradient(*state_->gradientVec,x,ftol);
132 state_->stepVec->set(x);
133 state_->stepVec->axpy(-one,state_->gradientVec->dual());
134 proj_->project(*state_->stepVec,outStream); state_->nproj++;
135 state_->stepVec->axpy(-one,x);
136 state_->gnorm = state_->stepVec->norm();
137 state_->snorm = ROL_INF<Real>();
140 alpha_ /= state_->gradientVec->norm();
143 if ( state_->searchSize <=
static_cast<Real
>(0) ) {
144 state_->searchSize = state_->gradientVec->norm();
148 rcon_ = makePtr<ReducedLinearConstraint<Real>>(proj_->getLinearConstraint(),
151 ns_ = makePtr<NullSpaceOperator<Real>>(rcon_,x,
152 *proj_->getResidual());
156template<
typename Real>
161 std::ostream &outStream ) {
163 Real tol0 = std::sqrt(ROL_EPSILON<Real>());
164 Real gfnorm(0), gfnormf(0), tol(0), stol(0), snorm(0);
165 Real ftrial(0), pRed(0), rho(1), q(0), delta(0);
166 int flagCG(0), iterCG(0), maxit(0);
168 initialize(x,g,obj,bnd,outStream);
169 Ptr<Vector<Real>> s = x.
clone();
170 Ptr<Vector<Real>> gmod = g.
clone(), gfree = g.
clone();
171 Ptr<Vector<Real>> pwa1 = x.
clone(), pwa2 = x.
clone(), pwa3 = x.
clone();
172 Ptr<Vector<Real>> dwa1 = g.
clone(), dwa2 = g.
clone(), dwa3 = g.
clone();
174 Real rhoNM(0), sigmac(0), sigmar(0);
175 Real fr(state_->value), fc(state_->value), fmin(state_->value);
180 if (verbosity_ > 0) writeOutput(outStream,
true);
182 while (status_->check(*state_)) {
184 model_->setData(obj,*state_->iterateVec,*state_->gradientVec);
188 snorm = dcauchy(*state_->stepVec,alpha_,q,*state_->iterateVec,
189 state_->gradientVec->dual(),state_->searchSize,
190 *model_,*dwa1,*dwa2,outStream);
191 x.
plus(*state_->stepVec);
192 state_->snorm = snorm;
193 delta = state_->searchSize - snorm;
198 gmod->plus(*state_->gradientVec);
201 pwa1->set(gfree->dual());
203 gfree->set(pwa1->dual());
205 applyFreePrecond(*pwa1,*gfree,x,*model_,bnd,tol0,*dwa1,*pwa2);
206 gfnorm = pwa1->norm();
209 gfnorm = gfree->norm();
211 SPiter_ = 0; SPflag_ = 0;
212 if (verbosity_ > 1) {
213 outStream <<
" Norm of free gradient components: " << gfnorm << std::endl;
214 outStream << std::endl;
219 for (
int i = 0; i < minit_; ++i) {
221 flagCG = 0; iterCG = 0;
223 tol = std::min(tol1_,tol2_*std::pow(gfnorm,spexp_));
225 if (gfnorm >
zero && delta >
zero) {
226 snorm = dtrpcg(*s,flagCG,iterCG,*gfree,x,
227 delta,*model_,bnd,tol,stol,maxit,
228 *pwa1,*dwa1,*pwa2,*dwa2,*pwa3,*dwa3,outStream);
231 if (verbosity_ > 1) {
232 outStream <<
" Computation of CG step" << std::endl;
233 outStream <<
" Current face (i): " << i << std::endl;
234 outStream <<
" CG step length: " << snorm << std::endl;
235 outStream <<
" Number of CG iterations: " << iterCG << std::endl;
236 outStream <<
" CG flag: " << flagCG << std::endl;
237 outStream <<
" Total number of iterations: " << SPiter_ << std::endl;
238 outStream << std::endl;
242 snorm = dprsrch(x,*s,q,gmod->
dual(),*model_,bnd,*pwa1,*dwa1,outStream);
246 state_->stepVec->plus(*s);
247 state_->snorm = state_->stepVec->norm();
248 delta = state_->searchSize - state_->snorm;
252 pwa1->set(gfree->dual());
254 gfree->set(pwa1->dual());
256 applyFreePrecond(*pwa1,*gfree,x,*model_,bnd,tol0,*dwa1,*pwa2);
257 gfnormf = pwa1->norm();
260 gfnormf = gfree->norm();
262 if (verbosity_ > 1) {
263 outStream <<
" Norm of free gradient components: " << gfnormf << std::endl;
264 outStream << std::endl;
269 if (gfnormf <= tol) {
273 else if (SPiter_ >= maxit_) {
277 else if (flagCG == 2) {
281 else if (delta <=
zero) {
293 ftrial = obj.
value(x,tol0);
298 TRUtils::analyzeRatio<Real>(rho,TRflag_,state_->value,ftrial,pRed,eps_,outStream,verbosity_>1);
300 TRUtils::analyzeRatio<Real>(rhoNM,TRflagNM,fr,ftrial,pRed+sigmar,eps_,outStream,verbosity_>1);
301 TRflag_ = (rho < rhoNM ? TRflagNM : TRflag_);
302 rho = (rho < rhoNM ? rhoNM : rho );
309 x.
set(*state_->iterateVec);
313 state_->searchSize = TRUtils::interpolateRadius<Real>(*state_->gradientVec,*state_->stepVec,
314 state_->snorm,pRed,state_->value,ftrial,state_->searchSize,gamma0_,gamma1_,eta2_,
315 outStream,verbosity_>1);
318 state_->searchSize = gamma1_*std::min(state_->snorm,state_->searchSize);
323 state_->value = ftrial;
326 sigmac += pRed; sigmar += pRed;
327 if (ftrial < fmin) { fmin = ftrial; fc = fmin; sigmac =
zero; L = 0; }
330 if (ftrial > fc) { fc = ftrial; sigmac =
zero; }
331 if (L == storageNM_) { fr = fc; sigmar = sigmac; }
335 if (rho >= eta2_) state_->searchSize = std::min(gamma2_*state_->searchSize, delMax_);
337 dwa1->set(*state_->gradientVec);
338 obj.
gradient(*state_->gradientVec,x,tol0);
341 state_->iterateVec->set(x);
343 model_->update(x,*state_->stepVec,*dwa1,*state_->gradientVec,
344 state_->snorm,state_->iter);
348 if (verbosity_ > 0) writeOutput(outStream,writeHeader_);
353template<
typename Real>
356 std::ostream &outStream)
const {
358 proj_->project(s,outStream); state_->nproj++;
359 s.
axpy(
static_cast<Real
>(-1),x);
363template<
typename Real>
372 std::ostream &outStream) {
373 const Real half(0.5);
375 Real tol = std::sqrt(ROL_EPSILON<Real>());
377 Real gs(0), snorm(0);
379 snorm = dgpstep(s,g,x,-alpha,outStream);
384 model.
hessVec(dwa,s,x,tol); nhess_++;
387 q = half * s.
apply(dwa) + gs;
388 interp = (q > mu0_*gs);
396 snorm = dgpstep(s,g,x,-alpha,outStream);
398 model.
hessVec(dwa,s,x,tol); nhess_++;
401 q = half * s.
apply(dwa) + gs;
402 search = (q > mu0_*gs) && (cnt < redlim_);
414 snorm = dgpstep(s,g,x,-alpha,outStream);
415 if (snorm <= del && cnt < explim_) {
416 model.
hessVec(dwa,s,x,tol); nhess_++;
419 q = half * s.
apply(dwa) + gs;
420 if (q <= mu0_*gs && std::abs(q-qs) > qtol_*std::abs(qs)) {
438 snorm = dgpstep(s,g,x,-alpha,outStream);
440 if (verbosity_ > 1) {
441 outStream <<
" Cauchy point" << std::endl;
442 outStream <<
" Step length (alpha): " << alpha << std::endl;
443 outStream <<
" Step length (alpha*g): " << snorm << std::endl;
444 outStream <<
" Model decrease (pRed): " << -q << std::endl;
446 outStream <<
" Number of extrapolation steps: " << cnt << std::endl;
452template<
typename Real>
458 std::ostream &outStream) {
459 const Real half(0.5);
460 Real tol = std::sqrt(ROL_EPSILON<Real>());
461 Real beta(1), snorm(0), gs(0);
467 snorm = dgpstep(pwa,s,x,beta,outStream);
468 model.
hessVec(dwa,pwa,x,tol); nhess_++;
471 q = half * pwa.
apply(dwa) + gs;
472 if (q <= mu0_*gs || nsteps > pslim_) {
481 if (verbosity_ > 1) {
482 outStream << std::endl;
483 outStream <<
" Projected search" << std::endl;
484 outStream <<
" Step length (beta): " << beta << std::endl;
485 outStream <<
" Step length (beta*s): " << snorm << std::endl;
486 outStream <<
" Model decrease (pRed): " << -q << std::endl;
487 outStream <<
" Number of steps: " << nsteps << std::endl;
492template<
typename Real>
496 const Real del)
const {
499 Real rad = ptx*ptx + ptp*(dsq-xtx);
500 rad = std::sqrt(std::max(rad,
zero));
503 sigma = (dsq-xtx)/(ptx+rad);
505 else if (rad >
zero) {
506 sigma = (rad-ptx)/ptp;
514template<
typename Real>
519 const Real tol,
const Real stol,
const int itermax,
522 std::ostream &outStream)
const {
527 Real tol0 = std::sqrt(ROL_EPSILON<Real>());
528 const Real
zero(0), one(1), two(2);
529 Real rho(0), kappa(0), beta(0), sigma(0), alpha(0);
530 Real rtr(0), tnorm(0), sMs(0), pMp(0), sMp(0);
537 applyFreePrecond(r,t,x,model,bnd,tol0,dwa,pwa);
542 pMp = (!hasEcon_ ? rho : p.
dot(p));
544 for (iter = 0; iter < itermax; ++iter) {
546 applyFreeHessian(q,p,x,model,bnd,tol0,pwa);
550 alpha = (kappa>
zero) ? rho/kappa :
zero;
551 sigma = dtrqsol(sMs,pMp,sMp,del);
553 if (kappa <= zero || alpha >= sigma) {
556 iflag = (kappa<=
zero) ? 2 : 3;
562 applyFreePrecond(r,t,x,model,bnd,tol0,dwa,pwa);
567 if (rtr <= stol*stol || tnorm <= tol) {
568 sMs = sMs + two*alpha*sMp + alpha*alpha*pMp;
580 sMs = sMs + two*alpha*sMp + alpha*alpha*pMp;
581 sMp = beta*(sMp + alpha*pMp);
582 pMp = (!hasEcon_ ? rho : p.
dot(p)) + beta*beta*pMp;
585 if (iter == itermax) {
591 return std::sqrt(sMs);
594template<
typename Real>
605 model.
hessVec(hv,pwa,x,tol); nhess_++;
615template<
typename Real>
638 rcon_->setX(makePtrFromRef(x));
641 ns_->apply(hv,pwa,tol);
645template<
typename Real>
647 std::stringstream hist;
648 if (verbosity_ > 1) {
649 hist << std::string(114,
'-') << std::endl;
650 hist <<
" Lin-More trust-region method status output definitions" << std::endl << std::endl;
651 hist <<
" iter - Number of iterates (steps taken)" << std::endl;
652 hist <<
" value - Objective function value" << std::endl;
653 hist <<
" gnorm - Norm of the gradient" << std::endl;
654 hist <<
" snorm - Norm of the step (update to optimization vector)" << std::endl;
655 hist <<
" delta - Trust-Region radius" << std::endl;
656 hist <<
" #fval - Number of times the objective function was evaluated" << std::endl;
657 hist <<
" #grad - Number of times the gradient was computed" << std::endl;
658 hist <<
" #hess - Number of times the Hessian was applied" << std::endl;
659 hist <<
" #proj - Number of times the projection was applied" << std::endl;
661 hist <<
" tr_flag - Trust-Region flag" << std::endl;
668 hist <<
" iterCG - Number of Truncated CG iterations" << std::endl << std::endl;
669 hist <<
" flagGC - Trust-Region Truncated CG flag" << std::endl;
675 hist << std::string(114,
'-') << std::endl;
678 hist << std::setw(6) << std::left <<
"iter";
679 hist << std::setw(15) << std::left <<
"value";
680 hist << std::setw(15) << std::left <<
"gnorm";
681 hist << std::setw(15) << std::left <<
"snorm";
682 hist << std::setw(15) << std::left <<
"delta";
683 hist << std::setw(10) << std::left <<
"#fval";
684 hist << std::setw(10) << std::left <<
"#grad";
685 hist << std::setw(10) << std::left <<
"#hess";
686 hist << std::setw(10) << std::left <<
"#proj";
687 hist << std::setw(10) << std::left <<
"tr_flag";
689 hist << std::setw(10) << std::left <<
"iterCG";
690 hist << std::setw(10) << std::left <<
"flagCG";
696template<
typename Real>
698 std::stringstream hist;
699 hist << std::endl <<
"Lin-More Trust-Region Method (Type B, Bound Constraints)" << std::endl;
703template<
typename Real>
705 std::stringstream hist;
706 hist << std::scientific << std::setprecision(6);
707 if ( state_->iter == 0 ) writeName(os);
708 if ( write_header ) writeHeader(os);
709 if ( state_->iter == 0 ) {
711 hist << std::setw(6) << std::left << state_->iter;
712 hist << std::setw(15) << std::left << state_->value;
713 hist << std::setw(15) << std::left << state_->gnorm;
714 hist << std::setw(15) << std::left <<
"---";
715 hist << std::setw(15) << std::left << state_->searchSize;
716 hist << std::setw(10) << std::left << state_->nfval;
717 hist << std::setw(10) << std::left << state_->ngrad;
718 hist << std::setw(10) << std::left << nhess_;
719 hist << std::setw(10) << std::left << state_->nproj;
720 hist << std::setw(10) << std::left <<
"---";
722 hist << std::setw(10) << std::left <<
"---";
723 hist << std::setw(10) << std::left <<
"---";
729 hist << std::setw(6) << std::left << state_->iter;
730 hist << std::setw(15) << std::left << state_->value;
731 hist << std::setw(15) << std::left << state_->gnorm;
732 hist << std::setw(15) << std::left << state_->snorm;
733 hist << std::setw(15) << std::left << state_->searchSize;
734 hist << std::setw(10) << std::left << state_->nfval;
735 hist << std::setw(10) << std::left << state_->ngrad;
736 hist << std::setw(10) << std::left << nhess_;
737 hist << std::setw(10) << std::left << state_->nproj;
738 hist << std::setw(10) << std::left << TRflag_;
740 hist << std::setw(10) << std::left << SPiter_;
741 hist << std::setw(10) << std::left << SPflag_;
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Provides the interface to apply upper and lower bound constraints.
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 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.
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.
void initialize(const Vector< Real > &x, const Vector< Real > &g)
Real optimalityCriterion(const Vector< Real > &x, const Vector< Real > &g, Vector< Real > &primal, std::ostream &outStream=std::cout) const
virtual void writeExitStatus(std::ostream &os) const
void initialize(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &bnd, std::ostream &outStream=std::cout)
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
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)
void writeOutput(std::ostream &os, bool write_header=false) const override
Print iterate status.
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...
Real dtrqsol(const Real xtx, const Real ptp, const Real ptx, const Real del) const
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
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
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)
void writeName(std::ostream &os) const override
Print step name.
Real dgpstep(Vector< Real > &s, const Vector< Real > &w, const Vector< Real > &x, const Real alpha, std::ostream &outStream=std::cout) const
LinMoreAlgorithm(ParameterList &list, const Ptr< Secant< Real > > &secant=nullPtr)
void writeHeader(std::ostream &os) const override
Print iterate header.
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 .
virtual Real dot(const Vector &x) const =0
Compute where .
std::string ETRFlagToString(ETRFlag trf)
std::string NumberToString(T Number)
ESecant StringToESecant(std::string s)
std::string ECGFlagToString(ECGFlag cgf)