Matrix-free Neural Preconditioner for the Dirac Operator in Lattice Gauge Theory

Phiala Shanahan; Prasanna Balaprakash; Srinivas Eswar; William Detmold; Xiaoye Li; Yang Liu; Yin Lin; Yixuan Sun

Matrix-free Neural Preconditioner for the Dirac Operator in Lattice Gauge Theory

T0 review · 2 major / 4 minor · reviewed 2026-08-04 · deepseek-v4-flash

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

Pith's one-line read A matrix-free neural network that rewrites the gauge field roughly halves Dirac solver iterations in the Schwinger model and transfers to new lattice sizes.

desk verdict A well-executed proof-of-concept for a matrix-free learned preconditioner, but the paper's own Section 4.3 shows a trivial analytic map does almost as well, which weakens the operator-learning claim. read the letter →

arxiv 2509.10378 v2 pith:GKJ4LBXK submitted 2025-09-12 hep-lat cs.LG

Yixuan Sun , Srinivas Eswar , Yin Lin , William Detmold , Phiala Shanahan , Xiaoye Li , Yang Liu , Prasanna Balaprakash This is my paper · ORCID

classification hep-latcs.LG MSC 65F0865F1068T0781T25 PACS 11.15.Ha12.38.Gc

keywords preconditionerneuraloperatorlatticegaugetheoryWilson-Diracconjugategradientmatrix-freezero-shottransferSchwingermodel

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

The paper proposes a way to precondition the large, ill-conditioned linear systems that dominate lattice gauge theory calculations, without ever forming or storing a matrix. The central idea is to train an operator-learning network to map a gauge field configuration U to a new configuration Ũ such that the Wilson-Dirac normal operator built from Ũ, M^{-1}=D†D[Ũ], acts as an approximate inverse of the original operator A=D†D[U]. Because M^{-1} has the same sparse stencil as A, applying it costs one matrix-vector product, and the whole pipeline is matrix-free from training through solver use. In the 2D Schwinger model the preconditioner reduces the condition number by roughly a factor of four and cuts CG iterations approximately in half, and a network trained on one lattice size works on other sizes and couplings without retraining. The learned mapping turns out to be nearly Ũ ≈ -½U, a simple analytic rule that itself already provides much of the acceleration.

What carries the argument

The load-bearing object is the map Gθ taking a U(1) gauge field U to a transformed field Ũ, together with the fixed rule that the preconditioner is M^{-1}=D†D[Ũ]. This ties the preconditioner to the same sparse stencil as A, so matrix-vector products are cheap and no matrix is ever assembled; two operator-learning architectures (a local convolutional network and a global spectral-mode network) are trained to produce Ũ. Training uses an unsupervised random-projection loss, Eq. (8): the squared norm of M^{-1}A v - v averaged over K random vectors v, which proxies for the Frobenius norm of M^{-1}A - I without computing spectra. A final piece is the identity A=D†D[U]=|1-2κH|², with H linear in U

What would settle it

Take the analytic rule Ũ = -½U, build M^{-1}=D†D[Ũ], and measure the condition number of M^{-1/2} A M^{-1/2} (or the CG iteration count) on SU(3) Wilson-Dirac operators at physical parameters in 4D. If the iteration reduction does not approach the roughly twofold gain seen in the Schwinger model—or if the spectrum of M^{-1}A develops clusters far from 1—then the learned mapping and the simple analytic form are artifacts of the Abelian 2D model, and the zero-shot transfer claim is restricted to that setting.

Watch

Extended reading notes

Core claim

The central claim is that an effective preconditioner for Wilson-Dirac normal equations can be learned as a function of the gauge field alone. The construction fixes M^{-1}=D†D[Ũ(x)] with Ũ(x)=Gθ(U(x)), so the preconditioner is the same discretized normal operator evaluated on a transformed gauge field and therefore shares the original matrix's sparse structure. Training is unsupervised and matrix-free: the loss is the average squared norm of (M^{-1}A - I) applied to K random vectors from an isotropic standard normal distribution. With this loss, a network trained on 1600 configurations of the two-flavor Schwinger model reduces the condition number of A by roughly a factor of four across L=8

Load-bearing premise

The load-bearing premise is that forcing the preconditioner to have exactly the same sparse structure as the original Wilson-Dirac normal operator—namely, evaluating D†D on a transformed gauge field—still leaves enough flexibility for M^{-1} to approximate A^{-1} well enough to accelerate conjugate gradient, and that the training loss, which checks the error only on random vectors, is a faithful enough proxy for that approximation.

Editorial extensions

If this is right

Preconditioners for lattice fermion solves can be generated on the fly from gauge fields with no matrix assembly or setup step, making the approach scalable to the very large operators (dimension greater than 10^8) encountered in lattice QCD.
A single training run on a small lattice can be amortized across many action parameters and lattice geometries: zero-shot volume transfer means no retraining for unseen configurations within the explored parameter ranges.
Since per CG iteration the neural preconditioner costs one matrix-vector product and no triangular solve, it avoids the per-iteration overhead of incomplete-factorization or even-odd preconditioners, at the cost of somewhat more iterations than the incomplete-factorization approach.
Using (M^{-1})^p with p=2, 3, or 4 as the preconditioner reduces iterations further (e.g., the FCN-based model at L=16 drops from 99 to 79–81 iterations), providing a tunable iteration-vs-per-step-cost knob.
The explicit analytic choice Ũ = -½U already gives most of the acceleration, implying that a useful preconditioner may be constructible without training at all in the Schwinger-model setting.

Reading between the lines

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

If the learned mapping generalizes to non-Abelian groups, the simple form Ũ ≈ -½U is unlikely to survive because link variables are unitary matrices rather than phases; a non-Abelian analog would more plausibly be a group-specific redefinition, so the preconditioner structure should be tailored to the gauge group.
The random-projection loss, which the paper acknowledges lacks theoretical grounding, could be replaced by a spectral proxy such as stochastic estimates of the condition number; the paper's own sensitivity study suggests K=128 is a pragmatic choice, and a principled objective might close the remaining gap to classical incomplete-factorization preconditioners.
The zero-shot volume transfer indicates the mapping depends mostly on local features of the gauge field, consistent with the local-network's superior transfer at L=8. A testable extension is that the preconditioner will continue to transfer to even larger lattices (e.g., L=128) with bounded loss of efficiency; if instead iteration counts grow as in the unpreconditioned case, the local-feature assu
One could apply the same construction to other fermion discretizations (staggered, domain-wall) by using the corresponding normal operator as the fixed stencil; this would test whether the learned mapping again reduces to a simple field redefinition or requires nontrivial correlations.

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

Read this if you want to see a careful attempt at a matrix-free neural preconditioner for Wilson-Dirac normal equations. The authors train FNO/FCN networks to map a gauge configuration U to a transformed configuration U~, then use D^†D[U~] as a preconditioner in CG. Experiments are on the 2D Schwinger model (U(1) gauge theory, two flavors). What is genuinely new: the matrix-free formulation (no explicit matrices anywhere, only operator applications), the random-projection loss that avoids condition-number computation, and zero-shot transfer across lattice sizes and coupling values. The paper is honest: it reports iteration counts with error bars, compares against IChol and even-odd preconditioners, and acknowledges the random-projection loss is an unproven proxy.

The main soft spot is the one the paper exposes in Section 4.3. The learned transformation is essentially U~ = -0.5U with a π phase shift. The authors test this simple analytic map explicitly and find it gives 43 vs 40 CG iterations for L=8, 104 vs 99 for L=16, 278 vs 267 for L=32. The neural network buys only a few percent over a parameter-free transformation. That doesn't invalidate the framework, but it does mean the 'operator learning' contribution is, for this test case, mostly learning a pointwise scaling. The paper notes the small gap suggests correlations across sites and directions matter, but it doesn't quantify how much or show a regime where the simple map fails. Without such evidence, the advantage of the network over the analytic baseline is not established.

Second soft spot: no wall-clock timings. The paper calls the preconditioner 'computationally efficient,' but applying M^-1 costs one D^†D application, the same complexity as applying A. If PCG iterations halve, time per iteration roughly doubles, so the wall-clock gain is unclear. A timing table would settle this. Also, the code is not yet public, and the higher-power results (Table 3) show modest gains at increased per-iteration cost.

Who should read this: people working on ML preconditioners for lattice field theory or sparse iterative solvers. It is a legitimate proof of concept with honest reporting, and the matrix-free idea may matter for 4D QCD if it transfers. The central empirical claim (iteration reduction and zero-shot transfer) holds up as an existence proof. But I would not let the paper's framing overstate the neural network's role before comparing to the analytic baseline in the target regime.

Recommendation: send for peer review. The method deserves discussion, and the flaws are addressable: add timings, compare against U~=-0.5U systematically across parameter ranges, and either show a case where the network outperforms the simple map or temper the claims accordingly.

Referee Report

2 major / 4 minor

Summary. The paper proposes a matrix-free neural preconditioner for Wilson-Dirac normal equations in lattice gauge theory. Given a gauge configuration U, a neural network Gθ (FNO or FCN) outputs a transformed configuration Ũ; the preconditioner M^{-1}=D†D[Ũ] is then applied as a linear operator inside preconditioned CG. Training uses a random-projection loss that minimizes ||M^{-1}D†D v - v||² without ever forming matrices explicitly. Experiments in the 2D U(1) Schwinger model with two Wilson fermions show that, for lattices L=8–64 at κ=0.276, β=2.0, the preconditioner reduces the condition number and roughly halves CG iterations; a model trained at L=16 transfers zero-shot to other volumes and to modified κ and β. Section 4.3 reports that the learned transformation is approximately Ũ = -1/2 e^{iπ} U, and Table 4 shows that the manually constructed analytic preconditioner D†D[-1/2 U] gives nearly identical CG iteration counts to the trained networks.

Significance. If the claimed gains hold, the matrix-free construction is a useful step toward ML preconditioners for large-scale lattice calculations, since it avoids explicit matrix storage and setup/decomposition costs and offers zero-shot volume transfer. The paper is transparent about its limitations: the loss is an acknowledged proxy, the simple learned map is identified, and higher powers of M^{-1} are proposed but not fully analyzed. However, the central added value of the learned operator over the analytic transformation is not established. The paper's own Table 4 shows the neural network outperforms the analytic pointwise mapping by only 3–11 CG iterations out of 78–548, which is a small fraction of the total solve. This weakens the 'operator learning' claim: the evidence suggests the network has essentially rediscovered a simple scaling/phase relation. The framework is plausible and the empirical methodology is sound, but the current manuscript does not yet demonstrate that the neural component is necessary or advantageous relative to a parameter-free analytic preconditioner.

major comments (2)

[Section 4.3, Table 4] This section reports that the learned maps satisfy |Ũ|≈1/2 and a phase difference of π, i.e., Ũ≈-1/2 U. Table 4 then shows that the manually constructed preconditioner D†D[-0.5U] gives CG iteration counts of 43±2, 104±7, 278±20 for L=8,16,32, while the trained FNO gives 40±2, 99±7, 267±21. The differences are only 3, 5, and 11 iterations, out of unpreconditioned counts of 78, 201, and 548. This is load-bearing: the paper's main contribution is that an operator-learning framework learns a useful preconditioner, but the evidence indicates the network has recovered a simple pointwise scaling/rotation that can be implemented without any training. The authors should (i) quantify how much Gθ deviates from -0.5U over the whole test set rather than a single L=64 instance (Fig. 7), (ii) test regimes where the small-κ expansion argument in §4.3 is expected to fail (e.g., κ further from criticality
[Section 4.1 and Remark after Table 2] The abstract and Section 4.1 claim computational efficiency and reduced complexity, and the Remark states that the NN preconditioner avoids IChol's setup cost and triangular solve. However, no wall-clock timings are reported anywhere. Figure 4 shows only iteration counts, even though the text in §4.1 says 'we compare the average number of iterations and the time required to reach convergence.' Since CG per-iteration cost differs between one sparse matrix-vector product with M^{-1} and an IChol triangular solve, iteration counts alone do not establish a computational advantage. The paper should report setup time, per-iteration cost, and total solve time on the same hardware (A100) for unpreconditioned, IChol, even-odd, and NN-preconditioned solves, at least for L=32 and L=64. Without these data, the 'computationally efficient' claim is unsubstantiated.

minor comments (4)

[Title page / affiliations] The affiliation 'Lawrance Berkeley National Laboratory' should be 'Lawrence Berkeley National Laboratory'.
[Eq. (2)] The fermion action term is rendered as '1X f=0' due to a formatting issue; it should be a sum over flavor index f=0,1.
[Eq. (5) and Fig. 6] The notation 'U(1)^{X×T×d} ⊂ C^{X×T×d}' is confusing: U(1) is the target space of each component, not a power of U(1). Clarify the spaces. In Fig. 6, the scatter plots would be more informative with a quantitative similarity metric such as a correlation coefficient or normalized distance.
[Code availability] The text says the code is available at the given URL, but the footnotes state the repository 'will be publicly accessible upon acceptance.' These statements should be made consistent.

Circularity Check

0 steps flagged · score 0.0 of 10

No significant circularity: training objective is a direct preconditioning objective evaluated on held-out/zero-shot data; the discovered simple mapping is a post-hoc control.

full rationale

The paper's derivation chain is not circular. The central object is a neural network Gθ mapping gauge configurations U to preconditioner-generating configurations U~, trained by minimizing the random-projection loss L(θ) in Eq. (8), which directly targets M^{-1}A ≈ I. All headline quantities (condition numbers, CG iterations) are measured on a held-out test set of 200 configurations, and the volume-transfer results apply the L=16 model to unseen sizes and parameters; these are not quantities fitted during training. The self-citation to Ref. [6] (Calì et al., which shares authors Lin and Shanahan) is contextual and methodological: the paper states 'Our work closely follows [6]' but the present contribution (matrix-free operator-learning preconditioner) is independently evaluated against IChol and even-odd baselines. No uniqueness theorem or ansatz is imported from the authors' prior work to force the conclusion. The post-hoc observation in Sec. 4.3 that the learned map approximates U~ = -0.5U, and the control in Table 4 showing the analytic preconditioner nearly matches the NN, is a critique of novelty/utility but not circularity: the simple map was not used to construct labels or training data, and the paper explicitly tests it as a check. The paper also openly notes the loss is a proxy and not theoretically grounded. Thus there is no step where a prediction reduces to an input by construction.

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

The central claim rests on the specific lattice discretization, the structural constraint that the preconditioner is itself a Dirac normal operator on a modified gauge field, and the empirical assumption that the random-projection loss trains a useful preconditioner. No new physical entities are introduced. The fitted parameters are the network weights and the hyperparameters K and p.

free parameters (3)

Network weights (FNO/FCN) = learned during training
All model parameters are fitted by minimizing the random-projection loss Eq. (8); the central claim depends on the specific learned mapping.
Number of random vectors K = 128
Hyperparameter controlling the loss estimate; Appendix B shows performance plateaus above K=64.
Power p of preconditioner = 1,2,3,4
Higher powers (M^{-1})^p densify the preconditioner; p=4 is fully dense for L=16 and reduces iterations further, at extra per-step cost.

assumptions (4)

domain assumption Wilson discretization of the Dirac operator and the 2D Schwinger action define the linear system.
The paper's entire setup is this lattice action; if this discretization is not representative of target applications, the method's value is limited.
standard math D†D is Hermitian positive definite and CG is an appropriate solver.
Standard result used implicitly.
domain assumption The random-projection loss is a sufficient proxy for the true objective of minimizing the condition number of the preconditioned system.
The paper acknowledges this is not theoretically grounded (Section 3.3) and uses an empirical proxy.
domain assumption A neural operator trained on one lattice size learns a mapping that generalizes to other sizes and parameters.
This is the zero-shot claim; tested empirically on L=8,16,32,64 with mixed success.

how reviews work

0 comments

Cite this review

Pith. "Pith review of Matrix-free Neural Preconditioner for the Dirac Operator in Lattice Gauge Theory." pith.science (2026). https://pith.science/paper/GKJ4LBXK

@misc{pith2026250910378,
  author       = {Pith},
  title        = {Pith review of: Matrix-free Neural Preconditioner for the Dirac Operator in Lattice Gauge Theory},
  year         = {2026},
  howpublished = {\url{https://pith.science/paper/GKJ4LBXK}},
  note         = {Machine review of arXiv:2509.10378}
}

read the original abstract

Linear systems arise in generating samples and in calculating observables in lattice quantum chromodynamics~(QCD). Solving the Hermitian positive definite systems, which are sparse but ill-conditioned, involves using iterative methods, such as Conjugate Gradient (CG), which are time-consuming and computationally expensive. Preconditioners can effectively accelerate this process, with the state-of-the-art being multigrid preconditioners. However, constructing useful preconditioners can be challenging, adding additional computational overhead, especially in large linear systems. We propose a framework, leveraging operator learning techniques, to construct linear maps as effective preconditioners. The method in this work does \emph{not} rely on explicit matrices from either the original linear systems or the produced preconditioners, allowing efficient model training and application in the CG solver. In the context of the Schwinger model U(1) gauge theory in 1+1 spacetime dimensions with two degenerate-mass fermions), this preconditioning scheme effectively decreases the condition number of the linear systems and approximately halves the number of iterations required for convergence in relevant parameter ranges. We further demonstrate the framework learns a general mapping dependent on the lattice structure which leads to zero-shot learning ability for the Dirac operators constructed from gauge field configurations of different sizes.