ROL
ROL_TypeG_AugmentedLagrangianAlgorithm_Def.hpp
Go to the documentation of this file.
1// @HEADER
2// ************************************************************************
3//
4// Rapid Optimization Library (ROL) Package
5// Copyright (2014) Sandia Corporation
6//
7// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8// license for use of this work by or on behalf of the U.S. Government.
9//
10// Redistribution and use in source and binary forms, with or without
11// modification, are permitted provided that the following conditions are
12// met:
13//
14// 1. Redistributions of source code must retain the above copyright
15// notice, this list of conditions and the following disclaimer.
16//
17// 2. Redistributions in binary form must reproduce the above copyright
18// notice, this list of conditions and the following disclaimer in the
19// documentation and/or other materials provided with the distribution.
20//
21// 3. Neither the name of the Corporation nor the names of the
22// contributors may be used to endorse or promote products derived from
23// this software without specific prior written permission.
24//
25// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36//
37// Questions? Contact lead developers:
38// Drew Kouri (dpkouri@sandia.gov) and
39// Denis Ridzal (dridzal@sandia.gov)
40//
41// ************************************************************************
42// @HEADER
43
44#ifndef ROL_TYPEG_AUGMENTEDLAGRANGIANALGORITHM_DEF_H
45#define ROL_TYPEG_AUGMENTEDLAGRANGIANALGORITHM_DEF_H
46
48
49namespace ROL {
50namespace TypeG {
51
52template<typename Real>
54 : TypeG::Algorithm<Real>::Algorithm(), list_(list), subproblemIter_(0) {
55 // Set status test
56 status_->reset();
57 status_->add(makePtr<ConstraintStatusTest<Real>>(list));
58
59 Real one(1), p1(0.1), p9(0.9), ten(1.e1), oe8(1.e8), oem8(1.e-8);
60 ParameterList& sublist = list.sublist("Step").sublist("Augmented Lagrangian");
61 useDefaultInitPen_ = sublist.get("Use Default Initial Penalty Parameter", true);
62 state_->searchSize = sublist.get("Initial Penalty Parameter", ten);
63 // Multiplier update parameters
64 scaleLagrangian_ = sublist.get("Use Scaled Augmented Lagrangian", false);
65 minPenaltyLowerBound_ = sublist.get("Penalty Parameter Reciprocal Lower Bound", p1);
66 penaltyUpdate_ = sublist.get("Penalty Parameter Growth Factor", ten);
67 maxPenaltyParam_ = sublist.get("Maximum Penalty Parameter", oe8);
69 // Optimality tolerance update
70 optIncreaseExponent_ = sublist.get("Optimality Tolerance Update Exponent", one);
71 optDecreaseExponent_ = sublist.get("Optimality Tolerance Decrease Exponent", one);
72 optToleranceInitial_ = sublist.get("Initial Optimality Tolerance", one);
73 // Feasibility tolerance update
74 feasIncreaseExponent_ = sublist.get("Feasibility Tolerance Update Exponent", p1);
75 feasDecreaseExponent_ = sublist.get("Feasibility Tolerance Decrease Exponent", p9);
76 feasToleranceInitial_ = sublist.get("Initial Feasibility Tolerance", one);
77 // Subproblem information
78 print_ = sublist.get("Print Intermediate Optimization History", false);
79 maxit_ = sublist.get("Subproblem Iteration Limit", 1000);
80 subStep_ = sublist.get("Subproblem Step Type", "Trust Region");
81 HessianApprox_ = sublist.get("Level of Hessian Approximation", 0);
82 list_.sublist("Step").set("Type",subStep_);
83 list_.sublist("Status Test").set("Iteration Limit",maxit_);
84 list_.sublist("Status Test").set("Use Relative Tolerances",false);
85 // Verbosity setting
86 verbosity_ = list.sublist("General").get("Output Level", 0);
88 print_ = (verbosity_ > 2 ? true : print_);
89 list_.sublist("General").set("Output Level",(print_ ? verbosity_ : 0));
90 // Outer iteration tolerances
91 outerFeasTolerance_ = list.sublist("Status Test").get("Constraint Tolerance", oem8);
92 outerOptTolerance_ = list.sublist("Status Test").get("Gradient Tolerance", oem8);
93 outerStepTolerance_ = list.sublist("Status Test").get("Step Tolerance", oem8);
94 useRelTol_ = list.sublist("Status Test").get("Use Relative Tolerances", false);
95 // Scaling
96 useDefaultScaling_ = sublist.get("Use Default Problem Scaling", true);
97 fscale_ = sublist.get("Objective Scaling", one);
98 cscale_ = sublist.get("Constraint Scaling", one);
99}
100
101template<typename Real>
103 const Vector<Real> &g,
104 const Vector<Real> &l,
105 const Vector<Real> &c,
108 Constraint<Real> &con,
109 std::ostream &outStream ) {
110 hasPolyProj_ = true;
111 if (proj_ == nullPtr) {
112 proj_ = makePtr<PolyhedralProjection<Real>>(makePtrFromRef(bnd));
113 hasPolyProj_ = false;
114 }
115 proj_->project(x,outStream);
116
117 const Real one(1), TOL(1.e-2);
118 Real tol = std::sqrt(ROL_EPSILON<Real>());
120
121 // Initialize the algorithm state
122 state_->nfval = 0;
123 state_->ncval = 0;
124 state_->ngrad = 0;
125
126 // Compute objective value
127 alobj.update(x,UpdateType::Initial,state_->iter);
128 state_->value = alobj.getObjectiveValue(x,tol);
129 alobj.gradient(*state_->gradientVec,x,tol);
130
131 // Compute constraint violation
132 state_->constraintVec->set(*alobj.getConstraintVec(x,tol));
133 state_->cnorm = state_->constraintVec->norm();
134
135 // Update evaluation counters
136 state_->ncval += alobj.getNumberConstraintEvaluations();
137 state_->nfval += alobj.getNumberFunctionEvaluations();
138 state_->ngrad += alobj.getNumberGradientEvaluations();
139
140 // Compute problem scaling
141 if (useDefaultScaling_) {
142 fscale_ = one/std::max(one,alobj.getObjectiveGradient(x,tol)->norm());
143 try {
144 Ptr<Vector<Real>> ji = x.clone();
145 Real maxji(0), normji(0);
146 for (int i = 0; i < c.dimension(); ++i) {
147 con.applyAdjointJacobian(*ji,*c.basis(i),x,tol);
148 normji = ji->norm();
149 maxji = std::max(normji,maxji);
150 }
151 cscale_ = one/std::max(one,maxji);
152 }
153 catch (std::exception &e) {
154 cscale_ = one;
155 }
156 }
157 alobj.setScaling(fscale_,cscale_);
158
159 // Compute gradient of the lagrangian
160 x.axpy(-one,state_->gradientVec->dual());
161 proj_->project(x,outStream);
162 x.axpy(-one/std::min(fscale_,cscale_),*state_->iterateVec);
163 state_->gnorm = x.norm();
164 x.set(*state_->iterateVec);
165
166 // Compute initial penalty parameter
167 if (useDefaultInitPen_) {
168 const Real oem8(1e-8), oem2(1e-2), two(2), ten(10);
169 state_->searchSize = std::max(oem8,
170 std::min(ten*std::max(one,std::abs(fscale_*state_->value))
171 / std::max(one,std::pow(cscale_*state_->cnorm,two)),oem2*maxPenaltyParam_));
172 }
173 // Initialize intermediate stopping tolerances
174 if (useRelTol_) outerOptTolerance_ *= state_->gnorm;
175 minPenaltyReciprocal_ = std::min(one/state_->searchSize,minPenaltyLowerBound_);
176 optTolerance_ = std::max<Real>(TOL*outerOptTolerance_,
177 optToleranceInitial_*std::pow(minPenaltyReciprocal_,optDecreaseExponent_));
178 optTolerance_ = std::min<Real>(optTolerance_,TOL*state_->gnorm);
179 feasTolerance_ = std::max<Real>(TOL*outerFeasTolerance_,
180 feasToleranceInitial_*std::pow(minPenaltyReciprocal_,feasDecreaseExponent_));
181
182 // Set data
183 alobj.reset(l,state_->searchSize);
184
185 if (verbosity_ > 1) {
186 outStream << std::endl;
187 outStream << "Augmented Lagrangian Initialize" << std::endl;
188 outStream << "Objective Scaling: " << fscale_ << std::endl;
189 outStream << "Constraint Scaling: " << cscale_ << std::endl;
190 outStream << "Penalty Parameter: " << state_->searchSize << std::endl;
191 outStream << std::endl;
192 }
193}
194
195template<typename Real>
197 const Vector<Real> &g,
198 Objective<Real> &obj,
200 Constraint<Real> &econ,
201 Vector<Real> &emul,
202 const Vector<Real> &eres,
203 std::ostream &outStream ) {
204 const Real one(1), oem2(1e-2);
205 Real tol(std::sqrt(ROL_EPSILON<Real>()));
206 // Initialize augmented Lagrangian data
207 AugmentedLagrangianObjective<Real> alobj(makePtrFromRef(obj),makePtrFromRef(econ),
208 state_->searchSize,g,eres,emul,
209 scaleLagrangian_,HessianApprox_);
210 initialize(x,g,emul,eres,alobj,bnd,econ,outStream);
211 Ptr<TypeB::Algorithm<Real>> algo;
212
213 // Output
214 if (verbosity_ > 0) writeOutput(outStream,true);
215
216 while (status_->check(*state_)) {
217 // Solve unconstrained augmented Lagrangian subproblem
218 list_.sublist("Status Test").set("Gradient Tolerance",optTolerance_);
219 list_.sublist("Status Test").set("Step Tolerance",1.e-6*optTolerance_);
220 algo = TypeB::AlgorithmFactory<Real>(list_);
221 if (hasPolyProj_) algo->run(x,g,alobj,bnd,*proj_->getLinearConstraint(),
222 *proj_->getMultiplier(),*proj_->getResidual(),
223 outStream);
224 else algo->run(x,g,alobj,bnd,outStream);
225 subproblemIter_ = algo->getState()->iter;
226
227 // Compute step
228 state_->stepVec->set(x);
229 state_->stepVec->axpy(-one,*state_->iterateVec);
230 state_->snorm = state_->stepVec->norm();
231
232 // Update iteration information
233 state_->iter++;
234 state_->iterateVec->set(x);
235 state_->value = alobj.getObjectiveValue(x,tol);
236 state_->constraintVec->set(*alobj.getConstraintVec(x,tol));
237 state_->cnorm = state_->constraintVec->norm();
238 alobj.gradient(*state_->gradientVec,x,tol);
239 if (scaleLagrangian_) {
240 state_->gradientVec->scale(state_->searchSize);
241 }
242 x.axpy(-one/std::min(fscale_,cscale_),state_->gradientVec->dual());
243 proj_->project(x,outStream);
244 x.axpy(-one,*state_->iterateVec);
245 state_->gnorm = x.norm();
246 x.set(*state_->iterateVec);
247 //alobj.update(x,UpdateType::Accept,state_->iter);
248
249 // Update evaluation counters
250 state_->nfval += alobj.getNumberFunctionEvaluations();
251 state_->ngrad += alobj.getNumberGradientEvaluations();
252 state_->ncval += alobj.getNumberConstraintEvaluations();
253
254 // Update multipliers
255 if ( algo->getState()->statusFlag == EXITSTATUS_CONVERGED ) {
256 minPenaltyReciprocal_ = std::min(one/state_->searchSize,minPenaltyLowerBound_);
257 if ( cscale_*state_->cnorm < feasTolerance_ ) {
258 emul.axpy(state_->searchSize*cscale_,state_->constraintVec->dual());
259 optTolerance_ = std::max(oem2*outerOptTolerance_,
260 optTolerance_*std::pow(minPenaltyReciprocal_,optIncreaseExponent_));
261 feasTolerance_ = std::max(oem2*outerFeasTolerance_,
262 feasTolerance_*std::pow(minPenaltyReciprocal_,feasIncreaseExponent_));
263 // Update Algorithm State
264 state_->snorm += state_->searchSize*cscale_*state_->cnorm;
265 state_->lagmultVec->set(emul);
266 }
267 else {
268 state_->searchSize = std::min(penaltyUpdate_*state_->searchSize,maxPenaltyParam_);
269 optTolerance_ = std::max(oem2*outerOptTolerance_,
270 optToleranceInitial_*std::pow(minPenaltyReciprocal_,optDecreaseExponent_));
271 feasTolerance_ = std::max(oem2*outerFeasTolerance_,
272 feasToleranceInitial_*std::pow(minPenaltyReciprocal_,feasDecreaseExponent_));
273 }
274 alobj.reset(emul,state_->searchSize);
275 }
276
277 // Update Output
278 if (verbosity_ > 0) writeOutput(outStream,printHeader_);
279 }
280 if (verbosity_ > 0) TypeG::Algorithm<Real>::writeExitStatus(outStream);
281}
282
283template<typename Real>
285 std::stringstream hist;
286 if(verbosity_>1) {
287 hist << std::string(114,'-') << std::endl;
288 hist << "Augmented Lagrangian status output definitions" << std::endl << std::endl;
289 hist << " iter - Number of iterates (steps taken)" << std::endl;
290 hist << " fval - Objective function value" << std::endl;
291 hist << " cnorm - Norm of the constraint violation" << std::endl;
292 hist << " gLnorm - Norm of the gradient of the Lagrangian" << std::endl;
293 hist << " snorm - Norm of the step" << std::endl;
294 hist << " penalty - Penalty parameter" << std::endl;
295 hist << " feasTol - Feasibility tolerance" << std::endl;
296 hist << " optTol - Optimality tolerance" << std::endl;
297 hist << " #fval - Number of times the objective was computed" << std::endl;
298 hist << " #grad - Number of times the gradient was computed" << std::endl;
299 hist << " #cval - Number of times the constraint was computed" << std::endl;
300 hist << " subIter - Number of iterations to solve subproblem" << std::endl;
301 hist << std::string(114,'-') << std::endl;
302 }
303 hist << " ";
304 hist << std::setw(6) << std::left << "iter";
305 hist << std::setw(15) << std::left << "fval";
306 hist << std::setw(15) << std::left << "cnorm";
307 hist << std::setw(15) << std::left << "gLnorm";
308 hist << std::setw(15) << std::left << "snorm";
309 hist << std::setw(10) << std::left << "penalty";
310 hist << std::setw(10) << std::left << "feasTol";
311 hist << std::setw(10) << std::left << "optTol";
312 hist << std::setw(8) << std::left << "#fval";
313 hist << std::setw(8) << std::left << "#grad";
314 hist << std::setw(8) << std::left << "#cval";
315 hist << std::setw(8) << std::left << "subIter";
316 hist << std::endl;
317 os << hist.str();
318}
319
320template<typename Real>
321void AugmentedLagrangianAlgorithm<Real>::writeName( std::ostream& os ) const {
322 std::stringstream hist;
323 hist << std::endl << "Augmented Lagrangian Solver (Type G, General Constraints)";
324 hist << std::endl;
325 hist << "Subproblem Solver: " << subStep_ << std::endl;
326 os << hist.str();
327}
328
329template<typename Real>
330void AugmentedLagrangianAlgorithm<Real>::writeOutput( std::ostream& os, const bool print_header ) const {
331 std::stringstream hist;
332 hist << std::scientific << std::setprecision(6);
333 if ( state_->iter == 0 ) writeName(os);
334 if ( print_header ) writeHeader(os);
335 if ( state_->iter == 0 ) {
336 hist << " ";
337 hist << std::setw(6) << std::left << state_->iter;
338 hist << std::setw(15) << std::left << state_->value;
339 hist << std::setw(15) << std::left << state_->cnorm;
340 hist << std::setw(15) << std::left << state_->gnorm;
341 hist << std::setw(15) << std::left << "---";
342 hist << std::scientific << std::setprecision(2);
343 hist << std::setw(10) << std::left << state_->searchSize;
344 hist << std::setw(10) << std::left << std::max(feasTolerance_,outerFeasTolerance_);
345 hist << std::setw(10) << std::left << std::max(optTolerance_,outerOptTolerance_);
346 hist << std::scientific << std::setprecision(6);
347 hist << std::setw(8) << std::left << state_->nfval;
348 hist << std::setw(8) << std::left << state_->ngrad;
349 hist << std::setw(8) << std::left << state_->ncval;
350 hist << std::setw(8) << std::left << "---";
351 hist << std::endl;
352 }
353 else {
354 hist << " ";
355 hist << std::setw(6) << std::left << state_->iter;
356 hist << std::setw(15) << std::left << state_->value;
357 hist << std::setw(15) << std::left << state_->cnorm;
358 hist << std::setw(15) << std::left << state_->gnorm;
359 hist << std::setw(15) << std::left << state_->snorm;
360 hist << std::scientific << std::setprecision(2);
361 hist << std::setw(10) << std::left << state_->searchSize;
362 hist << std::setw(10) << std::left << feasTolerance_;
363 hist << std::setw(10) << std::left << optTolerance_;
364 hist << std::scientific << std::setprecision(6);
365 hist << std::setw(8) << std::left << state_->nfval;
366 hist << std::setw(8) << std::left << state_->ngrad;
367 hist << std::setw(8) << std::left << state_->ncval;
368 hist << std::setw(8) << std::left << subproblemIter_;
369 hist << std::endl;
370 }
371 os << hist.str();
372}
373
374} // namespace TypeG
375} // namespace ROL
376
377#endif
Provides the interface to evaluate the augmented Lagrangian.
void reset(const Vector< Real > &multiplier, const Real penaltyParameter)
void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update objective function.
Real getObjectiveValue(const Vector< Real > &x, Real &tol)
void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
void setScaling(const Real fscale=1.0, const Real cscale=1.0)
const Ptr< const Vector< Real > > getObjectiveGradient(const Vector< Real > &x, Real &tol)
const Ptr< const Vector< Real > > getConstraintVec(const Vector< Real > &x, Real &tol)
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 objective functions.
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_
void writeOutput(std::ostream &os, const bool print_header=false) const override
Print iterate status.
void writeName(std::ostream &os) const override
Print step name.
void run(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &bnd, Constraint< Real > &econ, Vector< Real > &emul, const Vector< Real > &eres, 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, AugmentedLagrangianObjective< Real > &alobj, BoundConstraint< Real > &bnd, Constraint< Real > &con, std::ostream &outStream=std::cout)
void writeHeader(std::ostream &os) const override
Print iterate header.
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 .
@ EXITSTATUS_CONVERGED