finite_diff_hessian.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_MAT_FUNCTOR_FINITE_DIFF_HESSIAN_HPP
2 #define STAN_MATH_PRIM_MAT_FUNCTOR_FINITE_DIFF_HESSIAN_HPP
3
6
7 namespace stan {
8 namespace math {
9
10 template <typename F>
12  const F& f, const Eigen::Matrix<double, Eigen::Dynamic, 1>& x, int lambda,
13  double epsilon = 1e-03) {
14  using Eigen::Dynamic;
15  using Eigen::Matrix;
16
17  Matrix<double, Dynamic, 1> x_temp(x);
18
20  x_temp(lambda) = x(lambda) + 2.0 * epsilon;
22
23  x_temp(lambda) = x(lambda) + -2.0 * epsilon;
25
26  x_temp(lambda) = x(lambda) + epsilon;
27  grad += 8.0 * f(x_temp);
28
29  x_temp(lambda) = x(lambda) + -epsilon;
30  grad -= 8.0 * f(x_temp);
31
33 }
34
35 /**
36  * Calculate the value and the Hessian of the specified function
37  * at the specified argument using second-order finite difference.
38  *
39  * <p>The functor must implement
40  *
41  * <code>
42  * double
43  * operator()(const
44  * Eigen::Matrix<double, Eigen::Dynamic, 1>&)
45  * </code>
46  *
47  * Error should be on order of epsilon ^ 4, with 4
48  * calls to the function f.
49  *
50  * Reference:
51  * Eberly: Derivative Approximation by Finite Differences
52  * Page 6
53  *
54  * @tparam F Type of function
55  * @param[in] f Function
56  * @param[in] x Argument to function
57  * @param[out] fx Function applied to argument
59  * @param[out] hess_fx Hessian of function at argument
60  * @param[in] epsilon perturbation size
61  */
62 template <typename F>
63 void finite_diff_hessian(const F& f, const Eigen::Matrix<double, -1, 1>& x,
64  double& fx, Eigen::Matrix<double, -1, 1>& grad_fx,
65  Eigen::Matrix<double, -1, -1>& hess_fx,
66  double epsilon = 1e-03) {
67  using Eigen::Dynamic;
68  using Eigen::Matrix;
69
70  int d = x.size();
71
72  Matrix<double, Dynamic, 1> x_temp(x);
73  hess_fx.resize(d, d);
74
76  double f_diff(0.0);
77  for (int i = 0; i < d; ++i) {
78  for (int j = i; j < d; ++j) {
79  x_temp(i) += 2.0 * epsilon;
80  if (i != j) {
81  f_diff = -finite_diff_hess_helper(f, x_temp, j);
82  x_temp(i) = x(i) + -2.0 * epsilon;
83  f_diff += finite_diff_hess_helper(f, x_temp, j);
84  x_temp(i) = x(i) + epsilon;
85  f_diff += 8.0 * finite_diff_hess_helper(f, x_temp, j);
86  x_temp(i) = x(i) + -epsilon;
87  f_diff -= 8.0 * finite_diff_hess_helper(f, x_temp, j);
88  f_diff /= 12.0 * epsilon * 12.0 * epsilon;
89  } else {
90  f_diff = -f(x_temp);
91  f_diff -= 30 * fx;
92  x_temp(i) = x(i) + -2.0 * epsilon;
93  f_diff -= f(x_temp);
94  x_temp(i) = x(i) + epsilon;
95  f_diff += 16.0 * f(x_temp);
96  x_temp(i) = x(i) - epsilon;
97  f_diff += 16.0 * f(x_temp);
98  f_diff /= 12 * epsilon * epsilon;
99  }
100
101  x_temp(i) = x(i);
102
103  hess_fx(j, i) = f_diff;
104  hess_fx(i, j) = hess_fx(j, i);
105  }
106  }
107 }
108 } // namespace math
109 } // namespace stan
110 #endif
#define F(x, y, z)