Implementation of Multidimensional numerical integration in C++ and R

Question

I'm trying to perform a 4-dimensional numerical integration in R using a function I wrote in C++ code which is then sourced in R using the Rcpp package.

Below there is my code:

// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
#define EIGEN_PERMANENTLY_DISABLE_STUPID_WARNINGS
#include <Eigen/Eigen>
#include <Rcpp.h>
#include <math.h>
#include <iostream>
#include <cstdlib>
#include <boost/math/special_functions/bessel.hpp>
#include <boost/math/special_functions/gamma.hpp>

using namespace std;
using namespace Rcpp;
using namespace Numer;

// [[Rcpp::depends(BH)]]

// [[Rcpp::export]]
double gaussian_free_diffusion(double x,double x0, double sigma, double t) {
  double pi = 2 * acos(0.0);
  double a1 = (1/sqrt(2.0 * pi * sigma * t));
  double b1 = exp(-pow((x - x0), 2.0)/(2.0 * sigma * t));
  double res = a1 * b1;
  return res;
}

// [[Rcpp::export]]
double integral_CppFunctionCUBATURE(NumericVector Zt_pos, const double& xA0, const double& xB0, const double& yA0,
   const double& yB0, const double& t1, const double& sigma){

  double xAt_pos = Zt_pos[0];
  double xBt_pos = Zt_pos[1];
  double yAt_pos = Zt_pos[2];
  double yBt_pos = Zt_pos[3];

  double temp_pbxA = gaussian_free_diffusion(xAt_pos, xA0, sigma,t1);
  double temp_pbxB = gaussian_free_diffusion(xBt_pos, xB0, sigma, t1);
  double temp_pbyA = gaussian_free_diffusion(yAt_pos, yA0, sigma,t1);
  double temp_pbyB = gaussian_free_diffusion(yBt_pos, yB0, sigma, t1);

  return (temp_pbxB * temp_pbyB) * (temp_pbxA * temp_pbyA);

};

integral_CppFunctionCUBATURE is the function I use with cubintegrate. In R, I would then do:

sourceCpp('./myCppcode.cpp')

xA0<-3
xB0<-2
yA0<-5
yB0<-10
t<-500
sigma<-1


cubature::cubintegrate(integral_CppFunctionCUBATURE,lower=rep(-1000,4), upper=rep(1000,4),
                       xA0=xA0, xB0=xB0, yA0=yA0,yB0=yB0, t1=t, sigma=sigma,method='cuhre')$integral

Result: 0.9999978

The multidimensional integration takes about 4 sec, but I would like to speed it up as much as possible, since when I run the multidimensional integration on several xA0 values this can take quite some time. Do you have any suggestion on how I could improve the speed of the code?

And also, would you suggest an alternative way to perform fast 4-dimensional integration in C++/Rcpp?

Answer 1

Sources of slowness

The actual calculations you are doing are fine, I don't see anything you can optimize that the compiler won't already do for you itself. The main overhead will be caused by R repeatedly having to call integral_CppFunctionCUBATURE(). The first issue I see is that you are passing NumericVector Zt_pos by value, but all the doubles by reference. That's exactly the wrong way around. While Rccp::NumericVector itself already has reference semantics, it still needs to do reference counting, and that can be avoided by passing it as a reference. You are not trying to modify the double parameters, so passing it by reference just means you are going to pay to overhead of pointer indirection. Pass them by value instead.

The next issue is that R's cubature::cubintegrate() function is going to repeatedly have to call your function. You can avoid this by implementing the same algorithm as R uses in C++, so it can inline your integral_CppFunctionCubature(). It probably won't give a dramatic speedup though; the Cuhre method in R is implemented in C, and function call overhead is probably not that big compared to the actual cost of the calculations you are doing.

Choose the right integration algorithm

I am assuming that you don't just want to integrate the Guassian distribution function from effectively \$-\infty\$ to \$+\infty\$, as the result will by definition be 1. Assuming you want to integrate other functions, first check if you can integrate that function analytically. If not, then which algorithm works best to integrate it will depend on the shape of the function inside the range you want to integrate over. The Cuhre algorithm might not be the best, so try out others and benchmark their speed and accuracy if possible.

Don't #include headers you are not using

I see you are including a lot of headers, but almost none are actually used. This will just unnecessarily slow down compilation. The only headers you need are those for Rccp (for Rccp::NumericVector) and those for the math functions you use. For the latter, make sure you include the C++ version:

#include <cmath>