Variance reduction with probing and Multilevel Monte Carlo in Lattice QCD

Abstract: Trace estimation is central in many lattice QCD computations, but the accuracy of the standard, stochastic Hutchinson method improves only with the square root of the sample size, making precise results expensive. We investigate two complementary variance reduction strategies. First, multigrid multilevel Monte Carlo uses a multigrid hierarchy to construct an unbiased multilevel estimator via recursive coarse grid corrections available from the multigrid hierarchy of the solver. Second, stochastic probing uses distance-$d$ graph colorings; we propose a torus based coloring that requires substantially fewer colors than hierarchical probing at the same distance. We test these approaches on two representative problems: the connected pseudoscalar correlator and disconnected fermion loops. For the connected pseudoscalar two-point function, the multilevel decomposition yields a variance reduction of up to $\mathcal{O}(10^5)$ at large time separations and translates into a clear cost reduction at fixed accuracy, thus confirming earlier results of arXiv:2412.06347. For the disconnected loops, in contrast, the multilevel decomposition provides only moderate gains, whereas probing combined with dilution delivers a substantial cost reduction that improves as the number of probing vectors is increased. Overall, the results highlight a pronounced complementarity: deflation schemes are most effective for observables dominated by long distance propagation, while probing is most effective for localized quantities.

• The paper introduces MLMC, achieving up to 10^5 variance reduction for long-distance pseudoscalar correlators.
• The paper employs a new toroidal graph coloring strategy in stochastic probing to cut the number of required solves by up to a factor of four.
• The paper demonstrates how MLMC excels for global observables while probing targets localized fluctuations, suggesting a hybrid approach.

Variance Reduction Techniques for Trace Estimation in Lattice QCD

Introduction and Motivations

Estimating traces of matrix functions, particularly $tr[D^{-1}]$ for large sparse matrices such as the Wilson-Dirac operator in Lattice QCD, is a recurring computational bottleneck. Traditional stochastic estimators (e.g., Hutchinson’s method) suffer from slow convergence, requiring $O(1/\epsilon^2)$ solves to achieve a target relative precision $\epsilon$, leading to prohibitive compute costs. Frommer et al. systematically investigate two complementary variance reduction strategies—Multigrid Multilevel Monte Carlo (MLMC) and stochastic probing (including a new toroidal graph coloring)—for two classes of Lattice QCD observables: connected pseudoscalar correlators and disconnected fermion loops (2607.05157).

Multigrid Multilevel Monte Carlo: Hierarchical Deflation and Variance Localization

The MLMC estimator leverages the multigrid hierarchy of solvers to recursively decompose $D^{-1}$ via sequential coarse-grid corrections. The prolongators are constructed from near-null modes and encode long-range correlations, efficiently deflating the variance associated with low-lying eigenmodes. The resulting unbiased estimator,

$$\widehat{tr}(D^{-1}) = \sum_{l=1}^{L-1} \widehat{tr}(M_l) + \widehat{tr}(D_L^{-1}),$$

transfers residual variance to coarse levels where linear solves are significantly cheaper, while fine-level variance is dramatically suppressed.

Numerical results for the connected pseudoscalar correlator $G(t)$ show that MLMC achieves a variance reduction up to $\mathcal{O}(10^5)$ at large time separations. The decomposition concentrates variance onto coarse terms, enabling targeted sampling and cost savings commensurate with the expected overhead of multigrid setup.

Figure 1: MLMC delivers drastic variance suppression for the connected correlator, relegating most variance to computationally inexpensive coarse-level terms.

These results reinforce earlier findings for multigrid low-mode averaging and demonstrate that for observables dominated by propagation over large time distances (i.e., global modes), hierarchical deflation techniques are optimal and cost-effective.

Stochastic Probing: Exploiting Matrix Sparsity and Locality

Variance in trace estimation is often dominated by off-diagonal contributions. Probing methods attack this by constructing structurally orthogonal probing vectors via graph colorings, specifically exploiting the adjacency structure of the Wilson-Dirac operator. The paper introduces a toroidal coloring scheme, which significantly reduces the number of colors required relative to hierarchical probing (HP), especially for higher coloring distances. These vectors maintain stochastic characteristics (Rademacher structure) and yield improved trace estimators after combination with dilution schemes, which decouple internal degrees of freedom (spin/color components), further increasing estimator independence.

For disconnected fermion loop observables, probing exhibits substantial variance reduction and cost savings compared to MLMC. The new torus-based coloring achieves up to a factor of four fewer colors than HP for similar distances, leading to fewer solves and better scalability.

Figure 2: Estimated variances for disconnected loops; MLMC achieves only moderate reduction, failing to offset setup overhead for local observables.

Figure 3: Probing with full dilution and enhanced coloring reduces the number of solves required to reach target variance, with improvement scaling with the number of probing vectors.

This highlights that probing is maximally effective for quantifying local or short-range fluctuations—where the locality of the observable matches the locality encoded in the coloring and dilution.

Comparative Results and Complementary Regimes

The complementary efficacy of MLMC and probing is quantitatively established: MLMC yields up to $10^{5}$ variance reduction for long-distance correlators, while probing combined with dilution achieves up to $10^{4}$ reduction for localized observables. The latter is particularly relevant for disconnected diagrams, where MLMC's hierarchical structure does not adequately capture local variance sources. Both techniques are robust and can deliver reduction factors of one order of magnitude in computational cost for matched observables.

Implications and Outlook

These findings have immediate practical implications for Lattice QCD computations: deploying MLMC for long-range correlators and probing schemes for local quantities can systematically minimize solver overhead. The new torus-based coloring is generalizable and can be optimized offline for specific lattice sizes. Theoretically, the observed complementarity suggests hybrid schemes as a promising avenue—integrating multigrid MLMC with probing or dilution to simultaneously target global and local sources of stochastic variance.

Future directions include formal benchmarking of three-dimensional torus coloring against hierarchical probing at fixed computational cost, further optimizing transfer operator construction in MLMC (e.g., via additional exact modes), and exploring combined variance reduction strategies that seamlessly merge strengths of both methods.

Conclusion

This work establishes that variance reduction in stochastic trace estimation in Lattice QCD is highly sensitive to the structure of the observable. Multigrid MLMC is optimal for observables dominated by global, low-mode contributions, while stochastic probing (with advanced coloring and dilution) excels for localized fluctuations. Both deliver strong computational savings in their respective regimes, and their complementarity motivates integration for comprehensive variance control. The technical innovations in coloring and hierarchical sampling are readily applicable to broader classes of sparse matrix problems.