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