44 #ifndef ROL_TYPEG_STABILIZEDLCLALGORITHM_DEF_H 45 #define ROL_TYPEG_STABILIZEDLCLALGORITHM_DEF_H 53 template<
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);
100 template<
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;
192 template<
typename Real>
194 std::ostream &outStream ) {
197 problem.
finalize(
true,verbosity_>3,outStream);
213 template<
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_);
342 template<
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";
381 template<
typename Real>
383 std::stringstream hist;
384 hist << std::endl <<
"Stabilized LCL Solver (Type G, General Constraints)";
386 hist <<
"Subproblem Solver: " << subStep_ << std::endl;
390 template<
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 evaluate objective functions.
const Ptr< Constraint< Real > > & getConstraint()
Get the equality constraint.
Real optToleranceInitial_
Real feasIncreaseExponent_
const Ptr< Vector< Real > > & getPrimalOptimizationVector()
Get the primal optimization space vector.
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
Real feasDecreaseExponent_
void addBoundConstraint(const Ptr< BoundConstraint< Real >> &bnd)
Add a bound constraint.
const Ptr< const Vector< Real > > getConstraintVec(const Vector< Real > &x, Real &tol)
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
int getNumberGradientEvaluations(void) const
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.
Real optIncreaseExponent_
const Ptr< BoundConstraint< Real > > & getBoundConstraint()
Get the bound constraint.
virtual void writeHeader(std::ostream &os) const override
Print iterate header.
Real feasToleranceInitial_
virtual ROL::Ptr< Vector > basis(const int i) const
Return i-th basis vector.
Defines the linear algebra or vector space interface.
int getNumberConstraintEvaluations(void) const
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...
Real getObjectiveValue(const Vector< Real > &x, Real &tol)
StabilizedLCLAlgorithm(ParameterList &list)
Provides an interface to check status of optimization algorithms for problems with equality constrain...
Provides an interface to run general constrained optimization algorithms.
virtual void writeExitStatus(std::ostream &os) const
void setScaling(const Real fscale=1.0, const Real cscale=1.0)
int getNumberFunctionEvaluations(void) const
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...
Provides the interface to evaluate the elastic augmented Lagrangian.
const Ptr< Objective< Real > > & getObjective()
Get the objective function.
virtual void writeName(std::ostream &os) const override
Print step name.
const Ptr< AlgorithmState< Real > > state_
const Ptr< const Vector< Real > > getObjectiveGradient(const Vector< Real > &x, Real &tol)
virtual int dimension() const
Return dimension of the vector space.
const Ptr< Vector< Real > > & getResidualVector()
Get the primal constraint space vector.
const Ptr< Vector< Real > > & getMultiplierVector()
Get the dual constraint space vector.
void finalizeIteration()
Transform the optimization variables to the native parameterization after an optimization algorithm h...
const Ptr< AugmentedLagrangianObjective< Real > > getAugmentedLagrangian(void) const
Provides the interface to apply upper and lower bound constraints.
const Ptr< PolyhedralProjection< Real > > & getPolyhedralProjection()
Get the polyhedral projection object. This is a null pointer if no linear constraints and/or bounds a...
EProblem getProblemType()
Get the optimization problem type (U, B, E, or G).
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 .
virtual void writeOutput(std::ostream &os, const bool print_header=false) const override
Print iterate status.
virtual void set(const Vector &x)
Set where .
virtual Real norm() const =0
Returns where .
Real optDecreaseExponent_
void reset(const Vector< Real > &multiplier, Real penaltyParameter, Real sigma)
virtual void edit()
Resume editting optimization problem after finalize has been called.
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)
const Ptr< CombinedStatusTest< Real > > status_
const Ptr< Vector< Real > > & getDualOptimizationVector()
Get the dual optimization space vector.
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...
void initialize(const Vector< Real > &x, const Vector< Real > &g, const Vector< Real > &mul, const Vector< Real > &c)
Defines the general constraint operator interface.
void setProjectionAlgorithm(ParameterList &parlist)
Set polyhedral projection algorithm.