ROL
ROL_PrimalDualActiveSetStep.hpp
Go to the documentation of this file.
1 // @HEADER
2 // *****************************************************************************
3 // Rapid Optimization Library (ROL) Package
4 //
5 // Copyright 2014 NTESS and the ROL contributors.
6 // SPDX-License-Identifier: BSD-3-Clause
7 // *****************************************************************************
8 // @HEADER
9 #ifndef ROL_PRIMALDUALACTIVESETSTEP_H
10 #define ROL_PRIMALDUALACTIVESETSTEP_H
11 
12 #include "ROL_Step.hpp"
13 #include "ROL_Vector.hpp"
14 #include "ROL_KrylovFactory.hpp"
15 #include "ROL_Objective.hpp"
16 #include "ROL_BoundConstraint.hpp"
17 #include "ROL_Types.hpp"
18 #include "ROL_Secant.hpp"
19 #include "ROL_ParameterList.hpp"
20 
95 namespace ROL {
96 
97 template <class Real>
98 class PrimalDualActiveSetStep : public Step<Real> {
99 private:
100 
101  ROL::Ptr<Krylov<Real> > krylov_;
102 
103  // Krylov Parameters
104  int iterCR_;
105  int flagCR_;
106  Real itol_;
107 
108  // PDAS Parameters
109  int maxit_;
110  int iter_;
111  int flag_;
112  Real stol_;
113  Real gtol_;
114  Real scale_;
115  Real neps_;
116  bool feasible_;
117 
118  // Dual Variable
119  ROL::Ptr<Vector<Real> > lambda_;
120  ROL::Ptr<Vector<Real> > xlam_;
121  ROL::Ptr<Vector<Real> > x0_;
122  ROL::Ptr<Vector<Real> > xbnd_;
123  ROL::Ptr<Vector<Real> > As_;
124  ROL::Ptr<Vector<Real> > xtmp_;
125  ROL::Ptr<Vector<Real> > res_;
126  ROL::Ptr<Vector<Real> > Ag_;
127  ROL::Ptr<Vector<Real> > rtmp_;
128  ROL::Ptr<Vector<Real> > gtmp_;
129 
130  // Secant Information
131  ESecant esec_;
132  ROL::Ptr<Secant<Real> > secant_;
133  bool useSecantPrecond_;
134  bool useSecantHessVec_;
135 
136  class HessianPD : public LinearOperator<Real> {
137  private:
138  const ROL::Ptr<Objective<Real> > obj_;
139  const ROL::Ptr<BoundConstraint<Real> > bnd_;
140  const ROL::Ptr<Vector<Real> > x_;
141  const ROL::Ptr<Vector<Real> > xlam_;
142  ROL::Ptr<Vector<Real> > v_;
143  Real eps_;
144  const ROL::Ptr<Secant<Real> > secant_;
145  bool useSecant_;
146  public:
147  HessianPD(const ROL::Ptr<Objective<Real> > &obj,
148  const ROL::Ptr<BoundConstraint<Real> > &bnd,
149  const ROL::Ptr<Vector<Real> > &x,
150  const ROL::Ptr<Vector<Real> > &xlam,
151  const Real eps = 0,
152  const ROL::Ptr<Secant<Real> > &secant = ROL::nullPtr,
153  const bool useSecant = false )
154  : obj_(obj), bnd_(bnd), x_(x), xlam_(xlam),
155  eps_(eps), secant_(secant), useSecant_(useSecant) {
156  v_ = x_->clone();
157  if ( !useSecant || secant == ROL::nullPtr ) {
158  useSecant_ = false;
159  }
160  }
161  void apply( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
162  v_->set(v);
163  bnd_->pruneActive(*v_,*xlam_,eps_);
164  if ( useSecant_ ) {
165  secant_->applyB(Hv,*v_);
166  }
167  else {
168  obj_->hessVec(Hv,*v_,*x_,tol);
169  }
170  bnd_->pruneActive(Hv,*xlam_,eps_);
171  }
172  };
173 
174  class PrecondPD : public LinearOperator<Real> {
175  private:
176  const ROL::Ptr<Objective<Real> > obj_;
177  const ROL::Ptr<BoundConstraint<Real> > bnd_;
178  const ROL::Ptr<Vector<Real> > x_;
179  const ROL::Ptr<Vector<Real> > xlam_;
180  ROL::Ptr<Vector<Real> > v_;
181  Real eps_;
182  const ROL::Ptr<Secant<Real> > secant_;
183  bool useSecant_;
184  public:
185  PrecondPD(const ROL::Ptr<Objective<Real> > &obj,
186  const ROL::Ptr<BoundConstraint<Real> > &bnd,
187  const ROL::Ptr<Vector<Real> > &x,
188  const ROL::Ptr<Vector<Real> > &xlam,
189  const Real eps = 0,
190  const ROL::Ptr<Secant<Real> > &secant = ROL::nullPtr,
191  const bool useSecant = false )
192  : obj_(obj), bnd_(bnd), x_(x), xlam_(xlam),
193  eps_(eps), secant_(secant), useSecant_(useSecant) {
194  v_ = x_->dual().clone();
195  if ( !useSecant || secant == ROL::nullPtr ) {
196  useSecant_ = false;
197  }
198  }
199  void apply( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
200  Hv.set(v.dual());
201  }
202  void applyInverse( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
203  v_->set(v);
204  bnd_->pruneActive(*v_,*xlam_,eps_);
205  if ( useSecant_ ) {
206  secant_->applyH(Hv,*v_);
207  }
208  else {
209  obj_->precond(Hv,*v_,*x_,tol);
210  }
211  bnd_->pruneActive(Hv,*xlam_,eps_);
212  }
213  };
214 
227  Real computeCriticalityMeasure(Vector<Real> &x, Objective<Real> &obj, BoundConstraint<Real> &con, Real tol) {
228  Real one(1);
229  ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState();
230  obj.gradient(*(step_state->gradientVec),x,tol);
231  xtmp_->set(x);
232  xtmp_->axpy(-one,(step_state->gradientVec)->dual());
233  con.project(*xtmp_);
234  xtmp_->axpy(-one,x);
235  return xtmp_->norm();
236  }
237 
238 public:
245  PrimalDualActiveSetStep( ROL::ParameterList &parlist )
246  : Step<Real>::Step(), krylov_(ROL::nullPtr),
247  iterCR_(0), flagCR_(0), itol_(0),
248  maxit_(0), iter_(0), flag_(0), stol_(0), gtol_(0), scale_(0),
249  neps_(-ROL_EPSILON<Real>()), feasible_(false),
250  lambda_(ROL::nullPtr), xlam_(ROL::nullPtr), x0_(ROL::nullPtr),
251  xbnd_(ROL::nullPtr), As_(ROL::nullPtr), xtmp_(ROL::nullPtr),
252  res_(ROL::nullPtr), Ag_(ROL::nullPtr), rtmp_(ROL::nullPtr),
253  gtmp_(ROL::nullPtr),
254  esec_(SECANT_LBFGS), secant_(ROL::nullPtr), useSecantPrecond_(false),
255  useSecantHessVec_(false) {
256  Real one(1), oem6(1.e-6), oem8(1.e-8);
257  // Algorithmic parameters
258  maxit_ = parlist.sublist("Step").sublist("Primal Dual Active Set").get("Iteration Limit",10);
259  stol_ = parlist.sublist("Step").sublist("Primal Dual Active Set").get("Relative Step Tolerance",oem8);
260  gtol_ = parlist.sublist("Step").sublist("Primal Dual Active Set").get("Relative Gradient Tolerance",oem6);
261  scale_ = parlist.sublist("Step").sublist("Primal Dual Active Set").get("Dual Scaling", one);
262  // Build secant object
263  esec_ = StringToESecant(parlist.sublist("General").sublist("Secant").get("Type","Limited-Memory BFGS"));
264  useSecantHessVec_ = parlist.sublist("General").sublist("Secant").get("Use as Hessian", false);
265  useSecantPrecond_ = parlist.sublist("General").sublist("Secant").get("Use as Preconditioner", false);
266  if ( useSecantHessVec_ || useSecantPrecond_ ) {
267  secant_ = SecantFactory<Real>(parlist);
268  }
269  // Build Krylov object
270  krylov_ = KrylovFactory<Real>(parlist);
271  }
272 
284  using Step<Real>::initialize;
285  void initialize( Vector<Real> &x, const Vector<Real> &s, const Vector<Real> &g,
286  Objective<Real> &obj, BoundConstraint<Real> &con,
287  AlgorithmState<Real> &algo_state ) {
288  ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState();
289  Real zero(0), one(1);
290  // Initialize state descent direction and gradient storage
291  step_state->descentVec = s.clone();
292  step_state->gradientVec = g.clone();
293  step_state->searchSize = zero;
294  // Initialize additional storage
295  xlam_ = x.clone();
296  x0_ = x.clone();
297  xbnd_ = x.clone();
298  As_ = s.clone();
299  xtmp_ = x.clone();
300  res_ = g.clone();
301  Ag_ = g.clone();
302  rtmp_ = g.clone();
303  gtmp_ = g.clone();
304  // Project x onto constraint set
305  con.project(x);
306  // Update objective function, get value, and get gradient
307  Real tol = std::sqrt(ROL_EPSILON<Real>());
308  obj.update(x,true,algo_state.iter);
309  algo_state.value = obj.value(x,tol);
310  algo_state.nfval++;
311  algo_state.gnorm = computeCriticalityMeasure(x,obj,con,tol);
312  algo_state.ngrad++;
313  // Initialize dual variable
314  lambda_ = s.clone();
315  lambda_->set((step_state->gradientVec)->dual());
316  lambda_->scale(-one);
317  }
318 
344  using Step<Real>::compute;
345  void compute( Vector<Real> &s, const Vector<Real> &x, Objective<Real> &obj, BoundConstraint<Real> &con,
346  AlgorithmState<Real> &algo_state ) {
347  ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState();
348  Real zero(0), one(1);
349  s.zero();
350  x0_->set(x);
351  res_->set(*(step_state->gradientVec));
352  for ( iter_ = 0; iter_ < maxit_; iter_++ ) {
353  /********************************************************************/
354  // MODIFY ITERATE VECTOR TO CHECK ACTIVE SET
355  /********************************************************************/
356  xlam_->set(*x0_); // xlam = x0
357  xlam_->axpy(scale_,*(lambda_)); // xlam = x0 + c*lambda
358  /********************************************************************/
359  // PROJECT x ONTO PRIMAL DUAL FEASIBLE SET
360  /********************************************************************/
361  As_->zero(); // As = 0
362 
363  xbnd_->set(*con.getUpperBound()); // xbnd = u
364  xbnd_->axpy(-one,x); // xbnd = u - x
365  xtmp_->set(*xbnd_); // tmp = u - x
366  con.pruneUpperActive(*xtmp_,*xlam_,neps_); // tmp = I(u - x)
367  xbnd_->axpy(-one,*xtmp_); // xbnd = A(u - x)
368  As_->plus(*xbnd_); // As += A(u - x)
369 
370  xbnd_->set(*con.getLowerBound()); // xbnd = l
371  xbnd_->axpy(-one,x); // xbnd = l - x
372  xtmp_->set(*xbnd_); // tmp = l - x
373  con.pruneLowerActive(*xtmp_,*xlam_,neps_); // tmp = I(l - x)
374  xbnd_->axpy(-one,*xtmp_); // xbnd = A(l - x)
375  As_->plus(*xbnd_); // As += A(l - x)
376  /********************************************************************/
377  // APPLY HESSIAN TO ACTIVE COMPONENTS OF s AND REMOVE INACTIVE
378  /********************************************************************/
379  itol_ = std::sqrt(ROL_EPSILON<Real>());
380  if ( useSecantHessVec_ && secant_ != ROL::nullPtr ) { // IHAs = H*As
381  secant_->applyB(*gtmp_,*As_);
382  }
383  else {
384  obj.hessVec(*gtmp_,*As_,x,itol_);
385  }
386  con.pruneActive(*gtmp_,*xlam_,neps_); // IHAs = I(H*As)
387  /********************************************************************/
388  // SEPARATE ACTIVE AND INACTIVE COMPONENTS OF THE GRADIENT
389  /********************************************************************/
390  rtmp_->set(*(step_state->gradientVec)); // Inactive components
391  con.pruneActive(*rtmp_,*xlam_,neps_);
392 
393  Ag_->set(*(step_state->gradientVec)); // Active components
394  Ag_->axpy(-one,*rtmp_);
395  /********************************************************************/
396  // SOLVE REDUCED NEWTON SYSTEM
397  /********************************************************************/
398  rtmp_->plus(*gtmp_);
399  rtmp_->scale(-one); // rhs = -Ig - I(H*As)
400  s.zero();
401  if ( rtmp_->norm() > zero ) {
402  // Initialize Hessian and preconditioner
403  ROL::Ptr<Objective<Real> > obj_ptr = ROL::makePtrFromRef(obj);
404  ROL::Ptr<BoundConstraint<Real> > con_ptr = ROL::makePtrFromRef(con);
405  ROL::Ptr<LinearOperator<Real> > hessian
406  = ROL::makePtr<HessianPD>(obj_ptr,con_ptr,
407  algo_state.iterateVec,xlam_,neps_,secant_,useSecantHessVec_);
408  ROL::Ptr<LinearOperator<Real> > precond
409  = ROL::makePtr<PrecondPD>(obj_ptr,con_ptr,
410  algo_state.iterateVec,xlam_,neps_,secant_,useSecantPrecond_);
411  //solve(s,*rtmp_,*xlam_,x,obj,con); // Call conjugate residuals
412  krylov_->run(s,*hessian,*rtmp_,*precond,iterCR_,flagCR_);
413  con.pruneActive(s,*xlam_,neps_); // s <- Is
414  }
415  s.plus(*As_); // s = Is + As
416  /********************************************************************/
417  // UPDATE MULTIPLIER
418  /********************************************************************/
419  if ( useSecantHessVec_ && secant_ != ROL::nullPtr ) {
420  secant_->applyB(*rtmp_,s);
421  }
422  else {
423  obj.hessVec(*rtmp_,s,x,itol_);
424  }
425  gtmp_->set(*rtmp_);
426  con.pruneActive(*gtmp_,*xlam_,neps_);
427  lambda_->set(*rtmp_);
428  lambda_->axpy(-one,*gtmp_);
429  lambda_->plus(*Ag_);
430  lambda_->scale(-one);
431  /********************************************************************/
432  // UPDATE STEP
433  /********************************************************************/
434  x0_->set(x);
435  x0_->plus(s);
436  res_->set(*(step_state->gradientVec));
437  res_->plus(*rtmp_);
438  // Compute criticality measure
439  xtmp_->set(*x0_);
440  xtmp_->axpy(-one,res_->dual());
441  con.project(*xtmp_);
442  xtmp_->axpy(-one,*x0_);
443 // std::cout << s.norm() << " "
444 // << tmp->norm() << " "
445 // << res_->norm() << " "
446 // << lambda_->norm() << " "
447 // << flagCR_ << " "
448 // << iterCR_ << "\n";
449  if ( xtmp_->norm() < gtol_*algo_state.gnorm ) {
450  flag_ = 0;
451  break;
452  }
453  if ( s.norm() < stol_*x.norm() ) {
454  flag_ = 2;
455  break;
456  }
457  }
458  if ( iter_ == maxit_ ) {
459  flag_ = 1;
460  }
461  else {
462  iter_++;
463  }
464  }
465 
477  using Step<Real>::update;
478  void update( Vector<Real> &x, const Vector<Real> &s, Objective<Real> &obj, BoundConstraint<Real> &con,
479  AlgorithmState<Real> &algo_state ) {
480  ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState();
481  step_state->SPiter = (maxit_ > 1) ? iter_ : iterCR_;
482  step_state->SPflag = (maxit_ > 1) ? flag_ : flagCR_;
483 
484  x.plus(s);
485  feasible_ = con.isFeasible(x);
486  algo_state.snorm = s.norm();
487  algo_state.iter++;
488  Real tol = std::sqrt(ROL_EPSILON<Real>());
489  obj.update(x,true,algo_state.iter);
490  algo_state.value = obj.value(x,tol);
491  algo_state.nfval++;
492 
493  if ( secant_ != ROL::nullPtr ) {
494  gtmp_->set(*(step_state->gradientVec));
495  }
496  algo_state.gnorm = computeCriticalityMeasure(x,obj,con,tol);
497  algo_state.ngrad++;
498 
499  if ( secant_ != ROL::nullPtr ) {
500  secant_->updateStorage(x,*(step_state->gradientVec),*gtmp_,s,algo_state.snorm,algo_state.iter+1);
501  }
502  (algo_state.iterateVec)->set(x);
503  }
504 
510  std::string printHeader( void ) const {
511  std::stringstream hist;
512  hist << " ";
513  hist << std::setw(6) << std::left << "iter";
514  hist << std::setw(15) << std::left << "value";
515  hist << std::setw(15) << std::left << "gnorm";
516  hist << std::setw(15) << std::left << "snorm";
517  hist << std::setw(10) << std::left << "#fval";
518  hist << std::setw(10) << std::left << "#grad";
519  if ( maxit_ > 1 ) {
520  hist << std::setw(10) << std::left << "iterPDAS";
521  hist << std::setw(10) << std::left << "flagPDAS";
522  }
523  else {
524  hist << std::setw(10) << std::left << "iterCR";
525  hist << std::setw(10) << std::left << "flagCR";
526  }
527  hist << std::setw(10) << std::left << "feasible";
528  hist << "\n";
529  return hist.str();
530  }
531 
537  std::string printName( void ) const {
538  std::stringstream hist;
539  hist << "\nPrimal Dual Active Set Newton's Method\n";
540  return hist.str();
541  }
542 
550  virtual std::string print( AlgorithmState<Real> &algo_state, bool print_header = false ) const {
551  std::stringstream hist;
552  hist << std::scientific << std::setprecision(6);
553  if ( algo_state.iter == 0 ) {
554  hist << printName();
555  }
556  if ( print_header ) {
557  hist << printHeader();
558  }
559  if ( algo_state.iter == 0 ) {
560  hist << " ";
561  hist << std::setw(6) << std::left << algo_state.iter;
562  hist << std::setw(15) << std::left << algo_state.value;
563  hist << std::setw(15) << std::left << algo_state.gnorm;
564  hist << "\n";
565  }
566  else {
567  hist << " ";
568  hist << std::setw(6) << std::left << algo_state.iter;
569  hist << std::setw(15) << std::left << algo_state.value;
570  hist << std::setw(15) << std::left << algo_state.gnorm;
571  hist << std::setw(15) << std::left << algo_state.snorm;
572  hist << std::setw(10) << std::left << algo_state.nfval;
573  hist << std::setw(10) << std::left << algo_state.ngrad;
574  if ( maxit_ > 1 ) {
575  hist << std::setw(10) << std::left << iter_;
576  hist << std::setw(10) << std::left << flag_;
577  }
578  else {
579  hist << std::setw(10) << std::left << iterCR_;
580  hist << std::setw(10) << std::left << flagCR_;
581  }
582  if ( feasible_ ) {
583  hist << std::setw(10) << std::left << "YES";
584  }
585  else {
586  hist << std::setw(10) << std::left << "NO";
587  }
588  hist << "\n";
589  }
590  return hist.str();
591  }
592 
593 }; // class PrimalDualActiveSetStep
594 
595 } // namespace ROL
596 
597 #endif
598 
599 // void solve(Vector<Real> &sol, const Vector<Real> &rhs, const Vector<Real> &xlam, const Vector<Real> &x,
600 // Objective<Real> &obj, BoundConstraint<Real> &con) {
601 // Real rnorm = rhs.norm();
602 // Real rtol = std::min(tol1_,tol2_*rnorm);
603 // itol_ = std::sqrt(ROL_EPSILON<Real>());
604 // sol.zero();
605 //
606 // ROL::Ptr<Vector<Real> > res = rhs.clone();
607 // res->set(rhs);
608 //
609 // ROL::Ptr<Vector<Real> > v = x.clone();
610 // con.pruneActive(*res,xlam,neps_);
611 // obj.precond(*v,*res,x,itol_);
612 // con.pruneActive(*v,xlam,neps_);
613 //
614 // ROL::Ptr<Vector<Real> > p = x.clone();
615 // p->set(*v);
616 //
617 // ROL::Ptr<Vector<Real> > Hp = x.clone();
618 //
619 // iterCR_ = 0;
620 // flagCR_ = 0;
621 //
622 // Real kappa = 0.0, beta = 0.0, alpha = 0.0, tmp = 0.0, rv = v->dot(*res);
623 //
624 // for (iterCR_ = 0; iterCR_ < maxitCR_; iterCR_++) {
625 // if ( false ) {
626 // itol_ = rtol/(maxitCR_*rnorm);
627 // }
628 // con.pruneActive(*p,xlam,neps_);
629 // if ( secant_ == ROL::nullPtr ) {
630 // obj.hessVec(*Hp, *p, x, itol_);
631 // }
632 // else {
633 // secant_->applyB( *Hp, *p, x );
634 // }
635 // con.pruneActive(*Hp,xlam,neps_);
636 //
637 // kappa = p->dot(*Hp);
638 // if ( kappa <= 0.0 ) { flagCR_ = 2; break; }
639 // alpha = rv/kappa;
640 // sol.axpy(alpha,*p);
641 //
642 // res->axpy(-alpha,*Hp);
643 // rnorm = res->norm();
644 // if ( rnorm < rtol ) { break; }
645 //
646 // con.pruneActive(*res,xlam,neps_);
647 // obj.precond(*v,*res,x,itol_);
648 // con.pruneActive(*v,xlam,neps_);
649 // tmp = rv;
650 // rv = v->dot(*res);
651 // beta = rv/tmp;
652 //
653 // p->scale(beta);
654 // p->axpy(1.0,*v);
655 // }
656 // if ( iterCR_ == maxitCR_ ) {
657 // flagCR_ = 1;
658 // }
659 // else {
660 // iterCR_++;
661 // }
662 // }
663 
664 
665 // /** \brief Apply the inactive components of the Hessian operator.
666 //
667 // I.e., the components corresponding to \f$\mathcal{I}_k\f$.
668 //
669 // @param[out] hv is the result of applying the Hessian at @b x to
670 // @b v
671 // @param[in] v is the direction in which we apply the Hessian
672 // @param[in] x is the current iteration vector \f$x_k\f$
673 // @param[in] xlam is the vector \f$x_k + c\lambda_k\f$
674 // @param[in] obj is the objective function
675 // @param[in] con are the bound constraints
676 // */
677 // void applyInactiveHessian(Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x,
678 // const Vector<Real> &xlam, Objective<Real> &obj, BoundConstraint<Real> &con) {
679 // ROL::Ptr<Vector<Real> > tmp = v.clone();
680 // tmp->set(v);
681 // con.pruneActive(*tmp,xlam,neps_);
682 // if ( secant_ == ROL::nullPtr ) {
683 // obj.hessVec(hv,*tmp,x,itol_);
684 // }
685 // else {
686 // secant_->applyB(hv,*tmp,x);
687 // }
688 // con.pruneActive(hv,xlam,neps_);
689 // }
690 //
691 // /** \brief Apply the inactive components of the preconditioner operator.
692 //
693 // I.e., the components corresponding to \f$\mathcal{I}_k\f$.
694 //
695 // @param[out] hv is the result of applying the preconditioner at @b x to
696 // @b v
697 // @param[in] v is the direction in which we apply the preconditioner
698 // @param[in] x is the current iteration vector \f$x_k\f$
699 // @param[in] xlam is the vector \f$x_k + c\lambda_k\f$
700 // @param[in] obj is the objective function
701 // @param[in] con are the bound constraints
702 // */
703 // void applyInactivePrecond(Vector<Real> &pv, const Vector<Real> &v, const Vector<Real> &x,
704 // const Vector<Real> &xlam, Objective<Real> &obj, BoundConstraint<Real> &con) {
705 // ROL::Ptr<Vector<Real> > tmp = v.clone();
706 // tmp->set(v);
707 // con.pruneActive(*tmp,xlam,neps_);
708 // obj.precond(pv,*tmp,x,itol_);
709 // con.pruneActive(pv,xlam,neps_);
710 // }
711 //
712 // /** \brief Solve the inactive part of the PDAS optimality system.
713 //
714 // The inactive PDAS optimality system is
715 // \f[
716 // \nabla^2 f(x_k)_{\mathcal{I}_k,\mathcal{I}_k}s =
717 // -\nabla f(x_k)_{\mathcal{I}_k}
718 // -\nabla^2 f(x_k)_{\mathcal{I}_k,\mathcal{A}_k} (s_k)_{\mathcal{A}_k}.
719 // \f]
720 // Since the inactive part of the Hessian may not be positive definite, we solve
721 // using CR.
722 //
723 // @param[out] sol is the vector containing the solution
724 // @param[in] rhs is the right-hand side vector
725 // @param[in] xlam is the vector \f$x_k + c\lambda_k\f$
726 // @param[in] x is the current iteration vector \f$x_k\f$
727 // @param[in] obj is the objective function
728 // @param[in] con are the bound constraints
729 // */
730 // // Solve the inactive part of the optimality system using conjugate residuals
731 // void solve(Vector<Real> &sol, const Vector<Real> &rhs, const Vector<Real> &xlam, const Vector<Real> &x,
732 // Objective<Real> &obj, BoundConstraint<Real> &con) {
733 // // Initialize Residual
734 // ROL::Ptr<Vector<Real> > res = rhs.clone();
735 // res->set(rhs);
736 // Real rnorm = res->norm();
737 // Real rtol = std::min(tol1_,tol2_*rnorm);
738 // if ( false ) { itol_ = rtol/(maxitCR_*rnorm); }
739 // sol.zero();
740 //
741 // // Apply preconditioner to residual r = Mres
742 // ROL::Ptr<Vector<Real> > r = x.clone();
743 // applyInactivePrecond(*r,*res,x,xlam,obj,con);
744 //
745 // // Initialize direction p = v
746 // ROL::Ptr<Vector<Real> > p = x.clone();
747 // p->set(*r);
748 //
749 // // Apply Hessian to v
750 // ROL::Ptr<Vector<Real> > Hr = x.clone();
751 // applyInactiveHessian(*Hr,*r,x,xlam,obj,con);
752 //
753 // // Apply Hessian to p
754 // ROL::Ptr<Vector<Real> > Hp = x.clone();
755 // ROL::Ptr<Vector<Real> > MHp = x.clone();
756 // Hp->set(*Hr);
757 //
758 // iterCR_ = 0;
759 // flagCR_ = 0;
760 //
761 // Real kappa = 0.0, beta = 0.0, alpha = 0.0, tmp = 0.0, rHr = Hr->dot(*r);
762 //
763 // for (iterCR_ = 0; iterCR_ < maxitCR_; iterCR_++) {
764 // // Precondition Hp
765 // applyInactivePrecond(*MHp,*Hp,x,xlam,obj,con);
766 //
767 // kappa = Hp->dot(*MHp); // p' H M H p
768 // alpha = rHr/kappa; // r' M H M r
769 // sol.axpy(alpha,*p); // update step
770 // res->axpy(-alpha,*Hp); // residual
771 // r->axpy(-alpha,*MHp); // preconditioned residual
772 //
773 // // recompute rnorm and decide whether or not to exit
774 // rnorm = res->norm();
775 // if ( rnorm < rtol ) { break; }
776 //
777 // // Apply Hessian to v
778 // itol_ = rtol/(maxitCR_*rnorm);
779 // applyInactiveHessian(*Hr,*r,x,xlam,obj,con);
780 //
781 // tmp = rHr;
782 // rHr = Hr->dot(*r);
783 // beta = rHr/tmp;
784 // p->scale(beta);
785 // p->axpy(1.0,*r);
786 // Hp->scale(beta);
787 // Hp->axpy(1.0,*Hr);
788 // }
789 // if ( iterCR_ == maxitCR_ ) {
790 // flagCR_ = 1;
791 // }
792 // else {
793 // iterCR_++;
794 // }
795 // }
ROL::BoundConstraint::isFeasible
virtual bool isFeasible(const Vector< Real > &v)
Check if the vector, v, is feasible.
Definition: ROL_BoundConstraint_Def.hpp:98
ROL::PrimalDualActiveSetStep
Implements the computation of optimization steps with the Newton primal-dual active set method...
Definition: ROL_PrimalDualActiveSetStep.hpp:98
ROL::Objective
Provides the interface to evaluate objective functions.
Definition: ROL_Objective.hpp:44
ROL::PrimalDualActiveSetStep::printName
std::string printName(void) const
Print step name.
Definition: ROL_PrimalDualActiveSetStep.hpp:537
ROL::Vector::dual
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: ROL_Vector.hpp:192
ROL::Vector::clone
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
ROL::Vector::plus
virtual void plus(const Vector &x)=0
Compute , where .
ROL::PrimalDualActiveSetStep::printHeader
std::string printHeader(void) const
Print iterate header.
Definition: ROL_PrimalDualActiveSetStep.hpp:510
ROL::BoundConstraint::pruneUpperActive
virtual void pruneUpperActive(Vector< Real > &v, const Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the upper -active set.
Definition: ROL_BoundConstraint_Def.hpp:54
ROL::Objective::value
virtual Real value(const Vector< Real > &x, Real &tol)=0
Compute value.
ROL::Step
Provides the interface to compute optimization steps.
Definition: ROL_Step.hpp:34
ROL::PrimalDualActiveSetStep::PrecondPD::bnd_
const ROL::Ptr< BoundConstraint< Real > > bnd_
Definition: ROL_PrimalDualActiveSetStep.hpp:177
ROL::Objective::hessVec
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
Definition: ROL_ObjectiveDef.hpp:61
ROL_Types.hpp
Contains definitions of custom data types in ROL.
ROL::PrimalDualActiveSetStep::HessianPD::HessianPD
HessianPD(const ROL::Ptr< Objective< Real > > &obj, const ROL::Ptr< BoundConstraint< Real > > &bnd, const ROL::Ptr< Vector< Real > > &x, const ROL::Ptr< Vector< Real > > &xlam, const Real eps=0, const ROL::Ptr< Secant< Real > > &secant=ROL::nullPtr, const bool useSecant=false)
Definition: ROL_PrimalDualActiveSetStep.hpp:147
ROL::StringToESecant
ESecant StringToESecant(std::string s)
Definition: ROL_Types.hpp:509
ROL::PrimalDualActiveSetStep::secant_
ROL::Ptr< Secant< Real > > secant_
Secant object.
Definition: ROL_PrimalDualActiveSetStep.hpp:132
ROL::Vector::zero
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:133
ROL::PrimalDualActiveSetStep::scale_
Real scale_
Scale for dual variables in the active set, .
Definition: ROL_PrimalDualActiveSetStep.hpp:114
ROL::Vector
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
ROL::Objective::update
virtual void update(const Vector< Real > &x, UpdateType type, int iter=-1)
Update objective function.
Definition: ROL_Objective.hpp:62
ROL::PrimalDualActiveSetStep::gtmp_
ROL::Ptr< Vector< Real > > gtmp_
Container for temporary gradient storage.
Definition: ROL_PrimalDualActiveSetStep.hpp:128
zero
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Definition: ROL_Objective_SerialSimOpt.hpp:111
ROL::BoundConstraint::getLowerBound
virtual const Ptr< const Vector< Real > > getLowerBound(void) const
Return the ref count pointer to the lower bound vector.
Definition: ROL_BoundConstraint_Def.hpp:82
ROL::PrimalDualActiveSetStep::res_
ROL::Ptr< Vector< Real > > res_
Container for optimality system residual for quadratic model.
Definition: ROL_PrimalDualActiveSetStep.hpp:125
ROL::AlgorithmState
State for algorithm class. Will be used for restarts.
Definition: ROL_Types.hpp:109
ROL::PrimalDualActiveSetStep::x0_
ROL::Ptr< Vector< Real > > x0_
Container for initial priaml variables.
Definition: ROL_PrimalDualActiveSetStep.hpp:121
ROL::Objective::gradient
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
Definition: ROL_ObjectiveDef.hpp:40
ROL::PrimalDualActiveSetStep::PrimalDualActiveSetStep
PrimalDualActiveSetStep(ROL::ParameterList &parlist)
Constructor.
Definition: ROL_PrimalDualActiveSetStep.hpp:245
ROL::ESecant
ESecant
Enumeration of secant update algorithms.
Definition: ROL_Types.hpp:452
ROL::PrimalDualActiveSetStep::rtmp_
ROL::Ptr< Vector< Real > > rtmp_
Container for temporary right hand side storage.
Definition: ROL_PrimalDualActiveSetStep.hpp:127
ROL::PrimalDualActiveSetStep::Ag_
ROL::Ptr< Vector< Real > > Ag_
Container for gradient projected onto active set.
Definition: ROL_PrimalDualActiveSetStep.hpp:126
ROL::Step::getState
ROL::Ptr< StepState< Real > > getState(void)
Definition: ROL_Step.hpp:39
ROL::PrimalDualActiveSetStep::initialize
void initialize(Vector< Real > &x, const Vector< Real > &s, const Vector< Real > &g, Objective< Real > &obj, BoundConstraint< Real > &con, AlgorithmState< Real > &algo_state)
Initialize step with bound constraint.
Definition: ROL_PrimalDualActiveSetStep.hpp:285
ROL::PrimalDualActiveSetStep::HessianPD::bnd_
const ROL::Ptr< BoundConstraint< Real > > bnd_
Definition: ROL_PrimalDualActiveSetStep.hpp:139
ROL_BoundConstraint.hpp
ROL::PrimalDualActiveSetStep::feasible_
bool feasible_
Flag whether the current iterate is feasible or not.
Definition: ROL_PrimalDualActiveSetStep.hpp:116
ROL::BoundConstraint::pruneLowerActive
virtual void pruneLowerActive(Vector< Real > &v, const Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the lower -active set.
Definition: ROL_BoundConstraint_Def.hpp:68
ROL::Secant
Provides interface for and implements limited-memory secant operators.
Definition: ROL_Secant.hpp:45
ROL::PrimalDualActiveSetStep::HessianPD::apply
void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
Definition: ROL_PrimalDualActiveSetStep.hpp:161
ROL::PrimalDualActiveSetStep::computeCriticalityMeasure
Real computeCriticalityMeasure(Vector< Real > &x, Objective< Real > &obj, BoundConstraint< Real > &con, Real tol)
Compute the gradient-based criticality measure.
Definition: ROL_PrimalDualActiveSetStep.hpp:227
ROL::AlgorithmState::iterateVec
ROL::Ptr< Vector< Real > > iterateVec
Definition: ROL_Types.hpp:123
ROL::PrimalDualActiveSetStep::xbnd_
ROL::Ptr< Vector< Real > > xbnd_
Container for primal variable bounds.
Definition: ROL_PrimalDualActiveSetStep.hpp:122
ROL_Objective.hpp
ROL_Vector.hpp
ROL::BoundConstraint::project
virtual void project(Vector< Real > &x)
Project optimization variables onto the bounds.
Definition: ROL_BoundConstraint_Def.hpp:40
ROL::BoundConstraint::getUpperBound
virtual const Ptr< const Vector< Real > > getUpperBound(void) const
Return the ref count pointer to the upper bound vector.
Definition: ROL_BoundConstraint_Def.hpp:90
ROL::BoundConstraint::pruneActive
void pruneActive(Vector< Real > &v, const Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the -active set.
Definition: ROL_BoundConstraint_Def.hpp:169
ROL::LinearOperator
Provides the interface to apply a linear operator.
Definition: ROL_LinearOperator.hpp:37
ROL::BoundConstraint
Provides the interface to apply upper and lower bound constraints.
Definition: ROL_BoundConstraint.hpp:39
ROL::PrimalDualActiveSetStep::compute
void compute(Vector< Real > &s, const Vector< Real > &x, Objective< Real > &obj, BoundConstraint< Real > &con, AlgorithmState< Real > &algo_state)
Compute step.
Definition: ROL_PrimalDualActiveSetStep.hpp:345
ROL::PrimalDualActiveSetStep::gtol_
Real gtol_
PDAS gradient stopping tolerance.
Definition: ROL_PrimalDualActiveSetStep.hpp:113
ROL::PrimalDualActiveSetStep::PrecondPD::PrecondPD
PrecondPD(const ROL::Ptr< Objective< Real > > &obj, const ROL::Ptr< BoundConstraint< Real > > &bnd, const ROL::Ptr< Vector< Real > > &x, const ROL::Ptr< Vector< Real > > &xlam, const Real eps=0, const ROL::Ptr< Secant< Real > > &secant=ROL::nullPtr, const bool useSecant=false)
Definition: ROL_PrimalDualActiveSetStep.hpp:185
ROL::PrimalDualActiveSetStep::stol_
Real stol_
PDAS minimum step size stopping tolerance.
Definition: ROL_PrimalDualActiveSetStep.hpp:112
ROL::PrimalDualActiveSetStep::lambda_
ROL::Ptr< Vector< Real > > lambda_
Container for dual variables.
Definition: ROL_PrimalDualActiveSetStep.hpp:119
ROL::Vector::set
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:175
ROL_KrylovFactory.hpp
ROL::Vector::norm
virtual Real norm() const =0
Returns where .
ROL::ROL_EPSILON
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:57
ROL::PrimalDualActiveSetStep::PrecondPD::applyInverse
void applyInverse(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply inverse of linear operator.
Definition: ROL_PrimalDualActiveSetStep.hpp:202
ROL::PrimalDualActiveSetStep::xlam_
ROL::Ptr< Vector< Real > > xlam_
Container for primal plus dual variables.
Definition: ROL_PrimalDualActiveSetStep.hpp:120
ROL::PrimalDualActiveSetStep::esec_
ESecant esec_
Enum for secant type.
Definition: ROL_PrimalDualActiveSetStep.hpp:131
ROL_Step.hpp
ROL::PrimalDualActiveSetStep::update
void update(Vector< Real > &x, const Vector< Real > &s, Objective< Real > &obj, BoundConstraint< Real > &con, AlgorithmState< Real > &algo_state)
Update step, if successful.
Definition: ROL_PrimalDualActiveSetStep.hpp:478
ROL::PrimalDualActiveSetStep::print
virtual std::string print(AlgorithmState< Real > &algo_state, bool print_header=false) const
Print iterate status.
Definition: ROL_PrimalDualActiveSetStep.hpp:550
ROL::PrimalDualActiveSetStep::As_
ROL::Ptr< Vector< Real > > As_
Container for step projected onto active set.
Definition: ROL_PrimalDualActiveSetStep.hpp:123
ROL_Secant.hpp
ROL::PrimalDualActiveSetStep::xtmp_
ROL::Ptr< Vector< Real > > xtmp_
Container for temporary primal storage.
Definition: ROL_PrimalDualActiveSetStep.hpp:124
ROL::PrimalDualActiveSetStep::PrecondPD::apply
void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
Definition: ROL_PrimalDualActiveSetStep.hpp:199
ROL::PrimalDualActiveSetStep::maxit_
int maxit_
Maximum number of PDAS iterations.
Definition: ROL_PrimalDualActiveSetStep.hpp:109
Generated by   doxygen 1.8.5