grad_hessian.hpp

Go to the documentation of this file.

#ifndef STAN_MATH_MIX_FUNCTOR_GRAD_HESSIAN_HPP
#define STAN_MATH_MIX_FUNCTOR_GRAD_HESSIAN_HPP
 
#include <stan/math/rev/core.hpp>
#include <stan/math/fwd/core.hpp>
#include <stan/math/prim/fun/Eigen.hpp>
#include <stdexcept>
#include <vector>
 
namespace stan {
namespace math {
 
template <typename F>
void grad_hessian(
    const F& f, const Eigen::Matrix<double, Eigen::Dynamic, 1>& x, double& fx,
    Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic>& H,
    std::vector<Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> >&
        grad_H) {
  using Eigen::Dynamic;
  using Eigen::Matrix;
  fx = f(x);
  int d = x.size();
  H.resize(d, d);
  grad_H.resize(d, Matrix<double, Dynamic, Dynamic>(d, d));
  for (int i = 0; i < d; ++i) {
    for (int j = i; j < d; ++j) {
      // Run nested autodiff in this scope
      nested_rev_autodiff nested;
 
      Matrix<fvar<fvar<var> >, Dynamic, 1> x_ffvar(d);
      for (int k = 0; k < d; ++k) {
        x_ffvar(k)
            = fvar<fvar<var> >(fvar<var>(x(k), i == k), fvar<var>(j == k, 0));
      }
      fvar<fvar<var> > fx_ffvar = f(x_ffvar);
      H(i, j) = fx_ffvar.d_.d_.val();
      H(j, i) = H(i, j);
      grad(fx_ffvar.d_.d_.vi_);
      for (int k = 0; k < d; ++k) {
        grad_H[i](j, k) = x_ffvar(k).val_.val_.adj();
        grad_H[j](i, k) = grad_H[i](j, k);
      }
    }
  }
}
 
}  // namespace math
}  // namespace stan
#endif