QMCPy: A Python Software for Randomized Low-Discrepancy Sequences, Quasi-Monte Carlo, and Fast Kernel Methods

TL;DR

QMCPy addresses the need for a unified, Python-based toolset for randomized low-discrepancy sequences and fast kernel computations. It unifies rank-1 lattices, digital nets, and Halton sequences with versatile randomizations (shifts, LMS, NUS) and introduces higher-order scramblings and higher-order DSI kernels, enabling $O(n\log n)$ kernel operations via FFTBR/IFFTBR and FWHT. A key contribution is the LDData repository and the first Python interfaces supporting higher-order nets and scramblings, plus efficient eigenvalue-update procedures for doubling problem sizes. The work demonstrates rapid, scalable QMC and kernel-based computations, offering practical impact for high-dimensional integration, Bayesian cubature, and kernel interpolation in scientific computing.

Abstract

Low-discrepancy (LD) sequences have been extensively used as efficient experimental designs across many scientific disciplines. QMCPy (https://qmcsoftware.github.io/QMCSoftware/) is an accessible Python library which provides a unified implementation of randomized LD sequences, automatic variable transformations, adaptive Quasi-Monte Carlo error estimation algorithms, and fast kernel methods. This article focuses on recent updates to QMCPy which broaden support for randomized LD sequences and add new tools to enable fast kernel methods using LD sequences. Specifically, we give a unified description of the supported LD lattices, digital nets, and Halton point sets, along with randomization options including random permutations / shifts, linear matrix scrambling (LMS), and nested uniform scrambling (NUS). We also support higher-order digital nets, higher-order scrambling with LMS or NUS, and Halton scrambling with LMS or NUS. For fast kernel methods, we provide shift-invariant (SI) and digitally-shift-invariant (DSI) kernels, including a new set of higher-order smoothness DSI kernels. When SI and DSI kernels are respectively paired with n LD lattice and digital net points, the resulting Gram matrices permit multiplication and inversion at only O(n log n) cost. These fast operations utilize QMCPy's implementation of the fast Fourier transform in bit-reversed order (FFTBR), inverse FFTBR (IFFTBR), and fast Walsh--Hadamard transform (FWHT).

Figures (5)

• Figure 1: An independent identically distributed (IID) point set and low-discrepancy (LD) point sets of size $n=2^{13}=8192$. Notice the LD points more evenly fill the space than IID points, which leads to faster convergence of Quasi-Monte Carlo methods compared to IID Monte Carlo methods. Randomized rank-1 lattices, digital nets in base $2$, and Halton points are shown. Randomizations include shifts, linear matrix scramblings (LMS), digital shifts (DS), digital permutations (DP), and nested uniform scramblings (NUS). Digital interlacing of order $\alpha$ is used to create higher-order randomized digital nets in base $2$ (DN${}_\alpha$).
• Figure 2: Low-discrepancy randomization routines in dimension $d=1$. Each vertical line is a unit interval with $0$ at the bottom and $1$ at the top. A given interval is partitioned at the horizontal ticks extending to the right, and then rearranged following the arrows to create the partition of the right interval as shown by the ticks extending to the left. For the shifted rank-1 lattice and digitally shifted digital net, when the arrow hits a horizontal gray bar it is wrapped around to the next gray bar below. See for example the dotted line in the digitally shifted digital net. The lattice shift is $\Delta = 0.227$. All digital nets (DNs) use base $b=3$ and $t_\mathrm{max}=2$ digits of precision. The digital shift is ${\boldsymbol{\Delta}} = (2,1)^\intercal$. Dropping the $j$ subscript for dimension, the digital permutations in the third panel are $\pi_1 = (2,0,1)$ and $\pi_2 = (2,1,0)$. Notice $\pi_1$ is equivalent to a digital shift by $2$, but $\pi_2$ cannot be written as a digital shift. The nested uniform scramble (NUS) has digital permutations $\pi = (2,0,1)$, $\pi_0 = (2,1,0)$, $\pi_1 = (0,1,2)$, and $\pi_2 = (1,2,0)$. Notice the permuted digital net in the third panel has permutations depending only on the current digit in the base $b$ expansion. In contrast, the full NUS scrambling in the fourth panel has permutations which depend on all previous digits in the base $b$ expansion.
• Figure 3: Bit-reversed FFT and FWHT schemes. The bit-reversed FFT is performed via a decimation in time algorithm without the initial reordering of inputs. The bit-reversed IFFT may be performed by propagating $\{{\widetilde{y}}_i\}_{i=0}^{2^m-1}$ to $\{y_i\}_{i=0}^{2^m-1}$ in the other direction through the bit-reversed FFT algorithm.
• Figure 4: Comparison of time required to generate point sets and perform fast transforms. For IID and randomized LD generators, we vary the number of randomizations $R$, the number of dimensions $d$, and the sequence size $n$. Timings include initialization, and a new point set is generating for each plotted $(R,n,d)$. Fast transforms are applied to $r$ sequences of size $n$ in a vectorized fashion.
• Figure 5: The RMSE of the RQMC estimate for a few different integrands. Higher-order digital nets with linear matrix scrambling and digital shifts achieve higher-order convergence and significantly outperform both IID points and randomly shifted lattices.

Theorems & Definitions (8)

• Theorem 6.1
• proof
• Theorem 6.2
• proof
• proof : Proof of \ref{['thm:RKHS_DSI_contain']}
• lemma B.1: Walsh coefficients of low order monomials
• proof
• proof : Proof of \ref{['thm:explicit_DSI_low_order_forms']}