ROL
ROL_PrimalDualInteriorPointStep.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
45#ifndef ROL_PRIMALDUALINTERIORPOINTSTEP_H
46#define ROL_PRIMALDUALINTERIORPOINTSTEP_H
47
48#include "ROL_VectorsNorms.hpp"
52#include "ROL_Krylov.hpp"
53
54
112namespace ROL {
113namespace InteriorPoint {
114template <class Real>
115class PrimalDualInteriorPointStep : public Step<Real> {
116
117 typedef Vector<Real> V;
118 typedef PartitionedVector<Real> PV;
119 typedef Objective<Real> OBJ;
120 typedef BoundConstraint<Real> BND;
121 typedef Krylov<Real> KRYLOV;
122 typedef LinearOperator<Real> LINOP;
123 typedef LinearOperatorFromConstraint<Real> LOPEC;
124 typedef Constraint<Real> EQCON;
125 typedef StepState<Real> STATE;
126 typedef InteriorPointPenalty<Real> PENALTY;
127 typedef PrimalDualInteriorPointResidual<Real> RESIDUAL;
128
129private:
130
131 ROL::Ptr<KRYLOV> krylov_; // Krylov solver for the Primal Dual system
132 ROL::Ptr<LINOP> precond_; // Preconditioner for the Primal Dual system
133
134 ROL::Ptr<BND> pbnd_; // bound constraint for projecting x sufficiently away from the given bounds
135
136 ROL::Ptr<V> x_; // Optimization vector
137 ROL::Ptr<V> g_; // Gradient of the Lagrangian
138 ROL::Ptr<V> l_; // Lagrange multiplier
139
140 ROL::Ptr<V> xl_; // Lower bound vector
141 ROL::Ptr<V> xu_; // Upper bound vector
142
143 ROL::Ptr<V> zl_; // Lagrange multiplier for lower bound
144 ROL::Ptr<V> zu_; // Lagrange multiplier for upper bound
145
146 ROL::Ptr<V> xscratch_; // Scratch vector (size of x)
147 ROL::Ptr<V> lscratch_; // Scratch vector (size of l)
148
149 ROL::Ptr<V> zlscratch_; // Scratch vector (size of x)
150 ROL::Ptr<V> zuscratch_; // Scratch vector (size of x)
151
152 ROL::Ptr<V> maskL_; // Elements are 1 when xl>-INF, zero for xl = -INF
153 ROL::Ptr<V> maskU_; // Elements are 1 when xu< INF, zero for xu = INF
154
155 int iterKrylov_;
156 int flagKrylov_;
157
158 bool symmetrize_; // Symmetrize the Primal Dual system if true
159
160 Elementwise::Multiply<Real> mult_;
161
162 // Parameters used by the Primal-Dual Interior Point Method
163
164 Real mu_; // Barrier penalty parameter
165
166 Real eps_tol_; // Error function tolerance
167 Real tau_; // Current fraction-to-the-boundary parameter
168 Real tau_min_; // Minimum fraction-to-the-boundary parameter
169 Real kappa_eps_;
170 Real kappa_mu_;
171 Real kappa1_; // Feasibility projection parameter
172 Real kappa2_; // Feasibility projection parameter
173 Real kappa_eps_; //
174 Real lambda_max_; // multiplier maximum value
175 Real theta_mu_;
176 Real gamma_theta_;
177 Real gamma_phi_
178 Real sd_; // Lagragian gradient scaling parameter
179 Real scl_; // Lower bound complementarity scaling parameter
180 Real scu_; // Upper bound complementarity scaling parameter
181 Real smax_; // Maximum scaling parameter
182
183 Real diml_; // Dimension of constaint
184 Real dimx_; // Dimension of optimization vector
185
186
187
188 void updateState( const V& x, const V &l, OBJ &obj,
189 EQCON &con, BND &bnd, ALGO &algo_state ) {
190
191 Real tol = std::sqrt(ROL_EPSILON<Real>());
192 ROL::Ptr<STATE> state = Step<Real>::getState();
193
194 obj.update(x,true,algo_state.iter);
195 con.update(x,true,algo_state.iter);
196
197 algo_state.value = obj.value(tol);
198 con.value(*(state->constraintVec),x,tol);
199
200 obj.gradient(*(state->gradientVec),x,tol);
201 con.applyAdjointJacobian(*g_,l,x,tol);
202
203 state->gradientVec->plus(*g_); // \f$ \nabla f(x)-\nabla c(x)\lambda \f$
204
205 state->gradientVec->axpy(-1.0,*zl_);
206 state->gradientVec->axpy(-1.0,*zu_);
207
208 // Scaled Lagrangian gradient sup norm
209 algo_state.gnorm = normLinf(*(state->gradientVec))/sd_;
210
211 // Constraint sup norm
212 algo_state.cnorm = normLinf(*(state->constraintVec));
213
214 Elementwise::Multiply<Real> mult;
215
216 Real lowerViolation;
217 Real upperViolation;
218
219 // Deviation from complementarity
220 xscratch_->set(x);
221 xscratch_->applyBinary(mult,*zl_);
222
223 exactLowerViolation = normLinf(*xscratch_)/scl_;
224
225 xscratch_->set(x);
226 xscratch_->applyBunary(mult,*zu_);
227
228 exactUpperBound = normLinf(*xscratch_)/scu_;
229
230 // Measure ||xz||
231 algo_state.aggregateModelError = std::max(exactLowerViolation,
232 exactUpperViolation);
233
234 }
235
236 /* When the constraint Jacobians are ill-conditioned, we can compute
237 multiplier vectors with very large norms, making it difficult to
238 satisfy the primal-dual equations to a small tolerance. These
239 parameters allow us to rescale the constraint qualification component
240 of the convergence criteria
241 */
242 void updateScalingParameters(void) {
243
244 Real nl = normL1(*l_);
245 Real nzl = normL1(*zl_);
246 Real nzu = normL1(*zu_);
247
248 sd_ = (nl+nzl+nzu)/(diml_+2*dimx);
249 sd_ = std::max(smax_,sd_)
250
251 sd_ /= smax_;
252
253 scl_ = nzl/dimx_;
254 scl_ = std::max(smax_,scl_);
255 scl_ /= smax_;
256
257 scu_ = nzu/dimx_;
258 scu_ = std::max(smax_,scu_);
259 scu_ /= smax_;
260 }
261
262public:
263
264 using Step<Real>::initialize;
265 using Step<Real>::compute;
266 using Step<Real>::update;
267
268 PrimalDualInteriorPointStep( ROL::ParameterList &parlist,
269 const ROL::Ptr<Krylov<Real> > &krylov = ROL::nullPtr,
270 const ROL::Ptr<LinearOperator<Real> > &precond = ROL::nullPtr ) :
271 Step<Real>(), krylov_(krylov), precond_(precond), iterKrylov_(0), flagKrylov_(0) {
272
273 typedef ROL::ParameterList PL;
274
275 PL &iplist = parlist.sublist("Step").sublist("Primal Dual Interior Point");
276
277 kappa1_ = iplist.get("Bound Perturbation Coefficient 1", 1.e-2);
278 kappa2_ = iplist.get("Bound Perturbation Coefficient 2", 1.e-2);
279 lambda_max_ = iplist.get(" Multiplier Maximum Value", 1.e3 );
280 smax_ = iplist.get("Maximum Scaling Parameter", 1.e2 );
281 tau_min_ = iplist.get("Minimum Fraction-to-Boundary Parameter", 0.99 );
282 kappa_mu_ = iplist.get("Multiplicative Penalty Reduction Factor", 0.2 );
283 theta_mu_ = iplist.get("Penalty Update Power", 1.5 );
284 eps_tol_ = iplist.get("Error Tolerance", 1.e-8);
285 symmetrize_ = iplist.get("Symmetrize Primal Dual System", true );
286
287 PL &filter = iplist.sublist("Filter Parameters");
288
289 if(krylov_ == ROL::nullPtr) {
290 krylov_ = KrylovFactory<Real>(parlist);
291 }
292
293 if( precond_ == ROL::nullPtr) {
294 class IdentityOperator : public LINOP {
295 public:
296 apply( V& Hv, const V &v, Real tol ) const {
297 Hv.set(v);
298 }
299 }; // class IdentityOperator
300
301 precond_ = ROL::makePtr<IdentityOperator<Real>>();
302 }
303
304 }
305
306
307
308 ~PrimalDualInteriorPointStep() {
309
310
311 void initialize( Vector<Real> &x, const Vector<Real> &g, Vector<Real> &l, const Vector<Real> &c,
312 Objective<Real> &obj, Constraint<Real> &con, BoundConstraint<Real> &bnd,
313 AlgorithmState<Real> &algo_state ) {
314
315
316 using Elementwise::ValueSet;
317
318 ROL::Ptr<PENALTY> &ipPen = dynamic_cast<PENALTY&>(obj);
319
320 // Initialize step state
321 ROL::Ptr<STATE> state = Step<Real>::getState();
322 state->descentVec = x.clone();
323 state->gradientVec = g.clone();
324 state->constaintVec = c.clone();
325
326 diml_ = l.dimension();
327 dimx_ = x.dimension();
328
329 Real one(1.0);
330 Real zero(0.0);
331 Real tol = std::sqrt(ROL_EPSILON<Real>());
332
333 x_ = x.clone();
334 g_ = g.clone();
335 l_ = l.clone();
336 c_ = c.clone();
337
338 xscratch_ = x.clone();
339 lscratch_ = l.clone();
340
341 zlscratch_ = x.clone();
342 zuscratch_ = x.clone();
343
344 // Multipliers for lower and upper bounds
345 zl_ = x.clone();
346 zu_ = x.clone();
347
348 /*******************************************************************************************/
349 /* Identify (implicitly) the index sets of upper and lower bounds by creating mask vectors */
350 /*******************************************************************************************/
351
352 xl_ = bnd.getLowerBound();
353 xu_ = bnd.getUpperBound();
354
355 maskl_ = ipPen.getLowerMask();
356 masku_ = ipPen.getUpperMask();
357
358 // Initialize bound constraint multipliers to 1 one where the corresponding bounds are finite
359 zl_->set(maskl_);
360 zu_->set(masku_);
361
362 /*******************************************************************************************/
363 /* Create a new bound constraint with perturbed bounds */
364 /*******************************************************************************************/
365
366 ROL::Ptr<V> xdiff = xu_->clone();
367 xdiff->set(*xu_);
368 xdiff->axpy(-1.0,*xl_);
369 xdiff->scale(kappa2_);
370
371 class Max1X : public Elementwise::UnaryFunction<Real> {
372 public:
373 Real apply( const Real &x ) const {
374 return std::max(1.0,x);
375 }
376 };
377
378 Max1X max1x;
379 Elementwise::AbsoluteValue<Real> absval;
380 Elementwise::Min min;
381
382 // Lower perturbation vector
383 ROL::Ptr<V> pl = xl_->clone();
384 pl->applyUnary(absval);
385 pl->applyUnary(max1x); // pl_i = max(1,|xl_i|)
386 pl->scale(kappa1_);
387 pl->applyBinary(min,xdiff); // pl_i = min(kappa1*max(1,|xl_i|),kappa2*(xu_i-xl_i))
388
389 // Upper perturbation vector
390 ROL::Ptr<V> pu = xu_->clone();
391 pu->applyUnary(absval);
392 pu->applyUnary(max1x); // pu_i = max(1,|xu_i|)
393 pu->scale(kappa_1);
394 pu->applyBinary(min,xdiff); // pu_u = min(kappa1*max(1,|xu_i|,kappa2*(xu_i-xl_i)))
395
396 // Modified lower and upper bounds so that x in [xl+pl,xu-pu] using the above perturbation vectors
397 pl->plus(*xl_);
398 pu->scale(-1.0);
399 pu->plus(*xu_);
400
401 pbnd_ = ROL::makePtr<BoundConstraint<Real>>(pl,pu)
402
403 // Project the initial guess onto the perturbed bounds
404 pbnd_->project(x);
405
406
407 /*******************************************************************************************/
408 /* Solve least-squares problem for initial equality multiplier */
409 /*******************************************************************************************/
410
411 // \f$-(\nabla f-z_l+z_u) \f$
412 g_->set(*zl_);
413 g_->axpy(-1.0,g);
414 g_->axpy(-1.0,*zu_)
415
416 // We only need the updated multiplier
417 lscratch_->zero();
418 con.solveAugmentedSystem(*xscratch_,*l_,*g_,*lscratch_,x,tol);
419
420 // If the multiplier supremum is too large, zero the vector as this appears to avoid poor
421 // initial guesses of the multiplier in the case where the constraint Jacobian is
422 // ill-conditioned
423 if( normInf(l_) > lambda_max_ ) {
424 l_->zero();
425 }
426
427 // Initialize the algorithm state
428 algo_state.nfval = 0;
429 algo_state.ncval = 0;
430 algo_state.ngrad = 0;
431
432 updateState(x,l,obj,con,bnd,algo_state);
433
434 } // initialize()
435
436
437 void compute( Vector<Real> &s, const Vector<Real> &x, const Vector<Real> &l,
438 Objective<Real> &obj, Constraint<Real> &con,
439 BoundConstraint<Real> &bnd, AlgorithmState<Real> &algo_state ) {
440
441
442
443
444 Elementwise::Fill<Real> minus_mu(-mu_);
445 Elementwise::Divide<Real> div;
446 Elementwise::Multiply<Real> mult;
447
448 ROL::Ptr<STATE> state = Step<Real>::getState();
449
450 ROL::Ptr<OBJ> obj_ptr = ROL::makePtrFromRef(obj);
451 ROL::Ptr<EQCON> con_ptr = ROL::makePtrFromRef(con);
452 ROL::Ptr<BND> bnd_ptr = ROL::makePtrFromRef(bnd);
453
454
455 /*******************************************************************************************/
456 /* Form Primal-Dual system residual and operator then solve for direction vector */
457 /*******************************************************************************************/
458
459
460 ROL::Ptr<V> rhs = CreatePartitionedVector(state->gradientVec,
461 state->constraintVec,
462 resL_,
463 resU_);
464
465 ROL::Ptr<V> sysvec = CreatePartitionedVector( x_, l_, zl_, zu_ );
466
467
468 ROL::Ptr<RESIDUAL> residual = ROL::makePtr<RESIDUAL>(obj,con,bnd,*sol,maskL_,maskU_,w_,mu_,symmetrize_);
469
470 residual->value(*rhs,*sysvec,tol);
471
472 ROL::Ptr<V> sol = CreatePartitionedVector( xscratch_, lscratch_, zlscratch_, zuscratch_ );
473
474 LOPEC jacobian( sysvec, residual );
475
476
477 krylov_->run(*sol,jacobian,*residual,*precond_,iterKrylov_,flagKrylov_);
478
479 /*******************************************************************************************/
480 /* Perform line search */
481 /*******************************************************************************************/
482
483
484
485 } // compute()
486
487
488
489 void update( Vector<Real> &x, Vector<Real> &l, const Vector<Real> &s,
490 Objective<Real> &obj, Constraint<Real> &con,
491 BoundConstraint<Real> &bnd, AlgorithmState<Real> &algo_state ) {
492
493 // Check deviation from shifted complementarity
494 Elementwise::Shift<Real> minus_mu(-mu_);
495
496 xscratch_->set(x);
497 xscratch_->applyBinary(mult,*zl_);
498 xscratch_->applyUnary(minus_mu);
499
500 lowerViolation = normLinf(*xscratch_)/scl_; // \f$ \max_i xz_l^i - \mu \f$
501
502 xscratch_->set(x);
503 xscratch_->applyBinary(mult,*zu_);
504 xscratch_->applyUnary(minus_mu);
505
506 upperBound = normLinf(*xscratch_)/scu_;
507
508 // Evaluate \f$E_\mu(x,\lambda,z_l,z_u)\f$
509 Real Emu = algo_state.gnorm;
510 Emu = std::max(Emu,algo_state.cnorm);
511 Emu = std::max(Emu,upperBound);
512 Emu = std::max(Emu,lowerBound);
513
514 // If sufficiently converged for the current mu, update it
515 if(Emu < (kappa_epsilon_*mu_) ) {
516 Real mu_old = mu_;
517
518 /* \mu_{j+1} = \max\left{ \frac{\epsilon_\text{tol}}{10},
519 \min\{ \kappa_{\mu} \mu_j,
520 \mu_j^{\theta_\mu}\} \right\} */
521 mu_ = std::min(kappa_mu_*mu_old,std::pow(mu_old,theta_mu_));
522 mu_ = std::max(eps_tol_/10.0,mu_);
523
524 // Update fraction-to-boundary parameter
525 tau_ = std::max(tau_min_,1.0-mu_);
526
527
528
529 }
530
531 } // update()
532
533
534
535
536 // TODO: Implement header print out
537 std::string printHeader( void ) const {
538 std::string head("");
539 return head;
540 }
541
542 // TODO: Implement name print out
543 std::string printName( void ) const {
544 std::string name("");
545 return name;
546 }
547
548
549}; // class PrimalDualInteriorPointStep
550
551} // namespace InteriorPoint
552} // namespace ROL
553
554
555#endif // ROL_PRIMALDUALINTERIORPOINTSTEP_H
Vector< Real > V
PartitionedVector< Real > PV
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
ROL::BlockOperator2Diagonal BlockOperator2 apply(V &Hv, const V &v, Real &tol) const
Real normL1(const Vector< Real > &x)
ROL::Ptr< Vector< Real > > CreatePartitionedVector(const ROL::Ptr< Vector< Real > > &a)
virtual void update(const Vector< Real > &u, const Vector< Real > &z, bool flag=true, int iter=-1) override
Real normLinf(const Vector< Real > &x)