ROL
ROL_TruncatedCG_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_TRUNCATEDCG_U_H
11 #define ROL_TRUNCATEDCG_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 TruncatedCG_U : public TrustRegion_U<Real> {
24 private:
25  Ptr<Vector<Real>> s_, g_, v_, p_, Hp_;
26 
27  int maxit_;
28  Real tol1_;
29  Real tol2_;
30 
31 public:
32 
33  // Constructor
34  TruncatedCG_U( ParameterList &parlist ) {
35  // Unravel Parameter List
36  Real em4(1e-4), em2(1e-2);
37  maxit_ = parlist.sublist("General").sublist("Krylov").get("Iteration Limit",20);
38  tol1_ = parlist.sublist("General").sublist("Krylov").get("Absolute Tolerance",em4);
39  tol2_ = parlist.sublist("General").sublist("Krylov").get("Relative Tolerance",em2);
40  }
41 
42  void initialize(const Vector<Real> &x, const Vector<Real> &g) {
43  s_ = x.clone();
44  v_ = x.clone();
45  p_ = x.clone();
46  g_ = g.clone();
47  Hp_ = g.clone();
48  }
49 
50  void solve( Vector<Real> &s,
51  Real &snorm,
52  Real &pRed,
53  int &iflag,
54  int &iter,
55  const Real del,
56  TrustRegionModel_U<Real> &model ) {
57  Real tol = std::sqrt(ROL_EPSILON<Real>());
58  const Real zero(0), one(1), two(2), half(0.5);
59  // Initialize step
60  s.zero(); s_->zero();
61  snorm = zero;
62  Real snorm2(0), s1norm2(0);
63  // Compute (projected) gradient
64  g_->set(*model.getGradient());
65  Real gnorm = g_->norm(), normg = gnorm;
66  const Real gtol = std::min(tol1_,tol2_*gnorm);
67  // Preconditioned (projected) gradient vector
68  model.precond(*v_,*g_,s,tol);
69  // Initialize basis vector
70  p_->set(*v_); p_->scale(-one);
71  Real pnorm2 = v_->apply(*g_);
72  if ( pnorm2 <= zero ) {
73  iflag = 4;
74  iter = 0;
75  return;
76  }
77  // Initialize scalar storage
78  iter = 0; iflag = 0;
79  Real kappa(0), beta(0), sigma(0), alpha(0), tmp(0), sMp(0);
80  Real gv = pnorm2;
81  pRed = zero;
82  // Iterate CG
83  for (iter = 0; iter < maxit_; iter++) {
84  // Apply Hessian to direction p
85  model.hessVec(*Hp_,*p_,s,tol);
86  // Check positivity of Hessian
87  kappa = p_->apply(*Hp_);
88  if (kappa <= zero) {
89  sigma = (-sMp+sqrt(sMp*sMp+pnorm2*(del*del-snorm2)))/pnorm2;
90  s.axpy(sigma,*p_);
91  snorm = del;
92  iflag = 2;
93  break;
94  }
95  // Update step
96  alpha = gv/kappa;
97  s_->set(s);
98  s_->axpy(alpha,*p_);
99  s1norm2 = snorm2 + two*alpha*sMp + alpha*alpha*pnorm2;
100  // Check if step exceeds trust region radius
101  if (s1norm2 >= del*del) {
102  sigma = (-sMp+sqrt(sMp*sMp+pnorm2*(del*del-snorm2)))/pnorm2;
103  s.axpy(sigma,*p_);
104  snorm = del;
105  iflag = 3;
106  break;
107  }
108  // Update model predicted reduction
109  pRed += half*alpha*gv;
110  // Set step to temporary step and store norm
111  s.set(*s_);
112  snorm2 = s1norm2;
113  // Check for convergence
114  g_->axpy(alpha,*Hp_);
115  normg = g_->norm();
116  if (normg < gtol) break;
117  // Preconditioned updated (projected) gradient vector
118  model.precond(*v_,*g_,s,tol);
119  tmp = gv;
120  gv = v_->apply(*g_);
121  beta = gv/tmp;
122  // Update basis vector
123  p_->scale(beta);
124  p_->axpy(-one,*v_);
125  sMp = beta*(sMp+alpha*pnorm2);
126  pnorm2 = gv + beta*beta*pnorm2;
127  }
128  // Update model predicted reduction
129  if (iflag > 0) pRed += sigma*(gv-half*sigma*kappa);
130  else snorm = std::sqrt(snorm2);
131  // Check iteration count
132  if (iter == maxit_) iflag = 1;
133  if (iflag != 1) iter++;
134  }
135 };
136 
137 } // namespace ROL
138 
139 #endif
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
void initialize(const Vector< Real > &x, const Vector< Real > &g)
TruncatedCG_U(ParameterList &parlist)
Contains definitions of custom data types in ROL.
Ptr< Vector< Real > > Hp_
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
Ptr< Vector< Real > > v_
Objective_SerialSimOpt(const Ptr< Obj > &obj, const V &ui) z0_ zero()
Ptr< Vector< Real > > g_
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 > > s_
Provides interface for and implements trust-region subproblem solvers.
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:175
Provides interface for truncated CG trust-region subproblem solver.
virtual void precond(Vector< Real > &Pv, const Vector< Real > &v, const Vector< Real > &s, Real &tol) override
Apply preconditioner to vector.
Ptr< Vector< Real > > p_