Variance reduction with probing and Multilevel Monte Carlo in Lattice QCD

Andreas Frommer; Bruno Lang; Gustavo Ramirez-Hidalgo; Jose Jimenez-Merchan; Mario Papace

REVIEW 2 major objections 4 minor 16 references

Variance reduction with probing and Multilevel Monte Carlo in Lattice QCD

T0 review · 2 major / 4 minor · reviewed 2026-07-11 · grok-4.5

Challenge this review Re-run · record.json Download PDF Read on arXiv ↗ Review JSON

Pith's one-line read Multigrid multilevel Monte Carlo and torus-based probing cut variance in lattice QCD traces by complementary mechanisms: one for long-distance correlators, one for local loops.

desk verdict Solid numerical methods paper that cleanly shows multigrid MLMC and torus probing are complementary on two standard QCD observables; the O(10^5) variance claim is real on the given ensemble, but net-cost accounting is incomplete. read the letter →

arxiv 2607.05157 v1 pith:NJTCWQSE submitted 2026-07-06 hep-lat cs.NAmath.NA

Andreas Frommer , Jose Jimenez-Merchan , Bruno Lang , Mario Papace , Gustavo Ramirez-Hidalgo This is my paper · ORCID

classification hep-latcs.NAmath.NA

keywords latticeQCDtraceestimationmultilevelMonteCarlomultigridstochasticprobingvariancereductiondisconnectedloopspseudoscalarcorrelator

verification ladder T0 review T1 audit T2 compute T3 formal

The pith

A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.

The reading

Computing traces involving the inverse of the Wilson-Dirac operator is a bottleneck in lattice QCD, because the plain stochastic Hutchinson estimator only improves as the square root of the number of samples. This paper shows that two existing variance-reduction ideas become highly effective when adapted to the structure of that operator, and that each idea matches a different class of observables. Multigrid multilevel Monte Carlo recursively peels off coarse-grid corrections already available from the multigrid solver; for the connected pseudoscalar two-point function the fine-level variance collapses by as much as five orders of magnitude at large time separations, producing a clear net cost saving. Stochastic probing with a new torus coloring, combined with dilution, instead annihilates short-range off-diagonal contributions; for disconnected fermion loops it yields a substantial reduction in the number of linear solves that improves steadily with the number of probing vectors. The practical message is that deflation-style methods should be used when the observable is dominated by long-distance modes, while probing should be used when the dominant fluctuations are local.

What carries the argument

Recursive multigrid splitting of the inverse via the oblique projectors already present in the multigrid hierarchy (eqs. 8–11), together with a distance-d torus coloring (eq. 15) that produces far fewer colors than hierarchical probing at the same distance.

What would settle it

Repeat the connected-correlator measurement on a second ensemble with different volume or quark mass using the same multigrid hierarchy; if the fine-level variance reduction falls well below O(10^3–10^5), or if the torus coloring requires substantially more colors than hierarchical probing at equal distance, the claimed complementarity fails.

Watch

Extended reading notes

Core claim

The multigrid multilevel decomposition of the inverse yields a variance reduction of up to O(10^5) for the connected pseudoscalar correlator at large time separations and a clear cost reduction at fixed accuracy, while for disconnected loops the same multilevel scheme gives only moderate gains and probing plus dilution is the method that produces a substantial, scalable cost saving. The two techniques are therefore complementary rather than interchangeable.

Load-bearing premise

That transfer operators built from only 28 test vectors on one 64 imes32^{3} ensemble already capture enough low-mode content for the fine-level variance collapse to hold more generally, and that the torus-coloring coefficients found by exhaustive search stay near-optimal on other lattices.

Editorial extensions

If this is right

Connected long-distance correlators can be estimated at fixed accuracy with far fewer expensive fine-grid solves by shifting residual variance onto cheap coarse levels.
Disconnected local loops become cheaper by systematically increasing the number of torus-probing vectors plus dilution rather than by multilevel deflation.
The two techniques can be combined so that long-range and short-range variance contributions are attacked simultaneously.
Any observable whose variance is known to be dominated by either low modes or by local off-diagonal entries can now be assigned the matching reduction strategy a priori.

Reading between the lines

Editorial extensions of the paper, not claims the author makes directly.

The same multigrid splitting should apply, with only minor changes, to other connected Wick contractions that involve products of two or more inverses.
Once the torus-coloring coefficients are tabulated for standard lattice sizes they become a free, reusable library resource for any code that already uses probing.
If the moderate multilevel gain observed for disconnected loops can be amplified by adding a few exact low modes, the two methods may become synergistic rather than merely complementary.

Share X Bluesky LinkedIn Reddit HN

Editorial analysis

A structured set of objections, weighed in public.

Desk editor's note, referee report, and a circularity audit.

Desk Editor's Note

read the letter

The one thing worth knowing is that this paper gives a clean, side-by-side demonstration that multigrid multilevel Monte Carlo and a new torus coloring attack different parts of the variance in lattice-QCD traces. For the connected pseudoscalar correlator the multilevel split (eqs. 9–12, 17–20) collapses the expensive fine-level variance by up to O(10^5) at large t (Fig. 1); for disconnected loops the same multilevel only helps moderately, while probing + full dilution cuts the solve count substantially (Fig. 3). That complementarity is the real takeaway and matches the physics of long-range versus local contributions.

What is new is modest but useful. Multigrid MLMC for traces already existed (Frommer et al. 2022; Gruber et al. arXiv:2412.06347); the paper confirms and extends those numbers on the RQCD IV ensemble. The genuine increment is the torus coloring (eq. 15 and Table 1) that needs roughly four times fewer colors than hierarchical probing at the same distance. The algebra is clean: the recursive oblique projectors give an unbiased estimator with no free parameters, and the variance plots are consistent with the claimed reductions under fixed N=500.

The soft spot is exactly the one the stress-test flags, and it is real but not fatal. Figure 1 shows per-component variances, not the optimal sample allocation N_l from the MLMC formula (eq. 6), not the relative costs C_l of coarse versus fine solves, and not the setup cost of the 28 test vectors. Without those numbers the abstract’s “clear cost reduction at fixed accuracy” is only partially supported; the disconnected case already shows that moderate fine-level gains can be wiped out by overhead. Single-ensemble, fixed-N variance estimates also leave generalization open. These are ordinary limitations of a short conference paper, not hidden flaws.

Anyone running production correlators or disconnected loops will get value from the figures and the coloring formula. The math and citations are solid; the work is honest about where each method wins. I would send it to peer review without hesitation—referees can demand the missing cost tables—and I would cite the complementarity result and the torus coloring myself.

Referee Report

2 major / 4 minor

Summary. The manuscript studies two complementary variance-reduction techniques for stochastic trace estimation in lattice QCD: (i) multigrid multilevel Monte Carlo (MLMC), which builds an unbiased multilevel estimator from the recursive oblique-projector splitting induced by a Galerkin multigrid hierarchy (eqs. 8–12), and (ii) stochastic probing based on a new torus coloring (eq. 15) that uses substantially fewer colors than hierarchical probing at the same distance (Table 1). On a single RQCD IV configuration, the multilevel decomposition of the connected pseudoscalar two-point function yields fine-level variance reductions up to O(10^5) at large time separations (Fig. 1), while for disconnected fermion loops multilevel gains remain moderate (Fig. 2) and probing combined with full dilution produces a clear reduction in total solves versus target variance (Fig. 3). The authors conclude that deflation-type schemes are most effective for long-distance observables and probing for localized ones.

Significance. If the reported complementarity holds more generally, the paper supplies a practical guide for matching variance-reduction tools to the structure of lattice QCD observables, which is of direct use for precision disconnected and connected correlator campaigns. Strengths include an algebraically unbiased multilevel identity that reuses the existing multigrid solver hierarchy, an explicit and tableable torus coloring with fewer colors than hierarchical probing, and an honest demonstration that multilevel MLMC does not automatically improve cost for disconnected loops. The work also confirms the large fine-level variance collapse for connected correlators reported in arXiv:2412.06347. These are concrete, implementable contributions rather than purely formal ones.

major comments (2)

[Abstract, §4.1, Fig. 1, §5] The abstract and §4.1 assert that the multilevel decomposition 'translates into a clear cost reduction at fixed accuracy' for the connected pseudoscalar correlator, and §5 quantifies 'cost reductions of about one order of magnitude.' Figure 1 and the surrounding text report only per-component variances of the G_{i,j} terms under a fixed sample size N=500. Unlike the probing study (Fig. 3), the paper does not report the MLMC sample allocation N_l from eq. (6), the relative costs C_l of fine versus coarse solves, the multigrid setup cost for N_tv=28 test vectors, or a total-solve (or wall-clock) curve versus target variance. Because residual variance is concentrated in the coarsest term G_{3,3}, a net saving is plausible, but it is not demonstrated in this manuscript; the disconnected-loop discussion already shows that moderate fine-level variance reduction need not yield net cost improvem
[§4 setup, Table 2, §§4.1–4.2, §5] All numerical results (variance plots, cost curves, and the complementarity conclusion) are obtained on a single IV configuration (Table 2) with a fixed multigrid setup (three levels, N_tv=28). The weakest assumption underlying the fine-level variance collapse is that these transfer operators already capture enough low-mode content for the observed reductions to be representative. Without at least a second ensemble, a different volume, or a sensitivity check in N_tv, the strength of the claim that 'deflation schemes are most effective for observables dominated by long distance propagation, while probing is most effective for localized quantities' remains limited. A short additional data set or an explicit caveat that the complementarity is demonstrated on one ensemble would make the conclusion proportionate to the evidence.

minor comments (4)

[§3, Table 1, §5] The torus coloring (eq. 15) is a useful practical contribution, but the paper only compares color counts to hierarchical probing (Table 1), not variance reduction or cost at fixed accuracy. The outlook correctly flags a controlled benchmark; a sentence in §3 or §4.2 stating that such a comparison is left for future work would prevent readers from over-reading Table 1 as a performance comparison.
[§4, Fig. 1] Variance estimates themselves are obtained with fixed N=500. A brief remark on the statistical uncertainty of the reported variances (especially the O(10^5) ratio at large t) would strengthen confidence in Fig. 1.
[§3, §4.2, Figs. 1–3] Typographical issues: 'combinantion' should be 'combination' (end of §3); several run-on words appear in the compiled text (e.g. near the start of §4.2). Axis labels and legends in Figs. 1–3 should be checked for readability in the final PoS layout.
[§2, §4.1] The relation between the multilevel splitting used here and the Wick-contraction decomposition of Ref. [1] could be stated more explicitly in §2 or §4.1 so that the precise novelty relative to that work is clear.

Circularity Check

0 steps flagged · score 0.0 of 10

No significant circularity: multilevel identities are algebraic consequences of the Galerkin projector, probing vectors follow an explicit coloring formula, and reported variance reductions are independent numerical measurements.

full rationale

The paper's central claims are empirical variance and cost reductions measured on two lattice-QCD observables. The multigrid multilevel estimator rests on the algebraic identity tr(D^{-1}) = sum tr(M_l) + tr(D_L^{-1}) obtained by recursive insertion of the oblique projector Π_l = P_l D_{l+1}^{-1} P_l^† D_l (eqs. 8–11); this identity holds by construction of the Galerkin coarse operator and does not involve fitted parameters or data-dependent predictions. Stochastic probing vectors are defined directly from a distance-d torus coloring (eq. 15) whose coefficients are precomputed once by exhaustive search; the subsequent Hutchinson averages are ordinary Monte-Carlo estimates. Self-citations supply the multigrid hierarchy ([5]) and the general MLMC sample-allocation formula, but the numerical O(10^5) variance collapse for the connected correlator (Fig. 1) and the cost curves for probing+dilution (Fig. 3) are fresh measurements on a concrete ensemble; they are not forced by those citations. No uniqueness theorem, ansatz smuggling, or renaming of a known pattern appears. The derivation chain is therefore self-contained against external benchmarks.

Assumptions & free parameters 3 free parameters · 3 assumptions · 0 invented entities

The central claims rest on standard Monte-Carlo variance identities, the algebraic multigrid projector splitting already published by the authors, and the empirical observation that low modes dominate long-distance correlators while locality dominates disconnected loops. Free parameters are the usual multigrid setup choices (number of test vectors, number of levels) and the coloring coefficients obtained by search; no new physical constants or ad-hoc entities are introduced.

free parameters (3)

N_tv (number of multigrid test vectors) = 28
Chosen as 28 for the three-level hierarchy; controls how well low modes are captured and therefore the fine-level variance reduction.
torus coloring coefficients σ_j = lattice- and d-dependent (Table 1 example)
Determined by exhaustive search for given lattice dimensions and distance d; they fix the number of colors n_c and therefore the probing cost.
sample size N for variance estimation = 500
Fixed at 500 for all reported variance estimates; affects the reliability of the claimed reduction factors.

assumptions (3)

standard math Hutchinson estimator variance equals (1/2)||offdiag(D^{-1}+(D^{-1})^T)||_F^2 for Rademacher vectors
Standard result used to motivate both deflation and probing (eq. 4).
domain assumption Aggregation-based multigrid prolongators capture global near-null modes via local coherence
Invoked to justify that the recursive coarse-grid corrections shift variance to cheap levels (§2, citing Lüscher local coherence).
standard math Distance-d coloring annihilates contributions from matrix entries within graph distance d
Classical probing theory used to construct the stochastic probing vectors (eq. 13).

how reviews work

0 comments

Cite this review

Pith. "Pith review of Variance reduction with probing and Multilevel Monte Carlo in Lattice QCD." pith.science (2026). https://pith.science/paper/NJTCWQSE

@misc{pith2026260705157,
  author       = {Pith},
  title        = {Pith review of: Variance reduction with probing and Multilevel Monte Carlo in Lattice QCD},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/NJTCWQSE}},
  note         = {Machine review of arXiv:2607.05157}
}

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.

Figures

Figures reproduced from arXiv: 2607.05157 by the authors.

**Figure 1.** Comparison between the plain Hutchinson estimator and the multigrid Multilevel Monte Carlo estimator based on oblique projections. where 𝐺𝑖, 𝑗(𝑡) = 1 𝑇 ∑︁ 𝑇 𝑡 ′=1 tr 𝐵𝑖(𝑡 + 𝑡 ′ , 𝑡′ )Γ5 𝐵𝑗(𝑡 ′ , 𝑡 + 𝑡 ′ )Γ5 [PITH_FULL_IMAGE:figures/full_fig_p007_1.png] view at source ↗

**Figure 2.** Estimated variances for tr(Γ5 (𝑡)𝐷 −1 (𝑡, 𝑡)) on the IV configuration at 𝑡 = 1, comparing plain Hutchinson to the multilevel terms. 10-1 100 101 2 102 103 104 105 106 Ntotal [PITH_FULL_IMAGE:figures/full_fig_p008_2.png] view at source ↗

**Figure 3.** Total number of linear system solves (Ntotal) required to reach target variance 𝜀 2 for the disconnected correlator. Coloring distances (CD) from 1 to 3 are explored in combination with full dilution (Dil12) where all the 12 spin and color degrees of freedom are decoupled. where 𝐵𝑙(𝑡, 𝑡) is defined analogously to eqs. (17) and (18), now restricted to the (𝑡, 𝑡) time-slice block. The measurements were done for 𝑡 = 1.… view at source ↗

Discussion (0). Continue with ORCID to comment.

Reference graph

Works this paper leans on

16 extracted references · 4 linked inside Pith

[1]

Gruber, T

R. Gruber, T. Harris and M.K. Marinković,Multigrid low-mode averaging,Phys. Rev. D111 (2025) 7 074508

2025
[2]

Hutchinson,A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines,Commun

M.F. Hutchinson,A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines,Commun. Stat. Simul. Comput.19(1990) 2 433. 9 Variance reduction with probing and Multilevel Monte Carlo in Lattice QCDJose Jimenez-Merchan

1990
[3]

Bernardson, P

S. Bernardson, P. McCarty and C. Thron,Monte Carlo methods for estimating linear combinations of inverse matrix entries in lattice QCD,Comput. Phys. Commun.78(1994) 3 256

1994
[4]

Gambhir, A

A.S. Gambhir, A. Stathopoulos and K. Orginos,Deflation as a method of variance reduction for estimating the trace of a matrix inverse,SIAM J. Sci. Comput.39(2017) A532 [1603.05988]

arXiv 2017
[5]

Frommer, M.N

A. Frommer, M.N. Khalil and G. Ramirez-Hidalgo,A multilevel approach to variance reduction in the stochastic estimation of the trace of a matrix,SIAM J. Sci. Comput.44 (2022) A2536 [2108.11281]

arXiv 2022
[6]

Giles,Multilevel Monte Carlo methods,Acta Numerica24(2015) 259

M.B. Giles,Multilevel Monte Carlo methods,Acta Numerica24(2015) 259

2015
[7]

Giusti, T

L. Giusti, T. Harris, A. Nada and S. Schaefer,Frequency-splitting estimators of single-propagator traces,Eur. Phys. J. C79(2019) 7 586

2019
[8]

Whyte, A

T. Whyte, A. Stathopoulos, E. Romero and K. Orginos,Optimizing shift selection in multilevel Monte Carlo for disconnected diagrams in lattice QCD,Computer Physics Comm. 294(2024) 108928

2024

Show all 16 references

[9]

G.S. Bali, S. Collins and A. Schäfer,Effective noise reduction techniques for disconnected loops in lattice QCD,Comput. Phys. Commun.181(2010) 9 1570

2010
[10]

Frommer, K

A. Frommer, K. Kahl, S. Krieg, B. Leder and M. Rottmann,Adaptive aggregation-based domaindecompositionmultigridforthelatticeWilson–Diracoperator,SIAMJ.Sci.Comput. 36(2014) 4 A1581

2014
[11]

Lüscher,Local coherence and deflation of the low quark modes in lattice QCD,J

M. Lüscher,Local coherence and deflation of the low quark modes in lattice QCD,J. High Energy Phys.2007(2007) 07 081

2007
[12]

Wilcox,Noise methods for flavor singlet quantities, inInterdisciplinary Workshop on Numerical Challenges in Lattice QCD, pp

W. Wilcox,Noise methods for flavor singlet quantities, inInterdisciplinary Workshop on Numerical Challenges in Lattice QCD, pp. 127–141, 8, 1999 [hep-lat/9911013]

1999 arXiv
[13]

Frommer, M

A. Frommer, M. Rinelli and M. Schweitzer,Analysis of stochastic probing methods for estimating the trace of functions of sparse symmetric matrices,Math. Comp.94(2025) 801

2025
[14]

Stathopoulos, J

A. Stathopoulos, J. Laeuchli and K. Orginos,Hierarchical Probing for Estimating the Trace of the Matrix Inverse on Toroidal Lattices,SIAM J. Sci. Comput.35(2013) 5 S299 [1302.4018]

2013 arXiv
[15]

Foley, K.J

J. Foley, K.J. Juge, A.O. Cais, M. Peardon, S.M. Ryan and J.I. Skullerud,Practical all-to-all propagators for lattice QCD,Comput. Phys. Commun.172(2005) 3 145

2005
[16]

G.S. Bali, S. Collins, B. Gläßle, M. Göckeler, J. Najjar, R.H. Rödl et al.,The moment⟨𝑥⟩𝑢−𝑑 of the nucleon from𝑁𝑓 =2lattice QCD down to nearly physical quark masses,Phys. Rev. D 90(2014) 7 074510. 10

2014

Pith tools

Reviewed July 11, 2026 · model on record in the stance chip above.

[1] [1]

Gruber, T

R. Gruber, T. Harris and M.K. Marinković,Multigrid low-mode averaging,Phys. Rev. D111 (2025) 7 074508

2025

[2] [2]

Hutchinson,A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines,Commun

M.F. Hutchinson,A stochastic estimator of the trace of the influence matrix for Laplacian smoothing splines,Commun. Stat. Simul. Comput.19(1990) 2 433. 9 Variance reduction with probing and Multilevel Monte Carlo in Lattice QCDJose Jimenez-Merchan

1990

[3] [3]

Bernardson, P

S. Bernardson, P. McCarty and C. Thron,Monte Carlo methods for estimating linear combinations of inverse matrix entries in lattice QCD,Comput. Phys. Commun.78(1994) 3 256

1994

[4] [4]

Gambhir, A

A.S. Gambhir, A. Stathopoulos and K. Orginos,Deflation as a method of variance reduction for estimating the trace of a matrix inverse,SIAM J. Sci. Comput.39(2017) A532 [1603.05988]

arXiv 2017

[5] [5]

Frommer, M.N

A. Frommer, M.N. Khalil and G. Ramirez-Hidalgo,A multilevel approach to variance reduction in the stochastic estimation of the trace of a matrix,SIAM J. Sci. Comput.44 (2022) A2536 [2108.11281]

arXiv 2022

[6] [6]

Giles,Multilevel Monte Carlo methods,Acta Numerica24(2015) 259

M.B. Giles,Multilevel Monte Carlo methods,Acta Numerica24(2015) 259

2015

[7] [7]

Giusti, T

L. Giusti, T. Harris, A. Nada and S. Schaefer,Frequency-splitting estimators of single-propagator traces,Eur. Phys. J. C79(2019) 7 586

2019

[8] [8]

Whyte, A

T. Whyte, A. Stathopoulos, E. Romero and K. Orginos,Optimizing shift selection in multilevel Monte Carlo for disconnected diagrams in lattice QCD,Computer Physics Comm. 294(2024) 108928

2024