finite_diff_hessian_times_vector_auto.hpp

Go to the documentation of this file.

#ifndef STAN_MATH_REV_FUNCTOR_HESSIAN_TIMES_VECTOR_HPP
#define STAN_MATH_REV_FUNCTOR_HESSIAN_TIMES_VECTOR_HPP
 
#include <stan/math/prim/fun/constants.hpp>
#include <stan/math/rev/meta.hpp>
#include <stan/math/rev/core.hpp>
#include <stan/math/prim/fun/Eigen.hpp>
#include <stan/math/rev/functor.hpp>
#include <cmath>
 
namespace stan {
namespace math {
namespace internal {
 
template <typename F>
inline void finite_diff_hessian_times_vector_auto(const F& f,
                                                  const Eigen::VectorXd& x,
                                                  const Eigen::VectorXd& v,
                                                  double& fx,
                                                  Eigen::VectorXd& hvp) {
  fx = f(x);
 
  double epsilon = std::sqrt(EPSILON) * (1 + x.norm()) / v.norm();
 
  Eigen::VectorXd v_eps = epsilon * v;
 
  int d = x.size();
  double tmp;
  Eigen::VectorXd grad_forward(d);
  gradient(f, x + v_eps, tmp, grad_forward);
 
  Eigen::VectorXd grad_backward(d);
  gradient(f, x - v_eps, tmp, grad_backward);
 
  hvp.resize(d);
  hvp = (grad_forward - grad_backward) / (2 * epsilon);
}
}  // namespace internal
}  // namespace math
}  // namespace stan
#endif