AD with stan math, Rcpp

Algorithms
, ,
1

Hey Everyone,

I have been trying to use the stan math library and Rcpp to use autodiff in R. At this point, I’m just trying to figure out how to make it work. I haven’t used c++ before, so I’m sure that it is just a silly error on my part, but I’ve spent enough time trying to figure it out, that I figured I should ask. I am trying to get the gradients for:

\f(x,y, g) = exp{(x-y)^T * g^-1 * (x-y)} .

I understand this is easy to do symbolically, I’m just trying to figure out how to use autodiff.

The first code chunk below I am able to compile and the functions work correctly; when I don’t have the exponential part of the function:

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

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

// [[Rcpp::export]]
auto fun_test(Eigen::VectorXd x, Eigen::VectorXd y, Eigen::MatrixXd g) {
  using stan::math::transpose;
  using stan::math::inverse;
  return transpose( (x - y) ) * inverse( g ) * (x - y);
}

// [[Rcpp::export]]
auto grad_test(Eigen::VectorXd x, Eigen::VectorXd y, Eigen::MatrixXd g) {

  double fx;
  Eigen::VectorXd grad_fx;
  using stan::math::transpose;
  using stan::math::inverse;
  stan::math::gradient(
    [&](auto x) {
      return transpose( (x - y) ) * inverse(g) * (x - y);
    }, x, fx, grad_fx
  );
  return grad_fx;
}

But when I add the exponential component, it doesn’t compile (Rcpp::sourceCpp). Actually it works fine when adding the exponential to just the fun_test function, but it won’t compile when it is added to the grad_test function (the error if relevant is "Error 1 occurred building shared library.). Here is the code:

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

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

// [[Rcpp::export]]
auto fun_test(Eigen::VectorXd x, Eigen::VectorXd y, Eigen::MatrixXd g) {
  using stan::math::transpose;
  using stan::math::inverse;
  return exp( transpose( (x - y) ) * inverse( g ) * (x - y) );
}

// [[Rcpp::export]]
auto grad_test(Eigen::VectorXd x, Eigen::VectorXd y, Eigen::MatrixXd g) {

  double fx;
  Eigen::VectorXd grad_fx;
  using stan::math::transpose;
  using stan::math::inverse;
  stan::math::gradient(
    [&](auto x) {
      return exp( transpose( (x - y) ) * inverse(g) * (x - y) );
    }, x, fx, grad_fx
  );
  return grad_fx;
}

Thanks for any help you are willing to give!

Josh

1 Like
2

The error is happening because gradient() expects the function to return a single scalar value, but your function is returning a 1x1 matrix (I’m not entirely sure why the first didn’t error though).

To fix, just modify the functor in gradient() to return the first (only) element in the matrix:

  stan::math::gradient(
    [&](auto x) {
      auto res = exp(transpose( (x - y) ) * inverse(g) * (x - y));
      return res(0, 0);
    }, x, fx, grad_fx
  );
4 Likes
3

That did it, thank you very much!

1 Like