44#ifndef ROL_TYPEG_STABILIZEDLCLALGORITHM_DEF_H
45#define ROL_TYPEG_STABILIZEDLCLALGORITHM_DEF_H
53template<
typename Real>
60 Real one(1), p1(0.1), p9(0.9), ten(1.e1), oe8(1.e8), oem8(1.e-8);
61 ParameterList& sublist = list.sublist(
"Step").sublist(
"Stabilized LCL");
63 state_->searchSize = sublist.get(
"Initial Penalty Parameter", ten);
64 sigma_ = sublist.get(
"Initial Elastic Penalty Parameter", ten*ten);
67 penaltyUpdate_ = sublist.get(
"Penalty Parameter Growth Factor", ten);
69 sigmaMax_ = sublist.get(
"Maximum Elastic Penalty Parameter", oe8);
70 sigmaUpdate_ = sublist.get(
"Elastic Penalty Parameter Growth Rate", ten);
80 maxit_ = sublist.get(
"Subproblem Iteration Limit", 1000);
81 subStep_ = sublist.get(
"Subproblem Step Type",
"Trust Region");
84 list_.sublist(
"Status Test").set(
"Iteration Limit",
maxit_);
86 verbosity_ = list.sublist(
"General").get(
"Output Level", 0);
96 fscale_ = sublist.get(
"Objective Scaling", one);
97 cscale_ = sublist.get(
"Constraint Scaling", one);
100template<
typename Real>
108 std::ostream &outStream ) {
110 if (proj_ == nullPtr) {
111 proj_ = makePtr<PolyhedralProjection<Real>>(makePtrFromRef(bnd));
112 hasPolyProj_ =
false;
114 proj_->project(x,outStream);
116 const Real one(1), TOL(1.e-2);
117 Real tol = std::sqrt(ROL_EPSILON<Real>());
132 state_->cnorm = state_->constraintVec->norm();
140 if (useDefaultScaling_) {
143 Ptr<Vector<Real>> ji = x.
clone();
144 Real maxji(0), normji(0);
145 for (
int i = 0; i < c.
dimension(); ++i) {
148 maxji = std::max(normji,maxji);
150 cscale_ = one/std::max(one,maxji);
152 catch (std::exception &e) {
159 x.
axpy(-one,state_->gradientVec->dual());
160 proj_->project(x,outStream);
161 x.
axpy(-one/std::min(fscale_,cscale_),*state_->iterateVec);
162 state_->gnorm = x.
norm();
163 x.
set(*state_->iterateVec);
166 if (useDefaultInitPen_) {
167 const Real oem8(1e-8), oem2(1e-2), two(2), ten(10);
168 state_->searchSize = std::max(oem8,
169 std::min(ten*std::max(one,std::abs(fscale_*state_->value))
170 / std::max(one,std::pow(cscale_*state_->cnorm,two)),oem2*maxPenaltyParam_));
173 optTolerance_ = std::max<Real>(TOL*outerOptTolerance_,
174 optToleranceInitial_);
176 feasTolerance_ = std::max<Real>(TOL*outerFeasTolerance_,
177 feasToleranceInitial_);
180 alobj.
reset(l,state_->searchSize,sigma_);
182 if (verbosity_ > 1) {
183 outStream << std::endl;
184 outStream <<
"Stabilized LCL Initialize" << std::endl;
185 outStream <<
"Objective Scaling: " << fscale_ << std::endl;
186 outStream <<
"Constraint Scaling: " << cscale_ << std::endl;
187 outStream <<
"Penalty Parameter: " << state_->searchSize << std::endl;
188 outStream << std::endl;
192template<
typename Real>
194 std::ostream &outStream ) {
197 problem.
finalize(
true,verbosity_>3,outStream);
213template<
typename Real>
221 std::ostream &outStream ) {
222 const Real one(1), oem2(1e-2);
223 Real tol(std::sqrt(ROL_EPSILON<Real>())), cnorm(0), lnorm;;
226 state_->searchSize,sigma_,g,eres,emul,
227 scaleLagrangian_,HessianApprox_);
228 initialize(x,g,emul,eres,alobj,bnd,econ,outStream);
230 Ptr<Vector<Real>> u = eres.
clone(), v = eres.
clone(), c = eres.
clone();
231 Ptr<Vector<Real>> gu = emul.
clone(), gv = emul.
clone(), l = emul.
clone();
233 Ptr<ElasticLinearConstraint<Real>> lcon
234 = makePtr<ElasticLinearConstraint<Real>>(makePtrFromRef(x),
235 makePtrFromRef(econ),
236 makePtrFromRef(eres));
237 std::vector<Ptr<Vector<Real>>> vecList = {s,u,v};
238 Ptr<PartitionedVector<Real>> xp = makePtr<PartitionedVector<Real>>(vecList);
239 Ptr<PartitionedVector<Real>> gxp = makePtr<PartitionedVector<Real>>({gs,gu,gv});
240 Ptr<Vector<Real>> lb = u->clone(); lb->zero();
241 std::vector<Ptr<BoundConstraint<Real>>> bndList(3);
242 bndList[0] = makePtrFromRef(bnd);
243 bndList[1] = makePtr<Bounds<Real>>(*lb,
true);
244 bndList[2] = makePtr<Bounds<Real>>(*lb,
true);
245 Ptr<BoundConstraint<Real>> xbnd
246 = makePtr<BoundConstraint_Partitioned<Real>>(bndList,vecList);
247 ParameterList ppa_list;
248 if (c->dimension() == 1)
249 ppa_list.sublist(
"General").sublist(
"Polyhedral Projection").set(
"Type",
"Dai-Fletcher");
251 ppa_list.sublist(
"General").sublist(
"Polyhedral Projection").set(
"Type",
"Semismooth Newton");
256 elc.
finalize(
false,verbosity_>2,outStream);
259 Ptr<TypeB::Algorithm<Real>> algo;
262 if (verbosity_ > 0) writeOutput(outStream,
true);
264 while (status_->check(*state_)) {
265 lcon->setAnchor(state_->iterateVec);
266 if (verbosity_ > 3) elc.
check(
true,outStream);
269 list_.sublist(
"Status Test").set(
"Gradient Tolerance",optTolerance_);
270 list_.sublist(
"Status Test").set(
"Step Tolerance",1.e-6*optTolerance_);
271 algo = TypeB::AlgorithmFactory<Real>(list_);
272 algo->run(elc,outStream);
276 subproblemIter_ = algo->getState()->iter;
282 state_->stepVec->set(x);
283 state_->stepVec->axpy(-one,*state_->iterateVec);
284 state_->snorm = state_->stepVec->norm();
289 cnorm = cvec->norm();
290 if ( cscale_*cnorm < feasTolerance_ ) {
292 state_->iterateVec->set(x);
294 state_->constraintVec->set(*cvec);
295 state_->cnorm = cnorm;
299 emul.
axpy(state_->searchSize*cscale_,state_->constraintVec->dual());
302 if (scaleLagrangian_) state_->gradientVec->scale(state_->searchSize);
304 state_->gradientVec->axpy(-cscale_,*gs);
305 x.
axpy(-one/std::min(fscale_,cscale_),state_->gradientVec->dual());
306 proj_->project(x,outStream);
307 x.
axpy(-one,*state_->iterateVec);
308 state_->gnorm = x.
norm();
309 x.
set(*state_->iterateVec);
313 sigma_ = std::min(one+lnorm,sigmaMax_)/(one+state_->searchSize);
315 optTolerance_ = std::max(oem2*outerOptTolerance_,
316 optTolerance_/(one + std::pow(state_->searchSize,optIncreaseExponent_)));
318 feasTolerance_ = std::max(oem2*outerFeasTolerance_,
319 feasTolerance_/(one + std::pow(state_->searchSize,feasIncreaseExponent_)));
322 state_->snorm += lnorm + state_->searchSize*cscale_*state_->cnorm;
323 state_->lagmultVec->set(emul);
327 state_->searchSize = std::min(penaltyUpdate_*state_->searchSize,maxPenaltyParam_);
328 sigma_ /= sigmaUpdate_;
329 optTolerance_ = std::max(oem2*outerOptTolerance_,
330 optToleranceInitial_/(one + std::pow(state_->searchSize,optDecreaseExponent_)));
331 feasTolerance_ = std::max(oem2*outerFeasTolerance_,
332 feasToleranceInitial_/(one + std::pow(state_->searchSize,feasDecreaseExponent_)));
334 alobj.
reset(emul,state_->searchSize,sigma_);
337 if (verbosity_ > 0) writeOutput(outStream,printHeader_);
342template<
typename Real>
344 std::stringstream hist;
346 hist << std::string(114,
'-') << std::endl;
347 hist <<
"Stabilized LCL status output definitions" << std::endl << std::endl;
348 hist <<
" iter - Number of iterates (steps taken)" << std::endl;
349 hist <<
" fval - Objective function value" << std::endl;
350 hist <<
" cnorm - Norm of the constraint violation" << std::endl;
351 hist <<
" gLnorm - Norm of the gradient of the Lagrangian" << std::endl;
352 hist <<
" snorm - Norm of the step" << std::endl;
353 hist <<
" penalty - Penalty parameter" << std::endl;
354 hist <<
" sigma - Elastic Penalty parameter" << std::endl;
355 hist <<
" feasTol - Feasibility tolerance" << std::endl;
356 hist <<
" optTol - Optimality tolerance" << std::endl;
357 hist <<
" #fval - Number of times the objective was computed" << std::endl;
358 hist <<
" #grad - Number of times the gradient was computed" << std::endl;
359 hist <<
" #cval - Number of times the constraint was computed" << std::endl;
360 hist <<
" subIter - Number of iterations to solve subproblem" << std::endl;
361 hist << std::string(114,
'-') << std::endl;
364 hist << std::setw(6) << std::left <<
"iter";
365 hist << std::setw(15) << std::left <<
"fval";
366 hist << std::setw(15) << std::left <<
"cnorm";
367 hist << std::setw(15) << std::left <<
"gLnorm";
368 hist << std::setw(15) << std::left <<
"snorm";
369 hist << std::setw(10) << std::left <<
"penalty";
370 hist << std::setw(10) << std::left <<
"sigma";
371 hist << std::setw(10) << std::left <<
"feasTol";
372 hist << std::setw(10) << std::left <<
"optTol";
373 hist << std::setw(8) << std::left <<
"#fval";
374 hist << std::setw(8) << std::left <<
"#grad";
375 hist << std::setw(8) << std::left <<
"#cval";
376 hist << std::setw(8) << std::left <<
"subIter";
381template<
typename Real>
383 std::stringstream hist;
384 hist << std::endl <<
"Stabilized LCL Solver (Type G, General Constraints)";
386 hist <<
"Subproblem Solver: " << subStep_ << std::endl;
390template<
typename Real>
392 std::stringstream hist;
393 hist << std::scientific << std::setprecision(6);
394 if ( state_->iter == 0 ) writeName(os);
395 if ( print_header ) writeHeader(os);
396 if ( state_->iter == 0 ) {
398 hist << std::setw(6) << std::left << state_->iter;
399 hist << std::setw(15) << std::left << state_->value;
400 hist << std::setw(15) << std::left << state_->cnorm;
401 hist << std::setw(15) << std::left << state_->gnorm;
402 hist << std::setw(15) << std::left <<
"---";
403 hist << std::scientific << std::setprecision(2);
404 hist << std::setw(10) << std::left << state_->searchSize;
405 hist << std::setw(10) << std::left << sigma_;
406 hist << std::setw(10) << std::left << std::max(feasTolerance_,outerFeasTolerance_);
407 hist << std::setw(10) << std::left << std::max(optTolerance_,outerOptTolerance_);
408 hist << std::scientific << std::setprecision(6);
409 hist << std::setw(8) << std::left << state_->nfval;
410 hist << std::setw(8) << std::left << state_->ngrad;
411 hist << std::setw(8) << std::left << state_->ncval;
412 hist << std::setw(8) << std::left <<
"---";
417 hist << std::setw(6) << std::left << state_->iter;
418 hist << std::setw(15) << std::left << state_->value;
419 hist << std::setw(15) << std::left << state_->cnorm;
420 hist << std::setw(15) << std::left << state_->gnorm;
421 hist << std::setw(15) << std::left << state_->snorm;
422 hist << std::scientific << std::setprecision(2);
423 hist << std::setw(10) << std::left << state_->searchSize;
424 hist << std::setw(10) << std::left << sigma_;
425 hist << std::setw(10) << std::left << feasTolerance_;
426 hist << std::setw(10) << std::left << optTolerance_;
427 hist << std::scientific << std::setprecision(6);
428 hist << std::setw(8) << std::left << state_->nfval;
429 hist << std::setw(8) << std::left << state_->ngrad;
430 hist << std::setw(8) << std::left << state_->ncval;
431 hist << std::setw(8) << std::left << subproblemIter_;
Provides the interface to apply upper and lower bound constraints.
Provides an interface to check status of optimization algorithms for problems with equality constrain...
Defines the general constraint operator interface.
virtual void applyAdjointJacobian(Vector< Real > &ajv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply the adjoint of the the constraint Jacobian at , , to vector .
Provides the interface to evaluate the elastic augmented Lagrangian.
void reset(const Vector< Real > &multiplier, Real penaltyParameter, Real sigma)
const Ptr< const Vector< Real > > getConstraintVec(const Vector< Real > &x, Real &tol)
const Ptr< const Vector< Real > > getObjectiveGradient(const Vector< Real > &x, Real &tol)
const Ptr< AugmentedLagrangianObjective< Real > > getAugmentedLagrangian(void) const
int getNumberConstraintEvaluations(void) const
Real getObjectiveValue(const Vector< Real > &x, Real &tol)
void setScaling(const Real fscale=1.0, const Real cscale=1.0)
int getNumberFunctionEvaluations(void) const
int getNumberGradientEvaluations(void) const
Provides the interface to evaluate objective functions.
const Ptr< PolyhedralProjection< Real > > & getPolyhedralProjection()
Get the polyhedral projection object. This is a null pointer if no linear constraints and/or bounds a...
const Ptr< Vector< Real > > & getPrimalOptimizationVector()
Get the primal optimization space vector.
const Ptr< Vector< Real > > & getDualOptimizationVector()
Get the dual optimization space vector.
const Ptr< Vector< Real > > & getMultiplierVector()
Get the dual constraint space vector.
const Ptr< Constraint< Real > > & getConstraint()
Get the equality constraint.
EProblem getProblemType()
Get the optimization problem type (U, B, E, or G).
void addLinearConstraint(std::string name, const Ptr< Constraint< Real > > &linear_econ, const Ptr< Vector< Real > > &linear_emul, const Ptr< Vector< Real > > &linear_eres=nullPtr, bool reset=false)
Add a linear equality constraint.
void addBoundConstraint(const Ptr< BoundConstraint< Real > > &bnd)
Add a bound constraint.
void setProjectionAlgorithm(ParameterList &parlist)
Set polyhedral projection algorithm.
void finalizeIteration()
Transform the optimization variables to the native parameterization after an optimization algorithm h...
void check(bool printToStream=false, std::ostream &outStream=std::cout) const
Run vector, linearity and derivative checks for user-supplied vectors, objective function and constra...
const Ptr< Objective< Real > > & getObjective()
Get the objective function.
const Ptr< BoundConstraint< Real > > & getBoundConstraint()
Get the bound constraint.
const Ptr< Vector< Real > > & getResidualVector()
Get the primal constraint space vector.
virtual void finalize(bool lumpConstraints=false, bool printToStream=false, std::ostream &outStream=std::cout)
Tranform user-supplied constraints to consist of only bounds and equalities. Optimization problem can...
virtual void edit()
Resume editting optimization problem after finalize has been called.
Provides an interface to run general constrained optimization algorithms.
void initialize(const Vector< Real > &x, const Vector< Real > &g, const Vector< Real > &mul, const Vector< Real > &c)
virtual void writeExitStatus(std::ostream &os) const
const Ptr< CombinedStatusTest< Real > > status_
const Ptr< AlgorithmState< Real > > state_
Real optDecreaseExponent_
virtual void run(Problem< Real > &problem, std::ostream &outStream=std::cout) override
Run algorithm on general constrained problems (Type-G). This is the primary Type-G interface.
void initialize(Vector< Real > &x, const Vector< Real > &g, const Vector< Real > &l, const Vector< Real > &c, ElasticObjective< Real > &alobj, BoundConstraint< Real > &bnd, Constraint< Real > &con, std::ostream &outStream=std::cout)
Real optIncreaseExponent_
Real feasIncreaseExponent_
virtual void writeName(std::ostream &os) const override
Print step name.
StabilizedLCLAlgorithm(ParameterList &list)
Real feasDecreaseExponent_
Real optToleranceInitial_
virtual void writeHeader(std::ostream &os) const override
Print iterate header.
virtual void writeOutput(std::ostream &os, const bool print_header=false) const override
Print iterate status.
Real feasToleranceInitial_
Defines the linear algebra or vector space interface.
virtual Real norm() const =0
Returns where .
virtual void set(const Vector &x)
Set where .
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual int dimension() const
Return dimension of the vector space.
virtual ROL::Ptr< Vector > basis(const int i) const
Return i-th basis vector.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .