Conditioning versus preconditioning for gradient-based Dirac solvers in lattice QCD: net time-to-solution advantage at physical quark masses, production volumes, and many right-hand sides

Nearly every lattice-QCD computation is bounded by the solution of huge, sparse, extremely ill-conditioned Dirac linear systems, typically by conjugate-gradient-type (gradient/Krylov) iteration on the Hermitian positive-definite normal system. At (near-)physical quark masses and production four-dimensional volumes the condition number of the Dirac normal matrix reaches $10^6$-$10^8$, plain iteration stalls, and the recorded open question is the conditioning-versus-preconditioning tradeoff: whether deflation, multigrid, and gradient-trained (neural-network) preconditioners deliver a genuine net time-to-solution advantage - counting all setup, eigensolve, and training costs amortized over the many right-hand sides and quark masses of a production calculation. The feasibility of the deflation and multigrid branches is now established in production for Wilson-type and domain-wall fermions (multigrid is the acknowledged state of the art, with demonstrated order-of-magnitude wall-clock speedups at physical masses), but the cost-benefit residue recorded by the CLS source - whether a given setup pays off once all setup costs are amortized, and how solver parameters should be selected at large volumes - remains part of the open tradeoff. The surviving open core is the learned-preconditioner branch - no 4D production-scale demonstration of net advantage exists through early 2026 - together with this setup-amortization and parameter-selection residue, a principled quantitative understanding of the conditioning-preconditioning tradeoff (parameter scaling with volume, mass, and right-hand-side structure; robustness to degraded low-mode estimates), and the source-discussed conditioning bottleneck inside stochastic-gradient training of fermionic flow-based samplers.

Background

Lattice QCD is the first-principles discretization of quantum chromodynamics: space-time is replaced by a four-dimensional hypercubic grid, gluon fields by group-valued link variables, and quark fields by lattice spinors. Almost every lattice-QCD calculation reduces, at its computational core, to solving the lattice Dirac equation $D x = b$ many thousands of times - for each quark flavor, source point, and Monte Carlo gauge configuration - where $D$ is the discretized Dirac operator (a Wilson or Wilson-clover discretization - the standard fermion discretizations of production lattice QCD, the clover variant adding an improvement term - represented by a sparse matrix of dimension up to $\sim 10^{10}$ on production lattices). In practice one applies conjugate gradient (CG), a gradient-type Krylov-subspace iterative method (Krylov methods build their search space from repeated applications of the matrix to the starting residual), to the Hermitian positive-definite (HPD) normal system $D^\dagger D x = b$, or a non-Hermitian Krylov method such as the biconjugate gradient stabilized method (BiCGStab) to $D$ directly; the iteration count grows roughly like the square root of the condition number $\kappa$, the ratio of largest to smallest eigenvalue.

At physical quark masses and large volumes, $\kappa(D^\dagger D)$ reaches $10^6$-$10^8$: the smallest eigenvalues are set by the light-quark (pion) mass and become dense as the volume grows, so plain CG stalls - an instance of critical slowing down, the divergence of computational cost exactly at the physically interesting point. Two classical remedies dominate production codes. Deflation computes a set of low-lying eigenvectors of $D^\dagger D$ and projects them out of the Krylov space, removing precisely the modes that cause the stalling; Lüscher's local-coherence principle (the low quark modes vary slowly over blocks of neighboring sites) keeps the required number of modes growing only mildly with volume, making deflation scalable. Multigrid instead builds a hierarchy of coarsened lattices whose operators represent those same smooth low modes, correcting the error components that fine-grid iteration leaves behind; adaptive aggregation-based multigrid, which constructs coarse-grid interpolation from near-null-space test vectors (approximate null modes of the operator), is effective for Wilson-type fermions and is implemented in production solver libraries such as the GPU library QUDA. Both techniques cost setup time (eigensolves or hierarchy construction) that must be amortized over many solves, and both carry tunable parameters - deflation-subspace size, aggregate size, hierarchy depth - whose optimal scaling with volume, quark mass, and the right-hand-side structure has no closed-form theory and remains empirically tuned.

A third family is recent: neural-network preconditioners, linear maps trained by (stochastic) gradient descent - typically gauge-equivariant convolutional networks, since lattice physics is invariant under local color rotations - inserted into the CG iteration to approximate an effective inverse of $D^\dagger D$ or to supply restriction/prolongation operators inside a multigrid-like hierarchy. Their demonstrated successes are confined to the two-dimensional Schwinger model (two-dimensional quantum electrodynamics, a standard inexpensive testbed for lattice methods), to matrix/operator-learning variants, and to very small four-dimensional volumes; scaling to production 4D lattices with a genuine net wall-clock advantage is the explicitly stated open target of that literature. The same conditioning bottleneck also appears on the 'gradient side' of machine-learning lattice methods: flow-based samplers (neural normalizing-flow generative models trained to produce lattice field configurations) for theories with dynamical fermions represent the fermion determinant by pseudofermions - bosonic auxiliary fields whose quadratic action involves $D^\dagger D$ - so each stochastic-gradient training step requires ill-conditioned Dirac inversions, and standard preconditioners can conflict with the equivariance the flow architectures require.

The honest success metric throughout is net time-to-solution: total wall-clock time, including all setup, eigensolve, or training cost amortized over the full campaign of right-hand sides and quark masses, compared against the production-standard baseline - not iteration counts or operator-application counts alone, and not speedups measured with setup costs excluded.

Problem Statement

For the ill-conditioned Dirac systems of four-dimensional lattice QCD at (near-)physical quark masses, production-size volumes, and the many right-hand sides and quark masses of a production calculation, determine how the conditioning-versus-preconditioning tradeoff for gradient/Krylov Dirac solvers is optimally managed, and in particular whether gradient-trained (learned, typically gauge-equivariant neural-network) preconditioners - alone or combined with deflation and multigrid - deliver a genuine net time-to-solution advantage over the production-standard multigrid and deflation solvers, where 'genuine net advantage' counts all training, setup (eigensolves, hierarchy construction), and application costs, amortized over the full campaign of solves. Open aspects of the same recorded conditioning-versus-preconditioning tradeoff comprise: (a) a principled, quantitative characterization of how the tunable quantities of each preconditioning family - deflation-subspace size, multigrid hierarchy and aggregate parameters, learned-preconditioner capacity and training set - should scale with lattice volume, quark mass, and the right-hand-side/mass structure, replacing the current per-ensemble empirical tuning; (b) robustness of the resulting preconditioner when low-mode (eigenvector) estimates degrade, and when transferring across gauge configurations, topological sectors, volumes, and theories; and (c) the analogous conditioning bottleneck on the gradient side, namely whether the repeated $D^\dagger D$ inversions inside stochastic-gradient training of fermionic flow-based samplers (normalizing-flow generative models for lattice field configurations) can be preconditioned at QCD scale without breaking the equivariance the samplers require. Scope note: the 2007-2022 sources recorded this question jointly for deflation, multigrid, and neural-network preconditioners; production multigrid and deflation have since answered the feasibility/efficiency question for Wilson-type and domain-wall fermions and are the demonstrated state of the art at physical masses and large volumes (see References 5, 6, 8, 11, 16), while the cost-benefit residue recorded in the CLS source - whether a given setup pays off once all setup costs are amortized, and how solver parameters should be selected at large volumes - remains part of the open tradeoff. The open target therefore comprises the learned-preconditioner branch, this setup-amortization and parameter-selection residue for the classical families, and the regimes - master-field-scale volumes, sampler training - where no production-scale resolution yet exists. Completion boundary: the primary question is settled either way - by an affirmative demonstration of net time-to-solution advantage for gradient-trained preconditioners, or by a demonstrated absence of one under fair, fully-costed accounting at the stated production regime (near-physical quark masses, production-size four-dimensional volumes, many right-hand sides). Items (a)-(c) are open aspects of the same tradeoff rather than mandatory deliverables of every answer: an answer is judged on the fully-costed net time-to-solution question for the branch and regime it addresses, is not required to supply the scaling characterization of (a) or the robustness tests of (b) for preconditioning families it does not address, and need not resolve the sampler-side question (c) in order to settle the primary question. The problem is stated at the sources' natural generality: no fixed lattice size, operator, or answer form is assumed.

References

1. Abdou Abdel-Rehim, Ronald B. Morgan, Walter Wilcox, 'Deflated BiCGStab for linear equations in QCD problems' (2007), arXiv:0710.1988, https://arxiv.org/abs/0710.1988
2. Daniel Mohler, Stefan Schaefer, Jakob Simeth, 'CLS 2+1 flavor simulations at physical light- and strange-quark masses', EPJ Web of Conferences 175, 02010 (2018), doi:10.1051/epjconf/201817502010, arXiv:1712.04884, https://arxiv.org/abs/1712.04884
3. Salvatore Calì, Daniel C. Hackett, Yin Lin, Phiala E. Shanahan, Brian Xiao, 'Neural-network preconditioners for solving the Dirac equation in lattice gauge theory', Physical Review D 107, 034508 (2023), doi:10.1103/PhysRevD.107.034508, arXiv:2208.02728, https://arxiv.org/abs/2208.02728
4. Michael S. Albergo, Gurtej Kanwar, Sébastien Racanière, Danilo J. Rezende, Julian M. Urban, Denis Boyda, Kyle Cranmer, Daniel C. Hackett, Phiala E. Shanahan, 'Flow-based sampling for fermionic lattice field theories', Physical Review D 104, 114507 (2021), doi:10.1103/PhysRevD.104.114507, arXiv:2106.05934, https://arxiv.org/abs/2106.05934
5. Martin Lüscher, 'Deflation acceleration of lattice QCD simulations', Journal of High Energy Physics 0712, 011 (2007), doi:10.1088/1126-6708/2007/12/011, arXiv:0710.5417, https://arxiv.org/abs/0710.5417
6. Martin Lüscher, 'Local coherence and deflation of the low quark modes in lattice QCD', Journal of High Energy Physics 0707, 081 (2007), doi:10.1088/1126-6708/2007/07/081, arXiv:0706.2298, https://arxiv.org/abs/0706.2298
7. Walter Wilcox, 'Deflation Methods in Fermion Inverters', PoS(LATTICE 2007) 025 (2008), doi:10.22323/1.042.0025, arXiv:0710.1813, https://arxiv.org/abs/0710.1813
8. R. Babich, J. Brannick, R. C. Brower, M. A. Clark, T. A. Manteuffel, S. F. McCormick, J. C. Osborn, C. Rebbi, 'Adaptive multigrid algorithm for the lattice Wilson-Dirac operator', Physical Review Letters 105, 201602 (2010), doi:10.1103/PhysRevLett.105.201602, arXiv:1005.3043, https://arxiv.org/abs/1005.3043
9. Andreas Frommer, Karsten Kahl, Stefan Krieg, Björn Leder, Matthias Rottmann, 'Adaptive Aggregation Based Domain Decomposition Multigrid for the Lattice Wilson Dirac Operator' (2013), arXiv:1303.1377, https://arxiv.org/abs/1303.1377
10. J. C. Osborn, R. Babich, J. Brannick, R. C. Brower, M. A. Clark, S. D. Cohen, C. Rebbi, 'Multigrid solver for clover fermions' (2010), arXiv:1011.2775, https://arxiv.org/abs/1011.2775
11. M. A. Clark, Bálint Joó, Alexei Strelchenko, Michael Cheng, Arjun Gambhir, Richard Brower, 'Accelerating Lattice QCD Multigrid on GPUs Using Fine-Grained Parallelization' (2016), arXiv:1612.07873, https://arxiv.org/abs/1612.07873
12. James Brannick, Andreas Frommer, Karsten Kahl, Björn Leder, Matthias Rottmann, Artur Strebel, 'Multigrid Preconditioning for the Overlap Operator in Lattice QCD', Numerische Mathematik 132, 463 (2016), doi:10.1007/s00211-015-0725-6, arXiv:1410.7170, https://arxiv.org/abs/1410.7170
13. Azusa Yamaguchi, Peter Boyle, 'Hierarchically deflated conjugate residual', PoS LATTICE2016, 374 (2016), arXiv:1611.06944, https://arxiv.org/abs/1611.06944
14. Eloy Romero, Andreas Stathopoulos, Kostas Orginos, 'Multigrid deflation for lattice QCD', Journal of Computational Physics 409, 109356 (2020), doi:10.1016/j.jcp.2020.109356, arXiv:1909.12234, https://arxiv.org/abs/1909.12234
15. Peter A. Boyle, Dennis Bollweg, Christopher Kelly, Azusa Yamaguchi, 'Algorithms for Domain Wall Fermions' (2022), arXiv:2203.17119, https://arxiv.org/abs/2203.17119
16. Peter A. Boyle, 'Multiple right hand side multigrid for domain wall fermions with a multigrid preconditioned block conjugate gradient algorithm' (2024), arXiv:2409.03904, https://arxiv.org/abs/2409.03904
17. Anthony Francis, Patrick Fritzsch, Martin Lüscher, Antonio Rago, 'Master-field simulations of O(a)-improved lattice QCD: Algorithms, stability and exactness', Computer Physics Communications 256, 107355 (2020), doi:10.1016/j.cpc.2020.107355, arXiv:1911.04533, https://arxiv.org/abs/1911.04533
18. P. Fritzsch, J. Bulava, M. Cè, A. Francis, M. Lüscher, A. Rago, 'Master-field simulations of QCD', PoS(LATTICE2021) 465 (2022), doi:10.22323/1.396.0465, arXiv:2111.11544, https://arxiv.org/abs/2111.11544
19. Christoph Lehner, Tilo Wettig, 'Gauge-equivariant pooling layers for preconditioners in lattice QCD' (2023), arXiv:2304.10438, https://arxiv.org/abs/2304.10438
20. Yixuan Sun, Srinivas Eswar, Yin Lin, William Detmold, Phiala Shanahan, Xiaoye Li, Yang Liu, Prasanna Balaprakash, 'Matrix-free Neural Preconditioner for the Dirac Operator in Lattice Gauge Theory' (2025), arXiv:2509.10378, https://arxiv.org/abs/2509.10378
21. Simon Pfahler, Daniel Knüttel, Christoph Lehner, Tilo Wettig, 'A novel gauge-equivariant neural network architecture for preconditioners in lattice QCD', PoS(LATTICE2025) 027 (2026), doi:10.22323/1.518.0027, arXiv:2602.23840, https://arxiv.org/abs/2602.23840
22. Ryan Abbott, Michael S. Albergo, Denis Boyda, Kyle Cranmer, Daniel C. Hackett, Gurtej Kanwar, Sébastien Racanière, Danilo J. Rezende, Fernando Romero-López, Phiala E. Shanahan, Betsy Tian, Julian M. Urban, 'Gauge-equivariant flow models for sampling in lattice field theories with pseudofermions', Physical Review D 106, 074506 (2022), doi:10.1103/PhysRevD.106.074506, arXiv:2207.08945, https://arxiv.org/abs/2207.08945
23. Travis Whyte, Andreas Stathopoulos, Eloy Romero, 'Accelerating multigrid with streaming chiral SVD for Wilson fermions in lattice QCD' (2025), arXiv:2505.14399, https://arxiv.org/abs/2505.14399
24. Pauline Schauerte, Jaime Fabián Nieto Castellanos, Arnold Krechel, Marc Alexander Schweitzer, Stefan Krieg, 'A first look at Structured-Multiscale Algebraic Multigrid for Lattice Field Theory' (2026), arXiv:2608.03515, https://arxiv.org/abs/2608.03515