## Introduction

This R cubature package exposes both the hcubature and pcubature routines of the underlying C cubature library, including the vectorized interfaces.

Per the documentation, use of pcubature is advisable only for smooth integrands in dimensions up to three at most. In fact, the pcubature routines perform significantly worse than the vectorized hcubature in inappropriate cases. So when in doubt, you are better off using hcubature.

Version 2.0 of this package integrates the Cuba library as well, once again providing vectorized interfaces.

The main point of this note is to examine the difference vectorization makes. My recommendations are below in the summary section.

## A Timing Harness

Our harness will provide timing results for hcubature, pcubature (where appropriate) and Cuba cuhre calls. We begin by creating a harness for these calls.

library(benchr)
library(cubature)

harness <- function(which = NULL,
f, fv, lowerLimit, upperLimit, tol = 1e-3, times = 20, ...) {

fns <- c(hc = "Non-vectorized Hcubature",
hc.v = "Vectorized Hcubature",
pc = "Non-vectorized Pcubature",
pc.v = "Vectorized Pcubature",
cc = "Non-vectorized cubature::cuhre",
cc_v = "Vectorized cubature::cuhre")
cc <- function() cubature::cuhre(f = f,
lowerLimit = lowerLimit, upperLimit = upperLimit,
relTol = tol,
...)
cc_v <- function() cubature::cuhre(f = fv,
lowerLimit = lowerLimit, upperLimit = upperLimit,
relTol = tol,
nVec = 1024L,
...)

hc <- function() cubature::hcubature(f = f,
lowerLimit = lowerLimit,
upperLimit = upperLimit,
tol = tol,
...)

hc.v <- function() cubature::hcubature(f = fv,
lowerLimit = lowerLimit,
upperLimit = upperLimit,
tol = tol,
vectorInterface = TRUE,
...)

pc <- function() cubature::pcubature(f = f,
lowerLimit = lowerLimit,
upperLimit = upperLimit,
tol = tol,
...)

pc.v <- function() cubature::pcubature(f = fv,
lowerLimit = lowerLimit,
upperLimit = upperLimit,
tol = tol,
vectorInterface = TRUE,
...)

ndim = length(lowerLimit)

if (is.null(which)) {
fnIndices <- seq_along(fns)
} else {
fnIndices <- na.omit(match(which, names(fns)))
}
fnList <- lapply(names(fns)[fnIndices], function(x) call(x))

argList <- c(fnList, times = times, progress = FALSE)
result <- do.call(benchr::benchmark, args = argList)
d <- summary(result)[seq_along(fnIndices), ]
d$expr <- fns[fnIndices] d } We reel off the timing runs. ## Example 1. func <- function(x) sin(x) * cos(x) * exp(x) func.v <- function(x) { matrix(apply(x, 2, function(z) sin(z) * cos(z) * exp(z)), ncol = ncol(x)) } d <- harness(f = func, fv = func.v, lowerLimit = rep(0, 3), upperLimit = rep(1, 3), tol = 1e-5, times = 100) knitr::kable(d, digits = 3, row.names = FALSE) expr n.eval min lw.qu median mean up.qu max total relative Non-vectorized Hcubature 100 0.003 0.003 0.003 0.003 0.004 0.008 0.344 6.40 Vectorized Hcubature 100 0.000 0.000 0.000 0.001 0.001 0.003 0.055 1.00 Non-vectorized Pcubature 100 0.009 0.010 0.011 0.011 0.012 0.014 1.086 21.70 Vectorized Pcubature 100 0.001 0.001 0.001 0.002 0.001 0.004 0.151 2.79 Non-vectorized cubature::cuhre 100 0.005 0.005 0.005 0.006 0.006 0.008 0.579 11.00 Vectorized cubature::cuhre 100 0.001 0.001 0.001 0.001 0.001 0.003 0.073 1.39 ## Multivariate Normal Using cubature, we evaluate $\int_R\phi(x)dx$ where $$\phi(x)$$ is the three-dimensional multivariate normal density with mean 0, and variance $\Sigma = \left(\begin{array}{rrr} 1 &\frac{3}{5} &\frac{1}{3}\\ \frac{3}{5} &1 &\frac{11}{15}\\ \frac{1}{3} &\frac{11}{15} & 1 \end{array} \right)$ and $$R$$ is $$[-\frac{1}{2}, 1] \times [-\frac{1}{2}, 4] \times [-\frac{1}{2}, 2].$$ We construct a scalar function (my_dmvnorm) and a vector analog (my_dmvnorm_v). First the functions. m <- 3 sigma <- diag(3) sigma[2,1] <- sigma[1, 2] <- 3/5 ; sigma[3,1] <- sigma[1, 3] <- 1/3 sigma[3,2] <- sigma[2, 3] <- 11/15 logdet <- sum(log(eigen(sigma, symmetric = TRUE, only.values = TRUE)$values))
my_dmvnorm <- function (x, mean, sigma, logdet) {
x <- matrix(x, ncol = length(x))
distval <- stats::mahalanobis(x, center = mean, cov = sigma)
exp(-(3 * log(2 * pi) + logdet + distval)/2)
}

my_dmvnorm_v <- function (x, mean, sigma, logdet) {
distval <- stats::mahalanobis(t(x), center = mean, cov = sigma)
exp(matrix(-(3 * log(2 * pi) + logdet + distval)/2, ncol = ncol(x)))
}

Now the timing.

d <- harness(f = my_dmvnorm, fv = my_dmvnorm_v,
lowerLimit = rep(-0.5, 3),
upperLimit = c(1, 4, 2),
tol = 1e-5,
times = 10,
mean = rep(0, m), sigma = sigma, logdet = logdet)
knitr::kable(d, digits = 3)
expr n.eval min lw.qu median mean up.qu max total relative
Non-vectorized Hcubature 10 1.203 1.232 1.249 1.252 1.268 1.317 12.519 656.00
Vectorized Hcubature 10 0.002 0.002 0.003 0.003 0.003 0.003 0.026 1.32
Non-vectorized Pcubature 10 0.516 0.522 0.537 0.541 0.555 0.582 5.411 282.00
Vectorized Pcubature 10 0.002 0.002 0.002 0.002 0.002 0.002 0.019 1.00
Non-vectorized cubature::cuhre 10 0.499 0.509 0.514 0.515 0.518 0.540 5.147 270.00
Vectorized cubature::cuhre 10 0.004 0.005 0.005 0.005 0.006 0.007 0.053 2.67

The effect of vectorization is huge. So it makes sense for users to vectorize the integrands as much as possible for efficiency.

Furthermore, for this particular example, we expect mvtnorm::pmvnorm to do pretty well since it is specialized for the multivariate normal. The good news is that the vectorized versions of hcubature and pcubature are quite competitive if you compare the table above to the one below.

library(mvtnorm)
g1 <- function() mvtnorm::pmvnorm(lower = rep(-0.5, m),
upper = c(1, 4, 2), mean = rep(0, m), corr = sigma,
alg = Miwa(), abseps = 1e-5, releps = 1e-5)
g2 <- function() mvtnorm::pmvnorm(lower = rep(-0.5, m),
upper = c(1, 4, 2), mean = rep(0, m), corr = sigma,
alg = GenzBretz(), abseps = 1e-5, releps = 1e-5)
g3 <- function() mvtnorm::pmvnorm(lower = rep(-0.5, m),
upper = c(1, 4, 2), mean = rep(0, m), corr = sigma,
alg = TVPACK(), abseps = 1e-5, releps = 1e-5)

knitr::kable(summary(benchr::benchmark(g1(), g2(), g3(), times = 20, progress = FALSE)),
digits = 3, row.names = FALSE)
expr n.eval min lw.qu median mean up.qu max total relative
g1() 20 0.001 0.001 0.001 0.001 0.002 0.004 0.028 1.03
g2() 20 0.001 0.001 0.001 0.001 0.001 0.003 0.029 1.00
g3() 20 0.001 0.001 0.001 0.001 0.001 0.003 0.027 1.02

## Product of cosines

testFn0 <- function(x) prod(cos(x))
testFn0_v <- function(x) matrix(apply(x, 2, function(z) prod(cos(z))), ncol = ncol(x))

d <- harness(f = testFn0, fv = testFn0_v,
lowerLimit = rep(0, 2), upperLimit = rep(1, 2), times = 1000)
knitr::kable(d, digits = 3)
expr n.eval min lw.qu median mean up.qu max total relative
Non-vectorized Hcubature 1000 0.000 0.000 0.000 0.000 0.000 0.007 0.299 2.62
Vectorized Hcubature 1000 0.000 0.000 0.000 0.000 0.000 0.002 0.108 1.00
Non-vectorized Pcubature 1000 0.000 0.000 0.000 0.000 0.000 0.003 0.410 3.76
Vectorized Pcubature 1000 0.000 0.000 0.000 0.000 0.000 0.002 0.209 1.94
Non-vectorized cubature::cuhre 1000 0.002 0.003 0.003 0.003 0.003 0.042 3.065 28.00
Vectorized cubature::cuhre 1000 0.000 0.000 0.000 0.001 0.000 0.003 0.506 4.62

## Gaussian function

testFn1 <- function(x) {
val <- sum(((1 - x) / x)^2)
scale <- prod((2 / sqrt(pi)) / x^2)
exp(-val) * scale
}

testFn1_v <- function(x) {
val <- matrix(apply(x, 2, function(z) sum(((1 - z) / z)^2)), ncol(x))
scale <- matrix(apply(x, 2, function(z) prod((2 / sqrt(pi)) / z^2)), ncol(x))
exp(-val) * scale
}

d <- harness(f = testFn1, fv = testFn1_v,
lowerLimit = rep(0, 3), upperLimit = rep(1, 3), times = 10)

knitr::kable(d, digits = 3)
expr n.eval min lw.qu median mean up.qu max total relative
Non-vectorized Hcubature 10 0.020 0.021 0.022 0.022 0.023 0.024 0.221 98.6
Vectorized Hcubature 10 0.005 0.005 0.005 0.005 0.005 0.006 0.052 22.3
Non-vectorized Pcubature 10 0.000 0.000 0.000 0.000 0.000 0.001 0.005 2.0
Vectorized Pcubature 10 0.000 0.000 0.000 0.000 0.000 0.000 0.003 1.0
Non-vectorized cubature::cuhre 10 0.086 0.088 0.090 0.092 0.093 0.107 0.917 397.0
Vectorized cubature::cuhre 10 0.020 0.021 0.022 0.022 0.022 0.025 0.218 95.2

## Discontinuous function

testFn2 <- function(x) {
}

testFn2_v <- function(x) {
matrix(apply(x, 2, function(z) ifelse(sum(z * z) < radius * radius, 1, 0)), ncol = ncol(x))
}

d <- harness(which = c("hc", "hc.v", "cc", "cc_v"),
f = testFn2, fv = testFn2_v,
lowerLimit = rep(0, 2), upperLimit = rep(1, 2), times = 10)
knitr::kable(d, digits = 3)
expr n.eval min lw.qu median mean up.qu max total relative
Non-vectorized Hcubature 10 0.271 0.277 0.291 0.296 0.295 0.346 2.956 5.28
Vectorized Hcubature 10 0.051 0.054 0.055 0.056 0.057 0.065 0.564 1.00
Non-vectorized cubature::cuhre 10 4.662 4.686 4.719 4.724 4.752 4.805 47.241 85.60
Vectorized cubature::cuhre 10 1.058 1.080 1.099 1.104 1.104 1.232 11.039 19.90

## A Simple Polynomial (product of coordinates)

testFn3 <- function(x) prod(2 * x)
testFn3_v <- function(x) matrix(apply(x, 2, function(z) prod(2 * z)), ncol = ncol(x))

d <- harness(f = testFn3, fv = testFn3_v,
lowerLimit = rep(0, 3), upperLimit = rep(1, 3), times = 50)
knitr::kable(d, digits = 3)
expr n.eval min lw.qu median mean up.qu max total relative
Non-vectorized Hcubature 50 0.000 0.000 0.001 0.001 0.001 0.003 0.031 3.71
Vectorized Hcubature 50 0.000 0.000 0.000 0.000 0.000 0.000 0.008 1.03
Non-vectorized Pcubature 50 0.000 0.000 0.000 0.000 0.000 0.003 0.025 3.13
Vectorized Pcubature 50 0.000 0.000 0.000 0.000 0.000 0.000 0.007 1.00
Non-vectorized cubature::cuhre 50 0.005 0.005 0.006 0.006 0.007 0.009 0.303 42.60
Vectorized cubature::cuhre 50 0.001 0.001 0.001 0.001 0.001 0.003 0.045 5.83

## Gaussian centered at $$\frac{1}{2}$$

testFn4 <- function(x) {
a <- 0.1
s <- sum((x - 0.5)^2)
((2 / sqrt(pi)) / (2. * a))^length(x) * exp (-s / (a * a))
}

testFn4_v <- function(x) {
a <- 0.1
r <- apply(x, 2, function(z) {
s <- sum((z - 0.5)^2)
((2 / sqrt(pi)) / (2. * a))^length(z) * exp (-s / (a * a))
})
matrix(r, ncol = ncol(x))
}

d <- harness(f = testFn4, fv = testFn4_v,
lowerLimit = rep(0, 2), upperLimit = rep(1, 2), times = 20)
knitr::kable(d, digits = 3)
expr n.eval min lw.qu median mean up.qu max total relative
Non-vectorized Hcubature 20 0.010 0.011 0.011 0.011 0.012 0.014 0.225 5.71
Vectorized Hcubature 20 0.002 0.002 0.002 0.002 0.002 0.004 0.043 1.00
Non-vectorized Pcubature 20 0.015 0.017 0.018 0.020 0.019 0.058 0.392 9.41
Vectorized Pcubature 20 0.003 0.003 0.003 0.003 0.003 0.004 0.060 1.46
Non-vectorized cubature::cuhre 20 0.022 0.023 0.023 0.025 0.026 0.030 0.494 12.30
Vectorized cubature::cuhre 20 0.004 0.004 0.004 0.005 0.005 0.006 0.090 2.21

## Double Gaussian

testFn5 <- function(x) {
a <- 0.1
s1 <- sum((x - 1 / 3)^2)
s2 <- sum((x - 2 / 3)^2)
0.5 * ((2 / sqrt(pi)) / (2. * a))^length(x) * (exp(-s1 / (a * a)) + exp(-s2 / (a * a)))
}
testFn5_v <- function(x) {
a <- 0.1
r <- apply(x, 2, function(z) {
s1 <- sum((z - 1 / 3)^2)
s2 <- sum((z - 2 / 3)^2)
0.5 * ((2 / sqrt(pi)) / (2. * a))^length(z) * (exp(-s1 / (a * a)) + exp(-s2 / (a * a)))
})
matrix(r, ncol = ncol(x))
}

d <- harness(f = testFn5, fv = testFn5_v,
lowerLimit = rep(0, 2), upperLimit = rep(1, 2), times = 20)
knitr::kable(d, digits = 3)
expr n.eval min lw.qu median mean up.qu max total relative
Non-vectorized Hcubature 20 0.024 0.025 0.026 0.028 0.029 0.034 0.551 7.65
Vectorized Hcubature 20 0.005 0.005 0.005 0.005 0.006 0.010 0.108 1.40
Non-vectorized Pcubature 20 0.016 0.017 0.018 0.018 0.020 0.022 0.368 5.31
Vectorized Pcubature 20 0.003 0.003 0.003 0.004 0.004 0.005 0.073 1.00
Non-vectorized cubature::cuhre 20 0.042 0.043 0.044 0.045 0.046 0.055 0.893 12.60
Vectorized cubature::cuhre 20 0.009 0.009 0.009 0.010 0.010 0.014 0.198 2.68

## Tsuda’s Example

testFn6 <- function(x) {
a <- (1 + sqrt(10.0)) / 9.0
prod( a / (a + 1) * ((a + 1) / (a + x))^2)
}

testFn6_v <- function(x) {
a <- (1 + sqrt(10.0)) / 9.0
r <- apply(x, 2, function(z) prod( a / (a + 1) * ((a + 1) / (a + z))^2))
matrix(r, ncol = ncol(x))
}

d <- harness(f = testFn6, fv = testFn6_v,
lowerLimit = rep(0, 3), upperLimit = rep(1, 3), times = 20)
knitr::kable(d, digits = 3)
expr n.eval min lw.qu median mean up.qu max total relative
Non-vectorized Hcubature 20 0.013 0.014 0.016 0.016 0.017 0.025 0.323 6.67
Vectorized Hcubature 20 0.002 0.002 0.002 0.002 0.002 0.004 0.050 1.00
Non-vectorized Pcubature 20 0.076 0.079 0.081 0.085 0.089 0.116 1.693 34.30
Vectorized Pcubature 20 0.010 0.011 0.011 0.012 0.012 0.015 0.232 4.67
Non-vectorized cubature::cuhre 20 0.033 0.037 0.038 0.038 0.041 0.043 0.766 16.20
Vectorized cubature::cuhre 20 0.005 0.005 0.006 0.006 0.007 0.010 0.122 2.46

## Morokoff & Calflish Example

testFn7 <- function(x) {
n <- length(x)
p <- 1/n
(1 + p)^n * prod(x^p)
}
testFn7_v <- function(x) {
matrix(apply(x, 2, function(z) {
n <- length(z)
p <- 1/n
(1 + p)^n * prod(z^p)
}), ncol = ncol(x))
}

d <- harness(f = testFn7, fv = testFn7_v,
lowerLimit = rep(0, 3), upperLimit = rep(1, 3), times = 20)
knitr::kable(d, digits = 3)
expr n.eval min lw.qu median mean up.qu max total relative
Non-vectorized Hcubature 20 0.028 0.030 0.031 0.031 0.032 0.035 0.630 6.77
Vectorized Hcubature 20 0.004 0.005 0.005 0.005 0.005 0.007 0.100 1.00
Non-vectorized Pcubature 20 0.073 0.076 0.077 0.082 0.081 0.133 1.631 16.70
Vectorized Pcubature 20 0.010 0.010 0.010 0.011 0.012 0.013 0.217 2.25
Non-vectorized cubature::cuhre 20 0.310 0.315 0.321 0.331 0.333 0.414 6.627 69.50
Vectorized cubature::cuhre 20 0.047 0.048 0.050 0.051 0.053 0.063 1.020 10.70

## Wang-Landau Sampling 1d, 2d Examples

I.1d <- function(x) {
sin(4 * x) *
x * ((x * ( x * (x * x - 4) + 1) - 1))
}
I.1d_v <- function(x) {
matrix(apply(x, 2, function(z)
sin(4 * z) *
z * ((z * ( z * (z * z - 4) + 1) - 1))),
ncol = ncol(x))
}
d <- harness(f = I.1d, fv = I.1d_v,
lowerLimit = -2, upperLimit = 2, times = 100)
knitr::kable(d, digits = 3)
expr n.eval min lw.qu median mean up.qu max total relative
Non-vectorized Hcubature 100 0.001 0.001 0.001 0.002 0.001 0.008 0.157 5.17
Vectorized Hcubature 100 0.000 0.000 0.000 0.000 0.000 0.002 0.032 1.06
Non-vectorized Pcubature 100 0.000 0.000 0.000 0.001 0.000 0.002 0.051 1.75
Vectorized Pcubature 100 0.000 0.000 0.000 0.000 0.000 0.000 0.028 1.00
Non-vectorized cubature::cuhre 100 0.002 0.002 0.002 0.002 0.003 0.004 0.247 8.44
Vectorized cubature::cuhre 100 0.001 0.001 0.001 0.001 0.001 0.002 0.086 3.03
I.2d <- function(x) {
x1 <- x; x2 <- x
sin(4 * x1 + 1) * cos(4 * x2) * x1 * (x1 * (x1 * x1)^2 - x2 * (x2 * x2 - x1) +2)
}
I.2d_v <- function(x) {
matrix(apply(x, 2,
function(z) {
x1 <- z; x2 <- z
sin(4 * x1 + 1) * cos(4 * x2) * x1 * (x1 * (x1 * x1)^2 - x2 * (x2 * x2 - x1) +2)
}),
ncol = ncol(x))
}
d <- harness(f = I.2d, fv = I.2d_v,
lowerLimit = rep(-1, 2), upperLimit = rep(1, 2), times = 100)
knitr::kable(d, digits = 3)
expr n.eval min lw.qu median mean up.qu max total relative
Non-vectorized Hcubature 100 0.047 0.049 0.050 0.052 0.052 0.098 5.156 61.50
Vectorized Hcubature 100 0.005 0.006 0.006 0.006 0.007 0.013 0.625 7.16
Non-vectorized Pcubature 100 0.004 0.004 0.004 0.004 0.005 0.006 0.442 5.01
Vectorized Pcubature 100 0.001 0.001 0.001 0.001 0.001 0.003 0.090 1.00
Non-vectorized cubature::cuhre 100 0.010 0.010 0.011 0.011 0.012 0.018 1.125 13.40
Vectorized cubature::cuhre 100 0.002 0.002 0.002 0.002 0.002 0.004 0.182 2.00

## Implementation Notes

About the only real modification we have made to the underlying cubature-1.0.2 library is that we use M = 16 rather than the default M = 19 suggested by the original author for pcubature. This allows us to comply with CRAN package size limits and seems to work reasonably well for the above tests. Future versions will allow for such customization on demand.

The changes made to the Cuba-4.2 library are separated out so that the original source is left pristine and the changes (a small number) applied only during the building process.

## Summary

My recommendations are:

1. Vectorize your function. The time spent in so doing pays back enormously. This is easy to do and the examples above show how.

2. Vectorized hcubature seems to be a good starting point.

3. For smooth integrands in low dimensions ($$\leq 3$$), pcubature might be worth trying out. Experiment before using in a production package.

## Session Info

sessionInfo()
## R version 4.2.2 (2022-10-31)
## Platform: x86_64-apple-darwin21.6.0 (64-bit)
## Running under: macOS Ventura 13.1
##
## Matrix products: default
## BLAS:   /usr/local/Cellar/openblas/0.3.21/lib/libopenblasp-r0.3.21.dylib
## LAPACK: /usr/local/Cellar/r/4.2.2/lib/R/lib/libRlapack.dylib
##
## locale:
##  en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
##  stats     graphics  grDevices utils     datasets  methods   base
##
## other attached packages:
##  mvtnorm_1.1-3    cubature_2.0.4.6 benchr_0.2.5
##
## loaded via a namespace (and not attached):
##   Rcpp_1.0.9         knitr_1.41         magrittr_2.0.3     R6_2.5.1
##   rlang_1.0.6        fastmap_1.1.0      stringr_1.5.0      highr_0.10
##   RcppProgress_0.4.2 tools_4.2.2        xfun_0.36          cli_3.6.0
##  jquerylib_0.1.4    htmltools_0.5.4    yaml_2.3.6         digest_0.6.31
##  lifecycle_1.0.3    sass_0.4.4         vctrs_0.5.1        glue_1.6.2
##  cachem_1.0.6       evaluate_0.20      rmarkdown_2.18     stringi_1.7.12
##  compiler_4.2.2     bslib_0.4.2        jsonlite_1.8.4