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