Problem computing the Hessian


#1

Hi guys,

Im new to stan::math, im trying to build a stat model which requires me to calculate the hessian of some product of gaussian cdf with a dot product of param and data as a function argument.

Im testing with calculating the hessian of just one gaussian cdf e.g.: {\bigtriangledown}^2 \Phi(\vec{ \theta } \cdot \vec{x}), where \vec{\theta} is the param and \vec{x} is the data vector.

Im getting an error as shown in the pic:

It says that it has no matching for call to dot product. However, the dot product with the above arugments is well defined according to the documentation:
fvar stan::math::dot_product ( const Eigen::Matrix< fvar< T >, R1, C1 > & v1,const Eigen::Matrix<double, R2, C2 > & v2 )

The code for producing this error is shown below:

#include <Rcpp.h>
#include <stan/math.hpp>
#include <RcppEigen.h>

using namespace Rcpp;
using namespace Eigen;
using namespace std;

typedef Eigen::Matrix <double, Dynamic,Dynamic> MatrixXd;
typedef Eigen::Matrix <double,Dynamic,1> VectorXd;
typedef Eigen::Matrix <double,1,Dynamic> RowVectorXd;

struct f {
const Eigen::Matrix<double, Eigen::Dynamic, 1> y_;
f(const Eigen::Matrix<double, Eigen::Dynamic, 1>& y) : y_(y) { }
template
T operator()(const Eigen::Matrix<T, Dynamic, 1>& theta) const {
T ret = 0;
ret += stan::math::Phi(stan::math::dot_product(theta,y_));

return ret;

}
};

// [[Rcpp::export]]
NumericMatrix getHess(stan::math::vector_d v_){

stan::math::vector_d v = v_;
int len = v_.size();

MatrixXd zero = MatrixXd::Constant(len,2,0);
stan::math::vector_d theta = zero;

stan::math::vector_d grad(len);
stan::math::matrix_d hess = MatrixXd::Constant(len,len,0);

double fx;

stan::math::hessian(f(v),theta,fx,grad,hess);

stan::math::recover_memory();
stan::math::ChainableStack::memalloc_.free_all();
stan::math::set_zero_all_adjoints();

return wrap(hess);
}

I have no problem with using the same functor f to calculate the gradient using stan::math::gradient() nor the numerical hessian using stan::math::finite_diff_hessian. What is the problem here?

Thanks in advance!


#2

For me, it at least compiles if you start with

#include <Rcpp.h>
// [[Rcpp::depends(StanHeaders)]]
#include <stan/math.hpp>
#include <stan/math/mix/mat.hpp>
// [[Rcpp::depends(RcppEigen)]]
#include <RcppEigen.h>

and add <typename T> after template in your struct. See attachment.
johne.txt (1.2 KB)