I have been using Rcpp and RcppArmadillo in my R packages for a while now, but now I’d like to use some AD. I need some help/guidance as I try to understand how to use stanheaders.
I am able to reproduce the gradient example from the stanheaders vignette, and include it in a new and otherwise empty R package.
For those who are interested, this is my makevars file (works even with windows):
CXX14="$(BINPREF)g++ $(M_ARCH) -std=c++1y"
CXX14FLAGS="-O3 -Wall"
CXX14STD=-std=gnu1y
CXX11FLAGS=-O3-std=gnu++11
PKG_CXXFLAGS = -std=c++14
This is the linking to line in my DESCRIPTION file:
LinkingTo: Rcpp (>= 0.12.7), RcppEigen (>= 0.3.3.3.0), StanHeaders (>= 2.12.0), BH (>= 1.66.0)
I added a new .cpp file, and pasted the sample code from the vignette, calling the fuction gradExample:
#include <Rcpp.h>
#include <stan/math.hpp>
using namespace Rcpp;
// [[Rcpp::export]]
Eigen::VectorXd gradExample(Eigen::VectorXd x, Eigen::VectorXd a) { // gradient by AD using Stan
double fx;
Eigen::VectorXd grad_fx;
using stan::math::dot_self;
stan::math::gradient([&a](auto x) { return dot_self( (x - a).eval() ); },
x, fx, grad_fx);
return grad_fx;
}
It works fine, and just as expected.
As a next toy example, I wanted to move from scalars to vectors, so I decided to compute the partial derivative wrt x in Multivariate Normal distribution. I know, this is just a toy example to learn how stanheaders works. I call the function MultiNormalLpdfG.
// [[Rcpp::export]]
Eigen::VectorXd MultiNormalLpdfG(Eigen::VectorXd x, Eigen::VectorXd mu, Eigen::MatrixXd sig) {
double fx;
Eigen::VectorXd grad_fx;
using stan::math::multi_normal_lpdf;
stan::math::gradient([&mu,&sig](auto x) { return multi_normal_lpdf(x,mu,sig) ; },
x, fx, grad_fx);
return grad_fx;
}
Again, this works fine. It is a lot faster than using grad, and I also tested that results are consistent.
library(numDeriv)
library(microbenchmark)
library(mvtnorm)
microbenchmark(grad=grad(func = function(x){dmvnorm(x,c(0,0),diag(2)+.5,log=T)}, c(1:2))
, ad=MultiNormalLpdfG(1:2,c(0,0),diag(2)+.5)
)
So as long as I can put my target function into a lambda expression, using only stan math functions, I should be fine.
However, what if I need/want to create functions, and then call these functions? AD can be tricky, so I am curious about what’s possible here.
For my toy example, I create a function MultiNormalLpdftemp, which simply calls multi_normal_lpdf, but returns the result in a double.
// [[Rcpp::export]]
double MultiNormalLpdftemp(Eigen::VectorXd x, Eigen::VectorXd mu, Eigen::MatrixXd sig) {
double out;
using stan::math::multi_normal_lpdf;
out=multi_normal_lpdf(x,mu,sig);
return(out);
}
Now, what doesn’t work:
(1) lambda, calling custom function
MultiNormalLpdfG2 does not compile. I am calling MultiNormalLpdftemp instead of multi_normal_lpdf. I tried using to_var or not using to_var, but it does not compile.
Eigen::VectorXd MultiNormalLpdfG2(Eigen::VectorXd x, Eigen::VectorXd mu, Eigen::MatrixXd sig) {
double fx;
Eigen::VectorXd grad_fx;
using stan::math::to_var;
stan::math::gradient([&mu,&sig](auto x) { return (to_var(MultiNormalLpdftemp(x,mu,sig).eval())) ; },
x, fx, grad_fx);
return grad_fx;
}
(2) calling custom function
I also tried using .grad() and .adj(), however, Eigen::VectorXd does not support overloading / .adj() and stan::math::var does not seem to work with vectors.
// [[Rcpp::export]]
Eigen::VectorXd MultiNormalLpdfG3(Eigen::VectorXd x_, Eigen::VectorXd mu, Eigen::MatrixXd sig) {
using stan::math::to_vector;
// stan::math::var x = to_vector(x_);
Eigen::VectorXd x = to_vector(x_);
// function itself
stan::math::var lp = MultiNormalLpdftemp(x,mu,sig);
// Gradient
lp.grad();
//stan::math::recover_memory();
return x.adj();
}
I’d be really glad for hints/pointers/help. I intentionally kept the examples super minimal.