44#ifndef ROL_TYPEE_STABILIZEDLCLALGORITHM_DEF_H
45#define ROL_TYPEE_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>
107 std::ostream &outStream ) {
108 const Real one(1), TOL(1.e-2);
109 Real tol = std::sqrt(ROL_EPSILON<Real>());
124 state_->cnorm = state_->constraintVec->norm();
132 if (useDefaultScaling_) {
135 Ptr<Vector<Real>> ji = x.
clone();
136 Real maxji(0), normji(0);
137 for (
int i = 0; i < c.
dimension(); ++i) {
140 maxji = std::max(normji,maxji);
142 cscale_ = one/std::max(one,maxji);
144 catch (std::exception &e) {
151 state_->gnorm = state_->gradientVec->norm()/std::min(fscale_,cscale_);
154 if (useDefaultInitPen_) {
155 const Real oem8(1e-8), oem2(1e-2), two(2), ten(10);
156 state_->searchSize = std::max(oem8,
157 std::min(ten*std::max(one,std::abs(fscale_*state_->value))
158 / std::max(one,std::pow(cscale_*state_->cnorm,two)),oem2*maxPenaltyParam_));
161 optTolerance_ = std::max<Real>(TOL*outerOptTolerance_,
162 optToleranceInitial_);
164 feasTolerance_ = std::max<Real>(TOL*outerFeasTolerance_,
165 feasToleranceInitial_);
168 alobj.
reset(l,state_->searchSize,sigma_);
170 if (verbosity_ > 1) {
171 outStream << std::endl;
172 outStream <<
"Stabilized LCL Initialize" << std::endl;
173 outStream <<
"Objective Scaling: " << fscale_ << std::endl;
174 outStream <<
"Constraint Scaling: " << cscale_ << std::endl;
175 outStream <<
"Penalty Parameter: " << state_->searchSize << std::endl;
176 outStream << std::endl;
180template<
typename Real>
182 std::ostream &outStream ) {
185 problem.
finalize(
true,verbosity_>3,outStream);
200template<
typename Real>
207 std::ostream &outStream ) {
208 const Real one(1), oem2(1e-2);
209 Real tol(std::sqrt(ROL_EPSILON<Real>())), cnorm(0), lnorm;;
212 state_->searchSize,sigma_,g,eres,emul,
213 scaleLagrangian_,HessianApprox_);
214 initialize(x,g,emul,eres,alobj,econ,outStream);
216 Ptr<Vector<Real>> u = eres.
clone(), v = eres.
clone(), c = eres.
clone();
217 Ptr<Vector<Real>> gu = emul.
clone(), gv = emul.
clone(), l = emul.
clone();
219 Ptr<ElasticLinearConstraint<Real>> lcon
220 = makePtr<ElasticLinearConstraint<Real>>(makePtrFromRef(x),
221 makePtrFromRef(econ),
222 makePtrFromRef(eres));
223 std::vector<Ptr<Vector<Real>>> vecList = {s,u,v};
224 Ptr<PartitionedVector<Real>> xp = makePtr<PartitionedVector<Real>>(vecList);
225 Ptr<PartitionedVector<Real>> gxp = makePtr<PartitionedVector<Real>>({gs,gu,gv});
226 Ptr<Vector<Real>> lb = u->clone(); lb->zero();
227 std::vector<Ptr<BoundConstraint<Real>>> bndList(3);
228 bndList[0] = makePtr<BoundConstraint<Real>>(); bndList[0]->deactivate();
229 bndList[1] = makePtr<Bounds<Real>>(*lb,
true);
230 bndList[2] = makePtr<Bounds<Real>>(*lb,
true);
231 Ptr<BoundConstraint<Real>> xbnd
232 = makePtr<BoundConstraint_Partitioned<Real>>(bndList,vecList);
233 ParameterList ppa_list;
234 if (c->dimension() == 1)
235 ppa_list.sublist(
"General").sublist(
"Polyhedral Projection").set(
"Type",
"Dai-Fletcher");
237 ppa_list.sublist(
"General").sublist(
"Polyhedral Projection").set(
"Type",
"Semismooth Newton");
241 elc.
finalize(
false,verbosity_>2,outStream);
242 Ptr<Vector<Real>> b2 = eres.
clone(), xpwa = xp->clone(), mul = emul.
clone();
246 Ptr<TypeB::Algorithm<Real>> algo;
249 if (verbosity_ > 0) writeOutput(outStream,
true);
251 while (status_->check(*state_)) {
252 lcon->setAnchor(state_->iterateVec);
253 if (verbosity_ > 3) elc.
check(
true,outStream);
256 list_.sublist(
"Status Test").set(
"Gradient Tolerance",optTolerance_);
257 list_.sublist(
"Status Test").set(
"Step Tolerance",1.e-6*optTolerance_);
258 algo = TypeB::AlgorithmFactory<Real>(list_);
259 algo->run(elc,outStream);
263 subproblemIter_ = algo->getState()->iter;
269 state_->stepVec->set(x);
270 state_->stepVec->axpy(-one,*state_->iterateVec);
271 state_->snorm = state_->stepVec->norm();
276 cnorm = cvec->norm();
277 if ( cscale_*cnorm < feasTolerance_ ) {
279 state_->iterateVec->set(x);
281 state_->constraintVec->set(*cvec);
282 state_->cnorm = cnorm;
286 emul.
axpy(state_->searchSize*cscale_,state_->constraintVec->dual());
289 if (scaleLagrangian_) state_->gradientVec->scale(state_->searchSize);
291 state_->gradientVec->axpy(-cscale_,*gs);
292 state_->gnorm = state_->gradientVec->norm();
296 sigma_ = std::min(one+lnorm,sigmaMax_)/(one+state_->searchSize);
298 optTolerance_ = std::max(oem2*outerOptTolerance_,
299 optTolerance_/(one + std::pow(state_->searchSize,optIncreaseExponent_)));
301 feasTolerance_ = std::max(oem2*outerFeasTolerance_,
302 feasTolerance_/(one + std::pow(state_->searchSize,feasIncreaseExponent_)));
305 state_->snorm += lnorm + state_->searchSize*cscale_*state_->cnorm;
306 state_->lagmultVec->set(emul);
310 state_->searchSize = std::min(penaltyUpdate_*state_->searchSize,maxPenaltyParam_);
311 sigma_ /= sigmaUpdate_;
312 optTolerance_ = std::max(oem2*outerOptTolerance_,
313 optToleranceInitial_/(one + std::pow(state_->searchSize,optDecreaseExponent_)));
314 feasTolerance_ = std::max(oem2*outerFeasTolerance_,
315 feasToleranceInitial_/(one + std::pow(state_->searchSize,feasDecreaseExponent_)));
317 alobj.
reset(emul,state_->searchSize,sigma_);
320 if (verbosity_ > 0) writeOutput(outStream,printHeader_);
325template<
typename Real>
327 std::stringstream hist;
329 hist << std::string(114,
'-') << std::endl;
330 hist <<
"Stabilized LCL status output definitions" << std::endl << std::endl;
331 hist <<
" iter - Number of iterates (steps taken)" << std::endl;
332 hist <<
" fval - Objective function value" << std::endl;
333 hist <<
" cnorm - Norm of the constraint violation" << std::endl;
334 hist <<
" gLnorm - Norm of the gradient of the Lagrangian" << std::endl;
335 hist <<
" snorm - Norm of the step" << std::endl;
336 hist <<
" penalty - Penalty parameter" << std::endl;
337 hist <<
" sigma - Elastic Penalty parameter" << std::endl;
338 hist <<
" feasTol - Feasibility tolerance" << std::endl;
339 hist <<
" optTol - Optimality tolerance" << std::endl;
340 hist <<
" #fval - Number of times the objective was computed" << std::endl;
341 hist <<
" #grad - Number of times the gradient was computed" << std::endl;
342 hist <<
" #cval - Number of times the constraint was computed" << std::endl;
343 hist <<
" subIter - Number of iterations to solve subproblem" << std::endl;
344 hist << std::string(114,
'-') << std::endl;
347 hist << std::setw(6) << std::left <<
"iter";
348 hist << std::setw(15) << std::left <<
"fval";
349 hist << std::setw(15) << std::left <<
"cnorm";
350 hist << std::setw(15) << std::left <<
"gLnorm";
351 hist << std::setw(15) << std::left <<
"snorm";
352 hist << std::setw(10) << std::left <<
"penalty";
353 hist << std::setw(10) << std::left <<
"sigma";
354 hist << std::setw(10) << std::left <<
"feasTol";
355 hist << std::setw(10) << std::left <<
"optTol";
356 hist << std::setw(8) << std::left <<
"#fval";
357 hist << std::setw(8) << std::left <<
"#grad";
358 hist << std::setw(8) << std::left <<
"#cval";
359 hist << std::setw(8) << std::left <<
"subIter";
364template<
typename Real>
366 std::stringstream hist;
367 hist << std::endl <<
"Stabilized LCL Solver (Type E, Equality Constraints)";
369 hist <<
"Subproblem Solver: " << subStep_ << std::endl;
373template<
typename Real>
375 std::stringstream hist;
376 hist << std::scientific << std::setprecision(6);
377 if ( state_->iter == 0 ) writeName(os);
378 if ( print_header ) writeHeader(os);
379 if ( state_->iter == 0 ) {
381 hist << std::setw(6) << std::left << state_->iter;
382 hist << std::setw(15) << std::left << state_->value;
383 hist << std::setw(15) << std::left << state_->cnorm;
384 hist << std::setw(15) << std::left << state_->gnorm;
385 hist << std::setw(15) << std::left <<
"---";
386 hist << std::scientific << std::setprecision(2);
387 hist << std::setw(10) << std::left << state_->searchSize;
388 hist << std::setw(10) << std::left << sigma_;
389 hist << std::setw(10) << std::left << std::max(feasTolerance_,outerFeasTolerance_);
390 hist << std::setw(10) << std::left << std::max(optTolerance_,outerOptTolerance_);
391 hist << std::scientific << std::setprecision(6);
392 hist << std::setw(8) << std::left << state_->nfval;
393 hist << std::setw(8) << std::left << state_->ngrad;
394 hist << std::setw(8) << std::left << state_->ncval;
395 hist << std::setw(8) << std::left <<
"---";
400 hist << std::setw(6) << std::left << state_->iter;
401 hist << std::setw(15) << std::left << state_->value;
402 hist << std::setw(15) << std::left << state_->cnorm;
403 hist << std::setw(15) << std::left << state_->gnorm;
404 hist << std::setw(15) << std::left << state_->snorm;
405 hist << std::scientific << std::setprecision(2);
406 hist << std::setw(10) << std::left << state_->searchSize;
407 hist << std::setw(10) << std::left << sigma_;
408 hist << std::setw(10) << std::left << feasTolerance_;
409 hist << std::setw(10) << std::left << optTolerance_;
410 hist << std::scientific << std::setprecision(6);
411 hist << std::setw(8) << std::left << state_->nfval;
412 hist << std::setw(8) << std::left << state_->ngrad;
413 hist << std::setw(8) << std::left << state_->ncval;
414 hist << std::setw(8) << std::left << subproblemIter_;
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 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< 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.
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< AlgorithmState< Real > > state_
const Ptr< CombinedStatusTest< Real > > status_
StabilizedLCLAlgorithm(ParameterList &list)
Real optToleranceInitial_
void initialize(Vector< Real > &x, const Vector< Real > &g, const Vector< Real > &l, const Vector< Real > &c, ElasticObjective< Real > &alobj, Constraint< Real > &con, std::ostream &outStream=std::cout)
Real feasToleranceInitial_
virtual void writeOutput(std::ostream &os, const bool print_header=false) const override
Print iterate status.
virtual void writeHeader(std::ostream &os) const override
Print iterate header.
Real optDecreaseExponent_
Real feasDecreaseExponent_
virtual void run(Problem< Real > &problem, std::ostream &outStream=std::cout) override
Run algorithm on equality constrained problems (Type-E). This is the primary Type-E interface.
Real feasIncreaseExponent_
virtual void writeName(std::ostream &os) const override
Print step name.
Real optIncreaseExponent_
Defines the linear algebra or vector space interface.
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 .