finite_diff_gradient.hpp
Go to the documentation of this file.
1 #ifndef STAN_MATH_PRIM_MAT_FUNCTOR_FINITE_DIFF_GRADIENT_HPP
2 #define STAN_MATH_PRIM_MAT_FUNCTOR_FINITE_DIFF_GRADIENT_HPP
3 
5 
6 namespace stan {
7 namespace math {
8 
9 /**
10  * Calculate the value and the gradient of the specified function
11  * at the specified argument using finite difference.
12  *
13  * <p>The functor must implement
14  *
15  * <code>
16  * double
17  * operator()(const
18  * Eigen::Matrix<double, Eigen::Dynamic, 1>&)
19  * </code>
20  *
21  * Error should be on order of epsilon ^ 6.
22  * The reference for this algorithm is:
23  *
24  * De Levie: An improved numerical approximation
25  * for the first derivative, page 3
26  *
27  * This function involves 6 calls to f.
28  *
29  * @tparam F Type of function
30  * @param[in] f Function
31  * @param[in] x Argument to function
32  * @param[out] fx Function applied to argument
33  * @param[out] grad_fx Gradient of function at argument
34  * @param[in] epsilon perturbation size
35  */
36 template <typename F>
37 void finite_diff_gradient(const F& f, const Eigen::Matrix<double, -1, 1>& x,
38  double& fx, Eigen::Matrix<double, -1, 1>& grad_fx,
39  double epsilon = 1e-03) {
40  using Eigen::Dynamic;
41  using Eigen::Matrix;
42  Matrix<double, Dynamic, 1> x_temp(x);
43 
44  int d = x.size();
45  grad_fx.resize(d);
46 
47  fx = f(x);
48 
49  for (int i = 0; i < d; ++i) {
50  double delta_f = 0.0;
51 
52  x_temp(i) = x(i) + 3.0 * epsilon;
53  delta_f = f(x_temp);
54 
55  x_temp(i) = x(i) + 2.0 * epsilon;
56  delta_f -= 9.0 * f(x_temp);
57 
58  x_temp(i) = x(i) + epsilon;
59  delta_f += 45.0 * f(x_temp);
60 
61  x_temp(i) = x(i) + -3.0 * epsilon;
62  delta_f -= f(x_temp);
63 
64  x_temp(i) = x(i) + -2.0 * epsilon;
65  delta_f += 9.0 * f(x_temp);
66 
67  x_temp(i) = x(i) + -epsilon;
68  delta_f -= 45.0 * f(x_temp);
69 
70  delta_f /= 60 * epsilon;
71 
72  x_temp(i) = x(i);
73  grad_fx(i) = delta_f;
74  }
75 }
76 } // namespace math
77 } // namespace stan
78 #endif
#define F(x, y, z)
Float_t d
Definition: plot.C:236
void finite_diff_gradient(const F &f, const Eigen::Matrix< double,-1, 1 > &x, double &fx, Eigen::Matrix< double,-1, 1 > &grad_fx, double epsilon=1e-03)
double e()
Definition: constants.hpp:89
double epsilon