ROL
ROL_CauchyPoint.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_TRUSTREGIONFACTORY_H
12 #else
13 
14 #ifndef ROL_CAUCHYPOINT_H
15 #define ROL_CAUCHYPOINT_H
16 
21 #include "ROL_TrustRegion.hpp"
22 #include "ROL_Vector.hpp"
23 #include "ROL_Types.hpp"
24 #include "ROL_ParameterList.hpp"
25 
26 namespace ROL {
27 
28 template<class Real>
29 class CauchyPoint : public TrustRegion<Real> {
30 private:
31 
32  ROL::Ptr<Vector<Real> > g_;
33  ROL::Ptr<Vector<Real> > p_;
34  ROL::Ptr<Vector<Real> > Hp_;
35 
36  Real pRed_;
37  Real eps_;
38  Real alpha_;
39 
40  bool useCGTCP_;
41 
42 public:
43 
44  // Constructor
45  CauchyPoint( ROL::ParameterList &parlist )
46  : TrustRegion<Real>(parlist), pRed_(0), alpha_(-1), useCGTCP_(false) {
47  // Unravel Parameter List
48  Real oe2(100);
49  Real TRsafe = parlist.sublist("Step").sublist("Trust Region").get("Safeguard Size",oe2);
50  eps_ = TRsafe*ROL_EPSILON<Real>();
51  }
52 
53  void initialize( const Vector<Real> &x, const Vector<Real> &s, const Vector<Real> &g) {
55  Hp_ = g.clone();
56  p_ = s.clone();
57 // if ( useCGTCP_ ) {
58 // g_ = g.clone();
59 // p_ = s.clone();
60 // }
61  }
62 
63  void run( Vector<Real> &s,
64  Real &snorm,
65  int &iflag,
66  int &iter,
67  const Real del,
68  TrustRegionModel<Real> &model) {
69  //if ( pObj.isConActivated() ) {
70  // if ( useCGTCP_ ) {
71  // cauchypoint_CGT( s, snorm, iflag, iter, del, model );
72  // }
73  // else {
74  // cauchypoint_M( s, snorm, iflag, iter, del, model );
75  // }
76  //}
77  //else {
78  // cauchypoint_unc( s, snorm, iflag, iter, del, model );
79  //}
80  cauchypoint_unc( s, snorm, iflag, iter, del, model );
82  }
83 
84 private:
85  void cauchypoint_unc( Vector<Real> &s,
86  Real &snorm,
87  int &iflag,
88  int &iter,
89  const Real del,
90  TrustRegionModel<Real> &model) {
91  Real tol = std::sqrt(ROL_EPSILON<Real>());
92  // Set step to (projected) gradient
93  model.dualTransform(*Hp_,*model.getGradient());
94  s.set(Hp_->dual());
95  // Apply (reduced) Hessian to (projected) gradient
96  model.hessVec(*Hp_,s,s,tol);
97  Real gBg = Hp_->dot(s.dual());
98  Real gnorm = s.dual().norm();
99  Real gg = gnorm*gnorm;
100  Real alpha = del/gnorm;
101  if ( gBg > ROL_EPSILON<Real>() ) {
102  alpha = std::min(gg/gBg, del/gnorm);
103  }
104 
105  s.scale(-alpha);
106  model.primalTransform(*p_,s);
107  s.set(*p_);
108  snorm = s.norm(); //alpha*gnorm;
109  iflag = 0;
110  iter = 0;
111  pRed_ = alpha*(gg - static_cast<Real>(0.5)*alpha*gBg);
112  }
113 
114 // void cauchypoint_M( Vector<Real> &s,
115 // Real &snorm,
116 // int &iflag,
117 // int &iter,
118 // const Real del,
119 // const Vector<Real> &x,
120 // TrustRegionModel<Real> &model,
121 // BoundConstraint<Real> &bnd) {
122 // Real tol = std::sqrt(ROL_EPSILON<Real>()),
123 // const Real zero(0), half(0.5), oe4(1.e4), two(2);
124 // // Parameters
125 // Real mu0(1.e-2), mu1(1), beta1(0), beta2(0);
126 // bool decr = true;
127 // bool stat = true;
128 // // Initial step length
129 // Real alpha = (alpha_ > zero ? alpha_ : one);
130 // Real alpha0 = alpha;
131 // Real alphamax = oe4*alpha;
132 // // Set step to (projected) gradient
133 // s.zero();
134 // model.gradient(*Hp_,s,tol);
135 // s.set(Hp_->dual());
136 // // Initial model value
137 // s.scale(-alpha);
138 // bnd.computeProjectedStep(s,x);
139 // snorm = s.norm();
140 // Real val = model.value(s,tol);
141 // Real val0 = val;
142 //
143 // // Determine whether to increase or decrease alpha
144 // if ( val > mu0 * gs || snorm > mu1 * del ) {
145 // beta1 = half;
146 // beta2 = half;
147 // decr = true;
148 // }
149 // else {
150 // beta1 = two;
151 // beta2 = two;
152 // decr = false;
153 // }
154 //
155 // while ( stat ) {
156 // // Update step length
157 // alpha0 = alpha;
158 // val0 = val;
159 // alpha *= half*(beta1+beta2);
160 //
161 // // Update model value
162 // s.set(grad.dual());
163 // s.scale(-alpha);
164 // pObj.computeProjectedStep(s,x);
165 // snorm = s.norm();
166 // pObj.hessVec(*Hp_,s,x,tol);
167 // gs = s.dot(grad.dual());
168 // val = gs + half*s.dot(Hp_->dual());
169 //
170 // // Update termination criterion
171 // if ( decr ) {
172 // stat = ( val > mu0 * gs || snorm > mu1 * del );
173 // if ( std::abs(val) < eps_ && std::abs(mu0 *gs) < eps_ ) {
174 // stat = (snorm > mu1 * del);
175 // }
176 // }
177 // else {
178 // stat = !( val > mu0 * gs || snorm > mu1 * del );
179 // if ( std::abs(val) < eps_ && std::abs(mu0 *gs) < eps_ ) {
180 // stat = !(snorm > mu1 * del);
181 // }
182 // if ( alpha > alphamax ) {
183 // stat = false;
184 // }
185 // }
186 // }
187 // // Reset to last 'successful' step
188 // val = val0;
189 // alpha = alpha0;
190 // s.set(grad.dual());
191 // s.scale(-alpha);
192 // pObj.computeProjectedStep(s,x);
193 // snorm = s.norm();
194 //
195 // alpha_ = alpha;
196 // pRed_ = -val;
197 // }
198 //
199 // void cauchypoint_CGT( Vector<Real> &s, Real &snorm, Real &del, int &iflag, int &iter, const Vector<Real> &x,
200 // const Vector<Real> &grad, const Real &gnorm, ProjectedObjective<Real> &pObj ) {
201 // Real tol = std::sqrt(ROL_EPSILON<Real>()), one(1), half(0.5), two(2);
202 // bool tmax_flag = true;
203 // int maxit = 20;
204 // Real t = del/gnorm;
205 // Real tmax(1.e10), tmin(0), gs(0), pgnorm(0);
206 // Real c1(0.25), c2(0.75), c3(0.9), c4(0.25);
207 // for ( int i = 0; i < maxit; i++ ) {
208 // // Compute p = x + s = P(x - t*g)
209 // p_->set(x);
210 // p_->axpy(-t,grad.dual());
211 // pObj.project(*p_);
212 // // Compute s = p - x = P(x - t*g) - x
213 // s.set(*p_);
214 // s.axpy(-one,x);
215 // snorm = s.norm();
216 // // Evaluate Model
217 // pObj.hessVec(*Hp_,s,x,tol);
218 // gs = s.dot(grad.dual());
219 // pRed_ = -gs - half*s.dot(Hp_->dual());
220 //
221 // // Check Stopping Conditions
222 // g_->set(grad);
223 // pObj.pruneActive(*g_,grad,*p_); // Project gradient onto tangent cone at p
224 // pgnorm = g_->norm();
225 // if ( snorm > del || pRed_ < -c2*gs ) {
226 // tmax = t;
227 // tmax_flag = false;
228 // }
229 // else if ( snorm < c3*del && pRed_ > -c1*gs && pgnorm > c4*std::abs(gs)/del ) {
230 // tmin = t;
231 // }
232 // else {
233 // break;
234 // }
235 //
236 // // Update t
237 // if ( tmax_flag ) {
238 // t *= two;
239 // }
240 // else {
241 // t = half*(tmax + tmin);
242 // }
243 // }
244 // }
245 };
246 
247 }
248 
249 #endif
250 #endif
virtual void initialize(const Vector< Real > &x, const Vector< Real > &s, const Vector< Real > &g)
Contains definitions of custom data types in ROL.
void setPredictedReduction(const Real pRed)