ROL
ROL_PrimalDualActiveSetStep.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_PRIMALDUALACTIVESETSTEP_H
45#define ROL_PRIMALDUALACTIVESETSTEP_H
46
47#include "ROL_Step.hpp"
48#include "ROL_Vector.hpp"
49#include "ROL_KrylovFactory.hpp"
50#include "ROL_Objective.hpp"
52#include "ROL_Types.hpp"
53#include "ROL_Secant.hpp"
54#include "ROL_ParameterList.hpp"
55
130namespace ROL {
131
132template <class Real>
133class PrimalDualActiveSetStep : public Step<Real> {
134private:
135
136 ROL::Ptr<Krylov<Real> > krylov_;
137
138 // Krylov Parameters
141 Real itol_;
142
143 // PDAS Parameters
144 int maxit_;
145 int iter_;
146 int flag_;
147 Real stol_;
148 Real gtol_;
149 Real scale_;
150 Real neps_;
152
153 // Dual Variable
154 ROL::Ptr<Vector<Real> > lambda_;
155 ROL::Ptr<Vector<Real> > xlam_;
156 ROL::Ptr<Vector<Real> > x0_;
157 ROL::Ptr<Vector<Real> > xbnd_;
158 ROL::Ptr<Vector<Real> > As_;
159 ROL::Ptr<Vector<Real> > xtmp_;
160 ROL::Ptr<Vector<Real> > res_;
161 ROL::Ptr<Vector<Real> > Ag_;
162 ROL::Ptr<Vector<Real> > rtmp_;
163 ROL::Ptr<Vector<Real> > gtmp_;
164
165 // Secant Information
167 ROL::Ptr<Secant<Real> > secant_;
170
171 class HessianPD : public LinearOperator<Real> {
172 private:
173 const ROL::Ptr<Objective<Real> > obj_;
174 const ROL::Ptr<BoundConstraint<Real> > bnd_;
175 const ROL::Ptr<Vector<Real> > x_;
176 const ROL::Ptr<Vector<Real> > xlam_;
177 ROL::Ptr<Vector<Real> > v_;
178 Real eps_;
179 const ROL::Ptr<Secant<Real> > secant_;
181 public:
182 HessianPD(const ROL::Ptr<Objective<Real> > &obj,
183 const ROL::Ptr<BoundConstraint<Real> > &bnd,
184 const ROL::Ptr<Vector<Real> > &x,
185 const ROL::Ptr<Vector<Real> > &xlam,
186 const Real eps = 0,
187 const ROL::Ptr<Secant<Real> > &secant = ROL::nullPtr,
188 const bool useSecant = false )
189 : obj_(obj), bnd_(bnd), x_(x), xlam_(xlam),
190 eps_(eps), secant_(secant), useSecant_(useSecant) {
191 v_ = x_->clone();
192 if ( !useSecant || secant == ROL::nullPtr ) {
193 useSecant_ = false;
194 }
195 }
196 void apply( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
197 v_->set(v);
198 bnd_->pruneActive(*v_,*xlam_,eps_);
199 if ( useSecant_ ) {
200 secant_->applyB(Hv,*v_);
201 }
202 else {
203 obj_->hessVec(Hv,*v_,*x_,tol);
204 }
205 bnd_->pruneActive(Hv,*xlam_,eps_);
206 }
207 };
208
209 class PrecondPD : public LinearOperator<Real> {
210 private:
211 const ROL::Ptr<Objective<Real> > obj_;
212 const ROL::Ptr<BoundConstraint<Real> > bnd_;
213 const ROL::Ptr<Vector<Real> > x_;
214 const ROL::Ptr<Vector<Real> > xlam_;
215 ROL::Ptr<Vector<Real> > v_;
216 Real eps_;
217 const ROL::Ptr<Secant<Real> > secant_;
219 public:
220 PrecondPD(const ROL::Ptr<Objective<Real> > &obj,
221 const ROL::Ptr<BoundConstraint<Real> > &bnd,
222 const ROL::Ptr<Vector<Real> > &x,
223 const ROL::Ptr<Vector<Real> > &xlam,
224 const Real eps = 0,
225 const ROL::Ptr<Secant<Real> > &secant = ROL::nullPtr,
226 const bool useSecant = false )
227 : obj_(obj), bnd_(bnd), x_(x), xlam_(xlam),
228 eps_(eps), secant_(secant), useSecant_(useSecant) {
229 v_ = x_->dual().clone();
230 if ( !useSecant || secant == ROL::nullPtr ) {
231 useSecant_ = false;
232 }
233 }
234 void apply( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
235 Hv.set(v.dual());
236 }
237 void applyInverse( Vector<Real> &Hv, const Vector<Real> &v, Real &tol ) const {
238 v_->set(v);
239 bnd_->pruneActive(*v_,*xlam_,eps_);
240 if ( useSecant_ ) {
241 secant_->applyH(Hv,*v_);
242 }
243 else {
244 obj_->precond(Hv,*v_,*x_,tol);
245 }
246 bnd_->pruneActive(Hv,*xlam_,eps_);
247 }
248 };
249
263 Real one(1);
264 ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState();
265 obj.gradient(*(step_state->gradientVec),x,tol);
266 xtmp_->set(x);
267 xtmp_->axpy(-one,(step_state->gradientVec)->dual());
268 con.project(*xtmp_);
269 xtmp_->axpy(-one,x);
270 return xtmp_->norm();
271 }
272
273public:
280 PrimalDualActiveSetStep( ROL::ParameterList &parlist )
281 : Step<Real>::Step(), krylov_(ROL::nullPtr),
282 iterCR_(0), flagCR_(0), itol_(0),
283 maxit_(0), iter_(0), flag_(0), stol_(0), gtol_(0), scale_(0),
284 neps_(-ROL_EPSILON<Real>()), feasible_(false),
285 lambda_(ROL::nullPtr), xlam_(ROL::nullPtr), x0_(ROL::nullPtr),
286 xbnd_(ROL::nullPtr), As_(ROL::nullPtr), xtmp_(ROL::nullPtr),
287 res_(ROL::nullPtr), Ag_(ROL::nullPtr), rtmp_(ROL::nullPtr),
288 gtmp_(ROL::nullPtr),
289 esec_(SECANT_LBFGS), secant_(ROL::nullPtr), useSecantPrecond_(false),
290 useSecantHessVec_(false) {
291 Real one(1), oem6(1.e-6), oem8(1.e-8);
292 // Algorithmic parameters
293 maxit_ = parlist.sublist("Step").sublist("Primal Dual Active Set").get("Iteration Limit",10);
294 stol_ = parlist.sublist("Step").sublist("Primal Dual Active Set").get("Relative Step Tolerance",oem8);
295 gtol_ = parlist.sublist("Step").sublist("Primal Dual Active Set").get("Relative Gradient Tolerance",oem6);
296 scale_ = parlist.sublist("Step").sublist("Primal Dual Active Set").get("Dual Scaling", one);
297 // Build secant object
298 esec_ = StringToESecant(parlist.sublist("General").sublist("Secant").get("Type","Limited-Memory BFGS"));
299 useSecantHessVec_ = parlist.sublist("General").sublist("Secant").get("Use as Hessian", false);
300 useSecantPrecond_ = parlist.sublist("General").sublist("Secant").get("Use as Preconditioner", false);
302 secant_ = SecantFactory<Real>(parlist);
303 }
304 // Build Krylov object
305 krylov_ = KrylovFactory<Real>(parlist);
306 }
307
319 using Step<Real>::initialize;
320 void initialize( Vector<Real> &x, const Vector<Real> &s, const Vector<Real> &g,
322 AlgorithmState<Real> &algo_state ) {
323 ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState();
324 Real zero(0), one(1);
325 // Initialize state descent direction and gradient storage
326 step_state->descentVec = s.clone();
327 step_state->gradientVec = g.clone();
328 step_state->searchSize = zero;
329 // Initialize additional storage
330 xlam_ = x.clone();
331 x0_ = x.clone();
332 xbnd_ = x.clone();
333 As_ = s.clone();
334 xtmp_ = x.clone();
335 res_ = g.clone();
336 Ag_ = g.clone();
337 rtmp_ = g.clone();
338 gtmp_ = g.clone();
339 // Project x onto constraint set
340 con.project(x);
341 // Update objective function, get value, and get gradient
342 Real tol = std::sqrt(ROL_EPSILON<Real>());
343 obj.update(x,true,algo_state.iter);
344 algo_state.value = obj.value(x,tol);
345 algo_state.nfval++;
346 algo_state.gnorm = computeCriticalityMeasure(x,obj,con,tol);
347 algo_state.ngrad++;
348 // Initialize dual variable
349 lambda_ = s.clone();
350 lambda_->set((step_state->gradientVec)->dual());
351 lambda_->scale(-one);
352 }
353
379 using Step<Real>::compute;
381 AlgorithmState<Real> &algo_state ) {
382 ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState();
383 Real zero(0), one(1);
384 s.zero();
385 x0_->set(x);
386 res_->set(*(step_state->gradientVec));
387 for ( iter_ = 0; iter_ < maxit_; iter_++ ) {
388 /********************************************************************/
389 // MODIFY ITERATE VECTOR TO CHECK ACTIVE SET
390 /********************************************************************/
391 xlam_->set(*x0_); // xlam = x0
392 xlam_->axpy(scale_,*(lambda_)); // xlam = x0 + c*lambda
393 /********************************************************************/
394 // PROJECT x ONTO PRIMAL DUAL FEASIBLE SET
395 /********************************************************************/
396 As_->zero(); // As = 0
397
398 xbnd_->set(*con.getUpperBound()); // xbnd = u
399 xbnd_->axpy(-one,x); // xbnd = u - x
400 xtmp_->set(*xbnd_); // tmp = u - x
401 con.pruneUpperActive(*xtmp_,*xlam_,neps_); // tmp = I(u - x)
402 xbnd_->axpy(-one,*xtmp_); // xbnd = A(u - x)
403 As_->plus(*xbnd_); // As += A(u - x)
404
405 xbnd_->set(*con.getLowerBound()); // xbnd = l
406 xbnd_->axpy(-one,x); // xbnd = l - x
407 xtmp_->set(*xbnd_); // tmp = l - x
408 con.pruneLowerActive(*xtmp_,*xlam_,neps_); // tmp = I(l - x)
409 xbnd_->axpy(-one,*xtmp_); // xbnd = A(l - x)
410 As_->plus(*xbnd_); // As += A(l - x)
411 /********************************************************************/
412 // APPLY HESSIAN TO ACTIVE COMPONENTS OF s AND REMOVE INACTIVE
413 /********************************************************************/
414 itol_ = std::sqrt(ROL_EPSILON<Real>());
415 if ( useSecantHessVec_ && secant_ != ROL::nullPtr ) { // IHAs = H*As
416 secant_->applyB(*gtmp_,*As_);
417 }
418 else {
419 obj.hessVec(*gtmp_,*As_,x,itol_);
420 }
421 con.pruneActive(*gtmp_,*xlam_,neps_); // IHAs = I(H*As)
422 /********************************************************************/
423 // SEPARATE ACTIVE AND INACTIVE COMPONENTS OF THE GRADIENT
424 /********************************************************************/
425 rtmp_->set(*(step_state->gradientVec)); // Inactive components
427
428 Ag_->set(*(step_state->gradientVec)); // Active components
429 Ag_->axpy(-one,*rtmp_);
430 /********************************************************************/
431 // SOLVE REDUCED NEWTON SYSTEM
432 /********************************************************************/
433 rtmp_->plus(*gtmp_);
434 rtmp_->scale(-one); // rhs = -Ig - I(H*As)
435 s.zero();
436 if ( rtmp_->norm() > zero ) {
437 // Initialize Hessian and preconditioner
438 ROL::Ptr<Objective<Real> > obj_ptr = ROL::makePtrFromRef(obj);
439 ROL::Ptr<BoundConstraint<Real> > con_ptr = ROL::makePtrFromRef(con);
440 ROL::Ptr<LinearOperator<Real> > hessian
441 = ROL::makePtr<HessianPD>(obj_ptr,con_ptr,
443 ROL::Ptr<LinearOperator<Real> > precond
444 = ROL::makePtr<PrecondPD>(obj_ptr,con_ptr,
446 //solve(s,*rtmp_,*xlam_,x,obj,con); // Call conjugate residuals
447 krylov_->run(s,*hessian,*rtmp_,*precond,iterCR_,flagCR_);
448 con.pruneActive(s,*xlam_,neps_); // s <- Is
449 }
450 s.plus(*As_); // s = Is + As
451 /********************************************************************/
452 // UPDATE MULTIPLIER
453 /********************************************************************/
454 if ( useSecantHessVec_ && secant_ != ROL::nullPtr ) {
455 secant_->applyB(*rtmp_,s);
456 }
457 else {
458 obj.hessVec(*rtmp_,s,x,itol_);
459 }
460 gtmp_->set(*rtmp_);
462 lambda_->set(*rtmp_);
463 lambda_->axpy(-one,*gtmp_);
464 lambda_->plus(*Ag_);
465 lambda_->scale(-one);
466 /********************************************************************/
467 // UPDATE STEP
468 /********************************************************************/
469 x0_->set(x);
470 x0_->plus(s);
471 res_->set(*(step_state->gradientVec));
472 res_->plus(*rtmp_);
473 // Compute criticality measure
474 xtmp_->set(*x0_);
475 xtmp_->axpy(-one,res_->dual());
476 con.project(*xtmp_);
477 xtmp_->axpy(-one,*x0_);
478// std::cout << s.norm() << " "
479// << tmp->norm() << " "
480// << res_->norm() << " "
481// << lambda_->norm() << " "
482// << flagCR_ << " "
483// << iterCR_ << "\n";
484 if ( xtmp_->norm() < gtol_*algo_state.gnorm ) {
485 flag_ = 0;
486 break;
487 }
488 if ( s.norm() < stol_*x.norm() ) {
489 flag_ = 2;
490 break;
491 }
492 }
493 if ( iter_ == maxit_ ) {
494 flag_ = 1;
495 }
496 else {
497 iter_++;
498 }
499 }
500
512 using Step<Real>::update;
514 AlgorithmState<Real> &algo_state ) {
515 ROL::Ptr<StepState<Real> > step_state = Step<Real>::getState();
516 step_state->SPiter = (maxit_ > 1) ? iter_ : iterCR_;
517 step_state->SPflag = (maxit_ > 1) ? flag_ : flagCR_;
518
519 x.plus(s);
520 feasible_ = con.isFeasible(x);
521 algo_state.snorm = s.norm();
522 algo_state.iter++;
523 Real tol = std::sqrt(ROL_EPSILON<Real>());
524 obj.update(x,true,algo_state.iter);
525 algo_state.value = obj.value(x,tol);
526 algo_state.nfval++;
527
528 if ( secant_ != ROL::nullPtr ) {
529 gtmp_->set(*(step_state->gradientVec));
530 }
531 algo_state.gnorm = computeCriticalityMeasure(x,obj,con,tol);
532 algo_state.ngrad++;
533
534 if ( secant_ != ROL::nullPtr ) {
535 secant_->updateStorage(x,*(step_state->gradientVec),*gtmp_,s,algo_state.snorm,algo_state.iter+1);
536 }
537 (algo_state.iterateVec)->set(x);
538 }
539
545 std::string printHeader( void ) const {
546 std::stringstream hist;
547 hist << " ";
548 hist << std::setw(6) << std::left << "iter";
549 hist << std::setw(15) << std::left << "value";
550 hist << std::setw(15) << std::left << "gnorm";
551 hist << std::setw(15) << std::left << "snorm";
552 hist << std::setw(10) << std::left << "#fval";
553 hist << std::setw(10) << std::left << "#grad";
554 if ( maxit_ > 1 ) {
555 hist << std::setw(10) << std::left << "iterPDAS";
556 hist << std::setw(10) << std::left << "flagPDAS";
557 }
558 else {
559 hist << std::setw(10) << std::left << "iterCR";
560 hist << std::setw(10) << std::left << "flagCR";
561 }
562 hist << std::setw(10) << std::left << "feasible";
563 hist << "\n";
564 return hist.str();
565 }
566
572 std::string printName( void ) const {
573 std::stringstream hist;
574 hist << "\nPrimal Dual Active Set Newton's Method\n";
575 return hist.str();
576 }
577
585 virtual std::string print( AlgorithmState<Real> &algo_state, bool print_header = false ) const {
586 std::stringstream hist;
587 hist << std::scientific << std::setprecision(6);
588 if ( algo_state.iter == 0 ) {
589 hist << printName();
590 }
591 if ( print_header ) {
592 hist << printHeader();
593 }
594 if ( algo_state.iter == 0 ) {
595 hist << " ";
596 hist << std::setw(6) << std::left << algo_state.iter;
597 hist << std::setw(15) << std::left << algo_state.value;
598 hist << std::setw(15) << std::left << algo_state.gnorm;
599 hist << "\n";
600 }
601 else {
602 hist << " ";
603 hist << std::setw(6) << std::left << algo_state.iter;
604 hist << std::setw(15) << std::left << algo_state.value;
605 hist << std::setw(15) << std::left << algo_state.gnorm;
606 hist << std::setw(15) << std::left << algo_state.snorm;
607 hist << std::setw(10) << std::left << algo_state.nfval;
608 hist << std::setw(10) << std::left << algo_state.ngrad;
609 if ( maxit_ > 1 ) {
610 hist << std::setw(10) << std::left << iter_;
611 hist << std::setw(10) << std::left << flag_;
612 }
613 else {
614 hist << std::setw(10) << std::left << iterCR_;
615 hist << std::setw(10) << std::left << flagCR_;
616 }
617 if ( feasible_ ) {
618 hist << std::setw(10) << std::left << "YES";
619 }
620 else {
621 hist << std::setw(10) << std::left << "NO";
622 }
623 hist << "\n";
624 }
625 return hist.str();
626 }
627
628}; // class PrimalDualActiveSetStep
629
630} // namespace ROL
631
632#endif
633
634// void solve(Vector<Real> &sol, const Vector<Real> &rhs, const Vector<Real> &xlam, const Vector<Real> &x,
635// Objective<Real> &obj, BoundConstraint<Real> &con) {
636// Real rnorm = rhs.norm();
637// Real rtol = std::min(tol1_,tol2_*rnorm);
638// itol_ = std::sqrt(ROL_EPSILON<Real>());
639// sol.zero();
640//
641// ROL::Ptr<Vector<Real> > res = rhs.clone();
642// res->set(rhs);
643//
644// ROL::Ptr<Vector<Real> > v = x.clone();
645// con.pruneActive(*res,xlam,neps_);
646// obj.precond(*v,*res,x,itol_);
647// con.pruneActive(*v,xlam,neps_);
648//
649// ROL::Ptr<Vector<Real> > p = x.clone();
650// p->set(*v);
651//
652// ROL::Ptr<Vector<Real> > Hp = x.clone();
653//
654// iterCR_ = 0;
655// flagCR_ = 0;
656//
657// Real kappa = 0.0, beta = 0.0, alpha = 0.0, tmp = 0.0, rv = v->dot(*res);
658//
659// for (iterCR_ = 0; iterCR_ < maxitCR_; iterCR_++) {
660// if ( false ) {
661// itol_ = rtol/(maxitCR_*rnorm);
662// }
663// con.pruneActive(*p,xlam,neps_);
664// if ( secant_ == ROL::nullPtr ) {
665// obj.hessVec(*Hp, *p, x, itol_);
666// }
667// else {
668// secant_->applyB( *Hp, *p, x );
669// }
670// con.pruneActive(*Hp,xlam,neps_);
671//
672// kappa = p->dot(*Hp);
673// if ( kappa <= 0.0 ) { flagCR_ = 2; break; }
674// alpha = rv/kappa;
675// sol.axpy(alpha,*p);
676//
677// res->axpy(-alpha,*Hp);
678// rnorm = res->norm();
679// if ( rnorm < rtol ) { break; }
680//
681// con.pruneActive(*res,xlam,neps_);
682// obj.precond(*v,*res,x,itol_);
683// con.pruneActive(*v,xlam,neps_);
684// tmp = rv;
685// rv = v->dot(*res);
686// beta = rv/tmp;
687//
688// p->scale(beta);
689// p->axpy(1.0,*v);
690// }
691// if ( iterCR_ == maxitCR_ ) {
692// flagCR_ = 1;
693// }
694// else {
695// iterCR_++;
696// }
697// }
698
699
700// /** \brief Apply the inactive components of the Hessian operator.
701//
702// I.e., the components corresponding to \f$\mathcal{I}_k\f$.
703//
704// @param[out] hv is the result of applying the Hessian at @b x to
705// @b v
706// @param[in] v is the direction in which we apply the Hessian
707// @param[in] x is the current iteration vector \f$x_k\f$
708// @param[in] xlam is the vector \f$x_k + c\lambda_k\f$
709// @param[in] obj is the objective function
710// @param[in] con are the bound constraints
711// */
712// void applyInactiveHessian(Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x,
713// const Vector<Real> &xlam, Objective<Real> &obj, BoundConstraint<Real> &con) {
714// ROL::Ptr<Vector<Real> > tmp = v.clone();
715// tmp->set(v);
716// con.pruneActive(*tmp,xlam,neps_);
717// if ( secant_ == ROL::nullPtr ) {
718// obj.hessVec(hv,*tmp,x,itol_);
719// }
720// else {
721// secant_->applyB(hv,*tmp,x);
722// }
723// con.pruneActive(hv,xlam,neps_);
724// }
725//
726// /** \brief Apply the inactive components of the preconditioner operator.
727//
728// I.e., the components corresponding to \f$\mathcal{I}_k\f$.
729//
730// @param[out] hv is the result of applying the preconditioner at @b x to
731// @b v
732// @param[in] v is the direction in which we apply the preconditioner
733// @param[in] x is the current iteration vector \f$x_k\f$
734// @param[in] xlam is the vector \f$x_k + c\lambda_k\f$
735// @param[in] obj is the objective function
736// @param[in] con are the bound constraints
737// */
738// void applyInactivePrecond(Vector<Real> &pv, const Vector<Real> &v, const Vector<Real> &x,
739// const Vector<Real> &xlam, Objective<Real> &obj, BoundConstraint<Real> &con) {
740// ROL::Ptr<Vector<Real> > tmp = v.clone();
741// tmp->set(v);
742// con.pruneActive(*tmp,xlam,neps_);
743// obj.precond(pv,*tmp,x,itol_);
744// con.pruneActive(pv,xlam,neps_);
745// }
746//
747// /** \brief Solve the inactive part of the PDAS optimality system.
748//
749// The inactive PDAS optimality system is
750// \f[
751// \nabla^2 f(x_k)_{\mathcal{I}_k,\mathcal{I}_k}s =
752// -\nabla f(x_k)_{\mathcal{I}_k}
753// -\nabla^2 f(x_k)_{\mathcal{I}_k,\mathcal{A}_k} (s_k)_{\mathcal{A}_k}.
754// \f]
755// Since the inactive part of the Hessian may not be positive definite, we solve
756// using CR.
757//
758// @param[out] sol is the vector containing the solution
759// @param[in] rhs is the right-hand side vector
760// @param[in] xlam is the vector \f$x_k + c\lambda_k\f$
761// @param[in] x is the current iteration vector \f$x_k\f$
762// @param[in] obj is the objective function
763// @param[in] con are the bound constraints
764// */
765// // Solve the inactive part of the optimality system using conjugate residuals
766// void solve(Vector<Real> &sol, const Vector<Real> &rhs, const Vector<Real> &xlam, const Vector<Real> &x,
767// Objective<Real> &obj, BoundConstraint<Real> &con) {
768// // Initialize Residual
769// ROL::Ptr<Vector<Real> > res = rhs.clone();
770// res->set(rhs);
771// Real rnorm = res->norm();
772// Real rtol = std::min(tol1_,tol2_*rnorm);
773// if ( false ) { itol_ = rtol/(maxitCR_*rnorm); }
774// sol.zero();
775//
776// // Apply preconditioner to residual r = Mres
777// ROL::Ptr<Vector<Real> > r = x.clone();
778// applyInactivePrecond(*r,*res,x,xlam,obj,con);
779//
780// // Initialize direction p = v
781// ROL::Ptr<Vector<Real> > p = x.clone();
782// p->set(*r);
783//
784// // Apply Hessian to v
785// ROL::Ptr<Vector<Real> > Hr = x.clone();
786// applyInactiveHessian(*Hr,*r,x,xlam,obj,con);
787//
788// // Apply Hessian to p
789// ROL::Ptr<Vector<Real> > Hp = x.clone();
790// ROL::Ptr<Vector<Real> > MHp = x.clone();
791// Hp->set(*Hr);
792//
793// iterCR_ = 0;
794// flagCR_ = 0;
795//
796// Real kappa = 0.0, beta = 0.0, alpha = 0.0, tmp = 0.0, rHr = Hr->dot(*r);
797//
798// for (iterCR_ = 0; iterCR_ < maxitCR_; iterCR_++) {
799// // Precondition Hp
800// applyInactivePrecond(*MHp,*Hp,x,xlam,obj,con);
801//
802// kappa = Hp->dot(*MHp); // p' H M H p
803// alpha = rHr/kappa; // r' M H M r
804// sol.axpy(alpha,*p); // update step
805// res->axpy(-alpha,*Hp); // residual
806// r->axpy(-alpha,*MHp); // preconditioned residual
807//
808// // recompute rnorm and decide whether or not to exit
809// rnorm = res->norm();
810// if ( rnorm < rtol ) { break; }
811//
812// // Apply Hessian to v
813// itol_ = rtol/(maxitCR_*rnorm);
814// applyInactiveHessian(*Hr,*r,x,xlam,obj,con);
815//
816// tmp = rHr;
817// rHr = Hr->dot(*r);
818// beta = rHr/tmp;
819// p->scale(beta);
820// p->axpy(1.0,*r);
821// Hp->scale(beta);
822// Hp->axpy(1.0,*Hr);
823// }
824// if ( iterCR_ == maxitCR_ ) {
825// flagCR_ = 1;
826// }
827// else {
828// iterCR_++;
829// }
830// }
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Contains definitions of custom data types in ROL.
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 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 pruneUpperActive(Vector< Real > &v, const Vector< Real > &x, Real eps=Real(0))
Set variables to zero if they correspond to the upper -active set.
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.
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 apply a linear operator.
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.
const ROL::Ptr< BoundConstraint< Real > > bnd_
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)
void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
void apply(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply linear operator.
const ROL::Ptr< BoundConstraint< Real > > bnd_
void applyInverse(Vector< Real > &Hv, const Vector< Real > &v, Real &tol) const
Apply inverse of linear operator.
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)
Implements the computation of optimization steps with the Newton primal-dual active set method.
ROL::Ptr< Vector< Real > > xbnd_
Container for primal variable bounds.
Real scale_
Scale for dual variables in the active set, .
ROL::Ptr< Vector< Real > > Ag_
Container for gradient projected onto active set.
ROL::Ptr< Vector< Real > > res_
Container for optimality system residual for quadratic model.
ROL::Ptr< Vector< Real > > xlam_
Container for primal plus dual variables.
ROL::Ptr< Vector< Real > > x0_
Container for initial priaml variables.
ROL::Ptr< Vector< Real > > lambda_
Container for dual variables.
Real gtol_
PDAS gradient stopping tolerance.
ROL::Ptr< Vector< Real > > As_
Container for step projected onto active set.
std::string printHeader(void) const
Print iterate header.
Real stol_
PDAS minimum step size stopping tolerance.
ROL::Ptr< Vector< Real > > rtmp_
Container for temporary right hand side storage.
ROL::Ptr< Secant< Real > > secant_
Secant object.
ROL::Ptr< Vector< Real > > gtmp_
Container for temporary gradient storage.
void update(Vector< Real > &x, const Vector< Real > &s, Objective< Real > &obj, BoundConstraint< Real > &con, AlgorithmState< Real > &algo_state)
Update step, if successful.
ROL::Ptr< Vector< Real > > xtmp_
Container for temporary primal storage.
PrimalDualActiveSetStep(ROL::ParameterList &parlist)
Constructor.
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.
virtual std::string print(AlgorithmState< Real > &algo_state, bool print_header=false) const
Print iterate status.
Real computeCriticalityMeasure(Vector< Real > &x, Objective< Real > &obj, BoundConstraint< Real > &con, Real tol)
Compute the gradient-based criticality measure.
void compute(Vector< Real > &s, const Vector< Real > &x, Objective< Real > &obj, BoundConstraint< Real > &con, AlgorithmState< Real > &algo_state)
Compute step.
std::string printName(void) const
Print step name.
bool feasible_
Flag whether the current iterate is feasible or not.
int maxit_
Maximum number of PDAS iterations.
Provides interface for and implements limited-memory secant operators.
Definition: ROL_Secant.hpp:79
Provides the interface to compute optimization steps.
Definition: ROL_Step.hpp:68
ROL::Ptr< StepState< Real > > getState(void)
Definition: ROL_Step.hpp:73
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 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:226
virtual void plus(const Vector &x)=0
Compute , where .
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.
Real ROL_EPSILON(void)
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:91
ESecant StringToESecant(std::string s)
Definition: ROL_Types.hpp:543
@ SECANT_LBFGS
Definition: ROL_Types.hpp:487
State for algorithm class. Will be used for restarts.
Definition: ROL_Types.hpp:143
ROL::Ptr< Vector< Real > > iterateVec
Definition: ROL_Types.hpp:157