# How to get the Hessian Matrix in C++ using StanHeaders Package

Hi. I’m writing a MCMC code in C++ using the interface Rcpp.

I’m interested to get the gradient and Hessian for the full conditional posterior density of a block of parameters. I’m trying to use StanHeaders Package like showed in Using the Stan Math C++ Library.

The gradient is working fine to the latent Gaussian Process “f” of a log-Gaussian Cox Process as shown below:

``````
// [[Rcpp::export]]
Eigen::VectorXd lccGFast(Eigen::VectorXd f, Eigen::VectorXd y, Eigen::MatrixXd Km1){
double fx;
int n = y.rows();
stan::math::var fKm1f = -(1.0/2.0)*f.transpose()*Km1*f;
stan::math::vector_v efyf(n);
for(int i=0; i< n; i++){
efyf[i] = 0.0;
efyf[i] = (exp(f[i])-y[i]*f[i]);
}
return(fKm1f-sum(efyf));},
}

``````

but when I try to use the similar “function” for the same example to get the Hessian:

``````
// [[Rcpp::export]]
Eigen::MatrixXd lccHFast(Eigen::VectorXd f, Eigen::VectorXd y, Eigen::MatrixXd Km1){
double fx;
Eigen::MatrixXd H;
stan::math::hessian([&y,&Km1](auto f){
int n = y.rows();
stan::math::var fKm1f = -(1.0/2.0)*f.transpose()*Km1*f;
stan::math::vector_v efyf(n);
for(int i=0; i< n; i++){
efyf[i] = 0.0;
efyf[i] = (exp(f[i])-y[i]*f[i]);
}
return(fKm1f-sum(efyf));
},
return H;
}

``````

many errors are reported (there are hundred of lines of error in R, i don’t know which part is important to help me to realize what is not working).

What am I missing?

2 Likes

Hi. It’s me again. Sorry but I’m new in the Stan Forum.

Should I give more details or classify my request in other Section?

Hi, sorry for not getting to you earlier. This is a good place for this type of questions. Can you share the full error log you are seeing (if it’s too long, you can upload it as .txt file)

error-hessian.txt (1.4 MB)
Thank you to reply Martin! I’ve attachted the error message in “error-hessian.txt” file.

I’ll also show the full “.cpp” file just for you to confirm if I’m missing to include any package, etc.

``````// [[Rcpp::depends(BH)]]
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppParallel)]]

#include <stan/math/rev/core.hpp>
#include <stan/math/mix/mat/functor/hessian.hpp> // then stuff from mix/ must come next
#include <stan/math.hpp>                         // finally pull in everything from rev/ & prim/
#include <Rcpp.h>
#include <RcppEigen.h>

// [[Rcpp::plugins(cpp14)]]

// [[Rcpp::export]]
Eigen::MatrixXd lccHFast(Eigen::VectorXd f, Eigen::VectorXd y, Eigen::MatrixXd Km1){
double fx;
Eigen::MatrixXd H;
stan::math::hessian([&y,&Km1](auto f){
int n = y.rows();
stan::math::var fKm1f = -(1.0/2.0)*f.transpose()*Km1*f;
stan::math::vector_v efyf(n);
for(int i=0; i< n; i++){
efyf[i] = 0.0;
efyf[i] = (exp(f[i])-y[i]*f[i]);
}
return(fKm1f-sum(efyf));
},
return H;
}

``````

Hi Marcel,

When reproducing this I actually had a lot of trouble with the current StanHeaders release, since it’s a bit behind. So first I would recommend updating to the experimental preview of the next RStan release, so that you have access to the latest StanHeaders:

``````remove.packages(c("StanHeaders", "rstan"))
install.packages("rstan", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
``````

Onto the hessian!

There were a couple things to be fixed here (the hessian stuff is tricky!). First, you need to include `stan/math/fwd.hpp`, otherwise there are some type traits that the compiler misses.

Additionally, rather than creating intermediary variables in the hessian functor, try to define it using Stan’s built-in functions. This way they can handle the different autodiff types implicitly.

``````Eigen::MatrixXd lccHFast(Eigen::VectorXd f, Eigen::VectorXd y, Eigen::MatrixXd Km1){
double fx;
Eigen::MatrixXd H;
stan::math::hessian([&y,&Km1](auto f){
int n = y.rows();
stan::math::var fKm1f = -(1.0/2.0)*f.transpose()*Km1*f;
stan::math::vector_v efyf(n);
for(int i=0; i< n; i++){
efyf[i] = 0.0;
efyf[i] = (exp(f[i])-y[i]*f[i]);
}
return(fKm1f-sum(efyf));
},
return H;
}
``````

This could be equivalently written as:

``````Eigen::MatrixXd lccHFast(Eigen::VectorXd f, Eigen::VectorXd y, Eigen::MatrixXd Km1){
double fx;
Eigen::MatrixXd H;
stan::math::hessian([&y,&Km1](auto f){
- stan::math::sum(stan::math::exp(f) - stan::math::elt_multiply(y,f));
},
return H;
}
``````

Putting that all together, this looks like:

``````// [[Rcpp::depends(BH)]]
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppParallel)]]

#include <stan/math/rev.hpp>
#include <stan/math/fwd.hpp>
#include <stan/math/mix.hpp> // then stuff from mix/ must come next
#include <stan/math.hpp>                         // finally pull in everything from rev/ & prim/
#include <Rcpp.h>
#include <RcppEigen.h>

// [[Rcpp::plugins(cpp14)]]

// [[Rcpp::export]]
Eigen::MatrixXd lccHFast(Eigen::VectorXd f, Eigen::VectorXd y, Eigen::MatrixXd Km1){
double fx;
Eigen::MatrixXd H;
stan::math::hessian([&y,&Km1](auto f){
- stan::math::sum(stan::math::exp(f) - stan::math::elt_multiply(y,f));
},
return H;
}
``````
2 Likes

That `error-hessian.txt` has a ton of warnings but I think the important part is

``````C:/Users/marcel.quintana/Documents/R/win-library/4.0/StanHeaders/include/stan/math/mix/functor/hessian.hpp:50:8: error: cannot convert 'stan::math::fvar<stan::math::var_value<double> >' to 'double' in assignment
fx = f(x);
~~~^~~~~~
``````

which points here.

The hessian functor wants a “generic lambda” as the argument; no single type is appropriate. I think you could write some template metaprogramming to infer the type.

``````  stan::math::hessian([&y,&Km1](auto f){
using local_scalar = stan::scalar_type_t<decltype(f)>;
int n = y.rows();
local_scalar fKm1f = -(1.0/2.0)*f.transpose()*Km1*f;
Eigen::Matrix<local_scalar, -1, 1> efyf(n);
for(int i=0; i< n; i++){
efyf[i] = 0.0;
efyf[i] = (exp(f[i])-y[i]*f[i]);
}
return(fKm1f-sum(efyf));
},
Niko, I would like to understand better the functionality of `stan::scalar_type_t`. Does it “extract” the type infered for `f` so as to define others elements with the same type?
`decltype(f)` is whatever the type of `f` happens to be, I think it’s either `Eigen::Matrix<double, -1, 1>` or `Eigen::Matrix<stan::math::fvar<stan::math::var>, -1, 1>` depending on where in `hessian()` it’s being called, and `stan::scalar_type_t<...>` extracts the `double` or `stan::math::fvar<stan::math::var>` out of it so that `local_scalar` is an alias for the correct autodiff scalar type.