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