ROL
ROL_TypeE_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_TYPEE_AUGMENTEDLAGRANGIANALGORITHM_DEF_H
45#define ROL_TYPEE_AUGMENTEDLAGRANGIANALGORITHM_DEF_H
46
48
49namespace ROL {
50namespace TypeE {
51
52template<typename Real>
54 : TypeE::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,
107 Constraint<Real> &con,
108 std::ostream &outStream ) {
109 const Real one(1), TOL(1.e-2);
110 Real tol = std::sqrt(ROL_EPSILON<Real>());
112
113 // Initialize the algorithm state
114 state_->nfval = 0;
115 state_->ncval = 0;
116 state_->ngrad = 0;
117
118 // Compute objective value
119 alobj.update(x,UpdateType::Initial,state_->iter);
120 state_->value = alobj.getObjectiveValue(x,tol);
121 alobj.gradient(*state_->gradientVec,x,tol);
122
123 // Compute constraint violation
124 state_->constraintVec->set(*alobj.getConstraintVec(x,tol));
125 state_->cnorm = state_->constraintVec->norm();
126
127 // Update evaluation counters
128 state_->ncval += alobj.getNumberConstraintEvaluations();
129 state_->nfval += alobj.getNumberFunctionEvaluations();
130 state_->ngrad += alobj.getNumberGradientEvaluations();
131
132 // Compute problem scaling
133 if (useDefaultScaling_) {
134 fscale_ = one/std::max(one,alobj.getObjectiveGradient(x,tol)->norm());
135 try {
136 Ptr<Vector<Real>> ji = x.clone();
137 Real maxji(0), normji(0);
138 for (int i = 0; i < c.dimension(); ++i) {
139 con.applyAdjointJacobian(*ji,*c.basis(i),x,tol);
140 normji = ji->norm();
141 maxji = std::max(normji,maxji);
142 }
143 cscale_ = one/std::max(one,maxji);
144 }
145 catch (std::exception &e) {
146 cscale_ = one;
147 }
148 }
149 alobj.setScaling(fscale_,cscale_);
150
151 // Compute gradient of the lagrangian
152 state_->gnorm = state_->gradientVec->norm()/std::min(fscale_,cscale_);
153
154 // Compute initial penalty parameter
155 if (useRelTol_) outerOptTolerance_ *= state_->gnorm;
156 if (useDefaultInitPen_) {
157 const Real oem8(1e-8), oem2(1e-2), two(2), ten(10);
158 state_->searchSize = std::max(oem8,
159 std::min(ten*std::max(one,std::abs(fscale_*state_->value))
160 / std::max(one,std::pow(cscale_*state_->cnorm,two)),oem2*maxPenaltyParam_));
161 }
162 // Initialize intermediate stopping tolerances
163 minPenaltyReciprocal_ = std::min(one/state_->searchSize,minPenaltyLowerBound_);
164 optTolerance_ = std::max<Real>(TOL*outerOptTolerance_,
165 optToleranceInitial_*std::pow(minPenaltyReciprocal_,optDecreaseExponent_));
166 optTolerance_ = std::min<Real>(optTolerance_,TOL*state_->gnorm);
167 feasTolerance_ = std::max<Real>(TOL*outerFeasTolerance_,
168 feasToleranceInitial_*std::pow(minPenaltyReciprocal_,feasDecreaseExponent_));
169
170 // Set data
171 alobj.reset(l,state_->searchSize);
172
173 if (verbosity_ > 1) {
174 outStream << std::endl;
175 outStream << "Augmented Lagrangian Initialize" << std::endl;
176 outStream << "Objective Scaling: " << fscale_ << std::endl;
177 outStream << "Constraint Scaling: " << cscale_ << std::endl;
178 outStream << std::endl;
179 }
180}
181
182template<typename Real>
184 const Vector<Real> &g,
185 Objective<Real> &obj,
186 Constraint<Real> &econ,
187 Vector<Real> &emul,
188 const Vector<Real> &eres,
189 std::ostream &outStream ) {
190 const Real one(1), oem2(1e-2);
191 Real tol(std::sqrt(ROL_EPSILON<Real>()));
192 // Initialize augmented Lagrangian data
193 AugmentedLagrangianObjective<Real> alobj(makePtrFromRef(obj),makePtrFromRef(econ),
194 state_->searchSize,g,eres,emul,
195 scaleLagrangian_,HessianApprox_);
196 initialize(x,g,emul,eres,alobj,econ,outStream);
197 Ptr<TypeU::Algorithm<Real>> algo;
198
199 // Output
200 if (verbosity_ > 0) writeOutput(outStream,true);
201
202 while (status_->check(*state_)) {
203 // Solve unconstrained augmented Lagrangian subproblem
204 list_.sublist("Status Test").set("Gradient Tolerance",optTolerance_);
205 list_.sublist("Status Test").set("Step Tolerance",1.e-6*optTolerance_);
206 algo = TypeU::AlgorithmFactory<Real>(list_);
207 algo->run(x,g,alobj,outStream);
208 subproblemIter_ = algo->getState()->iter;
209
210 // Compute step
211 state_->stepVec->set(x);
212 state_->stepVec->axpy(-one,*state_->iterateVec);
213 state_->snorm = state_->stepVec->norm();
214
215 // Update iteration information
216 state_->iter++;
217 state_->iterateVec->set(x);
218 state_->value = alobj.getObjectiveValue(x,tol);
219 state_->constraintVec->set(*alobj.getConstraintVec(x,tol));
220 state_->cnorm = state_->constraintVec->norm();
221 alobj.gradient(*state_->gradientVec,x,tol);
222 if (scaleLagrangian_) {
223 state_->gradientVec->scale(state_->searchSize);
224 }
225 state_->gnorm = state_->gradientVec->norm()/std::min(fscale_,cscale_);
226 //alobj.update(x,UpdateType::Accept,state_->iter);
227
228 // Update evaluation counters
229 state_->nfval += alobj.getNumberFunctionEvaluations();
230 state_->ngrad += alobj.getNumberGradientEvaluations();
231 state_->ncval += alobj.getNumberConstraintEvaluations();
232
233 // Update multipliers
234 minPenaltyReciprocal_ = std::min(one/state_->searchSize,minPenaltyLowerBound_);
235 if ( cscale_*state_->cnorm < feasTolerance_ ) {
236 emul.axpy(state_->searchSize*cscale_,state_->constraintVec->dual());
237 if ( algo->getState()->statusFlag == EXITSTATUS_CONVERGED ) {
238 optTolerance_ = std::max(oem2*outerOptTolerance_,
239 optTolerance_*std::pow(minPenaltyReciprocal_,optIncreaseExponent_));
240 }
241 feasTolerance_ = std::max(oem2*outerFeasTolerance_,
242 feasTolerance_*std::pow(minPenaltyReciprocal_,feasIncreaseExponent_));
243 // Update Algorithm State
244 state_->snorm += state_->searchSize*cscale_*state_->cnorm;
245 state_->lagmultVec->set(emul);
246 }
247 else {
248 state_->searchSize = std::min(penaltyUpdate_*state_->searchSize,maxPenaltyParam_);
249 optTolerance_ = std::max(oem2*outerOptTolerance_,
250 optToleranceInitial_*std::pow(minPenaltyReciprocal_,optDecreaseExponent_));
251 feasTolerance_ = std::max(oem2*outerFeasTolerance_,
252 feasToleranceInitial_*std::pow(minPenaltyReciprocal_,feasDecreaseExponent_));
253 }
254 alobj.reset(emul,state_->searchSize);
255
256 // Update Output
257 if (verbosity_ > 0) writeOutput(outStream,printHeader_);
258 }
259 emul.scale(cscale_);
260 if (verbosity_ > 0) TypeE::Algorithm<Real>::writeExitStatus(outStream);
261}
262
263template<typename Real>
265 std::stringstream hist;
266 if(verbosity_>1) {
267 hist << std::string(114,'-') << std::endl;
268 hist << "Augmented Lagrangian status output definitions" << std::endl << std::endl;
269 hist << " iter - Number of iterates (steps taken)" << std::endl;
270 hist << " fval - Objective function value" << std::endl;
271 hist << " cnorm - Norm of the constraint violation" << std::endl;
272 hist << " gLnorm - Norm of the gradient of the Lagrangian" << std::endl;
273 hist << " snorm - Norm of the step" << std::endl;
274 hist << " penalty - Penalty parameter" << std::endl;
275 hist << " feasTol - Feasibility tolerance" << std::endl;
276 hist << " optTol - Optimality tolerance" << std::endl;
277 hist << " #fval - Number of times the objective was computed" << std::endl;
278 hist << " #grad - Number of times the gradient was computed" << std::endl;
279 hist << " #cval - Number of times the constraint was computed" << std::endl;
280 hist << " subIter - Number of iterations to solve subproblem" << std::endl;
281 hist << std::string(114,'-') << std::endl;
282 }
283 hist << " ";
284 hist << std::setw(6) << std::left << "iter";
285 hist << std::setw(15) << std::left << "fval";
286 hist << std::setw(15) << std::left << "cnorm";
287 hist << std::setw(15) << std::left << "gLnorm";
288 hist << std::setw(15) << std::left << "snorm";
289 hist << std::setw(10) << std::left << "penalty";
290 hist << std::setw(10) << std::left << "feasTol";
291 hist << std::setw(10) << std::left << "optTol";
292 hist << std::setw(8) << std::left << "#fval";
293 hist << std::setw(8) << std::left << "#grad";
294 hist << std::setw(8) << std::left << "#cval";
295 hist << std::setw(8) << std::left << "subIter";
296 hist << std::endl;
297 os << hist.str();
298}
299
300template<typename Real>
301void AugmentedLagrangianAlgorithm<Real>::writeName( std::ostream& os ) const {
302 std::stringstream hist;
303 hist << std::endl << "Augmented Lagrangian Solver (Type E, Equality Constraints)";
304 hist << std::endl;
305 hist << "Subproblem Solver: " << subStep_ << std::endl;
306 os << hist.str();
307}
308
309template<typename Real>
310void AugmentedLagrangianAlgorithm<Real>::writeOutput( std::ostream& os, const bool print_header ) const {
311 std::stringstream hist;
312 hist << std::scientific << std::setprecision(6);
313 if ( state_->iter == 0 ) writeName(os);
314 if ( print_header ) writeHeader(os);
315 if ( state_->iter == 0 ) {
316 hist << " ";
317 hist << std::setw(6) << std::left << state_->iter;
318 hist << std::setw(15) << std::left << state_->value;
319 hist << std::setw(15) << std::left << state_->cnorm;
320 hist << std::setw(15) << std::left << state_->gnorm;
321 hist << std::setw(15) << std::left << "---";
322 hist << std::scientific << std::setprecision(2);
323 hist << std::setw(10) << std::left << state_->searchSize;
324 hist << std::setw(10) << std::left << std::max(feasTolerance_,outerFeasTolerance_);
325 hist << std::setw(10) << std::left << std::max(optTolerance_,outerOptTolerance_);
326 hist << std::scientific << std::setprecision(6);
327 hist << std::setw(8) << std::left << state_->nfval;
328 hist << std::setw(8) << std::left << state_->ngrad;
329 hist << std::setw(8) << std::left << state_->ncval;
330 hist << std::setw(8) << std::left << "---";
331 hist << std::endl;
332 }
333 else {
334 hist << " ";
335 hist << std::setw(6) << std::left << state_->iter;
336 hist << std::setw(15) << std::left << state_->value;
337 hist << std::setw(15) << std::left << state_->cnorm;
338 hist << std::setw(15) << std::left << state_->gnorm;
339 hist << std::setw(15) << std::left << state_->snorm;
340 hist << std::scientific << std::setprecision(2);
341 hist << std::setw(10) << std::left << state_->searchSize;
342 hist << std::setw(10) << std::left << feasTolerance_;
343 hist << std::setw(10) << std::left << optTolerance_;
344 hist << std::scientific << std::setprecision(6);
345 hist << std::setw(8) << std::left << state_->nfval;
346 hist << std::setw(8) << std::left << state_->ngrad;
347 hist << std::setw(8) << std::left << state_->ncval;
348 hist << std::setw(8) << std::left << subproblemIter_;
349 hist << std::endl;
350 }
351 os << hist.str();
352}
353
354} // namespace TypeE
355} // namespace ROL
356
357#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 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.
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_
virtual void writeName(std::ostream &os) const override
Print step name.
virtual void writeHeader(std::ostream &os) const override
Print iterate header.
void initialize(Vector< Real > &x, const Vector< Real > &g, const Vector< Real > &l, const Vector< Real > &c, AugmentedLagrangianObjective< Real > &alobj, Constraint< Real > &con, std::ostream &outStream=std::cout)
virtual void writeOutput(std::ostream &os, const bool print_header=false) const override
Print iterate status.
virtual void run(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, Constraint< Real > &econ, Vector< Real > &emul, const Vector< Real > &eres, std::ostream &outStream=std::cout) override
Run algorithm on equality constrained problems (Type-E). This general interface supports the use of d...
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
virtual void scale(const Real alpha)=0
Compute 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.
Definition: ROL_Vector.hpp:196
virtual ROL::Ptr< Vector > basis(const int i) const
Return i-th basis vector.
Definition: ROL_Vector.hpp:182
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:153
@ EXITSTATUS_CONVERGED
Definition: ROL_Types.hpp:118