ROL
ROL_TypeB_PrimalDualActiveSetAlgorithm_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_TYPEB_PRIMALDUALACTIVESETALGORITHM_DEF_HPP
45#define ROL_TYPEB_PRIMALDUALACTIVESETALGORITHM_DEF_HPP
46
47namespace ROL {
48namespace TypeB {
49
50template<typename Real>
52 const Ptr<Secant<Real>> &secant)
53 : secant_(secant), esec_(SECANT_USERDEFINED),
54 neps_(-ROL_EPSILON<Real>()), itol_(std::sqrt(ROL_EPSILON<Real>())), hasPoly_(true) {
55 // Set status test
56 status_->reset();
57 status_->add(makePtr<StatusTest<Real>>(list));
58
59 if ( secant_ == nullPtr ) {
60 secantName_ = list.sublist("General").sublist("Secant").get("Type","Limited-Memory BFGS");
62 secant_ = SecantFactory<Real>(list);
63 }
64 else {
65 secantName_ = list.sublist("General").sublist("Secant").get("User Defined Secant Name",
66 "Unspecified User Defined Secant Method");
67 }
68 useSecantHessVec_ = list.sublist("General").sublist("Secant").get("Use as Hessian", false);
69 useSecantPrecond_ = list.sublist("General").sublist("Secant").get("Use as Preconditioner", false);
70
71 // Algorithmic parameters
72 maxit_ = list.sublist("Step").sublist("Primal Dual Active Set").get("Iteration Limit", 10);
73 stol_ = list.sublist("Step").sublist("Primal Dual Active Set").get("Relative Step Tolerance", 1e-8);
74 gtol_ = list.sublist("Step").sublist("Primal Dual Active Set").get("Relative Gradient Tolerance", 1e-6);
75 scale_ = list.sublist("Step").sublist("Primal Dual Active Set").get("Dual Scaling", 1.0);
76
77 // Krylov parameters
78 atolKrylov_ = list.sublist("General").sublist("Krylov").get("Absolute Tolerance", 1e-4);
79 rtolKrylov_ = list.sublist("General").sublist("Krylov").get("Relative Tolerance", 1e-2);
80 maxitKrylov_ = list.sublist("General").sublist("Krylov").get("Iteration Limit", 100);
81
82 verbosity_ = list.sublist("General").get("Output Level", 0);
84}
85
86template<typename Real>
88 const Vector<Real> &g,
89 Objective<Real> &obj,
91 std::ostream &outStream) {
92 // Initialize projection operator
93 if (proj_ == nullPtr) {
94 proj_ = makePtr<PolyhedralProjection<Real>>(makePtrFromRef(bnd));
95 hasPoly_ = false;
96 }
97 // Create Krylov solver
98 if (hasPoly_) {
99 ParameterList list;
100 list.sublist("General").sublist("Krylov").set("Absolute Tolerance", atolKrylov_);
101 list.sublist("General").sublist("Krylov").set("Relative Tolerance", rtolKrylov_);
102 list.sublist("General").sublist("Krylov").set("Iteration Limit", maxitKrylov_);
103 krylovName_ = "GMRES";
104 krylov_ = makePtr<GMRES<Real>>(list);
105 }
106 else {
107 krylovName_ = "CR";
108 krylov_ = makePtr<ConjugateResiduals<Real>>(atolKrylov_,rtolKrylov_,maxitKrylov_);
109 }
110 ekv_ = StringToEKrylov(krylovName_);
111 // Initialize data
112 const Real one(1);
114 // Update approximate gradient and approximate objective function.
115 Real ftol = std::sqrt(ROL_EPSILON<Real>());
116 proj_->project(x,outStream);
117 state_->iterateVec->set(x);
118 obj.update(x,UpdateType::Initial,state_->iter);
119 state_->value = obj.value(x,ftol); state_->nfval++;
120 obj.gradient(*state_->gradientVec,x,ftol); state_->ngrad++;
121 state_->stepVec->set(x);
122 state_->stepVec->axpy(-one,state_->gradientVec->dual());
123 proj_->project(*state_->stepVec,outStream);
124 state_->stepVec->axpy(-one,x);
125 state_->gnorm = state_->stepVec->norm();
126 state_->snorm = ROL_INF<Real>();
127}
128
129template<typename Real>
131 const Vector<Real> &g,
132 Objective<Real> &obj,
134 std::ostream &outStream ) {
135 const Real zero(0), one(1);
136 // Initialize PDAS data
137 initialize(x,g,obj,bnd,outStream);
138 Ptr<Vector<Real>> xlam = x.clone(), xtmp = x.clone(), As = x.clone();
139 Ptr<Vector<Real>> lambda = x.clone(), gtmp = g.clone();
140 Ptr<Vector<Real>> pwa = x.clone(), dwa = g.clone();
141 Ptr<Vector<Real>> mu, b, r;
142 if (hasPoly_) {
143 mu = proj_->getMultiplier()->clone(); mu->zero();
144 b = proj_->getResidual()->clone();
145 r = proj_->getResidual()->clone();
146 }
147 lambda->set(state_->gradientVec->dual());
148 lambda->scale(-one);
149 Real xnorm(0), snorm(0), rnorm(0), tol(std::sqrt(ROL_EPSILON<Real>()));
150
151 Ptr<LinearOperator<Real>> hessian, precond;
152
153 // Output
154 if (verbosity_ > 0) writeOutput(outStream,true);
155
156 while (status_->check(*state_)) {
157 totalKrylov_ = 0;
158 state_->stepVec->zero();
159 xnorm = x.norm();
160 if (hasPoly_) {
161 proj_->getLinearConstraint()->value(*r,x,tol);
162 }
163 for ( iter_ = 0; iter_ < maxit_; iter_++ ) {
164 /********************************************************************/
165 // MODIFY ITERATE VECTOR TO CHECK ACTIVE SET
166 /********************************************************************/
167 xlam->set(x); // xlam = x0
168 xlam->axpy(scale_,*lambda); // xlam = x0 + c*lambda
169 /********************************************************************/
170 // PROJECT x ONTO PRIMAL DUAL FEASIBLE SET
171 /********************************************************************/
172 As->zero(); // As = 0
173
174 xtmp->set(*bnd.getUpperBound()); // xtmp = u
175 xtmp->axpy(-one,*state_->iterateVec); // xtmp = u - x
176 bnd.pruneUpperInactive(*xtmp,*xlam,neps_); // xtmp = A(u - x)
177 As->plus(*xtmp); // As += A(u - x)
178
179 xtmp->set(*bnd.getLowerBound()); // xtmp = l
180 xtmp->axpy(-one,*state_->iterateVec); // xtmp = l - x
181 bnd.pruneLowerInactive(*xtmp,*xlam,neps_); // xtmp = A(l - x)
182 As->plus(*xtmp); // As += A(l - x)
183 /********************************************************************/
184 // APPLY HESSIAN TO ACTIVE COMPONENTS OF s AND REMOVE INACTIVE
185 /********************************************************************/
186 if ( useSecantHessVec_ && secant_ != nullPtr ) { // gtmp = H*As
187 secant_->applyB(*gtmp,*As);
188 }
189 else {
190 obj.hessVec(*gtmp,*As,*state_->iterateVec,itol_);
191 }
192 if (hasPoly_) {
193 proj_->getLinearConstraint()->applyJacobian(*b,*As,*state_->iterateVec,tol);
194 b->plus(*r);
195 b->scale(-one);
196 }
197 gtmp->plus(*state_->gradientVec); // gtmp = g + H*As + ...
198 gtmp->scale(-one); // gtmp = -(g + H*As + ...)
199 bnd.pruneActive(*gtmp,*xlam,neps_); // gtmp = -I(g + H*As + ...)
200 /********************************************************************/
201 // SOLVE REDUCED NEWTON SYSTEM
202 /********************************************************************/
203 if (hasPoly_) {
204 rnorm = std::sqrt(gtmp->dot(*gtmp)+b->dot(*b));
205 }
206 else {
207 rnorm = gtmp->norm();
208 }
209 if ( rnorm > zero ) {
210 if (hasPoly_) {
211 // Initialize Hessian and preconditioner
212 hessian = makePtr<HessianPDAS_Poly>(makePtrFromRef(obj),makePtrFromRef(bnd),
213 proj_->getLinearConstraint(),state_->iterateVec,xlam,neps_,secant_,
214 useSecantHessVec_,pwa,dwa);
215 precond = makePtr<PrecondPDAS_Poly>(makePtrFromRef(obj),makePtrFromRef(bnd),
216 state_->iterateVec,xlam,neps_,secant_,useSecantPrecond_,dwa);
217 PartitionedVector<Real> rhs(std::vector<Ptr<Vector<Real>>>({gtmp,b}));
218 PartitionedVector<Real> sol(std::vector<Ptr<Vector<Real>>>({state_->stepVec,mu}));
219 krylov_->run(sol,*hessian,rhs,*precond,iterKrylov_,flagKrylov_);
220 }
221 else {
222 // Initialize Hessian and preconditioner
223 hessian = makePtr<HessianPDAS>(makePtrFromRef(obj),makePtrFromRef(bnd),
224 state_->iterateVec,xlam,neps_,secant_,useSecantHessVec_,pwa);
225 precond = makePtr<PrecondPDAS>(makePtrFromRef(obj),makePtrFromRef(bnd),
226 state_->iterateVec,xlam,neps_,secant_,useSecantPrecond_,dwa);
227 krylov_->run(*state_->stepVec,*hessian,*gtmp,*precond,iterKrylov_,flagKrylov_);
228 }
229 totalKrylov_ += iterKrylov_;
230 bnd.pruneActive(*state_->stepVec,*xlam,neps_); // s <- Is
231 }
232 state_->stepVec->plus(*As); // s = Is + As
233 /********************************************************************/
234 // UPDATE STEP AND MULTIPLIER
235 /********************************************************************/
236 x.set(*state_->iterateVec);
237 x.plus(*state_->stepVec);
238 snorm = state_->stepVec->norm();
239 // Compute gradient of Lagrangian for QP
240 if ( useSecantHessVec_ && secant_ != nullPtr ) {
241 secant_->applyB(*gtmp,*state_->stepVec);
242 }
243 else {
244 obj.hessVec(*gtmp,*state_->stepVec,*state_->iterateVec,itol_);
245 }
246 if (hasPoly_) {
247 proj_->getLinearConstraint()->applyAdjointJacobian(*dwa,*mu,*state_->iterateVec,tol);
248 gtmp->plus(*dwa);
249 }
250 gtmp->plus(*state_->gradientVec);
251 gtmp->scale(-one);
252 lambda->set(gtmp->dual());
253 // Compute criticality measure
254 xtmp->set(x);
255 xtmp->plus(*lambda);
256 bnd.project(*xtmp);
257 xtmp->axpy(-one,x);
258 // Update multiplier
259 bnd.pruneInactive(*lambda,*xlam,neps_);
260 // Check stopping conditions
261 if ( xtmp->norm() < gtol_*state_->gnorm ) {
262 flag_ = 0;
263 break;
264 }
265 if ( snorm < stol_*xnorm ) {
266 flag_ = 2;
267 break;
268 }
269 }
270 if ( iter_ == maxit_ ) {
271 flag_ = 1;
272 }
273 else {
274 iter_++;
275 }
276
277 // Update algorithm state
278 state_->iter++;
279 state_->iterateVec->set(x);
280 feasible_ = bnd.isFeasible(x);
281 state_->snorm = snorm;
282 obj.update(x,UpdateType::Accept,state_->iter);
283 state_->value = obj.value(x,tol); state_->nfval++;
284
285 if ( secant_ != nullPtr ) {
286 gtmp->set(*state_->gradientVec);
287 }
288 obj.gradient(*state_->gradientVec,x,tol); state_->ngrad++;
289 xtmp->set(x); xtmp->axpy(-one,state_->gradientVec->dual());
290 proj_->project(*xtmp,outStream);
291 xtmp->axpy(-one,x);
292 state_->gnorm = xtmp->norm();
293 if ( secant_ != nullPtr ) {
294 secant_->updateStorage(x,*state_->gradientVec,*gtmp,*state_->stepVec,state_->snorm,state_->iter+1);
295 }
296
297 // Update Output
298 if (verbosity_ > 0) writeOutput(outStream,writeHeader_);
299 }
300 if (verbosity_ > 0) TypeB::Algorithm<Real>::writeExitStatus(outStream);
301}
302
303template<typename Real>
305 std::stringstream hist;
306 if (verbosity_ > 1) {
307 hist << std::string(114,'-') << std::endl;
308 if (!useSecantHessVec_) {
309 hist << "Primal Dual Active Set Newton's Method";
310 }
311 else {
312 hist << "Primal Dual Active Set Quasi-Newton Method with " << secantName_ << " Hessian approximation";
313 }
314 hist << " status output definitions" << std::endl << std::endl;
315 hist << " iter - Number of iterates (steps taken)" << std::endl;
316 hist << " value - Objective function value" << std::endl;
317 hist << " gnorm - Norm of the gradient" << std::endl;
318 hist << " snorm - Norm of the step (update to optimization vector)" << std::endl;
319 hist << " #fval - Cumulative number of times the objective function was evaluated" << std::endl;
320 hist << " #grad - Cumulative number of times the gradient was computed" << std::endl;
321 if (maxit_ > 1) {
322 hist << " iterPDAS - Number of Primal Dual Active Set iterations" << std::endl << std::endl;
323 hist << " flagPDAS - Primal Dual Active Set flag" << std::endl;
324 hist << " iterK - Number of Krylov iterations" << std::endl << std::endl;
325 }
326 else {
327 hist << " iterK - Number of Krylov iterations" << std::endl << std::endl;
328 hist << " flagK - Krylov flag" << std::endl;
329 for( int flag = CG_FLAG_SUCCESS; flag != CG_FLAG_UNDEFINED; ++flag ) {
330 hist << " " << NumberToString(flag) << " - "
331 << ECGFlagToString(static_cast<ECGFlag>(flag)) << std::endl;
332 }
333 }
334 hist << " feasible - Is iterate feasible?" << std::endl;
335 hist << std::string(114,'-') << std::endl;
336 }
337
338 hist << " ";
339 hist << std::setw(6) << std::left << "iter";
340 hist << std::setw(15) << std::left << "value";
341 hist << std::setw(15) << std::left << "gnorm";
342 hist << std::setw(15) << std::left << "snorm";
343 hist << std::setw(10) << std::left << "#fval";
344 hist << std::setw(10) << std::left << "#grad";
345 if (maxit_ > 1) {
346 hist << std::setw(10) << std::left << "iterPDAS";
347 hist << std::setw(10) << std::left << "flagPDAS";
348 hist << std::setw(10) << std::left << "iterK";
349 }
350 else {
351 hist << std::setw(10) << std::left << "iterK";
352 hist << std::setw(10) << std::left << "flagK";
353 }
354 hist << std::setw(10) << std::left << "feasible";
355 hist << std::endl;
356 os << hist.str();
357}
358
359template<typename Real>
360void PrimalDualActiveSetAlgorithm<Real>::writeName( std::ostream& os ) const {
361 std::stringstream hist;
362 if (!useSecantHessVec_) {
363 hist << std::endl << "Primal Dual Active Set Newton's Method (Type B, Bound Constraints)" << std::endl;
364 }
365 else {
366 hist << std::endl << "Primal Dual Active Set Quasi-Newton Method with "
367 << secantName_ << " Hessian approximation" << std::endl;
368 }
369 os << hist.str();
370}
371
372template<typename Real>
373void PrimalDualActiveSetAlgorithm<Real>::writeOutput( std::ostream& os, bool write_header ) const {
374 std::stringstream hist;
375 hist << std::scientific << std::setprecision(6);
376 if ( state_->iter == 0 ) writeName(os);
377 if ( write_header ) writeHeader(os);
378 if ( state_->iter == 0 ) {
379 hist << " ";
380 hist << std::setw(6) << std::left << state_->iter;
381 hist << std::setw(15) << std::left << state_->value;
382 hist << std::setw(15) << std::left << state_->gnorm;
383 hist << std::setw(15) << std::left << "---";
384 hist << std::setw(10) << std::left << state_->nfval;
385 hist << std::setw(10) << std::left << state_->ngrad;
386 hist << std::setw(10) << std::left << "---";
387 hist << std::setw(10) << std::left << "---";
388 if (maxit_ > 1) {
389 hist << std::setw(10) << std::left << "---";
390 }
391 if ( feasible_ ) {
392 hist << std::setw(10) << std::left << "YES";
393 }
394 else {
395 hist << std::setw(10) << std::left << "NO";
396 }
397 hist << std::endl;
398 }
399 else {
400 hist << " ";
401 hist << std::setw(6) << std::left << state_->iter;
402 hist << std::setw(15) << std::left << state_->value;
403 hist << std::setw(15) << std::left << state_->gnorm;
404 hist << std::setw(15) << std::left << state_->snorm;
405 hist << std::setw(10) << std::left << state_->nfval;
406 hist << std::setw(10) << std::left << state_->ngrad;
407 if (maxit_ > 1) {
408 hist << std::setw(10) << std::left << iter_;
409 hist << std::setw(10) << std::left << flag_;
410 hist << std::setw(10) << std::left << totalKrylov_;
411 }
412 else {
413 hist << std::setw(10) << std::left << iterKrylov_;
414 hist << std::setw(10) << std::left << flagKrylov_;
415 }
416 if ( feasible_ ) {
417 hist << std::setw(10) << std::left << "YES";
418 }
419 else {
420 hist << std::setw(10) << std::left << "NO";
421 }
422 hist << std::endl;
423 }
424 os << hist.str();
425}
426
427} // namespace TypeB
428} // namespace ROL
429
430#endif
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Provides the interface to apply upper and lower bound constraints.
virtual const Ptr< const Vector< Real > > getLowerBound(void) const
Return the ref count pointer to the lower bound vector.
virtual bool isFeasible(const Vector< Real > &v)
Check if the vector, v, is feasible.
void pruneLowerInactive(Vector< Real > &v, const Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the -inactive set.
void pruneInactive(Vector< Real > &v, const Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the -inactive set.
void pruneUpperInactive(Vector< Real > &v, const Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the -inactive set.
void pruneActive(Vector< Real > &v, const Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the -active set.
virtual void project(Vector< Real > &x)
Project optimization variables onto the bounds.
virtual const Ptr< const Vector< Real > > getUpperBound(void) const
Return the ref count pointer to the upper bound vector.
Provides the interface to evaluate objective functions.
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
virtual Real value(const Vector< Real > &x, Real &tol)=0
Compute value.
virtual void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update objective function.
Defines the linear algebra of vector space on a generic partitioned vector.
Provides interface for and implements limited-memory secant operators.
Definition: ROL_Secant.hpp:79
Provides an interface to check status of optimization algorithms.
void initialize(const Vector< Real > &x, const Vector< Real > &g)
virtual void writeExitStatus(std::ostream &os) const
const Ptr< CombinedStatusTest< Real > > status_
Real stol_
PDAS minimum step size stopping tolerance (default: 1e-8)
void writeName(std::ostream &os) const override
Print step name.
Real scale_
Scale for dual variables in the active set, (default: 1)
PrimalDualActiveSetAlgorithm(ParameterList &list, const Ptr< Secant< Real > > &secant=nullPtr)
Ptr< Secant< Real > > secant_
Secant object (used for quasi-Newton)
Real atolKrylov_
Absolute tolerance for Krylov solve (default: 1e-4)
Real gtol_
PDAS gradient stopping tolerance (default: 1e-6)
void initialize(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &bnd, std::ostream &outStream=std::cout)
void writeOutput(std::ostream &os, bool write_header=false) const override
Print iterate status.
void run(Vector< Real > &x, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &bnd, std::ostream &outStream=std::cout) override
Run algorithm on bound constrained problems (Type-B). This general interface supports the use of dual...
int maxit_
Maximum number of PDAS steps (default: 10)
void writeHeader(std::ostream &os) const override
Print iterate header.
int maxitKrylov_
Maximum number of Krylov iterations (default: 100)
bool useSecantHessVec_
Whether or not to use to a secant approximation as the Hessian.
Real rtolKrylov_
Relative tolerance for Krylov solve (default: 1e-2)
bool useSecantPrecond_
Whether or not to use a secant approximation to precondition inexact Newton.
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 plus(const Vector &x)=0
Compute , where .
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
EKrylov StringToEKrylov(std::string s)
std::string NumberToString(T Number)
Definition: ROL_Types.hpp:81
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:91
ESecant StringToESecant(std::string s)
Definition: ROL_Types.hpp:543
@ SECANT_USERDEFINED
Definition: ROL_Types.hpp:491
@ CG_FLAG_UNDEFINED
Definition: ROL_Types.hpp:827
@ CG_FLAG_SUCCESS
Definition: ROL_Types.hpp:822
std::string ECGFlagToString(ECGFlag cgf)
Definition: ROL_Types.hpp:831