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

General
1

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;
  Eigen::VectorXd grad_fx;
  stan::math::gradient([&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));},
    f, fx, grad_fx);
  return grad_fx;
}

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::VectorXd grad_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));
  },
  f, fx, grad_fx, H);
  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?

Thanks in advance.

2 Likes
2

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?

Thank’s in advance.

3

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)

4

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)]]
// [[Rcpp::depends(StanHeaders)]]

#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::VectorXd grad_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));
  },
  f, fx, grad_fx, H);
  return H;
}
5

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("StanHeaders", repos = c("https://mc-stan.org/r-packages/", getOption("repos")))
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.

Given your function:

Eigen::MatrixXd lccHFast(Eigen::VectorXd f, Eigen::VectorXd y, Eigen::MatrixXd Km1){
  double fx;
  Eigen::VectorXd grad_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));
  },
  f, fx, grad_fx, H);
  return H;
}

This could be equivalently written as:

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

Putting that all together, this looks like:

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

#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::VectorXd grad_fx;
  Eigen::MatrixXd H;
  stan::math::hessian([&y,&Km1](auto f){
    return -(1.0/2.0)*stan::math::quad_form(Km1,f)
             - stan::math::sum(stan::math::exp(f) - stan::math::elt_multiply(y,f));
  },
  f, fx, grad_fx, H);
  return H;
}
2 Likes
6

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));
  },
  f, fx, grad_fx, H);
7

Thank you so much Andrew and Niko for your help, It really improved my knowledge about the concepts involved, and both solved the problem in different ways.

The codes I’ve showed were a simplified example for the code I’m really interested in but now you’ve given me more tools to accomplish the job.

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?

8

Yes.
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.

2 Likes