RcppNumerical

Rcpp Integration for Numerical Computing Libraries <img src="https://statr.me/images/sticker-rcppnumerical.png" alt="RcppNumerical" height="150px" align="right" />

• Introduction
• Numerical Integration
  • One-dimensional
  • Multi-dimensional
  • Handling Infinite Limits
• Numerical Optimization
  • Unconstrained Minimization
  • Box-constrained Minimization
• A Practical Example

Introduction

Rcpp is a powerful tool to write fast C++ code to speed up R programs. However, it is not easy, or at least not straightforward, to compute numerical integration or do optimization using pure C++ code inside Rcpp.

RcppNumerical integrates a number of open source numerical computing libraries into Rcpp, so that users can call convenient functions to accomplish such tasks.

• To use RcppNumerical with Rcpp::sourceCpp(), add

// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]

in the C++ source file.

• To use RcppNumerical in your package, add Imports: RcppNumerical and LinkingTo: Rcpp, RcppEigen, RcppNumerical to the DESCRIPTION file, and import(RcppNumerical) to the NAMESPACE file.

Numerical Integration

One-dimensional

The one-dimensional numerical integration code contained in RcppNumerical is based on the NumericalIntegration library developed by Sreekumar Thaithara Balan, Mark Sauder, and Matt Beall.

To compute integration of a function, first define a functor derived from the Func class (under the namespace Numer):

class Func
{
public:
    virtual double operator()(const double& x) const = 0;
    virtual void eval(double* x, const int n) const
    {
        for(int i = 0; i < n; i++)
            x[i] = this->operator()(x[i]);
    }

    virtual ~Func() {}
};

The first function evaluates one point at a time, and the second version overwrites each point in the array by the corresponding function values. Only the second function will be used by the integration code, but usually it is easier to implement the first one.

RcppNumerical provides a wrapper function for the NumericalIntegration library with the following interface:

inline double integrate(
    const Func& f, const double& lower, const double& upper,
    double& err_est, int& err_code,
    const int subdiv = 100, const double& eps_abs = 1e-8, const double& eps_rel = 1e-6,
    const Integrator<double>::QuadratureRule rule = Integrator<double>::GaussKronrod41
)

• f: The functor of integrand.
• lower, upper: Limits of integral.
• err_est: Estimate of the error (output).
• err_code: Error code (output). See inst/include/integration/Integrator.h Line 676-704.
• subdiv: Maximum number of subintervals.
• eps_abs, eps_rel: Absolute and relative tolerance.
• rule: Integration rule. Possible values are GaussKronrod{15, 21, 31, 41, 51, 61, 71, 81, 91, 101, 121, 201}. Rules with larger values have better accuracy, but may involve more function calls.
• Return value: The final estimate of the integral.

See a full example below, which can be compiled using the Rcpp::sourceCpp function in Rcpp.

// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
#include <RcppNumerical.h>
using namespace Numer;

// P(0.3 < X < 0.8), X ~ Beta(a, b)
class BetaPDF: public Func
{
private:
    double a;
    double b;
public:
    BetaPDF(double a_, double b_) : a(a_), b(b_) {}

    double operator()(const double& x) const
    {
        return R::dbeta(x, a, b, 0);
    }
};

// [[Rcpp::export]]
Rcpp::List integrate_1d_test()
{
    const double a = 3, b = 10;
    const double lower = 0.3, upper = 0.8;
    const double true_val = R::pbeta(upper, a, b, 1, 0) -
                            R::pbeta(lower, a, b, 1, 0);

    BetaPDF f(a, b);
    double err_est;
    int err_code;
    const double res = integrate(f, lower, upper, err_est, err_code);
    return Rcpp::List::create(
        Rcpp::Named("true") = true_val,
        Rcpp::Named("approximate") = res,
        Rcpp::Named("error_estimate") = err_est,
        Rcpp::Named("error_code") = err_code
    );
}

Runing the integrate_1d_test() function in R gives

integrate_1d_test()
## $true
## [1] 0.2528108
##
## $approximate
## [1] 0.2528108
##
## $error_estimate
## [1] 2.806764e-15
##
## $error_code
## [1] 0

Multi-dimensional

Multi-dimensional integration in RcppNumerical is done by the Cuba library developed by Thomas Hahn.

To calculate the integration of a multivariate function, one needs to define a functor that inherits from the MFunc class:

class MFunc
{
public:
    virtual double operator()(Constvec& x) = 0;

    virtual ~MFunc() {}
};

Here Constvec represents a read-only vector with the definition

// Constant reference to a vector
using Constvec = const Eigen::Ref<const Eigen::VectorXd>;

(Basically you can treat Constvec as a const Eigen::VectorXd. Using Eigen::Ref is mainly to avoid memory copy. See the explanation here.)

The function provided by RcppNumerical for multi-dimensional integration is

inline double integrate(
    MFunc& f, Constvec& lower, Constvec& upper,
    double& err_est, int& err_code,
    const int maxeval = 1000,
    const double& eps_abs = 1e-6, const double& eps_rel = 1e-6
)

• f: The functor of integrand.
• lower, upper: Limits of integral. Both are vectors of the same dimension of f.
• err_est: Estimate of the error (output).
• err_code: Error code (output). Non-zero values indicate failure of convergence.
• maxeval: Maximum number of function evaluations.
• eps_abs, eps_rel: Absolute and relative tolerance.
• Return value: The final estimate of the integral.

See the example below:

// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
#include <RcppNumerical.h>
using namespace Numer;

// P(a1 < X1 < b1, a2 < X2 < b2), (X1, X2) ~ N([0], [1   rho])
//                                            ([0], [rho   1])
class BiNormal: public MFunc
{
private:
    const double rho;
    double const1;  // 2 * (1 - rho^2)
    double const2;  // 1 / (2 * PI) / sqrt(1 - rho^2)
public:
    BiNormal(const double& rho_) : rho(rho_)
    {
        const1 = 2.0 * (1.0 - rho * rho);
        const2 = 1.0 / (2 * M_PI) / std::sqrt(1.0 - rho * rho);
    }

    // PDF of bivariate normal
    double operator()(Constvec& x)
    {
        double z = x[0] * x[0] - 2 * rho * x[0] * x[1] + x[1] * x[1];
        return const2 * std::exp(-z / const1);
    }
};

// [[Rcpp::export]]
Rcpp::List integrate_md_test()
{
    BiNormal f(0.5);  // rho = 0.5
    Eigen::VectorXd lower(2);
    lower << -1, -1;
    Eigen::VectorXd upper(2);
    upper << 1, 1;
    double err_est;
    int err_code;
    const double res = integrate(f, lower, upper, err_est, err_code);
    return Rcpp::List::create(
        Rcpp::Named("approximate") = res,
        Rcpp::Named("error_estimate") = err_est,
        Rcpp::Named("error_code") = err_code
    );
}

We can test the result in R:

library(mvtnorm)
trueval = pmvnorm(c(-1, -1), c(1, 1), sigma = matrix(c(1, 0.5, 0.5, 1), 2))
integrate_md_test()
## $approximate
## [1] 0.4979718
##
## $error_estimate
## [1] 4.612333e-09
##
## $error_code
## [1] 0
as.numeric(trueval) - integrate_md_test()$approximate
## [1] 2.893336e-11

Handling Infinite Limits

Infinite intagral limits are also supported. In the case of one-dimensional integration:

// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
#include <RcppNumerical.h>
using namespace Numer;

class TestInf: public Func
{
public:
    double operator()(const double& x) const
    {
        return x * x * R::dnorm(x, 0.0, 1.0, 0);
    }
};

// [[Rcpp::export]]
Rcpp::List integrate_1d_inf_test(const double& lower, const double& upper)
{
    TestInf f;
    double err_est;
    int err_code;
    const double res = integrate(f, lower, upper, err_est, err_code);
    return Rcpp::List::create(
        Rcpp::Named("approximate") = res,
        Rcpp::Named("error_estimate") = err_est,
        Rcpp::Named("error_code") = err_code
    );
}

# integrate() in R
integrate(function(x) x^2 * dnorm(x), 0.5, Inf)
## 0.4845702 with absolute error < 3e-08
integrate_1d_inf_test(0.5, Inf)
## $approximate
## [1] 0.4845702
##
## $error_estimate
## [1] 1.633995e-08
##
## $error_code
## [1] 0

Similarly, for multi-dimensional integration, infinite limits are supported in each dimension by specifying std::numeric_limits<double>::infinity() or -std::numeric_limits<double>::infinity():

// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::depends(RcppNumerical)]]
#include <RcppNumerical.h>
using namespace Numer;

// Test 1: Semi-infinite [0, +Inf) x [0, 1]
// Integrate exp(-x) over x in [0, +Inf) and y in [0, 1]
// True value: 1.0
class SemiInfiniteTest: public MFunc
{
public:
    double operator()(Constvec& x)
    {
        return std::exp(-x[0]);
    }
};

// Test 2: Doubly-infinite (-Inf, +Inf) x [0, 1]
// Integrate exp(-x^2) over x in (-Inf, +Inf) and y in [0, 1]
// True value: sqrt(pi)
class DoublyInfiniteTest: public MFunc
{
public:
    double operator()(Constvec& x)
    {
        return std::exp(-x[0] * x[0]);
    }
};

// Test 3: All infinite bounds
// Integrate exp(-x^2 - y^2) over (-Inf, +Inf) x (-Inf, +Inf)
// Expected: pi
class Gaussian2D: public MFunc
{
public:
    doubl