ROL
ROL_DogLeg_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_DOGLEG_U_H
11 #define ROL_DOGLEG_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 DogLeg_U : public TrustRegion_U<Real> {
24 private:
25 
26  Ptr<Vector<Real>> primal_, dual_;
27 
28 public:
29 
30  // Constructor
31  DogLeg_U() {}
32 
33  void initialize(const Vector<Real> &x, const Vector<Real> &g) {
34  primal_ = x.clone();
35  dual_ = g.clone();
36  }
37 
38  void solve( Vector<Real> &s,
39  Real &snorm,
40  Real &pRed,
41  int &iflag,
42  int &iter,
43  const Real del,
44  TrustRegionModel_U<Real> &model ) {
45  Real tol = std::sqrt(ROL_EPSILON<Real>());
46  const Real zero(0), half(0.5), one(1), two(2);
47  iter = 0;
48  // Set s to be the gradient
49  s.set(model.getGradient()->dual());
50  // Compute (quasi-)Newton step
51  model.invHessVec(*primal_,*model.getGradient(),s,tol);
52  Real sNnorm = primal_->norm();
53  Real gsN = -primal_->dot(s);
54  // Check if (quasi-)Newton step is feasible
55  if ( gsN >= zero ) {
56  // Use the Cauchy point
57  model.hessVec(*dual_,s,s,tol);
58  Real gnorm = s.norm();
59  Real gnorm2 = gnorm*gnorm;
60  //Real gBg = dual_->dot(s.dual());
61  Real gBg = dual_->apply(s);
62  Real alpha = gnorm2/gBg;
63  if ( alpha*gnorm >= del || gBg <= zero ) {
64  alpha = del/gnorm;
65  }
66  s.scale(-alpha);
67  snorm = alpha*gnorm;
68  iflag = 2;
69  pRed = alpha*(gnorm2 - half*alpha*gBg);
70  }
71  else {
72  // Approximately solve trust region subproblem using double dogleg curve
73  if (sNnorm <= del) { // Use the (quasi-)Newton step
74  s.set(*primal_);
75  s.scale(-one);
76  snorm = sNnorm;
77  pRed = -half*gsN;
78  iflag = 0;
79  }
80  else { // The (quasi-)Newton step is outside of trust region
81  model.hessVec(*dual_,s,s,tol);
82  Real alpha = zero;
83  Real beta = zero;
84  Real gnorm = s.norm();
85  Real gnorm2 = gnorm*gnorm;
86  //Real gBg = dual_->dot(s.dual());
87  Real gBg = dual_->apply(s);
88  Real gamma = gnorm2/gBg;
89  if ( gamma*gnorm >= del || gBg <= zero ) {
90  // Use Cauchy point
91  alpha = zero;
92  beta = del/gnorm;
93  s.scale(-beta);
94  snorm = del;
95  iflag = 2;
96  }
97  else {
98  // Use a convex combination of Cauchy point and (quasi-)Newton step
99  Real a = sNnorm*sNnorm + two*gamma*gsN + gamma*gamma*gnorm2;
100  Real b = -gamma*gsN - gamma*gamma*gnorm2;
101  Real c = gamma*gamma*gnorm2 - del*del;
102  alpha = (-b + std::sqrt(b*b - a*c))/a;
103  beta = gamma*(one-alpha);
104  s.scale(-beta);
105  s.axpy(-alpha,*primal_);
106  snorm = del;
107  iflag = 1;
108  }
109  pRed = (alpha*(half*alpha-one)*gsN - half*beta*beta*gBg + beta*(one-alpha)*gnorm2);
110  }
111  }
112  }
113 };
114 
115 } // namespace ROL
116 
117 #endif
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
Provides interface for dog leg trust-region subproblem solver.
Contains definitions of custom data types in ROL.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
void solve(Vector< Real > &s, Real &snorm, Real &pRed, int &iflag, int &iter, const Real del, TrustRegionModel_U< Real > &model)
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &s, Real &tol) override
Apply Hessian approximation to vector.
Provides the interface to evaluate trust-region model functions.
virtual const Ptr< const Vector< Real > > getGradient(void) const
void initialize(const Vector< Real > &x, const Vector< Real > &g)
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 .
Ptr< Vector< Real > > primal_
Ptr< Vector< Real > > dual_