Multilevel RQMC Estimators

• Multilevel RQMC estimators are defined via telescoping decompositions that combine MLMC variance reduction with high-quality RQMC quadrature for smooth, high-dimensional integrands.
• They leverage randomized lattice and digital net rules to achieve superior convergence rates and reduce the computational complexity relative to standard MLMC methods.
• Practical applications include PDE-constrained optimization, Bayesian inverse problems, and financial risk, with adaptive sample allocations and CBC constructions enhancing efficiency.

Multilevel randomized quasi-Monte Carlo (MLRQMC) estimators constitute a class of high-performance numerical integration techniques that combine multilevel Monte Carlo (MLMC) methodology with quasi-Monte Carlo (QMC) or randomized QMC (RQMC) sampling on each level. They are designed to accelerate the convergence rate and reduce computational complexity for expectations arising in high-dimensional PDEs with random coefficients, Bayesian inverse problems, variational inference under simulation-based likelihoods, nested Monte Carlo integration, financial risk quantification, and PDE-constrained optimization under uncertainty. The primary advantage of MLRQMC estimators is the synergy between the variance reduction from multilevel telescoping decompositions and the superior convergence properties of (R)QMC quadrature for smooth, high-dimensional integrands.

1. Multilevel RQMC Framework and Telescoping Construction

The foundation of MLRQMC is the multilevel telescoping decomposition, which expresses an expectation of a target quantity computed at the finest discretization (or with highest computational cost) as a sum of expectations of differences between increasingly finer levels:

$$\mathbb{E}[Q_{L}] = \mathbb{E}[Q_{0}] + \sum_{\ell=1}^{L} \mathbb{E}[Q_{\ell} - Q_{\ell-1}],$$

where $Q_{\ell}$ denotes the approximation on level $\ell$ (for example, finite element mesh with mesh width $h_\ell$, or inner sample size $m_\ell$ in nested quadrature) (Cui et al., 2023, Kuo et al., 2015, Bartuska et al., 2024). Each level correction is then estimated by a (randomized) QMC quadrature rule, typically using high-quality lattice rules or digital nets with random shifts:

$$Q_{N_\ell}^{(\ell)}(f) = \frac{1}{N_\ell} \sum_{k=0}^{N_\ell-1} f(t_{\ell,k}),$$

with QMC points $t_{\ell,k}$ built from low-discrepancy sequences after mapping from the unit cube, ensuring unbiasedness in the randomized setting.

The overarching estimator is then:

$$Y^{MLRQMC} = \sum_{\ell=0}^L Q_{N_\ell}^{(\ell)}\left(Q_\ell - Q_{\ell-1}\right),$$

and variants may target ratios or include splitting-type formulations for Bayesian evidence estimation (Dick et al., 2016).

2. RQMC Rule Construction and Theoretical Error Bounds

On each level, the (R)QMC quadrature employs rules such as randomly shifted rank-1 lattice rules or high-order digital nets:

• For lattice rules in $s$ dimensions, each rule is determined by a generating vector $z_\ell\in\mathbb{Z}^s$ and a random shift $\Delta_\ell\sim \mathrm{Uniform}([0,1]^s)$:

$$t_{\ell, k} = \left\{ \frac{k z_\ell}{N_\ell} + \Delta_\ell \right\},\quad k=0,\dots, N_\ell-1.$$

• For higher-order QMC, interlaced polynomial lattice rules are used, allowing exploitation of higher smoothness (Dick et al., 2016). These rules yield

$$\mathbb{E}_\Delta \left|I(f) - Q_{N_\ell}^{(\ell)}(f)\right|^2 \lesssim N_\ell^{-1/\eta} V_{s,\gamma}(f),$$

for sufficiently smooth $f$, with $\eta \in (1/2, 1]$ and $V_{s,\gamma}$ a weighted variation norm controlling QMC tractability (Cui et al., 2023, Kuo et al., 2015, Guth et al., 2021).

Component-by-component (CBC) construction is common, with dimension-independent error bounds achievable under product-and-order-dependent weights, provided the integrand has sufficient mixed regularity, as established for PDE-based UQ problems (Kuo et al., 2015, Guth et al., 2021).

3. Error and Complexity Analysis

The global mean-squared error (MSE) of $Y^{MLRQMC}$ decomposes into bias (discretization) and multilevel sampling variance:

$$\mathrm{MSE} = \underbrace{|\mathbb{E}[Q] - \mathbb{E}[Q_{L}]|^2}_{\text{discretization}} + \underbrace{\sum_{\ell=0}^L \mathbb{V}\mathrm{ar}\big( Q_{N_\ell}^{(\ell)}(Q_\ell - Q_{\ell-1}) \big)}_{\text{sampling error}},$$

where, for geometrically refined mesh levels or increasing inner sample sizes, the bias decays as $O(h_L^\alpha)$ or $O(m_L^{-\alpha_\text{inner}})$, while the sampling error on each level satisfies

$$\mathrm{Var}_\Delta[Q_{N_\ell}^{(\ell)}(Q_\ell-Q_{\ell-1})] \lesssim N_\ell^{-1/\eta}\, h_\ell^\beta,$$

with exponents $\alpha,\beta$ and QMC parameter $\eta>1/2$ determined by the problem's spatial and stochastic regularity (Cui et al., 2023, Bartuska et al., 2024).

For fixed random shifts, the total work to achieve tolerance $\varepsilon^2$ is minimized via a Lagrange-multiplier argument, equidistributing the contribution of each level's variance per unit cost. The resulting complexity is

$$\mathrm{Cost}(\varepsilon) \lesssim \begin{cases} \varepsilon^{-2\eta}, & \text{if } \beta\eta > \gamma \ \varepsilon^{-2\eta} (\log \varepsilon^{-1})^{\eta+1}, & \text{if } \beta\eta = \gamma \ \varepsilon^{-2\eta - (\gamma-\beta\eta)/\alpha}, & \text{if } \beta\eta < \gamma \end{cases}$$

where $\gamma$ quantifies the cost per sample versus mesh refinement (Cui et al., 2023, Kuo et al., 2015).

In standard MLMC (Monte Carlo), $\eta = 1/2$, yielding $\mathrm{Cost} = O(\varepsilon^{-2})$ in the optimal regime, but with QMC or HoQMC, $\eta > 1/2$, improved exponents (e.g., $\theta \approx 1$ observed empirically for lognormal diffusion) are achieved (Kuo et al., 2015, Guth et al., 2021).

4. Application Classes and Methodological Variants

Elliptic and Convection-Diffusion PDEs: Multilevel (R)QMC estimators are extensively adopted for forward UQ in elliptic and convection-diffusion eigenvalue problems, as well as for gradient estimation in PDE-constrained optimization. The structure supports strategies like SUPG stabilization and mesh-wise homotopy in convection-dominated regimes (Cui et al., 2023, Guth et al., 2021).

Bayesian Inverse Problems: Ratio-type and splitting-type multilevel HoQMC estimators enable efficient computation of posterior expectations and normalization constants in operator Bayesian inversion, with provable algebraic convergence rates independent of parameter dimension under sparsity of the random field expansion (Dick et al., 2016).

Nested Monte Carlo and Variational Inference: MLRQMC methodology applies to nested expectations, e.g., expected information gain in Bayesian design or unbiased variational Bayes, leveraging Rhee–Glynn telescoping on inner sample sizes coupled with (R)QMC for improved variance decay, achieving complexities improving upon both MLMC and single-level rQMC (He et al., 2021, Bartuska et al., 2024).

Financial Mathematics: Fourier–RQMC multilevel estimators are developed for risk quantification problems such as multivariate shortfall risk, exploiting frequency-domain smoothness for RQMC integration at each iteration of SQP optimization (Hammouda et al., 6 Feb 2026).

Continuous-Level Extensions: In continuous-level unbiased MLMC, RQMC can be used to sample the random level variable, yielding variance reduction and sharper (e.g., optimal $O(\varepsilon^{-2})$) complexity in certain regimes (Beschle et al., 2023).

5. Parameter Choices, Regularity Requirements, and Optimization

The effectiveness of multilevel (R)QMC depends on several factors:

• Regularity: QMC efficiency requires high-order mixed derivative bounds in the stochastic parameters. This is often satisfied by analytic dependence (e.g., lognormal fields) or through regularity results for the PDE solution map (Kuo et al., 2015, Guth et al., 2021).
• Level Definition: Levels can correspond to mesh refinement, KL truncation dimension, or accuracy/order of quadrature (nested quadrature sample sizes for nested expectations) (Bartuska et al., 2024, He et al., 2021).
• Sample Allocation: Optimal $N_\ell$ are obtained by pilot estimation of per-level variances and costs, with adaptive algorithms iteratively reallocating effort as in (Guth et al., 2021) and (Cui et al., 2023). Fine levels, which are costly but contribute less variance, receive fewer sampling points.
• Dimensionality: Employing weight structures in the QMC construction (POD weights), plus dimension-reordering heuristics, helps maintain tractability in high or infinite dimensions (Kuo et al., 2015, Guth et al., 2021).
• Random Shifts: A small fixed number ($R\approx10$–20) of random QMC shifts is typically sufficient for variance estimation and unbiasedness without significant cost overhead.

6. Numerical Performance Evidence and Practical Recommendations

Empirical studies consistently support the theoretical predictions:

• For lognormal diffusion, MLQMC achieves cost $\sim O(\varepsilon^{-1})$–$O(\varepsilon^{-1.3})$, improving substantially over MLMC and single-level QMC which exhibit exponents $2$ or higher (Kuo et al., 2015, Guth et al., 2021).
• In nested integration for EIG, ML-rQMC reduces cost by up to $60\%$ compared to MLMC and achieves steep MSE-vs.-cost decay (Bartuska et al., 2024).
• In variational Bayes gradient estimation, RQMC boosts the variance-decay exponent, allows more aggressive weight decay, and reduces total cost factors by $\sim 16\%$ over standard MLMC in logistic mixed models (He et al., 2021).
• In PDE-constrained optimization and financial risk quantification, multilevel RQMC matches or surpasses state-of-the-art MLMC and stochastic optimization benchmarks (Guth et al., 2021, Hammouda et al., 6 Feb 2026).

Best practices include the use of nested FE grids, high-order QMC rules tailored via CBC to the effective stochastic dimension, adaptive per-level sample allocation, and stabilization approaches (e.g., SUPG, homotopy) when required by convection or other problem features. For ratio-type estimators in Bayesian applications, splitting formulations afford additional numerical stability in small noise settings (Dick et al., 2016).

7. Limitations, Open Directions, and Regime-Optimality

The benefits of multilevel RQMC are maximized when:

• The integrand is sufficiently smooth and the underlying random field expansion has fast eigenvalue decay.
• The problem structure allows for telescoping decompositions with rapidly decaying correction variances.
• Implementation of QMC rules is tractable at high dimensions, often via POD or tensor-product weight designs.

Limitations arise if regularity is poor, such as in problems with discontinuous integrands or extremely high-dimensional, non-sparse random expansions. In such settings, gains over MLMC may be moderated and classical rates may govern overall complexity (Kuo et al., 2015, Cui et al., 2023, He et al., 2021). In continuous-level MLMC, quasi–Monte Carlo sampling of the level variable never worsens, and often restores, optimal exponents versus the classical regime (Beschle et al., 2023).

The regime where $\min\{\beta,2\}>\alpha$ and $\gamma\le\alpha$ delivers the canonical $O(\varepsilon^{-2})$ scaling for both CLMC and QCLMC; in other situations, the QMC component improves the constant factors or achieves the exponent, but cannot outperform this benchmark (Beschle et al., 2023). Further research is ongoing in extending MLRQMC to more general stochastic processes, multi-index or continuous-index frameworks, and optimizing for modern, heterogeneous computational architectures.

---

References

• (Cui et al., 2023): Multilevel Monte Carlo methods for stochastic convection-diffusion eigenvalue problems
• (Kuo et al., 2015): Multilevel Quasi-Monte Carlo Methods for Lognormal Diffusion Problems
• (Bartuska et al., 2024): Multilevel randomized quasi-Monte Carlo estimator for nested integration
• (He et al., 2021): Unbiased MLMC-based variational Bayes for likelihood-free inference
• (Dick et al., 2016): Multilevel higher order Quasi-Monte Carlo Bayesian Estimation
• (Guth et al., 2021): Multilevel Quasi-Monte Carlo for Optimization under Uncertainty
• (Beschle et al., 2023): Complexity analysis of quasi continuous level Monte Carlo
• (Hammouda et al., 6 Feb 2026): Single- and Multi-Level Fourier-RQMC Methods for Multivariate Shortfall Risk