Cubature Vectorization Results

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(bench)
library(cubature)

harness <- function(which = NULL,
                    f, fv, lowerLimit, upperLimit, tol = 1e-3, iterations = 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, iterations = iterations, check = FALSE)
  result <- do.call(bench::mark, args = argList)
  d <- result[seq_along(fnIndices), 1:9]
  d$expression <- fns[fnIndices]
  d
}

We reel off the timing runs.

Example 1.

func <- function(x) sin(x[1]) * cos(x[2]) * exp(x[3])
func.v <- function(x) {
    matrix(apply(x, 2, function(z) sin(z[1]) * cos(z[2]) * exp(z[3])), ncol = ncol(x))
}

d <- harness(f = func, fv = func.v,
             lowerLimit = rep(0, 3),
             upperLimit = rep(1, 3),
             tol = 1e-5,
             iterations = 100)[, 1:9]

knitr::kable(d, digits = 3, row.names = FALSE)

expression	min	median	itr/sec	mem_alloc	gc/sec	n_itr	n_gc	total_time
Non-vectorized Hcubature	211µs	256.39µs	3818.791	59.3KB	38.574	99	1	25.9ms
Vectorized Hcubature	273µs	339.69µs	2685.312	54.5KB	27.124	99	1	36.9ms
Non-vectorized Pcubature	784µs	887.01µs	1065.711	31.2KB	21.749	98	2	92ms
Vectorized Pcubature	933µs	1.05ms	887.303	58.4KB	18.108	98	2	110.4ms
Non-vectorized cubature::cuhre	495µs	588.56µs	1590.889	40.5KB	32.467	98	2	61.6ms
Vectorized cubature::cuhre	442µs	515.58µs	1909.896	39.8KB	19.292	99	1	51.8ms

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,
             iterations = 10,
             mean = rep(0, m), sigma = sigma, logdet = logdet)
knitr::kable(d, digits = 3)

expression	min	median	itr/sec	mem_alloc	gc/sec	n_itr	n_gc	total_time
Non-vectorized Hcubature	554.08ms	621.3ms	1.566	139.68KB	22.078	10	141	6.39s
Vectorized Hcubature	1.99ms	2.46ms	409.183	1.89MB	0.000	10	0	24.44ms
Non-vectorized Pcubature	239.33ms	263.45ms	3.571	0B	21.427	10	60	2.8s
Vectorized Pcubature	885.57µs	940.04µs	1057.400	810.26KB	0.000	10	0	9.46ms
Non-vectorized cubature::cuhre	224.86ms	240.06ms	4.025	0B	23.347	10	58	2.48s
Vectorized cubature::cuhre	2.1ms	2.27ms	412.961	898.41KB	0.000	10	0	24.21ms

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 vectorized versions of hcubature and pcubature seem competitive and in some cases better, if you compare the table above to the one below.

library(mvtnorm)
g1 <- function() 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() 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() 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(bench::mark(g1(), g2(), g3(), iterations = 20, check = FALSE)[, 1:9],
             digits = 3, row.names = FALSE)

expression	min	median	itr/sec	mem_alloc	gc/sec	n_itr	n_gc	total_time
g1()	646µs	1.17ms	729.768	355.36KB	0.000	20	0	27.4ms
g2()	597µs	1.1ms	897.973	2.49KB	47.262	19	1	21.2ms
g3()	873µs	1.57ms	645.849	2.49KB	0.000	20	0	31ms

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), iterations = 1000)
knitr::kable(d, digits = 3)

expression	min	median	itr/sec	mem_alloc	gc/sec	n_itr	n_gc	total_time
Non-vectorized Hcubature	25.2µs	30µs	30561.378	7.66KB	30.592	999	1	32.7ms
Vectorized Hcubature	39.3µs	74.8µs	13836.817	1.16KB	13.851	999	1	72.2ms
Non-vectorized Pcubature	33.9µs	40.9µs	22470.192	0B	45.030	998	2	44.4ms
Vectorized Pcubature	70.3µs	84.9µs	10975.450	18.68KB	21.995	998	2	90.9ms
Non-vectorized cubature::cuhre	251.7µs	310.4µs	2972.349	0B	26.994	991	9	333.4ms
Vectorized cubature::cuhre	237.1µs	308.3µs	3023.119	16.38KB	24.380	992	8	328.1ms

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), iterations = 10)

knitr::kable(d, digits = 3)

expression	min	median	itr/sec	mem_alloc	gc/sec	n_itr	n_gc	total_time
Non-vectorized Hcubature	2.21ms	2.46ms	404.931	67.3KB	44.992	9	1	22.23ms
Vectorized Hcubature	3.75ms	4.29ms	214.061	290.5KB	23.785	9	1	42.04ms
Non-vectorized Pcubature	64.87µs	66.87µs	13712.064	0B	0.000	10	0	729.28µs
Vectorized Pcubature	110.62µs	127.32µs	7413.485	4.1KB	0.000	10	0	1.35ms
Non-vectorized cubature::cuhre	10.98ms	12.03ms	81.176	0B	34.790	7	3	86.23ms
Vectorized cubature::cuhre	15.14ms	15.95ms	60.670	971.5KB	60.670	5	5	82.41ms

Discontinuous function

testFn2 <- function(x) {
    radius <- 0.50124145262344534123412
    ifelse(sum(x * x) < radius * radius, 1, 0)
}

testFn2_v <- function(x) {
    radius <- 0.50124145262344534123412
    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), iterations = 10)
knitr::kable(d, digits = 3)

expression	min	median	itr/sec	mem_alloc	gc/sec	n_itr	n_gc	total_time
Non-vectorized Hcubature	38.5ms	47.5ms	20.960	17.7KB	37.727	10	18	477.11ms
Vectorized Hcubature	43.9ms	47.9ms	20.409	1011.8KB	32.655	10	16	489.97ms
Non-vectorized cubature::cuhre	790.2ms	829.9ms	1.173	0B	25.802	10	220	8.53s
Vectorized cubature::cuhre	840.4ms	948ms	1.047	21.2MB	21.471	10	205	9.55s

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), iterations = 50)
knitr::kable(d, digits = 3)

expression	min	median	itr/sec	mem_alloc	gc/sec	n_itr	n_gc	total_time
Non-vectorized Hcubature	40.8µs	45.7µs	20097.518	6.75KB	0.000	50	0	2.49ms
Vectorized Hcubature	58.2µs	64.1µs	14015.832	2.91KB	0.000	50	0	3.57ms
Non-vectorized Pcubature	36.8µs	40µs	21879.119	0B	0.000	50	0	2.29ms
Vectorized Pcubature	50.2µs	55.3µs	16447.019	19.12KB	0.000	50	0	3.04ms
Non-vectorized cubature::cuhre	509µs	564.2µs	1704.158	0B	34.779	49	1	28.75ms
Vectorized cubature::cuhre	461.5µs	539.7µs	1770.527	39.84KB	36.133	49	1	27.68ms

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), iterations = 20)
knitr::kable(d, digits = 3)

expression	min	median	itr/sec	mem_alloc	gc/sec	n_itr	n_gc	total_time
Non-vectorized Hcubature	1.15ms	1.34ms	719.993	85.2KB	37.894	19	1	26.4ms
Vectorized Hcubature	1.53ms	2.4ms	440.772	147.2KB	23.199	19	1	43.1ms
Non-vectorized Pcubature	1.9ms	2.43ms	411.197	0B	21.642	19	1	46.2ms
Vectorized Pcubature	2.13ms	2.36ms	394.728	68.9KB	20.775	19	1	48.1ms
Non-vectorized cubature::cuhre	2.57ms	2.89ms	320.457	0B	35.606	18	2	56.2ms
Vectorized cubature::cuhre	2.91ms	3.14ms	305.316	125.6KB	16.069	19	1	62.2ms

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), iterations = 20)
knitr::kable(d, digits = 3)

expression	min	median	itr/sec	mem_alloc	gc/sec	n_itr	n_gc	total_time
Non-vectorized Hcubature	2.96ms	4.73ms	227.902	133KB	25.322	18	2	79ms
Vectorized Hcubature	3.9ms	4.16ms	236.781	249KB	26.309	18	2	76ms
Non-vectorized Pcubature	2.12ms	2.27ms	405.944	0B	45.105	18	2	44.3ms
Vectorized Pcubature	2.74ms	2.85ms	347.183	69KB	18.273	19	1	54.7ms
Non-vectorized cubature::cuhre	5.73ms	5.93ms	161.207	0B	40.302	16	4	99.3ms
Vectorized cubature::cuhre	6.81ms	7.44ms	127.969	224KB	31.992	16	4	125ms

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), iterations = 20)
knitr::kable(d, digits = 3)

expression	min	median	itr/sec	mem_alloc	gc/sec	n_itr	n_gc	total_time
Non-vectorized Hcubature	1.43ms	1.68ms	541.147	64.6KB	28.481	19	1	35.1ms
Vectorized Hcubature	1.71ms	1.85ms	525.652	156KB	27.666	19	1	36.1ms
Non-vectorized Pcubature	6.79ms	7.36ms	126.335	0B	31.584	16	4	126.6ms
Vectorized Pcubature	8.73ms	9.51ms	102.459	386.2KB	34.153	15	5	146.4ms
Non-vectorized cubature::cuhre	3.6ms	3.98ms	240.702	0B	26.745	18	2	74.8ms
Vectorized cubature::cuhre	3.48ms	4.22ms	213.223	225.8KB	23.691	18	2	84.4ms

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), iterations = 20)
knitr::kable(d, digits = 3)

expression	min	median	itr/sec	mem_alloc	gc/sec	n_itr	n_gc	total_time
Non-vectorized Hcubature	3.04ms	3.81ms	234.419	32.89KB	23.442	20	2	85.3ms
Vectorized Hcubature	3.08ms	3.45ms	234.544	205.61KB	23.454	20	2	85.3ms
Non-vectorized Pcubature	7.1ms	8.76ms	103.834	0B	25.958	20	5	192.6ms
Vectorized Pcubature	7.46ms	8.61ms	102.982	386.24KB	25.745	20	5	194.2ms
Non-vectorized cubature::cuhre	38.83ms	42.51ms	22.694	0B	23.829	20	21	881.3ms
Vectorized cubature::cuhre	33.59ms	38.09ms	22.194	2.04MB	23.304	20	21	901.1ms

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, iterations = 100)
knitr::kable(d, digits = 3)

expression	min	median	itr/sec	mem_alloc	gc/sec	n_itr	n_gc	total_time
Non-vectorized Hcubature	95.5µs	106µs	9065.648	53.7KB	0.000	100	0	11.03ms
Vectorized Hcubature	140.4µs	159µs	6117.713	68.5KB	0.000	100	0	16.35ms
Non-vectorized Pcubature	38.9µs	44µs	20568.407	0B	207.762	99	1	4.81ms
Vectorized Pcubature	95.5µs	133µs	6856.702	0B	0.000	100	0	14.58ms
Non-vectorized cubature::cuhre	176.1µs	196µs	4963.207	0B	0.000	100	0	20.15ms
Vectorized cubature::cuhre	330.6µs	386µs	2358.810	0B	23.826	99	1	41.97ms

I.2d <- function(x) {
    x1 <- x[1]; x2 <- x[2]
    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[1]; x2 <- z[2]
                     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), iterations = 100)
knitr::kable(d, digits = 3)

expression	min	median	itr/sec	mem_alloc	gc/sec	n_itr	n_gc	total_time
Non-vectorized Hcubature	3.71ms	4.09ms	222.437	78.4KB	30.332	88	12	395.6ms
Vectorized Hcubature	4.21ms	4.63ms	211.668	304.7KB	28.864	88	12	415.7ms
Non-vectorized Pcubature	317.35µs	358.43µs	2715.608	0B	27.430	99	1	36.5ms
Vectorized Pcubature	417.82µs	496.99µs	2017.410	18.3KB	20.378	99	1	49.1ms
Non-vectorized cubature::cuhre	999.18µs	1.13ms	784.790	0B	24.272	97	3	123.6ms
Vectorized cubature::cuhre	1.01ms	1.17ms	802.998	60.1KB	24.835	97	3	120.8ms

Implementation Notes

About the only real modification we have made to the underlying cubature 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 library are managed in a Github repo branch: each time a new release is made, we update the main branch, and keep all changes for Unix platforms in a branch named R_pkg against the current main branch. Customization for windows is done in the package itself using the Makevars.win script.

Summary

The recommendations are:

Vectorize your function. The time spent in so doing pays back enormously. This is easy to do and the examples above show how.
Vectorized hcubature seems to be a good starting point.
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.5.1 (2025-06-13)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.3 LTS
## 
## Matrix products: default
## BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
## 
## locale:
##  [1] LC_CTYPE=C.UTF-8       LC_NUMERIC=C           LC_TIME=C.UTF-8       
##  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
##  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
## [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
## 
## time zone: Etc/UTC
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
## [1] mvtnorm_1.3-3    cubature_2.1.4-1 bench_1.1.4     
## 
## loaded via a namespace (and not attached):
##  [1] vctrs_0.6.5       cli_3.6.5         knitr_1.50        rlang_1.1.6      
##  [5] xfun_0.54         textshaping_1.0.4 jsonlite_2.0.0    glue_1.8.0       
##  [9] htmltools_0.5.8.1 ragg_1.5.0        sass_0.4.10       rmarkdown_2.30   
## [13] tibble_3.3.0      evaluate_1.0.5    jquerylib_0.1.4   fastmap_1.2.0    
## [17] profmem_0.7.0     yaml_2.3.10       lifecycle_1.0.4   compiler_4.5.1   
## [21] fs_1.6.6          pkgconfig_2.0.3   Rcpp_1.1.0        systemfonts_1.3.1
## [25] digest_0.6.39     R6_2.6.1          pillar_1.11.1     magrittr_2.0.4   
## [29] bslib_0.9.0       tools_4.5.1       pkgdown_2.2.0     cachem_1.1.0     
## [33] desc_1.4.3

Balasubramanian Narasimhan

2025-12-01

Introduction

A Timing Harness

Example 1.

Multivariate Normal

Product of cosines

Gaussian function

Discontinuous function

A Simple Polynomial (product of coordinates)

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

Double Gaussian

Tsuda’s Example

Morokoff & Calflish Example

Wang-Landau Sampling 1d, 2d Examples

Implementation Notes

Summary

Session Info

Cubature Vectorization Results

Balasubramanian Narasimhan

2025-12-01

Introduction

A Timing Harness

Example 1.

Multivariate Normal

Product of cosines

Gaussian function

Discontinuous function

A Simple Polynomial (product of coordinates)

Gaussian centered at 12\frac{1}{2}

Double Gaussian

Tsuda’s Example

Morokoff & Calflish Example

Wang-Landau Sampling 1d, 2d Examples

Implementation Notes

Summary

Session Info

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