finite_diff_grad_hessian.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_MIX_MAT_FUNCTOR_FINITE_DIFF_GRAD_HESSIAN_HPP
2 #define STAN_MATH_MIX_MAT_FUNCTOR_FINITE_DIFF_GRAD_HESSIAN_HPP
3 
5 #include <stan/math/rev/core.hpp>
7 #include <vector>
8 
9 namespace stan {
10 namespace math {
11 
12 /**
13  * Calculate the value and the gradient of the hessian of the specified
14  * function at the specified argument using second-order autodiff and
15  * first-order finite difference.
16  *
17  * <p>The functor must implement
18  *
19  * <code>
20  * double
21  * operator()(const
22  * Eigen::Matrix<double, Eigen::Dynamic, 1>&)
23  * </code>
24  *
25  * Reference:
26  *
27  * De Levie: An improved numerical approximation
28  * for the first derivative, page 3
29  *
30  * 4 calls to the function, f.
31  *
32  * @tparam F Type of function
33  * @param[in] f Function
34  * @param[in] x Argument to function
35  * @param[out] fx Function applied to argument
36  * @param[out] hess Hessian matrix
37  * @param[out] grad_hess_fx gradient of Hessian of function at argument
38  * @param[in] epsilon perturbation size
39  */
40 template <typename F>
42  const F& f, const Eigen::Matrix<double, -1, 1>& x, double& fx,
43  Eigen::Matrix<double, -1, -1>& hess,
44  std::vector<Eigen::Matrix<double, -1, -1> >& grad_hess_fx,
45  double epsilon = 1e-04) {
46  using Eigen::Dynamic;
47  using Eigen::Matrix;
48 
49  int d = x.size();
50  double dummy_fx_eval;
51 
52  Matrix<double, Dynamic, 1> x_temp(x);
53  Matrix<double, Dynamic, 1> grad_auto(d);
54  Matrix<double, Dynamic, Dynamic> hess_auto(d, d);
55  Matrix<double, Dynamic, Dynamic> hess_diff(d, d);
56 
57  hessian(f, x, fx, grad_auto, hess);
58  for (int i = 0; i < d; ++i) {
59  hess_diff.setZero();
60 
61  x_temp(i) = x(i) + 2.0 * epsilon;
62  hessian(f, x_temp, dummy_fx_eval, grad_auto, hess_auto);
63  hess_diff = -hess_auto;
64 
65  x_temp(i) = x(i) + -2.0 * epsilon;
66  hessian(f, x_temp, dummy_fx_eval, grad_auto, hess_auto);
67  hess_diff += hess_auto;
68 
69  x_temp(i) = x(i) + epsilon;
70  hessian(f, x_temp, dummy_fx_eval, grad_auto, hess_auto);
71  hess_diff += 8.0 * hess_auto;
72 
73  x_temp(i) = x(i) + -epsilon;
74  hessian(f, x_temp, dummy_fx_eval, grad_auto, hess_auto);
75  hess_diff -= 8.0 * hess_auto;
76 
77  x_temp(i) = x(i);
78  hess_diff /= 12.0 * epsilon;
79 
80  grad_hess_fx.push_back(hess_diff);
81  }
82  fx = f(x);
83 }
84 
85 } // namespace math
86 } // namespace stan
87 #endif
#define F(x, y, z)
void finite_diff_grad_hessian(const F &f, const Eigen::Matrix< double,-1, 1 > &x, double &fx, Eigen::Matrix< double,-1,-1 > &hess, std::vector< Eigen::Matrix< double,-1,-1 > > &grad_hess_fx, double epsilon=1e-04)
void hessian(const F &f, const Eigen::Matrix< T, Eigen::Dynamic, 1 > &x, T &fx, Eigen::Matrix< T, Eigen::Dynamic, 1 > &grad, Eigen::Matrix< T, Eigen::Dynamic, Eigen::Dynamic > &H)
Definition: hessian.hpp:42
Float_t d
Definition: plot.C:236
double e()
Definition: constants.hpp:89
double epsilon