ROL
ROL_TruncatedCG.hpp
Go to the documentation of this file.
1 // @HEADER
2 // ************************************************************************
3 //
4 // Rapid Optimization Library (ROL) Package
5 // Copyright (2014) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact lead developers:
38 // Drew Kouri (dpkouri@sandia.gov) and
39 // Denis Ridzal (dridzal@sandia.gov)
40 //
41 // ************************************************************************
42 // @HEADER
43 
44 #ifndef ROL_TRUNCATEDCG_H
45 #define ROL_TRUNCATEDCG_H
46 
51 #include "ROL_TrustRegion.hpp"
52 #include "ROL_Types.hpp"
53 #include "ROL_HelperFunctions.hpp"
54 
55 namespace ROL {
56 
57 template<class Real>
58 class TruncatedCG : public TrustRegion<Real> {
59 private:
60  Teuchos::RCP<Vector<Real> > primalVector_;
61 
62  Teuchos::RCP<Vector<Real> > s_;
63  Teuchos::RCP<Vector<Real> > g_;
64  Teuchos::RCP<Vector<Real> > v_;
65  Teuchos::RCP<Vector<Real> > p_;
66  Teuchos::RCP<Vector<Real> > Hp_;
67 
68  int maxit_;
69  Real tol1_;
70  Real tol2_;
71 
72  Real pRed_;
73 
74 public:
75 
76  // Constructor
77  TruncatedCG( Teuchos::ParameterList &parlist ) : TrustRegion<Real>(parlist), pRed_(0) {
78  // Unravel Parameter List
79  Real em4(1e-4), em2(1e-2);
80  maxit_ = parlist.sublist("General").sublist("Krylov").get("Iteration Limit",20);
81  tol1_ = parlist.sublist("General").sublist("Krylov").get("Absolute Tolerance",em4);
82  tol2_ = parlist.sublist("General").sublist("Krylov").get("Relative Tolerance",em2);
83  }
84 
85  void initialize( const Vector<Real> &x, const Vector<Real> &s, const Vector<Real> &g) {
87 
88  primalVector_ = x.clone();
89 
90  s_ = s.clone();
91  g_ = g.clone();
92  v_ = s.clone();
93  p_ = s.clone();
94  Hp_ = g.clone();
95  }
96 
97  void run( Vector<Real> &s,
98  Real &snorm,
99  int &iflag,
100  int &iter,
101  const Real del,
102  TrustRegionModel<Real> &model ) {
103  Real tol = std::sqrt(ROL_EPSILON<Real>());
104  const Real zero(0), one(1), two(2), half(0.5);
105  // Initialize step
106  s.zero(); s_->zero();
107  snorm = zero;
108  Real snorm2(0), s1norm2(0);
109  // Compute (projected) gradient
110  model.dualTransform(*g_,*model.getGradient());
111  Real gnorm = g_->norm(), normg = gnorm;
112  const Real gtol = std::min(tol1_,tol2_*gnorm);
113  // Preconditioned (projected) gradient vector
114  model.precond(*v_,*g_,s,tol);
115  // Initialize basis vector
116  p_->set(*v_); p_->scale(-one);
117  Real pnorm2 = v_->dot(g_->dual());
118  // Initialize scalar storage
119  iter = 0; iflag = 0;
120  Real kappa(0), beta(0), sigma(0), alpha(0), tmp(0), sMp(0);
121  Real gv = v_->dot(g_->dual());
122  pRed_ = zero;
123  // Iterate CG
124  for (iter = 0; iter < maxit_; iter++) {
125  // Apply Hessian to direction p
126  model.hessVec(*Hp_,*p_,s,tol);
127  // Check positivity of Hessian
128  kappa = p_->dot(Hp_->dual());
129  if (kappa <= zero) {
130  sigma = (-sMp+sqrt(sMp*sMp+pnorm2*(del*del-snorm2)))/pnorm2;
131  s.axpy(sigma,*p_);
132  iflag = 2;
133  break;
134  }
135  // Update step
136  alpha = gv/kappa;
137  s_->set(s);
138  s_->axpy(alpha,*p_);
139  s1norm2 = snorm2 + two*alpha*sMp + alpha*alpha*pnorm2;
140  // Check if step exceeds trust region radius
141  if (s1norm2 >= del*del) {
142  sigma = (-sMp+sqrt(sMp*sMp+pnorm2*(del*del-snorm2)))/pnorm2;
143  s.axpy(sigma,*p_);
144  iflag = 3;
145  break;
146  }
147  // Update model predicted reduction
148  pRed_ += half*alpha*gv;
149  // Set step to temporary step and store norm
150  s.set(*s_);
151  snorm2 = s1norm2;
152  // Check for convergence
153  g_->axpy(alpha,*Hp_);
154  normg = g_->norm();
155  if (normg < gtol) {
156  break;
157  }
158  // Preconditioned updated (projected) gradient vector
159  model.precond(*v_,*g_,s,tol);
160  tmp = gv;
161  gv = v_->dot(g_->dual());
162  beta = gv/tmp;
163  // Update basis vector
164  p_->scale(beta);
165  p_->axpy(-one,*v_);
166  sMp = beta*(sMp+alpha*pnorm2);
167  pnorm2 = gv + beta*beta*pnorm2;
168  }
169  // Update model predicted reduction
170  if (iflag > 0) {
171  pRed_ += sigma*(gv-half*sigma*kappa);
172  }
173  // Check iteration count
174  if (iter == maxit_) {
175  iflag = 1;
176  }
177  if (iflag != 1) {
178  iter++;
179  }
180  // Update norm of step and update model predicted reduction
181  model.primalTransform(*s_,s);
182  s.set(*s_);
183  snorm = s.norm();
185  }
186 
187 #if 0
188  void truncatedCG_proj( Vector<Real> &s, Real &snorm, Real &del, int &iflag, int &iter, const Vector<Real> &x,
189  const Vector<Real> &grad, const Real &gnorm, ProjectedObjective<Real> &pObj ) {
190  Real tol = std::sqrt(ROL_EPSILON<Real>());
191 
192  const Real gtol = std::min(tol1_,tol2_*gnorm);
193 
194  // Compute Cauchy Point
195  Real scnorm = 0.0;
196  Teuchos::RCP<Vector<Real> > sc = x.clone();
197  cauchypoint(*sc,scnorm,del,iflag,iter,x,grad,gnorm,pObj);
198  Teuchos::RCP<Vector<Real> > xc = x.clone();
199  xc->set(x);
200  xc->plus(*sc);
201 
202  // Old and New Step Vectors
203  s.set(*sc);
204  snorm = s.norm();
205  Real snorm2 = snorm*snorm;
206  Teuchos::RCP<Vector<Real> > s1 = x.clone();
207  s1->zero();
208  Real s1norm2 = 0.0;
209 
210  // Gradient Vector
211  Teuchos::RCP<Vector<Real> > g = x.clone();
212  g->set(grad);
213  Teuchos::RCP<Vector<Real> > Hs = x.clone();
214  pObj.reducedHessVec(*Hs,s,*xc,x,tol);
215  g->plus(*Hs);
216  Real normg = g->norm();
217 
218  // Preconditioned Gradient Vector
219  Teuchos::RCP<Vector<Real> > v = x.clone();
220  pObj.reducedPrecond(*v,*g,*xc,x,tol);
221 
222  // Basis Vector
223  Teuchos::RCP<Vector<Real> > p = x.clone();
224  p->set(*v);
225  p->scale(-1.0);
226  Real pnorm2 = v->dot(*g);
227 
228  // Hessian Times Basis Vector
229  Teuchos::RCP<Vector<Real> > Hp = x.clone();
230 
231  iter = 0;
232  iflag = 0;
233  Real kappa = 0.0;
234  Real beta = 0.0;
235  Real sigma = 0.0;
236  Real alpha = 0.0;
237  Real tmp = 0.0;
238  Real gv = v->dot(*g);
239  Real sMp = 0.0;
240  pRed_ = 0.0;
241 
242  for (iter = 0; iter < maxit_; iter++) {
243  pObj.reducedHessVec(*Hp,*p,*xc,x,tol);
244 
245  kappa = p->dot(*Hp);
246  if (kappa <= 0) {
247  sigma = (-sMp+sqrt(sMp*sMp+pnorm2*(del*del-snorm2)))/pnorm2;
248  s.axpy(sigma,*p);
249  iflag = 2;
250  break;
251  }
252 
253  alpha = gv/kappa;
254  s1->set(s);
255  s1->axpy(alpha,*p);
256  s1norm2 = snorm2 + 2.0*alpha*sMp + alpha*alpha*pnorm2;
257 
258  if (s1norm2 >= del*del) {
259  sigma = (-sMp+sqrt(sMp*sMp+pnorm2*(del*del-snorm2)))/pnorm2;
260  s.axpy(sigma,*p);
261  iflag = 3;
262  break;
263  }
264 
265  pRed_ += 0.5*alpha*gv;
266 
267  s.set(*s1);
268  snorm2 = s1norm2;
269 
270  g->axpy(alpha,*Hp);
271  normg = g->norm();
272  if (normg < gtol) {
273  break;
274  }
275 
276  pObj.reducedPrecond(*v,*g,*xc,x,tol);
277  tmp = gv;
278  gv = v->dot(*g);
279  beta = gv/tmp;
280 
281  p->scale(beta);
282  p->axpy(-1.0,*v);
283  sMp = beta*(sMp+alpha*pnorm2);
284  pnorm2 = gv + beta*beta*pnorm2;
285  }
286  if (iflag > 0) {
287  pRed_ += sigma*(gv-0.5*sigma*kappa);
288  }
289 
290  if (iter == maxit_) {
291  iflag = 1;
292  }
293  if (iflag != 1) {
294  iter++;
295  }
296 
297  snorm = s.norm();
298  }
299 #endif
300 };
301 
302 }
303 
304 #endif
void initialize(const Vector< Real > &x, const Vector< Real > &s, const Vector< Real > &g)
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:143
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.
Teuchos::RCP< Vector< Real > > primalVector_
Teuchos::RCP< Vector< Real > > p_
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
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:157
Contains definitions for helper functions in ROL.
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
virtual void dualTransform(Vector< Real > &tv, const Vector< Real > &v)
Teuchos::RCP< Vector< Real > > v_
virtual const Teuchos::RCP< const Vector< Real > > getGradient(void) const
void setPredictedReduction(const Real pRed)
Teuchos::RCP< Vector< Real > > Hp_
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)
Teuchos::RCP< Vector< Real > > g_
Provides interface for truncated CG trust-region subproblem solver.
virtual void set(const Vector &x)
Set where .
Definition: ROL_Vector.hpp:196
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 ...
virtual void precond(Vector< Real > &Pv, const Vector< Real > &v, const Vector< Real > &s, Real &tol)
Apply preconditioner to vector.
Teuchos::RCP< Vector< Real > > s_
TruncatedCG(Teuchos::ParameterList &parlist)
virtual void primalTransform(Vector< Real > &tv, const Vector< Real > &v)