ROL
ROL_TypeG_MoreauYosidaAlgorithm_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_MOREAUYOSIDAALGORITHM_DEF_H
45#define ROL_TYPEG_MOREAUYOSIDAALGORITHM_DEF_H
46
48
49namespace ROL {
50namespace TypeG {
51
52template<typename Real>
54 : TypeG::Algorithm<Real>::Algorithm(),
55 tau_(10), print_(false), list_(list), subproblemIter_(0) {
56 // Set status test
57 status_->reset();
58 status_->add(makePtr<ConstraintStatusTest<Real>>(list));
59
60 // Parse parameters
61 Real ten(10), oem6(1.e-6), oem8(1.e-8), oe8(1e8);
62 ParameterList& steplist = list.sublist("Step").sublist("Moreau-Yosida Penalty");
63 state_->searchSize = steplist.get("Initial Penalty Parameter", ten);
64 maxPenalty_ = steplist.get("Maximum Penalty Parameter", oe8);
65 tau_ = steplist.get("Penalty Parameter Growth Factor", ten);
66 updatePenalty_ = steplist.get("Update Penalty", true);
67 updateMultiplier_ = steplist.get("Update Multiplier", true);
68 print_ = steplist.sublist("Subproblem").get("Print History", false);
69 // Set parameters for step subproblem
70 Real gtol = steplist.sublist("Subproblem").get("Optimality Tolerance", oem8);
71 Real ctol = steplist.sublist("Subproblem").get("Feasibility Tolerance", oem8);
72 int maxit = steplist.sublist("Subproblem").get("Iteration Limit", 1000);
73 Real stol = oem6*std::min(gtol,ctol);
74 list_.sublist("Status Test").set("Gradient Tolerance", gtol);
75 list_.sublist("Status Test").set("Constraint Tolerance", ctol);
76 list_.sublist("Status Test").set("Step Tolerance", stol);
77 list_.sublist("Status Test").set("Iteration Limit", maxit);
78 // Get step name from parameterlist
79 stepname_ = steplist.sublist("Subproblem").get("Step Type","Augmented Lagrangian");
80
81 // Output settings
82 verbosity_ = list.sublist("General").get("Output Level", 0);
84 print_ = (verbosity_ > 2 ? true : print_);
85 list_.sublist("General").set("Output Level",(print_ ? verbosity_ : 0));
86}
87
88template<typename Real>
90 const Vector<Real> &g,
91 const Vector<Real> &l,
92 const Vector<Real> &c,
96 Vector<Real> &pwa,
97 Vector<Real> &dwa,
98 std::ostream &outStream) {
99 hasPolyProj_ = true;
100 if (proj_ == nullPtr) {
101 proj_ = makePtr<PolyhedralProjection<Real>>(makePtrFromRef(bnd));
102 hasPolyProj_ = false;
103 }
104 proj_->project(x,outStream);
105 // Initialize data
107 // Initialize the algorithm state
108 state_->nfval = 0;
109 state_->ngrad = 0;
110 state_->ncval = 0;
111 updateState(x,l,myobj,bnd,con,pwa,dwa,outStream);
112}
113
114
115template<typename Real>
117 const Vector<Real> &l,
120 Constraint<Real> &con,
121 Vector<Real> &pwa,
122 Vector<Real> &dwa,
123 std::ostream &outStream) {
124 const Real one(1);
125 Real zerotol = std::sqrt(ROL_EPSILON<Real>());
126 // Update objective and constraint
127 if (state_->iter == 0) {
128 myobj.update(x,UpdateType::Initial,state_->iter);
129 con.update(x,UpdateType::Initial,state_->iter);
130 }
131 //else {
132 // myobj.update(x,UpdateType::Accept,state_->iter);
133 // con.update(x,UpdateType::Accept,state_->iter);
134 //}
135 // Compute norm of the gradient of the Lagrangian
136 state_->value = myobj.getObjectiveValue(x, zerotol);
137 myobj.getObjectiveGradient(*state_->gradientVec, x, zerotol);
138 //myobj.gradient(*state_->gradientVec, x, zerotol);
139 con.applyAdjointJacobian(dwa, l, x, zerotol);
140 state_->gradientVec->plus(dwa);
141 //gnorm_ = state_->gradientVec->norm();
142 pwa.set(x);
143 pwa.axpy(-one,state_->gradientVec->dual());
144 proj_->project(pwa,outStream);
145 pwa.axpy(-one,x);
146 gnorm_ = pwa.norm();
147 // Compute constraint violation
148 con.value(*state_->constraintVec, x, zerotol);
149 state_->cnorm = state_->constraintVec->norm();
150 compViolation_ = myobj.testComplementarity(x);
151 state_->gnorm = std::max(gnorm_,compViolation_);
152 // Update state
153 state_->nfval++;
154 state_->ngrad++;
155 state_->ncval++;
156}
157
158template<typename Real>
160 const Vector<Real> &g,
161 Objective<Real> &obj,
163 Constraint<Real> &econ,
164 Vector<Real> &emul,
165 const Vector<Real> &eres,
166 std::ostream &outStream ) {
167 const Real one(1);
168 Ptr<Vector<Real>> pwa = x.clone(), dwa = g.clone();
169 // Initialize Moreau-Yosida data
170 MoreauYosidaObjective<Real> myobj(makePtrFromRef(obj),makePtrFromRef(bnd),
171 x,g,state_->searchSize,updateMultiplier_,
172 updatePenalty_);
173 initialize(x,g,emul,eres,myobj,bnd,econ,*pwa,*dwa,outStream);
174 Ptr<TypeE::Algorithm<Real>> algo;
175
176 // Output
177 if (verbosity_ > 0) writeOutput(outStream,true);
178
179 while (status_->check(*state_)) {
180 // Solve augmented Lagrangian subproblem
181 algo = TypeE::AlgorithmFactory<Real>(list_);
182 emul.zero();
183 if (hasPolyProj_) algo->run(x,g,myobj,econ,emul,eres,
184 *proj_->getLinearConstraint(),
185 *proj_->getMultiplier(),
186 *proj_->getResidual(),outStream);
187 else algo->run(x,g,myobj,econ,emul,eres,outStream);
188 subproblemIter_ = algo->getState()->iter;
189 state_->nfval += algo->getState()->nfval;
190 state_->ngrad += algo->getState()->ngrad;
191 state_->ncval += algo->getState()->ncval;
192
193 // Compute step
194 state_->stepVec->set(x);
195 state_->stepVec->axpy(-one,*state_->iterateVec);
196 state_->snorm = state_->stepVec->norm();
197 state_->lagmultVec->axpy(-one,emul);
198 state_->snorm += state_->lagmultVec->norm();
199
200 // Update iterate and Lagrange multiplier
201 state_->iterateVec->set(x);
202 state_->lagmultVec->set(emul);
203
204 // Update objective and constraint
205 state_->iter++;
206
207 // Update state
208 updateState(x,emul,myobj,bnd,econ,*pwa,*dwa);
209
210 // Update multipliers
211 if (updatePenalty_) {
212 state_->searchSize *= tau_;
213 state_->searchSize = std::min(state_->searchSize,maxPenalty_);
214 }
215 myobj.updateMultipliers(state_->searchSize,x);
216
217 // Update Output
218 if (verbosity_ > 0) writeOutput(outStream,printHeader_);
219 }
220 if (verbosity_ > 0) TypeG::Algorithm<Real>::writeExitStatus(outStream);
221}
222
223template<typename Real>
224void MoreauYosidaAlgorithm<Real>::writeHeader( std::ostream& os ) const {
225 std::stringstream hist;
226 if (verbosity_ > 1) {
227 hist << std::string(109,'-') << std::endl;
228 hist << "Moreau-Yosida Penalty Solver";
229 hist << " status output definitions" << std::endl << std::endl;
230 hist << " iter - Number of iterates (steps taken)" << std::endl;
231 hist << " fval - Objective function value" << std::endl;
232 hist << " cnorm - Norm of the constraint" << std::endl;
233 hist << " gLnorm - Norm of the gradient of the Lagrangian" << std::endl;
234 hist << " ifeas - Infeasibility metric" << std::endl;
235 hist << " snorm - Norm of the step (update to optimization vector)" << std::endl;
236 hist << " penalty - Penalty parameter for bound constraints" << std::endl;
237 hist << " #fval - Cumulative number of times the objective function was evaluated" << std::endl;
238 hist << " #grad - Cumulative number of times the gradient was computed" << std::endl;
239 hist << " #cval - Cumulative number of times the constraint was evaluated" << std::endl;
240 hist << " subiter - Number of subproblem iterations" << std::endl;
241 hist << std::string(109,'-') << std::endl;
242 }
243
244 hist << " ";
245 hist << std::setw(6) << std::left << "iter";
246 hist << std::setw(15) << std::left << "fval";
247 hist << std::setw(15) << std::left << "cnorm";
248 hist << std::setw(15) << std::left << "gLnorm";
249 hist << std::setw(15) << std::left << "ifeas";
250 hist << std::setw(15) << std::left << "snorm";
251 hist << std::setw(10) << std::left << "penalty";
252 hist << std::setw(8) << std::left << "#fval";
253 hist << std::setw(8) << std::left << "#grad";
254 hist << std::setw(8) << std::left << "#cval";
255 hist << std::setw(8) << std::left << "subIter";
256 hist << std::endl;
257 os << hist.str();
258}
259
260template<typename Real>
261void MoreauYosidaAlgorithm<Real>::writeName( std::ostream& os ) const {
262 std::stringstream hist;
263 hist << std::endl << "Moreau-Yosida Penalty Solver (Type G, General Constraints)";
264 hist << std::endl;
265 hist << "Subproblem Solver: " << stepname_ << std::endl;
266 os << hist.str();
267}
268
269template<typename Real>
270void MoreauYosidaAlgorithm<Real>::writeOutput( std::ostream& os, const bool print_header ) const {
271 std::stringstream hist;
272 hist << std::scientific << std::setprecision(6);
273 if ( state_->iter == 0 ) writeName(os);
274 if ( print_header ) writeHeader(os);
275 if ( state_->iter == 0 ) {
276 hist << " ";
277 hist << std::setw(6) << std::left << state_->iter;
278 hist << std::setw(15) << std::left << state_->value;
279 hist << std::setw(15) << std::left << state_->cnorm;
280 hist << std::setw(15) << std::left << gnorm_;
281 hist << std::setw(15) << std::left << compViolation_;
282 hist << std::setw(15) << std::left << "---";
283 hist << std::scientific << std::setprecision(2);
284 hist << std::setw(10) << std::left << state_->searchSize;
285 hist << std::setw(8) << std::left << state_->nfval;
286 hist << std::setw(8) << std::left << state_->ngrad;
287 hist << std::setw(8) << std::left << state_->ncval;
288 hist << std::setw(8) << std::left << "---";
289 hist << std::endl;
290 }
291 else {
292 hist << " ";
293 hist << std::setw(6) << std::left << state_->iter;
294 hist << std::setw(15) << std::left << state_->value;
295 hist << std::setw(15) << std::left << state_->cnorm;
296 hist << std::setw(15) << std::left << gnorm_;
297 hist << std::setw(15) << std::left << compViolation_;
298 hist << std::setw(15) << std::left << state_->snorm;
299 hist << std::scientific << std::setprecision(2);
300 hist << std::setw(10) << std::left << state_->searchSize;
301 hist << std::scientific << std::setprecision(6);
302 hist << std::setw(8) << std::left << state_->nfval;
303 hist << std::setw(8) << std::left << state_->ngrad;
304 hist << std::setw(8) << std::left << state_->ncval;
305 hist << std::setw(8) << std::left << subproblemIter_;
306 hist << std::endl;
307 }
308 os << hist.str();
309}
310
311} // namespace TypeG
312} // namespace ROL
313
314#endif
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 value(Vector< Real > &c, const Vector< Real > &x, Real &tol)=0
Evaluate the constraint operator at .
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 update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update constraint function.
Provides the interface to evaluate the Moreau-Yosida penalty function.
void getObjectiveGradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Real getObjectiveValue(const Vector< Real > &x, Real &tol)
Real testComplementarity(const Vector< Real > &x)
void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update Moreau-Yosida penalty function.
void updateMultipliers(Real mu, const Vector< Real > &x)
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 writeName(std::ostream &os) const override
Print step name.
void initialize(Vector< Real > &x, const Vector< Real > &g, const Vector< Real > &l, const Vector< Real > &c, MoreauYosidaObjective< Real > &myobj, BoundConstraint< Real > &bnd, Constraint< Real > &con, Vector< Real > &pwa, Vector< Real > &dwa, std::ostream &outStream=std::cout)
void updateState(const Vector< Real > &x, const Vector< Real > &l, MoreauYosidaObjective< Real > &myobj, BoundConstraint< Real > &bnd, Constraint< Real > &con, Vector< Real > &pwa, Vector< Real > &dwa, std::ostream &outStream=std::cout)
void writeOutput(std::ostream &os, const bool print_header=false) const override
Print iterate status.
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 writeHeader(std::ostream &os) const override
Print iterate header.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:84
virtual Real norm() const =0
Returns where .
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:209
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:167
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:153