ROL
ROL_ObjectiveDef.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_OBJECTIVE_DEF_H
45 #define ROL_OBJECTIVE_DEF_H
46 
51 namespace ROL {
52 
53 template <class Real>
54 Real Objective<Real>::dirDeriv( const Vector<Real> &x, const Vector<Real> &d, Real &tol) {
55  Real ftol = std::sqrt(ROL_EPSILON);
56 
57  Teuchos::RCP<Vector<Real> > xd = d.clone();
58  xd->set(x);
59  xd->axpy(tol, d);
60  return (this->value(*xd,ftol) - this->value(x,ftol)) / tol;
61 }
62 
63 template <class Real>
64 void Objective<Real>::gradient( Vector<Real> &g, const Vector<Real> &x, Real &tol ) {
65  g.zero();
66  Real deriv = 0.0, h = 0.0, xi = 0.0;
67  for (int i = 0; i < g.dimension(); i++) {
68  xi = std::abs(x.dot(*x.basis(i)));
69  h = ((xi < ROL_EPSILON) ? 1. : xi)*tol;
70  deriv = this->dirDeriv(x,*x.basis(i),h);
71  g.axpy(deriv,*g.basis(i));
72  }
73 }
74 
75 template <class Real>
76 void Objective<Real>::hessVec( Vector<Real> &hv, const Vector<Real> &v, const Vector<Real> &x, Real &tol ) {
77  // Get Step Length
78  if ( v.norm() == 0. ) {
79  hv.zero();
80  }
81  else {
82  Real gtol = std::sqrt(ROL_EPSILON);
83 
84  Real h = std::max(1.0,x.norm()/v.norm())*tol;
85  //Real h = 2.0/(v.norm()*v.norm())*tol;
86 
87  // Compute Gradient at x
88  Teuchos::RCP<Vector<Real> > g = hv.clone();
89  this->gradient(*g,x,gtol);
90 
91  // Compute New Step x + h*v
92  Teuchos::RCP<Vector<Real> > xnew = x.clone();
93  xnew->set(x);
94  xnew->axpy(h,v);
95  this->update(*xnew);
96 
97  // Compute Gradient at x + h*v
98  hv.zero();
99  this->gradient(hv,*xnew,gtol);
100 
101  // Compute Newton Quotient
102  hv.axpy(-1.0,*g);
103  hv.scale(1.0/h);
104  }
105 }
106 
107 
108 
109 template <class Real>
110 std::vector<std::vector<Real> > Objective<Real>::checkGradient( const Vector<Real> &x,
111  const Vector<Real> &g,
112  const Vector<Real> &d,
113  const bool printToStream,
114  std::ostream & outStream,
115  const int numSteps,
116  const int order ) {
117 
118  std::vector<Real> steps(numSteps);
119  for(int i=0;i<numSteps;++i) {
120  steps[i] = pow(10,-i);
121  }
122 
123  return checkGradient(x,g,d,steps,printToStream,outStream,order);
124 
125 } // checkGradient
126 
127 
128 
129 template <class Real>
130 std::vector<std::vector<Real> > Objective<Real>::checkGradient( const Vector<Real> &x,
131  const Vector<Real> &g,
132  const Vector<Real> &d,
133  const std::vector<Real> &steps,
134  const bool printToStream,
135  std::ostream & outStream,
136  const int order ) {
137 
138  TEUCHOS_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
139  "Error: finite difference order must be 1,2,3, or 4" );
140 
143 
144  Real tol = std::sqrt(ROL_EPSILON);
145 
146  int numSteps = steps.size();
147  int numVals = 4;
148  std::vector<Real> tmp(numVals);
149  std::vector<std::vector<Real> > gCheck(numSteps, tmp);
150 
151  // Save the format state of the original outStream.
152  Teuchos::oblackholestream oldFormatState;
153  oldFormatState.copyfmt(outStream);
154 
155  // Evaluate objective value at x.
156  this->update(x);
157  Real val = this->value(x,tol);
158 
159  // Compute gradient at x.
160  Teuchos::RCP<Vector<Real> > gtmp = g.clone();
161  this->gradient(*gtmp, x, tol);
162  Real dtg = d.dot(gtmp->dual());
163 
164  // Temporary vectors.
165  Teuchos::RCP<Vector<Real> > xnew = x.clone();
166 
167  for (int i=0; i<numSteps; i++) {
168 
169  Real eta = steps[i];
170 
171  xnew->set(x);
172 
173  // Compute gradient, finite-difference gradient, and absolute error.
174  gCheck[i][0] = eta;
175  gCheck[i][1] = dtg;
176 
177  gCheck[i][2] = weights[order-1][0] * val;
178 
179  for(int j=0; j<order; ++j) {
180  // Evaluate at x <- x+eta*c_i*d.
181  xnew->axpy(eta*shifts[order-1][j], d);
182 
183  // Only evaluate at shifts where the weight is nonzero
184  if( weights[order-1][j+1] != 0 ) {
185  this->update(*xnew);
186  gCheck[i][2] += weights[order-1][j+1] * this->value(*xnew,tol);
187  }
188  }
189 
190  gCheck[i][2] /= eta;
191 
192  gCheck[i][3] = std::abs(gCheck[i][2] - gCheck[i][1]);
193 
194  if (printToStream) {
195  if (i==0) {
196  outStream << std::right
197  << std::setw(20) << "Step size"
198  << std::setw(20) << "grad'*dir"
199  << std::setw(20) << "FD approx"
200  << std::setw(20) << "abs error"
201  << "\n";
202  }
203  outStream << std::scientific << std::setprecision(11) << std::right
204  << std::setw(20) << gCheck[i][0]
205  << std::setw(20) << gCheck[i][1]
206  << std::setw(20) << gCheck[i][2]
207  << std::setw(20) << gCheck[i][3]
208  << "\n";
209  }
210 
211  }
212 
213  // Reset format state of outStream.
214  outStream.copyfmt(oldFormatState);
215 
216  return gCheck;
217 } // checkGradient
218 
219 
220 
221 
222 
223 
224 
225 
226 template <class Real>
227 std::vector<std::vector<Real> > Objective<Real>::checkHessVec( const Vector<Real> &x,
228  const Vector<Real> &hv,
229  const Vector<Real> &v,
230  const bool printToStream,
231  std::ostream & outStream,
232  const int numSteps,
233  const int order ) {
234  std::vector<Real> steps(numSteps);
235  for(int i=0;i<numSteps;++i) {
236  steps[i] = pow(10,-i);
237  }
238 
239  return checkHessVec(x,hv,v,steps,printToStream,outStream,order);
240 } // checkHessVec
241 
242 
243 
244 template <class Real>
245 std::vector<std::vector<Real> > Objective<Real>::checkHessVec( const Vector<Real> &x,
246  const Vector<Real> &hv,
247  const Vector<Real> &v,
248  const std::vector<Real> &steps,
249  const bool printToStream,
250  std::ostream & outStream,
251  const int order ) {
252 
253  TEUCHOS_TEST_FOR_EXCEPTION( order<1 || order>4, std::invalid_argument,
254  "Error: finite difference order must be 1,2,3, or 4" );
255 
258 
259 
260  Real tol = std::sqrt(ROL_EPSILON);
261 
262  int numSteps = steps.size();
263  int numVals = 4;
264  std::vector<Real> tmp(numVals);
265  std::vector<std::vector<Real> > hvCheck(numSteps, tmp);
266 
267  // Save the format state of the original outStream.
268  Teuchos::oblackholestream oldFormatState;
269  oldFormatState.copyfmt(outStream);
270 
271  // Compute gradient at x.
272  Teuchos::RCP<Vector<Real> > g = hv.clone();
273  this->update(x);
274  this->gradient(*g, x, tol);
275 
276  // Compute (Hessian at x) times (vector v).
277  Teuchos::RCP<Vector<Real> > Hv = hv.clone();
278  this->hessVec(*Hv, v, x, tol);
279  Real normHv = Hv->norm();
280 
281  // Temporary vectors.
282  Teuchos::RCP<Vector<Real> > gdif = hv.clone();
283  Teuchos::RCP<Vector<Real> > gnew = hv.clone();
284  Teuchos::RCP<Vector<Real> > xnew = x.clone();
285 
286  for (int i=0; i<numSteps; i++) {
287 
288  Real eta = steps[i];
289 
290  // Evaluate objective value at x+eta*d.
291  xnew->set(x);
292 
293  gdif->set(*g);
294  gdif->scale(weights[order-1][0]);
295 
296  for(int j=0; j<order; ++j) {
297 
298  // Evaluate at x <- x+eta*c_i*d.
299  xnew->axpy(eta*shifts[order-1][j], v);
300 
301  // Only evaluate at shifts where the weight is nonzero
302  if( weights[order-1][j+1] != 0 ) {
303  this->update(*xnew);
304  this->gradient(*gnew, *xnew, tol);
305  gdif->axpy(weights[order-1][j+1],*gnew);
306  }
307 
308  }
309 
310  gdif->scale(1.0/eta);
311 
312  // Compute norms of hessvec, finite-difference hessvec, and error.
313  hvCheck[i][0] = eta;
314  hvCheck[i][1] = normHv;
315  hvCheck[i][2] = gdif->norm();
316  gdif->axpy(-1.0, *Hv);
317  hvCheck[i][3] = gdif->norm();
318 
319  if (printToStream) {
320  if (i==0) {
321  outStream << std::right
322  << std::setw(20) << "Step size"
323  << std::setw(20) << "norm(Hess*vec)"
324  << std::setw(20) << "norm(FD approx)"
325  << std::setw(20) << "norm(abs error)"
326  << "\n";
327  }
328  outStream << std::scientific << std::setprecision(11) << std::right
329  << std::setw(20) << hvCheck[i][0]
330  << std::setw(20) << hvCheck[i][1]
331  << std::setw(20) << hvCheck[i][2]
332  << std::setw(20) << hvCheck[i][3]
333  << "\n";
334  }
335 
336  }
337 
338  // Reset format state of outStream.
339  outStream.copyfmt(oldFormatState);
340 
341  return hvCheck;
342 } // checkHessVec
343 
344 
345 
346 template<class Real>
347 std::vector<Real> Objective<Real>::checkHessSym( const Vector<Real> &x,
348  const Vector<Real> &hv,
349  const Vector<Real> &v,
350  const Vector<Real> &w,
351  const bool printToStream,
352  std::ostream & outStream ) {
353 
354  Real tol = std::sqrt(ROL_EPSILON);
355 
356  // Compute (Hessian at x) times (vector v).
357  Teuchos::RCP<Vector<Real> > h = hv.clone();
358  this->hessVec(*h, v, x, tol);
359  Real wHv = w.dot(h->dual());
360 
361  this->hessVec(*h, w, x, tol);
362  Real vHw = v.dot(h->dual());
363 
364  std::vector<Real> hsymCheck(3, 0);
365 
366  hsymCheck[0] = wHv;
367  hsymCheck[1] = vHw;
368  hsymCheck[2] = std::abs(vHw-wHv);
369 
370  // Save the format state of the original outStream.
371  Teuchos::oblackholestream oldFormatState;
372  oldFormatState.copyfmt(outStream);
373 
374  if (printToStream) {
375  outStream << std::right
376  << std::setw(20) << "<w, H(x)v>"
377  << std::setw(20) << "<v, H(x)w>"
378  << std::setw(20) << "abs error"
379  << "\n";
380  outStream << std::scientific << std::setprecision(11) << std::right
381  << std::setw(20) << hsymCheck[0]
382  << std::setw(20) << hsymCheck[1]
383  << std::setw(20) << hsymCheck[2]
384  << "\n";
385  }
386 
387  // Reset format state of outStream.
388  outStream.copyfmt(oldFormatState);
389 
390  return hsymCheck;
391 
392 } // checkHessSym
393 
394 
395 
396 } // namespace ROL
397 
398 #endif
virtual void scale(const Real alpha)=0
Compute where .
virtual int dimension() const
Return dimension of the vector space.
Definition: ROL_Vector.hpp:183
const double weights[4][5]
Definition: ROL_Types.hpp:1101
virtual void axpy(const Real alpha, const Vector &x)
Compute where .
Definition: ROL_Vector.hpp:143
virtual Real dirDeriv(const Vector< Real > &x, const Vector< Real > &d, Real &tol)
Compute directional derivative.
virtual void hessVec(Vector< Real > &hv, const Vector< Real > &v, const Vector< Real > &x, Real &tol)
Apply Hessian approximation to vector.
virtual Teuchos::RCP< Vector > clone() const =0
Clone to make a new (uninitialized) vector.
virtual void zero()
Set to zero vector.
Definition: ROL_Vector.hpp:157
Defines the linear algebra or vector space interface.
Definition: ROL_Vector.hpp:74
virtual Real dot(const Vector &x) const =0
Compute where .
virtual std::vector< std::vector< Real > > checkGradient(const Vector< Real > &x, const Vector< Real > &d, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Finite-difference gradient check.
virtual void gradient(Vector< Real > &g, const Vector< Real > &x, Real &tol)
Compute gradient.
virtual std::vector< std::vector< Real > > checkHessVec(const Vector< Real > &x, const Vector< Real > &v, const bool printToStream=true, std::ostream &outStream=std::cout, const int numSteps=ROL_NUM_CHECKDERIV_STEPS, const int order=1)
Finite-difference Hessian-applied-to-vector check.
virtual Teuchos::RCP< Vector > basis(const int i) const
Return i-th basis vector.
Definition: ROL_Vector.hpp:172
virtual Real norm() const =0
Returns where .
virtual std::vector< Real > checkHessSym(const Vector< Real > &x, const Vector< Real > &v, const Vector< Real > &w, const bool printToStream=true, std::ostream &outStream=std::cout)
Hessian symmetry check.
static const double ROL_EPSILON
Platform-dependent machine epsilon.
Definition: ROL_Types.hpp:118