ROL
ROL_TruncatedCG.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_TRUNCATEDCG_H
11 #define ROL_TRUNCATEDCG_H
12 
17 #include "ROL_TrustRegion.hpp"
18 #include "ROL_Types.hpp"
19 
20 namespace ROL {
21 
22 template<class Real>
23 class TruncatedCG : public TrustRegion<Real> {
24 private:
25  ROL::Ptr<Vector<Real> > primalVector_;
26 
27  ROL::Ptr<Vector<Real> > s_;
28  ROL::Ptr<Vector<Real> > g_;
29  ROL::Ptr<Vector<Real> > v_;
30  ROL::Ptr<Vector<Real> > p_;
31  ROL::Ptr<Vector<Real> > Hp_;
32 
33  int maxit_;
34  Real tol1_;
35  Real tol2_;
36 
37  Real pRed_;
38 
39 public:
40 
41  // Constructor
42  TruncatedCG( ROL::ParameterList &parlist ) : TrustRegion<Real>(parlist), pRed_(0) {
43  // Unravel Parameter List
44  Real em4(1e-4), em2(1e-2);
45  maxit_ = parlist.sublist("General").sublist("Krylov").get("Iteration Limit",20);
46  tol1_ = parlist.sublist("General").sublist("Krylov").get("Absolute Tolerance",em4);
47  tol2_ = parlist.sublist("General").sublist("Krylov").get("Relative Tolerance",em2);
48  }
49 
50  void initialize( const Vector<Real> &x, const Vector<Real> &s, const Vector<Real> &g) {
52 
53  primalVector_ = x.clone();
54 
55  s_ = s.clone();
56  g_ = g.clone();
57  v_ = s.clone();
58  p_ = s.clone();
59  Hp_ = g.clone();
60  }
61 
62  void run( Vector<Real> &s,
63  Real &snorm,
64  int &iflag,
65  int &iter,
66  const Real del,
67  TrustRegionModel<Real> &model ) {
68  Real tol = std::sqrt(ROL_EPSILON<Real>());
69  const Real zero(0), one(1), two(2), half(0.5);
70  // Initialize step
71  s.zero(); s_->zero();
72  snorm = zero;
73  Real snorm2(0), s1norm2(0);
74  // Compute (projected) gradient
75  model.dualTransform(*g_,*model.getGradient());
76  Real gnorm = g_->norm(), normg = gnorm;
77  const Real gtol = std::min(tol1_,tol2_*gnorm);
78  // Preconditioned (projected) gradient vector
79  model.precond(*v_,*g_,s,tol);
80  // Initialize basis vector
81  p_->set(*v_); p_->scale(-one);
82  Real pnorm2 = v_->dot(g_->dual());
83  if ( pnorm2 <= zero ) {
84  iflag = 4;
85  iter = 0;
86  return;
87  }
88  // Initialize scalar storage
89  iter = 0; iflag = 0;
90  Real kappa(0), beta(0), sigma(0), alpha(0), tmp(0), sMp(0);
91  Real gv = v_->dot(g_->dual());
92  pRed_ = zero;
93  // Iterate CG
94  for (iter = 0; iter < maxit_; iter++) {
95  // Apply Hessian to direction p
96  model.hessVec(*Hp_,*p_,s,tol);
97  // Check positivity of Hessian
98  kappa = p_->dot(Hp_->dual());
99  if (kappa <= zero) {
100  sigma = (-sMp+sqrt(sMp*sMp+pnorm2*(del*del-snorm2)))/pnorm2;
101  s.axpy(sigma,*p_);
102  iflag = 2;
103  break;
104  }
105  // Update step
106  alpha = gv/kappa;
107  s_->set(s);
108  s_->axpy(alpha,*p_);
109  s1norm2 = snorm2 + two*alpha*sMp + alpha*alpha*pnorm2;
110  // Check if step exceeds trust region radius
111  if (s1norm2 >= del*del) {
112  sigma = (-sMp+sqrt(sMp*sMp+pnorm2*(del*del-snorm2)))/pnorm2;
113  s.axpy(sigma,*p_);
114  iflag = 3;
115  break;
116  }
117  // Update model predicted reduction
118  pRed_ += half*alpha*gv;
119  // Set step to temporary step and store norm
120  s.set(*s_);
121  snorm2 = s1norm2;
122  // Check for convergence
123  g_->axpy(alpha,*Hp_);
124  normg = g_->norm();
125  if (normg < gtol) {
126  break;
127  }
128  // Preconditioned updated (projected) gradient vector
129  model.precond(*v_,*g_,s,tol);
130  tmp = gv;
131  gv = v_->dot(g_->dual());
132  beta = gv/tmp;
133  // Update basis vector
134  p_->scale(beta);
135  p_->axpy(-one,*v_);
136  sMp = beta*(sMp+alpha*pnorm2);
137  pnorm2 = gv + beta*beta*pnorm2;
138  }
139  // Update model predicted reduction
140  if (iflag > 0) {
141  pRed_ += sigma*(gv-half*sigma*kappa);
142  }
143  // Check iteration count
144  if (iter == maxit_) {
145  iflag = 1;
146  }
147  if (iflag != 1) {
148  iter++;
149  }
150  // Update norm of step and update model predicted reduction
151  model.primalTransform(*s_,s);
152  s.set(*s_);
153  snorm = s.norm();
155  }
156 
157 #if 0
158  void truncatedCG_proj( Vector<Real> &s, Real &snorm, Real &del, int &iflag, int &iter, const Vector<Real> &x,
159  const Vector<Real> &grad, const Real &gnorm, ProjectedObjective<Real> &pObj ) {
160  Real tol = std::sqrt(ROL_EPSILON<Real>());
161 
162  const Real gtol = std::min(tol1_,tol2_*gnorm);
163 
164  // Compute Cauchy Point
165  Real scnorm = 0.0;
166  ROL::Ptr<Vector<Real> > sc = x.clone();
167  cauchypoint(*sc,scnorm,del,iflag,iter,x,grad,gnorm,pObj);
168  ROL::Ptr<Vector<Real> > xc = x.clone();
169  xc->set(x);
170  xc->plus(*sc);
171 
172  // Old and New Step Vectors
173  s.set(*sc);
174  snorm = s.norm();
175  Real snorm2 = snorm*snorm;
176  ROL::Ptr<Vector<Real> > s1 = x.clone();
177  s1->zero();
178  Real s1norm2 = 0.0;
179 
180  // Gradient Vector
181  ROL::Ptr<Vector<Real> > g = x.clone();
182  g->set(grad);
183  ROL::Ptr<Vector<Real> > Hs = x.clone();
184  pObj.reducedHessVec(*Hs,s,*xc,x,tol);
185  g->plus(*Hs);
186  Real normg = g->norm();
187 
188  // Preconditioned Gradient Vector
189  ROL::Ptr<Vector<Real> > v = x.clone();
190  pObj.reducedPrecond(*v,*g,*xc,x,tol);
191 
192  // Basis Vector
193  ROL::Ptr<Vector<Real> > p = x.clone();
194  p->set(*v);
195  p->scale(-1.0);
196  Real pnorm2 = v->dot(*g);
197 
198  // Hessian Times Basis Vector
199  ROL::Ptr<Vector<Real> > Hp = x.clone();
200 
201  iter = 0;
202  iflag = 0;
203  Real kappa = 0.0;
204  Real beta = 0.0;
205  Real sigma = 0.0;
206  Real alpha = 0.0;
207  Real tmp = 0.0;
208  Real gv = v->dot(*g);
209  Real sMp = 0.0;
210  pRed_ = 0.0;
211 
212  for (iter = 0; iter < maxit_; iter++) {
213  pObj.reducedHessVec(*Hp,*p,*xc,x,tol);
214 
215  kappa = p->dot(*Hp);
216  if (kappa <= 0) {
217  sigma = (-sMp+sqrt(sMp*sMp+pnorm2*(del*del-snorm2)))/pnorm2;
218  s.axpy(sigma,*p);
219  iflag = 2;
220  break;
221  }
222 
223  alpha = gv/kappa;
224  s1->set(s);
225  s1->axpy(alpha,*p);
226  s1norm2 = snorm2 + 2.0*alpha*sMp + alpha*alpha*pnorm2;
227 
228  if (s1norm2 >= del*del) {
229  sigma = (-sMp+sqrt(sMp*sMp+pnorm2*(del*del-snorm2)))/pnorm2;
230  s.axpy(sigma,*p);
231  iflag = 3;
232  break;
233  }
234 
235  pRed_ += 0.5*alpha*gv;
236 
237  s.set(*s1);
238  snorm2 = s1norm2;
239 
240  g->axpy(alpha,*Hp);
241  normg = g->norm();
242  if (normg < gtol) {
243  break;
244  }
245 
246  pObj.reducedPrecond(*v,*g,*xc,x,tol);
247  tmp = gv;
248  gv = v->dot(*g);
249  beta = gv/tmp;
250 
251  p->scale(beta);
252  p->axpy(-1.0,*v);
253  sMp = beta*(sMp+alpha*pnorm2);
254  pnorm2 = gv + beta*beta*pnorm2;
255  }
256  if (iflag > 0) {
257  pRed_ += sigma*(gv-0.5*sigma*kappa);
258  }
259 
260  if (iter == maxit_) {
261  iflag = 1;
262  }
263  if (iflag != 1) {
264  iter++;
265  }
266 
267  snorm = s.norm();
268  }
269 #endif
270 };
271 
272 }
273 
274 #endif
ROL::Ptr< Vector< Real > > v_
void initialize(const Vector< Real > &x, const Vector< Real > &s, const Vector< Real > &g)
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)
void reducedHessVec(Vector< Real > &Hv, const Vector< Real > &v, const Vector< Real > &p, const Vector< Real > &d, const Vector< Real > &x, Real &tol)
Apply the reduced Hessian to a vector, v. The reduced Hessian first removes elements of v correspondi...
Contains definitions of custom data types in ROL.
ROL::Ptr< Vector< Real > > s_
Provides interface for and implements trust-region subproblem solvers.
Provides the interface to evaluate trust-region model functions.
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:133
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:46
virtual const Ptr< const Vector< Real > > getGradient(void) const
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
virtual void dualTransform(Vector< Real > &tv, const Vector< Real > &v)
void setPredictedReduction(const Real pRed)
TruncatedCG(ROL::ParameterList &parlist)
ROL::Ptr< Vector< Real > > g_
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &s, Real &tol)
Apply Hessian approximation to vector.
void run(Vector< Real > &s, Real &snorm, int &iflag, int &iter, const Real del, TrustRegionModel< Real > &model)
ROL::Ptr< Vector< Real > > Hp_
Provides interface for truncated CG trust-region subproblem solver.
ROL::Ptr< Vector< Real > > primalVector_
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:175
virtual Real norm() const =0
Returns where .
void reducedPrecond(Vector< Real > &Mv, const Vector< Real > &v, const Vector< Real > &p, const Vector< Real > &d, const Vector< Real > &x, Real &tol)
Apply the reduced preconditioner to a vector, v. The reduced preconditioner first removes elements of...
ROL::Ptr< Vector< Real > > p_
virtual void precond(Vector< Real > &Pv, const Vector< Real > &v, const Vector< Real > &s, Real &tol)
Apply preconditioner to vector.
virtual void primalTransform(Vector< Real > &tv, const Vector< Real > &v)