Doubly blocked variant of Genz’s method for multivariate normal integration

TL;DR

This paper presents a doubly blocked (DB) variant of Genz’s method for multivariate normal integration, implemented in Julia with multithreading, cache blocking, and SIMD/BLAS-friendly design. In the tested settings it is often close to or above 100× faster than existing solvers such as mvtnorm and tlrmvnmvt for large dimensions, with optional FP32 and sparse variants.

Key findings

• The DB algorithm was fastest for all tested matrix sizes; for \(n > 2^9\) it was often close to or above 100× faster than the other implementations that could handle larger matrices. For example, on the M2 Ultra at \(n=4096\), \(m=245{,}760\), DB took 5.1061 s versus 1081.6290 s for MvNormalCDF.jl (211.83×) and 515.2860 s for tlrmvnmvt::GenzBretz (100.92×).
• DB-FP32 was typically about 1.5–2× faster than DB, with observed speedups from approximately 1.3 to 2.4; DB-Sparse-FP32 was up to approximately 3× faster than DB-FP32. FP32 machine epsilon is approximately \(1.19\times10^{-7}\) versus \(2.22\times10^{-16}\) for FP64.
• The tile-low-rank (TLR) approach was slower than DB but substantially faster than the classical Genz implementation in tlrmvnmvt without low-rank approximations. In an independent equicorrelated accuracy check, FP32 variants closely tracked their FP64 counterparts, and TLR had errors of broadly similar order but did not uniformly dominate.
• In a real-data Gaussian-copula example using S&P 500 and FTSE All-World data, the DB algorithm estimated high-dimensional portfolio probabilities in a few minutes rather than a few hours; more than 90% of the FTSE All-World index value was captured.

Methods

• The paper implements Genz’s algorithm for \(P(a < X \le b)\) with \(X \sim N(0,\Sigma)\), using Cholesky factorization, variable reordering, unit-cube transformation, and randomized QMC (Sobol with Matoušek affine scrambling or Richtmyer with Cranley–Patterson rotations), repeated \(n_r\) times. The DB variant reorganizes the same operations into blocks over QMC points and dimensions, exposing large matrix multiplications and using BLAS, multithreading, SIMD, and cache-aware data layout.
• Benchmarks used an Apple M2 Ultra (16 threads) and a dual-socket Intel Xeon Gold 5420+ (56 cores), with Julia 1.12.2/1.12.4 and R 4.5.2/4.4, using Apple Accelerate or Intel MKL. Simulations used covariance constructions \(\Sigma_1\)–\(\Sigma_4\), repeated at least 100 times, reporting minimum execution time; comparisons included mvtnorm, tlrmvnmvt, MvNormalCDF.jl, and TLR.
• The real-data example used Yahoo Finance data for S&P 500 and FTSE All-World stocks from 2000 to 2025, transformed log-returns to normal scores by rank, estimated the covariance matrix, and computed \(p(b)=P(R_i \le b, i \in P)\).

Limitations

• The paper states it focuses on positive definite covariance matrices; singular-covariance extensions such as Genz–Kwong are not covered. FP32 variants trade numerical accuracy for speed, and sparse variants use a standard non-pivoted Cholesky decomposition, which may be less numerically robust for ill-conditioned problems.
• The real-data example is not a full production risk-monitoring workflow, and the independent accuracy check uses the special equicorrelated structure \(\Sigma_2 = I + 11^T\), so it does not establish accuracy for arbitrary covariance matrices.
• The text notes that performance comparisons were performed within each system and are not controlled cross-system comparisons; mvtnorm can handle covariance matrices only up to size 1000, so larger sizes compare only the other implementations. The performance analysis is described as empirical rather than a formal proof of optimality.