ROL
ROL_DoubleDogLeg.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 
10 #ifndef ROL_DOUBLEDOGLEG_H
11 #define ROL_DOUBLEDOGLEG_H
12 
17 #include "ROL_TrustRegion.hpp"
18 #include "ROL_Types.hpp"
19 
20 namespace ROL {
21 
22 template<class Real>
23 class DoubleDogLeg : public TrustRegion<Real> {
24 private:
25 
26  ROL::Ptr<CauchyPoint<Real> > cpt_;
27 
28  ROL::Ptr<Vector<Real> > s_;
29  ROL::Ptr<Vector<Real> > v_;
30  ROL::Ptr<Vector<Real> > Hp_;
31 
32  Real pRed_;
33 
34 public:
35 
36  // Constructor
37  DoubleDogLeg( ROL::ParameterList &parlist ) : TrustRegion<Real>(parlist), pRed_(0) {
38  cpt_ = ROL::makePtr<CauchyPoint<Real>>(parlist);
39  }
40 
41  void initialize( const Vector<Real> &x, const Vector<Real> &s, const Vector<Real> &g) {
43  cpt_->initialize(x,s,g);
44  s_ = s.clone();
45  v_ = s.clone();
46  Hp_ = g.clone();
47  }
48 
49  void run( Vector<Real> &s,
50  Real &snorm,
51  int &iflag,
52  int &iter,
53  const Real del,
54  TrustRegionModel<Real> &model ) {
55  Real tol = std::sqrt(ROL_EPSILON<Real>());
56  const Real one(1), zero(0), half(0.5), p2(0.2), p8(0.8), two(2);
57  // Set s to be the (projected) gradient
58  model.dualTransform(*Hp_,*model.getGradient());
59  s.set(Hp_->dual());
60  // Compute (quasi-)Newton step
61  model.invHessVec(*s_,*Hp_,s,tol);
62  Real sNnorm = s_->norm();
63  Real tmp = -s_->dot(s);
64  bool negCurv = (tmp > zero ? true : false);
65  Real gsN = std::abs(tmp);
66  // Check if (quasi-)Newton step is feasible
67  if ( negCurv ) {
68  // Use Cauchy point
69  cpt_->run(s,snorm,iflag,iter,del,model);
70  pRed_ = cpt_->getPredictedReduction();
71  iflag = 2;
72  }
73  else {
74  // Approximately solve trust region subproblem using double dogleg curve
75  if (sNnorm <= del) { // Use the (quasi-)Newton step
76  s.set(*s_);
77  s.scale(-one);
78  snorm = sNnorm;
79  pRed_ = half*gsN;
80  iflag = 0;
81  }
82  else { // The (quasi-)Newton step is outside of trust region
83  model.hessVec(*Hp_,s,s,tol);
84  Real alpha = zero;
85  Real beta = zero;
86  Real gnorm = s.norm();
87  Real gnorm2 = gnorm*gnorm;
88  Real gBg = Hp_->dot(s.dual());
89  Real gamma1 = gnorm/gBg;
90  Real gamma2 = gnorm/gsN;
91  Real eta = p8*gamma1*gamma2 + p2;
92  if (eta*sNnorm <= del || gBg <= zero) { // Dogleg Point is inside trust region
93  alpha = del/sNnorm;
94  beta = zero;
95  s.set(*s_);
96  s.scale(-alpha);
97  snorm = del;
98  iflag = 1;
99  }
100  else {
101  if (gnorm2*gamma1 >= del) { // Cauchy Point is outside trust region
102  alpha = zero;
103  beta = -del/gnorm;
104  s.scale(beta);
105  snorm = del;
106  iflag = 2;
107  }
108  else { // Find convex combination of Cauchy and Dogleg point
109  s.scale(-gamma1*gnorm);
110  v_->set(s);
111  v_->axpy(eta,*s_);
112  v_->scale(-one);
113  Real wNorm = v_->dot(*v_);
114  Real sigma = del*del-std::pow(gamma1*gnorm,two);
115  Real phi = s.dot(*v_);
116  Real theta = (-phi + std::sqrt(phi*phi+wNorm*sigma))/wNorm;
117  s.axpy(theta,*v_);
118  snorm = del;
119  alpha = theta*eta;
120  beta = (one-theta)*(-gamma1*gnorm);
121  iflag = 3;
122  }
123  }
124  pRed_ = -(alpha*(half*alpha-one)*gsN + half*beta*beta*gBg + beta*(one-alpha)*gnorm2);
125  }
126  }
127  model.primalTransform(*s_,s);
128  s.set(*s_);
129  snorm = s.norm();
131  }
132 };
133 
134 }
135 
136 #endif
ROL::Ptr< Vector< Real > > s_
void initialize(const Vector< Real > &x, const Vector< Real > &s, const Vector< Real > &g)
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
virtual void scale(const Real alpha)=0
Compute where .
virtual ROL::Ptr< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:119
virtual void initialize(const Vector< Real > &x, const Vector< Real > &s, const Vector< Real > &g)
Contains definitions of custom data types in ROL.
Provides interface for and implements trust-region subproblem solvers.
Provides the interface to evaluate trust-region model functions.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
virtual const Ptr< const Vector< Real > > getGradient(void) const
virtual Real dot(const Vector &x) const =0
Compute where .
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
virtual void dualTransform(Vector< Real > &tv, const Vector< Real > &v)
Provides interface for the double dog leg trust-region subproblem solver.
void setPredictedReduction(const Real pRed)
ROL::Ptr< Vector< Real > > v_
DoubleDogLeg(ROL::ParameterList &parlist)
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &s, Real &tol)
Apply Hessian approximation to vector.
virtual void invHessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &s, Real &tol)
Apply inverse Hessian approximation to vector.
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:175
virtual Real norm() const =0
Returns where .
void run(Vector< Real > &s, Real &snorm, int &iflag, int &iter, const Real del, TrustRegionModel< Real > &model)
ROL::Ptr< Vector< Real > > Hp_
ROL::Ptr< CauchyPoint< Real > > cpt_
virtual void primalTransform(Vector< Real > &tv, const Vector< Real > &v)