44#ifndef ROL_TYPEB_LSecantBALGORITHM_DEF_HPP
45#define ROL_TYPEB_LSecantBALGORITHM_DEF_HPP
51template<
typename Real>
58 ParameterList &trlist = list.sublist(
"Step").sublist(
"Line Search");
59 state_->searchSize =
static_cast<Real
>(1);
61 maxit_ = list.sublist(
"General").sublist(
"Krylov").get(
"Iteration Limit", 20);
62 tol1_ = list.sublist(
"General").sublist(
"Krylov").get(
"Absolute Tolerance", 1e-4);
63 tol2_ = list.sublist(
"General").sublist(
"Krylov").get(
"Relative Tolerance", 1e-2);
65 ROL::ParameterList &lmlist = trlist.sublist(
"Quasi-Newton").sublist(
"L-Secant-B");
66 mu0_ = lmlist.get(
"Sufficient Decrease Parameter", 1e-2);
67 spexp_ = lmlist.get(
"Relative Tolerance Exponent", 1.0);
68 spexp_ = std::max(
static_cast<Real
>(1),std::min(spexp_,
static_cast<Real
>(2)));
69 redlim_ = lmlist.sublist(
"Cauchy Point").get(
"Maximum Number of Reduction Steps", 10);
70 explim_ = lmlist.sublist(
"Cauchy Point").get(
"Maximum Number of Expansion Steps", 10);
71 alpha_ = lmlist.sublist(
"Cauchy Point").get(
"Initial Step Size", 1.0);
72 normAlpha_ = lmlist.sublist(
"Cauchy Point").get(
"Normalize Initial Step Size",
false);
73 interpf_ = lmlist.sublist(
"Cauchy Point").get(
"Reduction Rate", 0.1);
74 extrapf_ = lmlist.sublist(
"Cauchy Point").get(
"Expansion Rate", 10.0);
75 qtol_ = lmlist.sublist(
"Cauchy Point").get(
"Decrease Tolerance", 1e-8);
76 interpfPS_ = trlist.sublist(
"Line-Search Method").get(
"Backtracking Rate", 0.5);
78 verbosity_ = list.sublist(
"General").get(
"Output Level",0);
79 writeHeader_ = verbosity_ > 2;
81 useSecantPrecond_ =
true;
82 useSecantHessVec_ =
true;
84 if (secant == nullPtr) {
85 esec_ =
StringToESecant(list.sublist(
"General").sublist(
"Secant").get(
"Type",
"Limited-Memory Secant"));
86 secant_ = SecantFactory<Real>(list,mode);
90template<
typename Real>
95 std::ostream &outStream) {
98 if (proj_ == nullPtr) {
99 proj_ = makePtr<PolyhedralProjection<Real>>(makePtrFromRef(bnd));
105 Real ftol =
static_cast<Real
>(0.1)*ROL_OVERFLOW<Real>();
106 proj_->project(x,outStream); state_->nproj++;
107 state_->iterateVec->set(x);
109 state_->value = obj.
value(x,ftol); state_->nfval++;
110 obj.
gradient(*state_->gradientVec,x,ftol); state_->ngrad++;
111 state_->stepVec->set(x);
112 state_->stepVec->axpy(-one,state_->gradientVec->dual());
113 proj_->project(*state_->stepVec,outStream); state_->nproj++;
114 state_->stepVec->axpy(-one,x);
115 state_->gnorm = state_->stepVec->norm();
116 state_->snorm = ROL_INF<Real>();
119 alpha_ /= state_->gradientVec->norm();
123 rcon_ = makePtr<ReducedLinearConstraint<Real>>(proj_->getLinearConstraint(),
126 ns_ = makePtr<NullSpaceOperator<Real>>(rcon_,x,
127 *proj_->getResidual());
131template<
typename Real>
136 std::ostream &outStream ) {
137 const Real
zero(0), one(1);
138 Real tol0 = ROL_EPSILON<Real>(), ftol = ROL_EPSILON<Real>();
139 Real gfnorm(0), gs(0), ftrial(0), q(0), tol(0), stol(0), snorm(0);
141 initialize(x,g,obj,bnd,outStream);
142 Ptr<Vector<Real>> s = x.
clone();
143 Ptr<Vector<Real>> gmod = g.
clone(), gfree = g.
clone();
144 Ptr<Vector<Real>> pwa1 = x.
clone(), pwa2 = x.
clone(), pwa3 = x.
clone();
145 Ptr<Vector<Real>> dwa1 = g.
clone(), dwa2 = g.
clone(), dwa3 = g.
clone();
148 if (verbosity_ > 0) writeOutput(outStream,
true);
150 while (status_->check(*state_)) {
157 state_->gradientVec->
dual(),
166 gmod->plus(*state_->gradientVec);
169 pwa1->set(gfree->dual());
171 gfree->set(pwa1->dual());
173 applyFreePrecond(*pwa1,*gfree,x,*secant_,bnd,tol0,*dwa1,*pwa2);
174 gfnorm = pwa1->norm();
177 gfnorm = gfree->norm();
179 SPiter_ = 0; SPflag_ = 0;
180 if (verbosity_ > 1) {
181 outStream <<
" Norm of free gradient components: " << gfnorm << std::endl;
182 outStream << std::endl;
187 tol = std::min(tol1_,tol2_*std::pow(gfnorm,spexp_));
190 snorm = dpcg(*s, SPflag_, SPiter_, *gfree, x, *secant_, bnd,
192 *pwa1, *dwa1, *pwa2, *dwa2, *pwa3, *dwa3,
194 if (verbosity_ > 1) {
195 outStream <<
" Computation of CG step" << std::endl;
196 outStream <<
" CG step length: " << snorm << std::endl;
197 outStream <<
" Number of CG iterations: " << SPiter_ << std::endl;
198 outStream <<
" CG flag: " << SPflag_ << std::endl;
199 outStream << std::endl;
204 proj_->project(x,outStream); state_->nproj++;
206 state_->stepVec->set(x);
207 state_->stepVec->axpy(-one,*state_->iterateVec);
208 gs = state_->gradientVec->apply(*state_->stepVec);
211 state_->snorm = dsrch(*state_->iterateVec,*state_->stepVec,ftrial,state_->searchSize,
212 state_->value,gs,obj,*pwa1,outStream);
213 x.
set(*state_->iterateVec);
217 state_->value = ftrial;
219 dwa1->set(*state_->gradientVec);
220 obj.
gradient(*state_->gradientVec,x,ftol); state_->ngrad++;
224 secant_->updateStorage(x,*state_->gradientVec,*dwa1,*state_->stepVec,
225 state_->snorm,state_->iter);
228 if (verbosity_ > 0) writeOutput(outStream,writeHeader_);
233template<
typename Real>
236 std::ostream &outStream)
const {
239 proj_->project(s,outStream); state_->nproj++;
244template<
typename Real>
252 std::ostream &outStream) {
253 const Real half(0.5);
255 Real gs(0), snorm(0);
257 snorm = dgpstep(s,g,x,-alpha,outStream);
260 q = half * s.
apply(dwa) + gs;
261 interp = (q > mu0_*gs);
268 snorm = dgpstep(s,g,x,-alpha,outStream);
271 q = half * s.
apply(dwa) + gs;
272 search = (q > mu0_*gs) && (cnt < redlim_);
283 snorm = dgpstep(s,g,x,-alpha,outStream);
287 q = half * s.
apply(dwa) + gs;
288 if (q <= mu0_*gs && std::abs(q-qs) > qtol_*std::abs(qs)) {
306 snorm = dgpstep(s,g,x,-alpha,outStream);
308 if (verbosity_ > 1) {
309 outStream <<
" Cauchy point" << std::endl;
310 outStream <<
" Step length (alpha): " << alpha << std::endl;
311 outStream <<
" Step length (alpha*g): " << snorm << std::endl;
312 outStream <<
" Model decrease: " << -q << std::endl;
314 outStream <<
" Number of extrapolation steps: " << cnt << std::endl;
320template<
typename Real>
325 Real tol = std::sqrt(ROL_EPSILON<Real>());
335 fnew = obj.
value(pwa,tol); state_->nfval++;
336 if (fnew <= fold + mu0_*beta*gs) search =
false;
337 else beta *= interpfPS_;
342 if (verbosity_ > 1) {
343 outStream << std::endl;
344 outStream <<
" Line search" << std::endl;
345 outStream <<
" Step length (beta): " << beta << std::endl;
346 outStream <<
" Step length (beta*s): " << snorm << std::endl;
347 outStream <<
" New objective value: " << fnew << std::endl;
348 outStream <<
" Old objective value: " << fold << std::endl;
349 outStream <<
" Descent verification (gs): " << gs << std::endl;
350 outStream <<
" Number of steps: " << nsteps << std::endl;
355template<
typename Real>
360 const Real tol,
const Real stol,
const int itermax,
363 std::ostream &outStream)
const {
368 Real tol0 = ROL_EPSILON<Real>(), tolB(0);
369 const Real minIntSize = ROL_EPSILON<Real>();
370 const Real
zero(0), half(0.5), one(1), two(2);
371 Real rho(0), kappa(0), beta(0), sigma(0), alpha(0), alphatmp(0);
372 Real rtr(0), tnorm(0), sMs(0), pMp(0), sMp(0);
379 applyFreePrecond(r,t,x,secant,bnd,tol0,dwa,pwa);
384 pMp = (!hasEcon_ ? rho : p.
dot(p));
386 for (iter = 0; iter < itermax; ++iter) {
388 applyFreeHessian(q,p,x,secant,bnd,tol0,pwa);
396 tolB = minIntSize * alpha;
397 while (alpha-sigma > tolB) {
398 alphatmp = sigma + half * (alpha - sigma);
401 else sigma = alphatmp;
405 sMs = sMs + two*sigma*sMp + sigma*sigma*pMp;
412 applyFreePrecond(r,t,x,secant,bnd,tol0,dwa,pwa);
417 if (rtr <= stol*stol || tnorm <= tol) {
418 sMs = sMs + two*alpha*sMp + alpha*alpha*pMp;
430 sMs = sMs + two*alpha*sMp + alpha*alpha*pMp;
431 sMp = beta*(sMp + alpha*pMp);
432 pMp = (!hasEcon_ ? rho : p.
dot(p)) + beta*beta*pMp;
435 if (iter == itermax) iflag = 1;
436 if (iflag != 1) iter++;
437 return std::sqrt(sMs);
440template<
typename Real>
457template<
typename Real>
476 rcon_->setX(makePtrFromRef(x));
479 ns_->apply(hv,pwa,tol);
483template<
typename Real>
485 std::stringstream hist;
486 if (verbosity_ > 1) {
487 hist << std::string(114,
'-') << std::endl;
488 hist <<
" L-Secant-B line search method status output definitions" << std::endl << std::endl;
489 hist <<
" iter - Number of iterates (steps taken)" << std::endl;
490 hist <<
" value - Objective function value" << std::endl;
491 hist <<
" gnorm - Norm of the gradient" << std::endl;
492 hist <<
" snorm - Norm of the step (update to optimization vector)" << std::endl;
493 hist <<
" LSpar - Line-Search parameter" << std::endl;
494 hist <<
" #fval - Number of times the objective function was evaluated" << std::endl;
495 hist <<
" #grad - Number of times the gradient was computed" << std::endl;
496 hist <<
" #proj - Number of times the projection was applied" << std::endl;
497 hist <<
" iterCG - Number of Truncated CG iterations" << std::endl << std::endl;
498 hist <<
" flagGC - Trust-Region Truncated CG flag" << std::endl;
499 hist <<
" 0 - Converged" << std::endl;
500 hist <<
" 1 - Iteration Limit Exceeded" << std::endl;
501 hist <<
" 2 - Bounds Exceeded" << std::endl;
502 hist << std::string(114,
'-') << std::endl;
505 hist << std::setw(6) << std::left <<
"iter";
506 hist << std::setw(15) << std::left <<
"value";
507 hist << std::setw(15) << std::left <<
"gnorm";
508 hist << std::setw(15) << std::left <<
"snorm";
509 hist << std::setw(15) << std::left <<
"LSpar";
510 hist << std::setw(10) << std::left <<
"#fval";
511 hist << std::setw(10) << std::left <<
"#grad";
512 hist << std::setw(10) << std::left <<
"#proj";
513 hist << std::setw(10) << std::left <<
"iterCG";
514 hist << std::setw(10) << std::left <<
"flagCG";
519template<
typename Real>
521 std::stringstream hist;
522 hist << std::endl <<
"L-Secant-B Line-Search Method (Type B, Bound Constraints)" << std::endl;
526template<
typename Real>
528 std::stringstream hist;
529 hist << std::scientific << std::setprecision(6);
530 if ( state_->iter == 0 ) writeName(os);
531 if ( write_header ) writeHeader(os);
532 if ( state_->iter == 0 ) {
534 hist << std::setw(6) << std::left << state_->iter;
535 hist << std::setw(15) << std::left << state_->value;
536 hist << std::setw(15) << std::left << state_->gnorm;
537 hist << std::setw(15) << std::left <<
"---";
538 hist << std::setw(15) << std::left << state_->searchSize;
539 hist << std::setw(10) << std::left << state_->nfval;
540 hist << std::setw(10) << std::left << state_->ngrad;
541 hist << std::setw(10) << std::left << state_->nproj;
542 hist << std::setw(10) << std::left <<
"---";
543 hist << std::setw(10) << std::left <<
"---";
548 hist << std::setw(6) << std::left << state_->iter;
549 hist << std::setw(15) << std::left << state_->value;
550 hist << std::setw(15) << std::left << state_->gnorm;
551 hist << std::setw(15) << std::left << state_->snorm;
552 hist << std::setw(15) << std::left << state_->searchSize;
553 hist << std::setw(10) << std::left << state_->nfval;
554 hist << std::setw(10) << std::left << state_->ngrad;
555 hist << std::setw(10) << std::left << state_->nproj;
556 hist << std::setw(10) << std::left << SPiter_;
557 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.
virtual bool isFeasible(const Vector< Real > &v)
Check if the vector, v, is feasible.
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.
virtual void applyB(Vector< Real > &Bv, const Vector< Real > &v) const =0
virtual void applyH(Vector< Real > &Hv, const Vector< Real > &v) const =0
Provides an interface to check status of optimization algorithms.
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 applyFreeHessian(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Secant< Real > &secant, BoundConstraint< Real > &bnd, Real &tol, Vector< Real > &pwa) const
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 dcauchy(Vector< Real > &s, Real &alpha, Real &q, const Vector< Real > &x, const Vector< Real > &g, Secant< Real > &secant, Vector< Real > &dwa, Vector< Real > &dwa1, std::ostream &outStream=std::cout)
Real dpcg(Vector< Real > &w, int &iflag, int &iter, const Vector< Real > &g, const Vector< Real > &x, Secant< Real > &secant, 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, Secant< Real > &secant, BoundConstraint< Real > &bnd, Real &tol, Vector< Real > &dwa, Vector< Real > &pwa) const
Real dgpstep(Vector< Real > &s, const Vector< Real > &w, const Vector< Real > &x, const Real alpha, std::ostream &outStream=std::cout) const
void writeHeader(std::ostream &os) const override
Print iterate header.
void initialize(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &bnd, std::ostream &outStream=std::cout)
LSecantBAlgorithm(ParameterList &list, const Ptr< Secant< Real > > &secant=nullPtr)
Real dsrch(Vector< Real > &x, Vector< Real > &s, Real &fnew, Real &beta, Real fold, Real gs, Objective< Real > &obj, Vector< Real > &pwa, std::ostream &outStream=std::cout)
void writeName(std::ostream &os) const override
Print step name.
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 .
ESecant StringToESecant(std::string s)