Peterlin computes a 4096-dimensional Gaussian rectangle probability in 5.1061 seconds using 245,760 quasi-Monte Carlo points. On a desktop Mac, MvNormalCDF.jl needs 1081.629 seconds for the same problem, while tlrmvnmvt's Genz-Bretz routine takes 515.286 seconds.
The mathematics is unchanged.
The target is P(a < X <= b) for X ~ N(0, Sigma). Beyond trivial dimensions it has no closed form, yet it appears throughout Gaussian-copula credit and portfolio models, multivariate probit and simultaneous inference. Genz's method factors Sigma, reorders variables by interval width, maps the region onto the unit cube, then averages a product of conditional increments over randomized QMC points.
Most of the expense hides in the bound calculation. At step i, the conditional bounds depend on z_1 through z_{i-1}. Each step therefore requires a dot product of length i-1 for every QMC point and every replication. Across all steps, the operation becomes a rectangular matrix multiplied by a triangular one. Peterlin's complexity table gives the full algorithm as (1/3)n^3 + (1/2)n^2 m nr mul-adds. For large n, the pivoted Cholesky factorization and the bound update dominate. Since m exceeds n in any serious integration problem, Peterlin identifies the multiplication step as often the largest opportunity for optimization.
What does double blocking buy?
Peterlin leaves Genz-Bretz intact and changes how memory is accessed. His doubly blocked variant, DB, uses nr replications, usually 12 to 20. Parallel threads process the m QMC points in blocks of BJ, from 2^6 to 2^7. The dimension n is divided at two nested levels, with BI around 2^9 to 2^10 and BI2 around 2^6.
Once DB has computed a set of contributions inside a block, a BLAS matrix-matrix multiply prepares the bound updates owed to every later point. Phi and Phi inverse pass through erf and erfinv on contiguous vectors. Each thread receives preallocated buffers, while BLAS stays pinned to one thread inside the loop to prevent oversubscription. Parallel work is assigned across QMC points. The conditioning recursion within each point remains sequential, and Cao et al. already observed that it resists parallelization.
DB leads every reported configuration: n from 2^4 to 2^12, point budgets of 24,576 and 245,760, on both machines. The timings are minimums over 100 runs. They include Cholesky, reordering and the main loop, while package loading and first-use compilation are excluded.
The accuracy claim is parity. At n = 4096 with 245,760 points, using random s.p.d. matrices on the M2 Ultra, DB reports 0.2514 with SD 1.28e-05 and mean error 4.07e-05. MvNormalCDF.jl records SD 1.32e-05 and error 3.85e-05; tlrmvnmvt records SD 1.92e-05 and error 3.72e-05. In that setting on that machine, DB produces the tightest spread and the largest mean error among the three.
Where the ratio comes from
The M2 Ultra used 16 threads. Its 100.92x advantage over tlrmvnmvt works out to about 6x per core. The Xeon used 56 threads, making 471.04x about 8x. The paper gives core counts only for its own runs, so we do not know how the comparators were threaded.
Peterlin also ran the R and Julia comparators with Apple Accelerate and MKL. He reports that MKL substantially sped up MvNormalCDF.jl's dot product on the Xeon, and argues that the competing implementations used those libraries far less effectively. We read the result the same way: the block layout presents BLAS with enough work to matter.
Peterlin qualifies the headline himself. The abstract reports speedups "often close to or above 100x" over mvtnorm and tlrmvnmvt for large dimensions "in the tested settings". The conclusion limits the claim to some tested configurations on two machines, with block sizes tuned by hand. Its remaining defence is borne out by the tables: those settings show no systematic loss of precision, with further gains when lower precision is acceptable.
The conclusion gives priority to the claim we think will travel, namely that runtime approaches the cost of the dominant linear-algebra operations. We largely agree with Peterlin's framing. The lingering question is why the abstract foregrounds the multiple.
Hardware changes the comparator ranking. At n = 4096, MvNormalCDF.jl is slowest on the M2 Ultra at 211.83x, while tlrmvnmvt is slowest on the Xeon at 471.04x. DB takes 3.0676 s on the Xeon and 5.1061 s on the Mac. Peterlin states that every comparison was made within its own system and that these timings are not intended as controlled cross-system comparisons.
Block choices also vary by machine. The Mac uses Cholesky levels 2^5/2^7/2^10/2^11, against 2^9/2^9/2^10/2^11 on the Xeon, while BI changes from 2^10 to 2^9. Peterlin tuned them by hand. The Mac also benefits from Accelerate reaching the AMX coprocessor. At vector length 2^20, close to the simulations here, the M2 Ultra retained more than 75% of measured peak bandwidth; the Xeon managed under 25%.
The tlrmvnmvt ratios require a code check before quotation because the budget conventions differ. Peterlin discloses that tlrmvnmvt interprets N as points per repetition and performs 20 repetitions. By contrast, mvtnorm and MvNormalCDF.jl interpret qmc_pts as the total. Figure captions describe the plotted budgets as total QMC points including repetitions. The repository script contains the per-call argument supplied to tlrmvnmvt, leaving this as a code-level question that the text does not settle.
The BLAS comparison travels
The benchmark against linear algebra says more than the speedup table. When a rectangular-by-triangular multiply dominates cost, an implementation running near that multiply has captured most of the available performance. Peterlin compares the DB loop with in-place BLAS triangular multiplication at 7 * 2^11 QMC points, a count divisible by both 16 and 56 cores. For large n, the timings are close and DB is occasionally faster.
On the Xeon, the DB loop can beat the OpenBLAS triangular multiply. Peterlin attributes this to more direct AVX-512 use and describes the result as "an empirical performance comparison rather than a formal proof of optimality." His hand-blocked Cholesky variants run close to, and sometimes below, Accelerate, MKL and OpenBLAS at large n. Those variants still do extra work by scaling Sigma and the bounds, and they incur the copy cost.
Accuracy rests on a narrower base than timing. The independent reference is a one-dimensional integral evaluated to relative tolerance 1e-12, and it applies only to the equicorrelated Sigma2 = I + 11^T. Peterlin explicitly says this cannot establish accuracy for arbitrary covariance matrices. Individual estimates were not retained. Only the median and empirical standard deviation remain, leaving the estimator's tail behaviour uncharacterized.
The FP32 route is typically 1.5 to 2x faster than FP64, with observed gains ranging from about 1.3 to 2.4x. Sparsity pushes further. For Sigma4 at n = 4096 with 245,760 points on the M2 Ultra, DB-Sparse-FP32 finishes in 0.9216 s, against 4.5619 s for DB.
The sparse variants preserve sparsity by dropping pivoting and reordering. Peterlin warns that this may be less reliable for ill-conditioned problems, and it forfeits the variance reduction from reordering. Sigma3 shows the cost at the same n and budget: DB has mean error 1.10e-03 with SD 4.08e-04, compared with 4.55e-03 and SD 1.53e-03 for DB-Sparse-FP32. Sigma4 reverses the ordering, with 3.93e-07 for DB-Sparse-FP32 and 6.97e-07 for DB.
A copula stopwatch, carefully qualified
Peterlin presents the finance section as a stopwatch exercise. It uses daily Yahoo Finance prices from 1 January 2000 through 31 December 2025, covering roughly 500 S&P 500 names and more than 3300 of about 4200 FTSE All-World constituents, representing over 90% of index value. Each asset's log returns are rank-transformed into normal scores, whose covariance is then estimated. The calculation asks for the probability that the largest one-day gain remains below a threshold.
Runtime falls from hours to minutes. Peterlin says the example benchmarks only the probability-evaluation step and does not represent a full production risk-monitoring workflow. He also cites Embrechts et al. on criticism of the Gaussian copula while observing that industry still uses it.
A qualification absent from his discussion matters for any tail number built from the same matrix. Membership is fixed as of 2025-12-31, with history extending back to 2000, so the fitted correlation matrix contains survivors. That choice is adequate for a timing demonstration. About 900 index constituents are missing from the FTSE All-World pull, and no accuracy check or alternative solver is run at that scale.
We could not run this on our own data because the paper contains no strategy to test. Its text does not provide the 2025-12-31 constituent lists, identify which 900 FTSE names failed to download, or print the threshold grid used for the curves. Those details remain in the public repository scripts, which the author says reproduce the figures.
A third machine, using a different BLAS and untuned block sizes, would change our view if it landed far from triangular-multiply time. For now, the transferable result is the proximity to that multiply, the same claim Peterlin places first in his conclusion. The 100x belongs to two well-provisioned boxes whose block sizes were selected by hand.