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.0) {
78  // Unravel Parameter List
79  maxit_ = parlist.sublist("General").sublist("Krylov").get("Iteration Limit",20);
80  tol1_ = parlist.sublist("General").sublist("Krylov").get("Absolute Tolerance",1.e-4);
81  tol2_ = parlist.sublist("General").sublist("Krylov").get("Relative Tolerance",1.e-2);
82  }
83 
84  void initialize( const Vector<Real> &x, const Vector<Real> &s, const Vector<Real> &g) {
86 
87  primalVector_ = x.clone();
88 
89  s_ = s.clone();
90  g_ = g.clone();
91  v_ = s.clone();
92  p_ = s.clone();
93  Hp_ = g.clone();
94  }
95 
96  void run( Vector<Real> &s, Real &snorm, Real &del, int &iflag, int &iter, const Vector<Real> &x,
97  const Vector<Real> &grad, const Real &gnorm, ProjectedObjective<Real> &pObj ) {
98  Real tol = std::sqrt(ROL_EPSILON);
99  const Real gtol = std::min(tol1_,tol2_*gnorm);
100 
101  // Gradient Vector
102  g_->set(grad);
103  Real normg = gnorm;
104  if ( pObj.isConActivated() ) {
105  primalVector_->set(grad.dual());
106  pObj.pruneActive(*primalVector_,grad.dual(),x);
107  g_->set(primalVector_->dual());
108  normg = g_->norm();
109  }
110 
111  // Old and New Step Vectors
112  s.zero(); s_->zero();
113  snorm = 0.0;
114  Real snorm2 = 0.0, s1norm2 = 0.0;
115 
116  // Preconditioned Gradient Vector
117  //pObj.precond(*v,*g,x,tol);
118  pObj.reducedPrecond(*v_,*g_,x,grad.dual(),x,tol);
119 
120  // Basis Vector
121  p_->set(*v_);
122  p_->scale(-1.0);
123  Real pnorm2 = v_->dot(g_->dual());
124 
125  iter = 0;
126  iflag = 0;
127  Real kappa = 0.0;
128  Real beta = 0.0;
129  Real sigma = 0.0;
130  Real alpha = 0.0;
131  Real tmp = 0.0;
132  Real gv = v_->dot(g_->dual());
133  Real sMp = 0.0;
134  pRed_ = 0.0;
135 
136  for (iter = 0; iter < maxit_; iter++) {
137  //pObj.hessVec(*Hp,*p,x,tol);
138  pObj.reducedHessVec(*Hp_,*p_,x,grad.dual(),x,tol);
139 
140  kappa = p_->dot(Hp_->dual());
141  if (kappa <= 0.0) {
142  sigma = (-sMp+sqrt(sMp*sMp+pnorm2*(del*del-snorm2)))/pnorm2;
143  s.axpy(sigma,*p_);
144  iflag = 2;
145  break;
146  }
147 
148  alpha = gv/kappa;
149  s_->set(s);
150  s_->axpy(alpha,*p_);
151  s1norm2 = snorm2 + 2.0*alpha*sMp + alpha*alpha*pnorm2;
152 
153  if (s1norm2 >= del*del) {
154  sigma = (-sMp+sqrt(sMp*sMp+pnorm2*(del*del-snorm2)))/pnorm2;
155  s.axpy(sigma,*p_);
156  iflag = 3;
157  break;
158  }
159 
160  pRed_ += 0.5*alpha*gv;
161 
162  s.set(*s_);
163  snorm2 = s1norm2;
164 
165  g_->axpy(alpha,*Hp_);
166  normg = g_->norm();
167  if (normg < gtol) {
168  break;
169  }
170 
171  //pObj.precond(*v,*g,x,tol);
172  pObj.reducedPrecond(*v_,*g_,x,grad.dual(),x,tol);
173  tmp = gv;
174  gv = v_->dot(g_->dual());
175  beta = gv/tmp;
176 
177  p_->scale(beta);
178  p_->axpy(-1.0,*v_);
179  sMp = beta*(sMp+alpha*pnorm2);
180  pnorm2 = gv + beta*beta*pnorm2;
181  }
182  if (iflag > 0) {
183  pRed_ += sigma*(gv-0.5*sigma*kappa);
184  }
185 
186  if (iter == maxit_) {
187  iflag = 1;
188  }
189  if (iflag != 1) {
190  iter++;
191  }
192 
193  snorm = s.norm();
195  }
196 
197 #if 0
198  void truncatedCG_proj( Vector<Real> &s, Real &snorm, Real &del, int &iflag, int &iter, const Vector<Real> &x,
199  const Vector<Real> &grad, const Real &gnorm, ProjectedObjective<Real> &pObj ) {
200  Real tol = std::sqrt(ROL_EPSILON);
201 
202  const Real gtol = std::min(tol1_,tol2_*gnorm);
203 
204  // Compute Cauchy Point
205  Real scnorm = 0.0;
206  Teuchos::RCP<Vector<Real> > sc = x.clone();
207  cauchypoint(*sc,scnorm,del,iflag,iter,x,grad,gnorm,pObj);
208  Teuchos::RCP<Vector<Real> > xc = x.clone();
209  xc->set(x);
210  xc->plus(*sc);
211 
212  // Old and New Step Vectors
213  s.set(*sc);
214  snorm = s.norm();
215  Real snorm2 = snorm*snorm;
216  Teuchos::RCP<Vector<Real> > s1 = x.clone();
217  s1->zero();
218  Real s1norm2 = 0.0;
219 
220  // Gradient Vector
221  Teuchos::RCP<Vector<Real> > g = x.clone();
222  g->set(grad);
223  Teuchos::RCP<Vector<Real> > Hs = x.clone();
224  pObj.reducedHessVec(*Hs,s,*xc,x,tol);
225  g->plus(*Hs);
226  Real normg = g->norm();
227 
228  // Preconditioned Gradient Vector
229  Teuchos::RCP<Vector<Real> > v = x.clone();
230  pObj.reducedPrecond(*v,*g,*xc,x,tol);
231 
232  // Basis Vector
233  Teuchos::RCP<Vector<Real> > p = x.clone();
234  p->set(*v);
235  p->scale(-1.0);
236  Real pnorm2 = v->dot(*g);
237 
238  // Hessian Times Basis Vector
239  Teuchos::RCP<Vector<Real> > Hp = x.clone();
240 
241  iter = 0;
242  iflag = 0;
243  Real kappa = 0.0;
244  Real beta = 0.0;
245  Real sigma = 0.0;
246  Real alpha = 0.0;
247  Real tmp = 0.0;
248  Real gv = v->dot(*g);
249  Real sMp = 0.0;
250  pRed_ = 0.0;
251 
252  for (iter = 0; iter < maxit_; iter++) {
253  pObj.reducedHessVec(*Hp,*p,*xc,x,tol);
254 
255  kappa = p->dot(*Hp);
256  if (kappa <= 0) {
257  sigma = (-sMp+sqrt(sMp*sMp+pnorm2*(del*del-snorm2)))/pnorm2;
258  s.axpy(sigma,*p);
259  iflag = 2;
260  break;
261  }
262 
263  alpha = gv/kappa;
264  s1->set(s);
265  s1->axpy(alpha,*p);
266  s1norm2 = snorm2 + 2.0*alpha*sMp + alpha*alpha*pnorm2;
267 
268  if (s1norm2 >= del*del) {
269  sigma = (-sMp+sqrt(sMp*sMp+pnorm2*(del*del-snorm2)))/pnorm2;
270  s.axpy(sigma,*p);
271  iflag = 3;
272  break;
273  }
274 
275  pRed_ += 0.5*alpha*gv;
276 
277  s.set(*s1);
278  snorm2 = s1norm2;
279 
280  g->axpy(alpha,*Hp);
281  normg = g->norm();
282  if (normg < gtol) {
283  break;
284  }
285 
286  pObj.reducedPrecond(*v,*g,*xc,x,tol);
287  tmp = gv;
288  gv = v->dot(*g);
289  beta = gv/tmp;
290 
291  p->scale(beta);
292  p->axpy(-1.0,*v);
293  sMp = beta*(sMp+alpha*pnorm2);
294  pnorm2 = gv + beta*beta*pnorm2;
295  }
296  if (iflag > 0) {
297  pRed_ += sigma*(gv-0.5*sigma*kappa);
298  }
299 
300  if (iter == maxit_) {
301  iflag = 1;
302  }
303  if (iflag != 1) {
304  iter++;
305  }
306 
307  snorm = s.norm();
308  }
309 #endif
310 };
311 
312 }
313 
314 #endif
void initialize(const Vector< Real > &x, const Vector< Real > &s, const Vector< Real > &g)
virtual const Vector & dual() const
Return dual representation of , for example, the result of applying a Riesz map, or change of basis...
Definition: ROL_Vector.hpp:213
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.
void pruneActive(Vector< Real > &v, const Vector< Real > &g, const Vector< Real > &x)
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
Teuchos::RCP< Vector< Real > > v_
void setPredictedReduction(const Real pRed)
Teuchos::RCP< Vector< Real > > Hp_
Teuchos::RCP< Vector< Real > > g_
Provides interface for truncated CG trust-region subproblem solver.
void run(Vector< Real > &s, Real &snorm, Real &del, int &iflag, int &iter, const Vector< Real > &x, const Vector< Real > &grad, const Real &gnorm, ProjectedObjective< Real > &pObj)
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 of...
Teuchos::RCP< Vector< Real > > s_
TruncatedCG(Teuchos::ParameterList &parlist)
static const double ROL_EPSILON
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:118