44 #ifndef ROL_TYPEB_LSecantBALGORITHM_DEF_HPP
45 #define ROL_TYPEB_LSecantBALGORITHM_DEF_HPP
50 template<
typename Real>
57 ParameterList &trlist = list.sublist(
"Step").sublist(
"Line Search");
58 state_->searchSize =
static_cast<Real
>(1);
60 maxit_ = list.sublist(
"General").sublist(
"Krylov").get(
"Iteration Limit", 20);
61 tol1_ = list.sublist(
"General").sublist(
"Krylov").get(
"Absolute Tolerance", 1e-4);
62 tol2_ = list.sublist(
"General").sublist(
"Krylov").get(
"Relative Tolerance", 1e-2);
64 ROL::ParameterList &lmlist = trlist.sublist(
"Quasi-Newton").sublist(
"L-Secant-B");
65 mu0_ = lmlist.get(
"Sufficient Decrease Parameter", 1e-2);
66 spexp_ = lmlist.get(
"Relative Tolerance Exponent", 1.0);
67 spexp_ = std::max(static_cast<Real>(1),std::min(spexp_,static_cast<Real>(2)));
68 redlim_ = lmlist.sublist(
"Cauchy Point").get(
"Maximum Number of Reduction Steps", 10);
69 explim_ = lmlist.sublist(
"Cauchy Point").get(
"Maximum Number of Expansion Steps", 10);
70 alpha_ = lmlist.sublist(
"Cauchy Point").get(
"Initial Step Size", 1.0);
71 normAlpha_ = lmlist.sublist(
"Cauchy Point").get(
"Normalize Initial Step Size",
false);
72 interpf_ = lmlist.sublist(
"Cauchy Point").get(
"Reduction Rate", 0.1);
73 extrapf_ = lmlist.sublist(
"Cauchy Point").get(
"Expansion Rate", 10.0);
74 qtol_ = lmlist.sublist(
"Cauchy Point").get(
"Decrease Tolerance", 1e-8);
75 interpfPS_ = trlist.sublist(
"Line-Search Method").get(
"Backtracking Rate", 0.5);
77 verbosity_ = list.sublist(
"General").get(
"Output Level",0);
78 writeHeader_ = verbosity_ > 2;
80 useSecantPrecond_ =
true;
81 useSecantHessVec_ =
true;
83 if (secant == nullPtr) {
84 esec_ =
StringToESecant(list.sublist(
"General").sublist(
"Secant").get(
"Type",
"Limited-Memory Secant"));
85 secant_ = SecantFactory<Real>(list,mode);
89 template<
typename Real>
94 std::ostream &outStream) {
97 if (proj_ == nullPtr) {
98 proj_ = makePtr<PolyhedralProjection<Real>>(makePtrFromRef(bnd));
104 Real ftol =
static_cast<Real
>(0.1)*ROL_OVERFLOW<Real>();
105 proj_->project(x,outStream); state_->nproj++;
106 state_->iterateVec->set(x);
108 state_->value = obj.
value(x,ftol); state_->nfval++;
109 obj.
gradient(*state_->gradientVec,x,ftol); state_->ngrad++;
110 state_->stepVec->set(x);
111 state_->stepVec->axpy(-one,state_->gradientVec->dual());
112 proj_->project(*state_->stepVec,outStream); state_->nproj++;
113 state_->stepVec->axpy(-one,x);
114 state_->gnorm = state_->stepVec->norm();
115 state_->snorm = ROL_INF<Real>();
118 alpha_ /= state_->gradientVec->norm();
122 rcon_ = makePtr<ReducedLinearConstraint<Real>>(proj_->getLinearConstraint(),
125 ns_ = makePtr<NullSpaceOperator<Real>>(rcon_,x,
126 *proj_->getResidual());
130 template<
typename Real>
135 std::ostream &outStream ) {
136 const Real
zero(0), one(1);
137 Real tol0 = ROL_EPSILON<Real>(), ftol = ROL_EPSILON<Real>();
138 Real gfnorm(0), gs(0), ftrial(0), q(0), tol(0), stol(0), snorm(0);
140 initialize(x,g,obj,bnd,outStream);
141 Ptr<Vector<Real>> s = x.
clone();
142 Ptr<Vector<Real>> gmod = g.
clone(), gfree = g.
clone();
143 Ptr<Vector<Real>> pwa1 = x.
clone(), pwa2 = x.
clone(), pwa3 = x.
clone();
144 Ptr<Vector<Real>> dwa1 = g.
clone(), dwa2 = g.
clone(), dwa3 = g.
clone();
147 if (verbosity_ > 0) writeOutput(outStream,
true);
149 while (status_->check(*state_)) {
156 state_->gradientVec->
dual(),
165 gmod->plus(*state_->gradientVec);
168 pwa1->set(gfree->dual());
170 gfree->set(pwa1->dual());
172 applyFreePrecond(*pwa1,*gfree,x,*secant_,bnd,tol0,*dwa1,*pwa2);
173 gfnorm = pwa1->norm();
176 gfnorm = gfree->norm();
178 SPiter_ = 0; SPflag_ = 0;
179 if (verbosity_ > 1) {
180 outStream <<
" Norm of free gradient components: " << gfnorm << std::endl;
181 outStream << std::endl;
186 tol = std::min(tol1_,tol2_*std::pow(gfnorm,spexp_));
189 snorm = dpcg(*s, SPflag_, SPiter_, *gfree, x, *secant_, bnd,
191 *pwa1, *dwa1, *pwa2, *dwa2, *pwa3, *dwa3,
193 if (verbosity_ > 1) {
194 outStream <<
" Computation of CG step" << std::endl;
195 outStream <<
" CG step length: " << snorm << std::endl;
196 outStream <<
" Number of CG iterations: " << SPiter_ << std::endl;
197 outStream <<
" CG flag: " << SPflag_ << std::endl;
198 outStream << std::endl;
203 proj_->project(x,outStream); state_->nproj++;
205 state_->stepVec->set(x);
206 state_->stepVec->axpy(-one,*state_->iterateVec);
207 gs = state_->gradientVec->apply(*state_->stepVec);
210 state_->snorm = dsrch(*state_->iterateVec,*state_->stepVec,ftrial,state_->searchSize,
211 state_->value,gs,obj,*pwa1,outStream);
212 x.
set(*state_->iterateVec);
216 state_->value = ftrial;
218 dwa1->set(*state_->gradientVec);
219 obj.
gradient(*state_->gradientVec,x,ftol); state_->ngrad++;
223 secant_->updateStorage(x,*state_->gradientVec,*dwa1,*state_->stepVec,
224 state_->snorm,state_->iter);
227 if (verbosity_ > 0) writeOutput(outStream,writeHeader_);
232 template<
typename Real>
235 std::ostream &outStream)
const {
238 proj_->project(s,outStream); state_->nproj++;
243 template<
typename Real>
251 std::ostream &outStream) {
252 const Real half(0.5);
254 Real gs(0), snorm(0);
256 snorm = dgpstep(s,g,x,-alpha,outStream);
259 q = half * s.
apply(dwa) + gs;
260 interp = (q > mu0_*gs);
267 snorm = dgpstep(s,g,x,-alpha,outStream);
270 q = half * s.
apply(dwa) + gs;
271 search = (q > mu0_*gs) && (cnt < redlim_);
282 snorm = dgpstep(s,g,x,-alpha,outStream);
286 q = half * s.
apply(dwa) + gs;
287 if (q <= mu0_*gs && std::abs(q-qs) > qtol_*std::abs(qs)) {
305 snorm = dgpstep(s,g,x,-alpha,outStream);
307 if (verbosity_ > 1) {
308 outStream <<
" Cauchy point" << std::endl;
309 outStream <<
" Step length (alpha): " << alpha << std::endl;
310 outStream <<
" Step length (alpha*g): " << snorm << std::endl;
311 outStream <<
" Model decrease: " << -q << std::endl;
313 outStream <<
" Number of extrapolation steps: " << cnt << std::endl;
319 template<
typename Real>
324 Real tol = std::sqrt(ROL_EPSILON<Real>());
334 fnew = obj.
value(pwa,tol); state_->nfval++;
335 if (fnew <= fold + mu0_*beta*gs) search =
false;
336 else beta *= interpfPS_;
341 if (verbosity_ > 1) {
342 outStream << std::endl;
343 outStream <<
" Line search" << std::endl;
344 outStream <<
" Step length (beta): " << beta << std::endl;
345 outStream <<
" Step length (beta*s): " << snorm << std::endl;
346 outStream <<
" New objective value: " << fnew << std::endl;
347 outStream <<
" Old objective value: " << fold << std::endl;
348 outStream <<
" Descent verification (gs): " << gs << std::endl;
349 outStream <<
" Number of steps: " << nsteps << std::endl;
354 template<
typename Real>
359 const Real tol,
const Real stol,
const int itermax,
362 std::ostream &outStream)
const {
367 Real tol0 = ROL_EPSILON<Real>(), tolB(0);
368 const Real minIntSize = ROL_EPSILON<Real>();
369 const Real
zero(0), half(0.5), one(1), two(2);
370 Real rho(0), kappa(0), beta(0), sigma(0), alpha(0), alphatmp(0);
371 Real rtr(0), tnorm(0), sMs(0), pMp(0), sMp(0);
378 applyFreePrecond(r,t,x,secant,bnd,tol0,dwa,pwa);
383 pMp = (!hasEcon_ ? rho : p.
dot(p));
385 for (iter = 0; iter < itermax; ++iter) {
387 applyFreeHessian(q,p,x,secant,bnd,tol0,pwa);
395 tolB = minIntSize * alpha;
396 while (alpha-sigma > tolB) {
397 alphatmp = sigma + half * (alpha - sigma);
400 else sigma = alphatmp;
404 sMs = sMs + two*sigma*sMp + sigma*sigma*pMp;
411 applyFreePrecond(r,t,x,secant,bnd,tol0,dwa,pwa);
416 if (rtr <= stol*stol || tnorm <= tol) {
417 sMs = sMs + two*alpha*sMp + alpha*alpha*pMp;
429 sMs = sMs + two*alpha*sMp + alpha*alpha*pMp;
430 sMp = beta*(sMp + alpha*pMp);
431 pMp = (!hasEcon_ ? rho : p.
dot(p)) + beta*beta*pMp;
434 if (iter == itermax) iflag = 1;
435 if (iflag != 1) iter++;
436 return std::sqrt(sMs);
439 template<
typename Real>
456 template<
typename Real>
475 rcon_->setX(makePtrFromRef(x));
478 ns_->apply(hv,pwa,tol);
482 template<
typename Real>
484 std::ios_base::fmtflags osFlags(os.flags());
485 if (verbosity_ > 1) {
486 os << std::string(114,
'-') << std::endl;
487 os <<
" L-Secant-B line search method status output definitions" << std::endl << std::endl;
488 os <<
" iter - Number of iterates (steps taken)" << std::endl;
489 os <<
" value - Objective function value" << std::endl;
490 os <<
" gnorm - Norm of the gradient" << std::endl;
491 os <<
" snorm - Norm of the step (update to optimization vector)" << std::endl;
492 os <<
" LSpar - Line-Search parameter" << std::endl;
493 os <<
" #fval - Number of times the objective function was evaluated" << std::endl;
494 os <<
" #grad - Number of times the gradient was computed" << std::endl;
495 os <<
" #proj - Number of times the projection was applied" << std::endl;
496 os <<
" iterCG - Number of Truncated CG iterations" << std::endl << std::endl;
497 os <<
" flagGC - Trust-Region Truncated CG flag" << std::endl;
498 os <<
" 0 - Converged" << std::endl;
499 os <<
" 1 - Iteration Limit Exceeded" << std::endl;
500 os <<
" 2 - Bounds Exceeded" << std::endl;
501 os << std::string(114,
'-') << std::endl;
504 os << std::setw(6) << std::left <<
"iter";
505 os << std::setw(15) << std::left <<
"value";
506 os << std::setw(15) << std::left <<
"gnorm";
507 os << std::setw(15) << std::left <<
"snorm";
508 os << std::setw(15) << std::left <<
"LSpar";
509 os << std::setw(10) << std::left <<
"#fval";
510 os << std::setw(10) << std::left <<
"#grad";
511 os << std::setw(10) << std::left <<
"#proj";
512 os << std::setw(10) << std::left <<
"iterCG";
513 os << std::setw(10) << std::left <<
"flagCG";
518 template<
typename Real>
520 std::ios_base::fmtflags osFlags(os.flags());
521 os << std::endl <<
"L-Secant-B Line-Search Method (Type B, Bound Constraints)" << std::endl;
525 template<
typename Real>
527 std::ios_base::fmtflags osFlags(os.flags());
528 os << std::scientific << std::setprecision(6);
529 if ( state_->iter == 0 ) writeName(os);
530 if ( write_header ) writeHeader(os);
531 if ( state_->iter == 0 ) {
533 os << std::setw(6) << std::left << state_->iter;
534 os << std::setw(15) << std::left << state_->value;
535 os << std::setw(15) << std::left << state_->gnorm;
536 os << std::setw(15) << std::left <<
"---";
537 os << std::setw(15) << std::left << state_->searchSize;
538 os << std::setw(10) << std::left << state_->nfval;
539 os << std::setw(10) << std::left << state_->ngrad;
540 os << std::setw(10) << std::left << state_->nproj;
541 os << std::setw(10) << std::left <<
"---";
542 os << std::setw(10) << std::left <<
"---";
547 os << std::setw(6) << std::left << state_->iter;
548 os << std::setw(15) << std::left << state_->value;
549 os << std::setw(15) << std::left << state_->gnorm;
550 os << std::setw(15) << std::left << state_->snorm;
551 os << std::setw(15) << std::left << state_->searchSize;
552 os << std::setw(10) << std::left << state_->nfval;
553 os << std::setw(10) << std::left << state_->ngrad;
554 os << std::setw(10) << std::left << state_->nproj;
555 os << std::setw(10) << std::left << SPiter_;
556 os << std::setw(10) << std::left << SPflag_;
virtual bool isFeasible(const Vector< Real > &v)
Check if the vector, v, is feasible.
Provides the interface to evaluate objective functions.
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
virtual void scale(const Real alpha)=0
Compute where .
LSecantBAlgorithm(ParameterList &list, const Ptr< Secant< Real >> &secant=nullPtr)
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 .
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.
virtual void applyH(Vector< Real > &Hv, const Vector< Real > &v) const =0
virtual void writeExitStatus(std::ostream &os) 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()
Real dgpstep(Vector< Real > &s, const Vector< Real > &w, const Vector< Real > &x, const Real alpha, std::ostream &outStream=std::cout) const
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
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
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)
void writeName(std::ostream &os) const override
Print step name.
Provides interface for and implements limited-memory secant operators.
Provides an interface to check status of optimization algorithms.
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...
virtual void applyB(Vector< Real > &Bv, const Vector< Real > &v) const =0
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.
void writeOutput(std::ostream &os, const bool write_header=false) const override
Print iterate status.
Real optimalityCriterion(const Vector< Real > &x, const Vector< Real > &g, Vector< Real > &primal, 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
void initialize(const Vector< Real > &x, const Vector< Real > &g)
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 writeHeader(std::ostream &os) const override
Print iterate header.
virtual void set(const Vector &x)
Set where .
virtual Real norm() const =0
Returns where .
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, Secant< Real > &secant, BoundConstraint< Real > &bnd, Real &tol, Vector< Real > &pwa) const