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