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