Machine learning for lattice QCD

LA-UR-26-29844 MIT-CTP/6147

Machine learning for lattice QCD Gurtej Kanwar1 , Scott Lawrence2 , Alessandro Nada3 , and Phiala Shanahan4,5 1 Higgs Centre for Theoretical Physics, The University of Edinburgh, James Clerk Maxwell

Building, Peter Guthrie Tait Road, Edinburgh, EH9 3FD, UK 2 Theoretical Division, Los Alamos National Laboratory, Los Alamos, NM 87545, USA

arXiv:2611.09817v1 [hep-lat] 6 Nov 2026

3 Dipartimento di Fisica, Università degli Studi di Torino and INFN, Sezione di Torino, Via Pietro

Giuria 1, 10125 Turin, Italy 4 Center for Theoretical Physics – a Leinweber Institute, Massachusetts Institute of Technology,

Cambridge, MA 02139, USA 5 The NSF AI Institute for Artificial Intelligence and Fundamental Interactions

November 9, 2026 Abstract Machine learning (ML) methods have been explored as a means of optimizing lattice field theory calculations, with increasing interest in recent years. This chapter gives an introduction and summary to the use of ML methods to optimize lattice field theory without introducing any uncontrolled approximations. We focus on those directions which are clearly established as of late 2025, marking the 50th anniversary of the beginning of the study of lattice gauge theory. This chapter is a contribution to the online book “Lattice QCD at 50 years,” (LQCD@50), edited by Tanmoy Bhattacharya, Maarten Golterman, Rajan Gupta, Laurent Lellouch, and Steve Sharpe.

Contents 1 Introduction

2

2 Brief introduction to ML/AI

4

3 Monte Carlo acceleration 3.1 Markov Chain Monte Carlo . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 New sampling algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Generative models . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Estimating the partition function . . . . . . . . . . . . . . . . . . . . . . . .

7 8 9 13 20

4 Observable estimators 4.1 Motivation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Control variates . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

21 22 24

1

4.3 4.4

Contour deformations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Correlated samples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

25 29

5 Observable acceleration 5.1 Surrogate observables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.2 Preconditioners . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

30 30 31

1

Introduction

Machine learning (ML) is a branch of artificial intelligence that enables pattern identification and the development of data-driven predictions or decision-making. Many applications of ML are based on neural networks, which are highly flexible models, typically with huge numbers of free parameters, that are capable of fitting high-dimensional, nonlinear functions and are particularly well-suited for tasks such as interpolation, regression, and the emulation of computationally expensive operations. In scientific computing, ML techniques can be leveraged to optimize algorithms, accelerate simulations, and uncover structures in complex data; see Refs. [1–3] for accessible introductions to ML. In Lattice Quantum Field Theory (LQFT), ML has now been applied (with varying levels of maturity and/or success) to accelerate every part of the calculational workflow, including: defining lattice actions, sampling gauge field configurations, accelerating solvers, observable estimation, and aspects of analysis [4]. A separation can be made between applications of ML in LQFT that are model-based (e.g., efforts to define order parameters, the use of neural networks in fits and parameterisations of data, and applications to inverse problems), and those that demand exactness guarantees in the sense of asymptotic correctness. This is an important distinction; unlike many areas of applied ML where approximations and heuristic or ‘black box’ methods are acceptable, LQFT is inherently a first-principles approach: therefore, integration of ML techniques into many aspects of the LQFT workflow must be undertaken in such a way that any mathematical guarantees afforded by conventional algorithms are preserved. This category of applications—those where unbiasedness must be guaranteed—is the focus of this chapter. Methods of accelerating computations through ML while maintaining the mathematical rigour of traditional approaches can be broadly grouped into three classes based on the role of ML within the computational workflow: • Optimise existing algorithms: One straightforward strategy is to accelerate existing computations by tuning algorithmic parameters using ML. This approach preserves the structure and logic of the original algorithm while enhancing its performance. For instance, free parameters in algebraic multigrid methods for solving linear systems (e.g., for propagator inversions) can be optimized using ML techniques. Similarly, replacing stochastic optimization approaches in various contexts with automatic differentiation can lead to performance improvements without sacrificing exactness. • Transform the problem: A second strategy involves using ML to transform the problem into an equivalent formulation that is computationally more tractable to solve than the 2

original. Standard algorithmic approaches, with the associated guarantees, are then applied to the transformed problem. This includes techniques such as preconditioning, which can significantly improve the convergence rate of numerical solvers, and change of variables—for example, deforming the integration contour—which can reduce the variance of numerical integration via Monte Carlo without altering the central value. • Solve a related problem and correct: A third approach leverages ML to solve a related but more accessible problem and then apply a known correction to achieve the exact solution. One example includes learning a mapping between a low-precision and a highprecision solution in linear systems and subsequently correcting for the introduced bias. Another example involves sampling from a nearby probability distribution followed by a reweighting or resampling step to obtain the correct distribution. This method allows for flexibility in computational methodology while ensuring that the final results remain provably exact. Table 1 collects a (non-exhaustive) grid of these three classes of applications of asymptoticallycorrect ML-accelerated approaches to the different steps of the LQFT workflow. Some applications are in the early stages of investigation at the current time, while others are approaching maturity. None are, as of late 2025, yet the default approach in state-of-the-art LQFT calculations at scale. An important reason that ML approaches are still far from being default in LQFT applications is the wide gap between the standard ML/AI workflow and the specialized workflows and infrastructure that have developed within LQFT. Most ML research and engineering is carried out in Python, using a small set of powerful toolkits—most prominently TensorFlow [5], PyTorch [6], and JAX [7]. These libraries provide highly optimized routines for defining neural network architectures and performing gradient-based optimization (e.g., through backpropagation). They have been designed to lower the barrier to entry, allowing rapid prototyping and flexible experimentation for applications close to typical ML use cases. By contrast, the computational infrastructure in LQFT has been shaped by decades of largescale scientific computing. LQFT calculations are typically performed using custom-built software packages such as Chroma [8], CPS [9], or QUDA [10], which have been engineered to exploit the architectures of high-performance computing (HPC) platforms. These codes are tuned for massively parallel execution on distributed-memory supercomputers, optimized for efficient communication over interconnects, and adapted to exploit GPU acceleration at scale. Every layer of the software stack has been refined to address the specific numerical challenges of simulating gauge field configurations, propagators, and correlators in gauge theories. Bringing ML methods into this environment introduces substantial technical challenges. Standard neural network architectures—such as convolutional or transformer-based models—are not naturally aligned with the data structures or symmetry properties of lattice gauge theory. In typical applications of ML, the conventional wisdom (the “bitter lesson” [11]) is that incorporating symmetries is often not necessary or even helpful. The (gauge) symmetries of lattice QCD, however, are extremely high-dimensional, and anecdotal 3

evidence to date strongly indicates that it is worth including manifest gauge (in/co)variance in the design of networks being trained [12]. Symmetry-respecting neural networks have been developed to address this challenge [13]: by explicitly encoding gauge invariance or rotational symmetries into the architecture, these models can capture the essential structure of the problem in a way that black-box networks cannot. These applications, however, are specialized, and scaling custom models efficiently across tens of thousands of nodes or GPUs, with careful management of memory hierarchies and communication bottlenecks, requires significant re-engineering. That is, the high-throughput, distributed training strategies common in commercial ML do not map easily onto the specialized HPC workflows that define state-of-the-art lattice calculations. Thus, even when an ML model shows promise in smallscale exploratory studies, extending it to production-level simulations demands overcoming both algorithmic and systems-level hurdles. The long-term challenge, then, is twofold: first, to continue developing specialized ML architectures that respect the mathematical structure of gauge theories, and second, to integrate them into HPC environments at a level of scalability and efficiency comparable to traditional lattice methods. Meeting both demands is essential if ML is to become a true complement—or even a replacement—in the established LQFT computational toolkit. Despite the magnitude of these challenges, the rapid pace of innovation at the intersection of ML and physics gives strong reason for optimism that ML will ultimately become an integral part of the future of LQFT. The remainder of this chapter is divided into four sections, with a brief introduction to ML followed by the ML approaches to LQFT categorized depending on the manner in which the lattice calculation is accelerated or improved. If the main computational obstacle is performing importance sampling—that is, if the observable of interest is fast to evaluate on any given sample—then using machine learning methods to accelerate the Monte Carlo sampling itself is of greatest value. This class of approaches, exemplified by normalizing flows, is discussed in Section 3. If the total number of samples on which an observable is to be estimated is held fixed, then ML-based acceleration must improve the signal-to-noise ratio by modifying the observable itself (e.g. via control variates); this is discussed in Section 4. Finally, if there are a large number of samples available but evaluating the observable is itself computationally challenging, then ML may be used to accelerate (or approximate) the evaluation of the observable; this is discussed in Section 5.

2

Brief introduction to ML/AI

Here we include a brief introduction to machine learning methods, with a bias towards the points of view that have been most productive in computational physics. This introduction is not intended to be comprehensive for a reader entirely new to the field. We recommend such readers refer to the introductory material for PyTorch, which provides a good practical introduction to ML, or Refs. [1–3]. In great generality, many calculations may be accelerated or improved if we have numerical access to some function, which has some desired mathematical property, and which can 4

Approach

ML application

Correlated sampling Optimize existing algorithm

Tune MG params. Learned interpolators Perfect actions

Transform the problem

Learned-manifold HMC Learned preconditioners Contour deformation

Solve related problem (and correct)

Generative models Surrogate observables

Area

Obs. estimation Config. generation, Obs. estimation Obs. estimation Action

Config. generation Solvers & inversion Obs. estimation, Action Config. generation Obs. estimation

Table 1: Classification of some ML applications to LQFT based on the part of the computational pipeline that is targeted, and the approach to maintaining exactness. Table is not exhaustive.

5

be evaluated sufficiently quickly. “Machine learning” is the name given to the process of searching for such a function. To have a concrete example, we might try to learn a function f (U ) which does a good job of approximating the action S(U ). This would be useful when the action is slow to evaluate, for example, and is not too different from some examples we will see later in this text. We want our function f (U ) to be as close as possible to S(U ), and we want it to be fast to evaluate. The paradigm of machine learning proceeds in the following steps: 1. Identify a family of functions, fθ (U ), defined by a set of real parameters θ. 2. Define a measure of ‘loss’; a functional L[fθ ] which is small when fθ (·) has our desired property. 3. Perform some search algorithm—typically an elaboration on gradient descent—to find parameters θ such that L(θ) ≡ L[fθ ] is as small as possible. This optimization step is often referred to as “training”. 4. Employ the resulting fθ (·) in a calculation. In the first step we must identify a family of functions. The reader will not be surprised that this family is often defined by a neural network: the weights constitute the parameters θ, and the gauge configuration U shall be passed as an input. The most straightforward neural network is the multi-layer perceptron: an alternating composition of element-wise nonlinear transformations with arbitrary linear transformations1 . In the case of one hidden layer (meaning one nonlinear transformation, sandwiched by arbitrary matrix multiplications), this may be written:   (1) (2) (1) NN(x)i = Mij σ Mjk xk , where the function σ : R → R is a nonlinear function, termed an “activation function”.2 Why a neural network? Historically, neural networks were justified by somewhat stretched analogies to the structure of the nervous systems of living things. Today it is probably more accurate to simply observe that, because they rely so heavily on linear algebra, neural networks are an effective mechanism for packing a large number of free parameters into a relatively small amount of computation. Moreover, the computations that are involved map nicely onto modern computing architectures, such as GPU-based machines. In the second step we must define a loss functional L[fθ ], which is also implicitly a loss function L(θ) of the parameters. This loss function will typically be evaluated stochastically. In the case where we are seeking some function which approximates the action, a suitable loss functional might be Z L[fθ ] = dU [S(U ) − fθ (U )]2 , 1

Sparser architectures are possible and in fact common, but the fundamental idea—that element-wise nonlinear transformations and fairly arbitrary linear transformations are both fast—recurs. 2 Typical examples are σ(x) = tanh(x) or σ(x) = max(0, x).

6

which may be approximated by Monte Carlo sampling from the measure dU . It is in the third step that the computational part of the learning actually happens: we search for parameters θ that achieve a small loss L(θ). If L(·) has been chosen wisely, then achieving a small-but-not-minimal loss will still be helpful. Because L(·) is typically expressed as an integral and can only be evaluated by Monte Carlo, we perform stochastic gradient descent or one of its relatives. All3 follow the same basic pattern:   ∂L(θ) . (2) θn+1 = θn − ϵ ∂θ In stochastic gradient descent, ϵ is a fixed number (termed the “learning rate”) and the partial derivative is evaluated by Monte Carlo as discussed above. In more elaborate methods, including the widely used Adam [14], the estimate of the gradient and/or the learning rate is modified, taking into account the history of the gradients (e.g. to add momentum, or estimate local curvature). Much of the philosophy behind modern machine learning and AI is reflected in the following heuristic: having more parameters, and spending more time training, is what ultimately drives a successful minimization of the loss. These factors have historically overwhelmed advanced ansätze or clever minimization methods. This observation is sometimes called the “bitter lesson” [11]. This does not imply that no thought is required in machine learning: there is no escaping the need to select a meaningful loss function, for example. But it does provide a guide. The purpose of the methods in this chapter is, by and large, to find new ways of bringing a great deal of computational power to bear on specific aspects of LQCD.

3

Monte Carlo acceleration

The fundamental problem in lattice QCD is calculating thermal expectation values via the Euclidean path integral Z 1 DU e−S(U ) O(U ) (3) ⟨O⟩ = Z with O an operator to be measured and Z being the partition function of the theory, Z Z = DU e−S(U ) . (4) This problem can be numerically solved by importance sampling according to the Boltzmann distribution e−S(U ) /Z defined by Eq. (3). Given an ensemble of field configurations {U (1) , . . . , U (N ) } sampled in this way, the sample mean provides a Monte Carlo estimate of the expectation value, N 1 X ⟨O⟩ ≃ O(U (i) ). (5) N i=1 Increasing the number of samples N systematically lowers the uncertainties on this estimate. 3

Of course there are approximate minimization algorithms which depart from this paradigm. The most extreme examples are gradient-free methods. These lie outside the scope of this chapter.

7

Cost

Hybrid Monte Carlo Flow-based sampling

N∗

Generated samples

Figure 1: Schematic comparison of the sampling cost using a traditional Monte Carlo sampler, such as Hybrid Monte Carlo, versus an optimized flow-based sampler. The up-front cost of training a flow-based sampler is typically larger than the cost of thermalizing a Markov Chain sampler, while the cost per uncorrelated sample is typically smaller for a flow-based sampler. If the number of required samples exceeds the threshold N ∗ , the overall cost of sampling is lower when using the flow-based sampling strategy.

3.1

Markov Chain Monte Carlo

Markov Chain Monte Carlo algorithms provide a practical way to perform the importance sampling above by applying a memoryless sequence of stochastic updates U (1) → U (2) → U (3) → . . .

(6)

that are constructed to asymptotically converge to the Boltzmann distribution. The overall cost of producing N independent, thermalized field configurations over which the average of Eq. (5) is calculated can be expressed as the cost Cth of the burn-in (or thermalization) of the Markov Chain plus the cost N Cs of sampling the N decorrelated configurations, such that the total cost is Ctot (N ) = Cth + N Cs ,

(7)

as shown in Figure 1. The cost per sample Cs is determined by a combination of the computational cost of individual Markov chain updates and the number of updates required to obtain an independent sample. Hybrid Monte Carlo [15] updates are the workhorse of lattice QCD simulations, based on integrating fictitious Hamiltonian dynamics to propose the next configuration in the Markov chain, U (i−1) −−−−−−→ U ′ . (8) Ham. dyn.

8

In Hybrid Monte Carlo, as with many Markov chain methods, this proposal must be corrected with a Metropolis-Hastings step which accepts the proposal for the next state in the chain with probability
A(U (i−1) , U ′ ) = min 1, e−∆H . (9) Here the integrator error determines the violation ∆H of energy conservation in the integrated dynamics and therefore the acceptance rate and efficiency of the calculation. 3.1.1

Critical slowing down

For a given observable O and at fixed cost per Markov chain update, the cost Cs per independent (thermalized) sample of O(U ) typically scales as Cs ∝ 2τint (O) ∝ a−z .

(10)

The scaling of this cost is set by the integrated autocorrelation time τint (O), which measures the correlations in the measured observables O between sequential field configurations in the Markov chain [16]. It is well known that τint (O) diverges as the lattice spacing a approaches zero, or equivalently as the parameters are tuned towards a critical point where the correlation length diverges. In particular, when taking the continuum limit it diverges at least polynomially (as in Eq. (10)), with an exponent z depending on the observable O and on the specific MCMC algorithm used during the simulation [17]. When sampling topological observables, the slowing down is particularly severe, with the exponent z being particularly large, in a phenomenon denoted as topological freezing [18–22]. Let us stress that this is an issue of Monte Carlo algorithms based on local updates and not of the observables themselves.

3.2

New sampling algorithms

Sampling algorithms that avoid or modify traditional MCMC approaches have been developed with the explicit goal of mitigating the rapid increase in costs indicated in Eq. (10). We expect these methods to be computationally cheaper than traditional MCMC methods only in a certain regime determined by the total start-up time (e.g., training time for ML methods) and the cost Cs per independent sample. For complex ML architectures, for example, the training costs may easily exceed the thermalization costs, providing an advantage only when sufficiently many configurations need to be drawn. A qualitative picture of this situation is shown in Figure 1, and identifying the associated cross-over in advantage may be useful as methods move towards practical applications. In this section we provide a general description of some of these approaches, in particular those in which ML tools have the potential to be highly disruptive. We will describe the implementation of ML tools for flow-based approaches in section 3.3.

9

3.2.1

Learned MCMC dynamics

One approach to circumventing the challenges of critical slowing down is to explicitly modify the MCMC dynamics to avoid slow-mixing modes as criticality is approached. For example, Fourier acceleration [23] was proposed to mitigate critical slowing down for theories behaving similarly to a free scalar theory; however, obstacles remain in implementing this approach for gauge theories [24, 25]. A generalization of this approach is to learn more arbitrary reparameterizations: implementations for Hybrid Monte Carlo dynamics have been investigated [26–29]. 3.2.2

Learned actions

Critical slowing down in MCMC simulations is generally related to the use of local updates. In fact, in the case of statistical models for which a global update is known, critical slowing down in the proximity of phase transitions is successfully mitigated [30, 31]. In general, though, the development of a global MCMC update is a highly non-trivial problem. One way to design such an update is to introduce an effective action (or Hamiltonian) which can be efficiently simulated using global updates. In the Self-Learning Monte Carlo method [32–41] the couplings of this effective action are learned through a training process using configurations generated with the original action. Then, a new MCMC is run by proposing new configurations using the effective action and accepting/rejecting them using the original one. In this way, a properly-trained effective long-range action can generate uncorrelated samples, without compromising exactness. A machine learning approach can be naturally applied in this setup employing an effective action parametrized by neural networks. In the case of lattice gauge theory, this has led to the construction of gaugecovariant neural networks [42]; for further applications see refs. [43, 44]. Another intriguing approach to circumvent the issue of critical slowing down in lattice gauge theory is to use a new action characterized by strongly-suppressed lattice artifacts: this is the goal of the so-called fixed point (or classically-perfect) actions [45, 46], which are based on renormalization group transformations, and would allow better control of continuum extrapolations from simulations at coarser lattice spacings. A typical challenge when building this class of actions is finding a compact and accurate parametrization. This problem can be tackled in a new way using neural network-based architectures, which are much more expressive than the original approaches. In particular, recent efforts [47, 48] show how lattice gauge-equivariant convolutional neural networks (L-CNN) [13, 49] hold great promise in designing fixed-point actions in an efficient way. 3.2.3

Trivializing flows

An alternate approach is to transform the probability density to be sampled. Taking the path integral of Eq. (3), if we perform a change of variables V = f −1 (U ), we obtain Z 1 DV O(f (V )) e−S(f (V ))+log | det Jf (V )| (11) ⟨O⟩ = Z 10

with det J being the determinant of the Jacobian of the transformation induced by f . For this transformation to be well-defined, we assume f is invertible and differentiable. Now, the path integral can be computed by sampling from the new effective action Seff (V ) = S(f (V )) − log | det Jf (V )|. If Seff is constant, i.e., it corresponds to a uniform distribution, then f is called a trivializing map [50]. Naturally, the challenge in this approach is to find a general non-perturbative method to find a true trivializing map f . However, we stress that even an approximate trivializing map can effectively mitigate critical slowing down [51, 52]. The strategy of the simulation would follow the sketch in Figure 2: run a Monte Carlo using the new effective action and then flow back to the target distribution using f . If the approximate trivializing procedure is successful, smaller integrated autocorrelation times may be observed when running an MCMC. More precisely, this alternative sampling procedure may exhibit an improved scaling of the integrated autocorrelation time in the continuum, see Eq. (10); however, we also stress that the scaling of the computational cost of an MCMC update using the new action has to be monitored as well. In section 3.3 we review some machine learning architectures that can be used to build the map f . 3.2.4

Flow-based importance sampling

A different but similar approach can be followed if the map f is applied to field configurations V sampled from a known, tractable distribution q0 . In this case, it is possible to explicitly write the probability distribution that the transformed configurations U = f (V ) follow as q(U ) = q0 (f −1 (U )) |det J|−1 ,

(12)

and expectation values ⟨·⟩q can easily be obtained by building an ensemble of samples U as above. Going back to the path integral of Eq. (3), one can then insert it i h −S(U ) R DU q(U ) e q(U ) O(U ) h i (13) ⟨O⟩ = R −S(U ) DU q(U ) e q(U ) and use q(U ) for sampling. The factor w(U ) =

e−S(U ) q(U )

(14)

plays a crucial role, as its expectation value calculated with configurations sampled from q(U ) is the partition function of Eq. (4). A sketch of the structure of this sampling approach is shown in Figure 2. Provided that q and the target distribution p are similar enough, there are essentially two alternatives in order to compute exactly the expectation values in Eq. (3). • The first possibility is to insert w into the observable and compute the path integral with a reweighting procedure simply as ⟨O⟩ =

⟨O(U )w(U )⟩q . ⟨w(U )⟩q 11

(15)

Model density (q)

t=0

t = 0.5

t=1

1.0

0.5

t

0.0

Prior density (q0 )

Compute forces

Model density (q)

t=0

t = 0.5

t=1

1.0

0.5

t

0.0

Prior density (q0 )

Flow

Figure 2: Comparison of applying a trivializing flow to improve HMC mixing time (top panel) versus using the same trivializing flow to perform direct importance sampling (bottom panel). The target model density is selected as a simple toy model with two widely separated modes, for which traditional sampling methods suffer from freezing.

12

The variance of this estimator clearly depends on the degree of overlap between q(U ) and p(U ). This can be monitored by looking at the Effective Sample Size (ESS), which is defined as ⟨w(U )⟩2q (16) ESS = ⟨w(U )2 ⟩q and goes from 0 (vanishing overlap) to 1 (perfect overlap). It can be seen as the asymptotic ratio N/M relating the number of raw samples M from q(U ) that must be drawn to achieve precision equivalent to N independent samples from p(U ). Consequently, the average sampling cost of each independent sample in this case can then be (approximately) expressed as Cs(flow) ∝

1 . ESS

(17)

• The second possibility is to build a new Markov Chain using samples from q as proposals, with a Metropolis-Hastings acceptance probability   w(U ′ ) (i−1) ′ (18) A(U , U ) = min 1, w(U (i−1) ) to guarantee asymptotic convergence to the distribution p(U ). The sampling costs in this case still follow Eq. (10), but the integrated autocorrelation time now becomes purely a measure of the overlap of q(U ) with the target. Generally speaking, finding a suitable transformation f such that the overlap between q(U ) and the target p(U ) is sufficient represents a very difficult task. However, machine learning methods have begun to make promising inroads, as outlined in section 3.3.

3.3

Generative models

The alternative approaches discussed in the previous section that rely on a “flow” transformation require an essential ingredient, i.e., how can we move from one easy or trivial distribution to another which is more difficult to sample from? Generative models provide a powerful tool to do so: they are designed to generate samples from a target probability distribution, which in the case of lattice field theories is analytically known, generally starting from an analytically tractable prior distribution. The generative model induces on the target space a probability density qθ that approximates the target density p. The need for exactness typically requires that qθ must be analytically tractable, although other methods of unbiased sampling such as path weights are also feasible, as discussed below. Crucially, the scaling of the costs changes completely when compared to the one of a standard MCMC, see Eq. (7): (flow)

(flow)

Ctot (N ) = Ctrain + N Cs(flow) 13

(19)

(flow)

(flow)

with Ctrain being the fixed training costs and Cs given in the case of importance sampling by Eq. (17). A necessary condition to mitigate critical slowing down is that the sampling cost now scales in the continuum better than a traditional Monte Carlo. However, the costs of training an architecture that features such an improved scaling might be very large compared to the thermalization time or the sampling costs for systems not yet affected by severe autocorrelations, see the cartoon in Figure 1 for a simplified graphical comparison. We refer to [53] for a thorough discussion of the trade-off between training and sampling times and how to perform an honest comparison with MCMC costs. Furthermore, we mention the ability of some generative flow architectures to be generalizable across parameter space, i.e., being trained simultaneously for different values of the model parameters; see [54–57] for some examples. The remainder of this section explores the three broad strategies commonly used to define the flow transformation itself. The distinction between these three approaches is whether and how the flow transformation is constrained to follow specific intermediate distributions between the prior and target distribution: 1. General normalizing flow models follow no defined path. 2. Non-equilibrium approaches follow an explicitly defined path, specified by a changing action. 3. Diffusion models and related approaches follow an implicitly defined path, specified by a transformation of the field configuration samples. These approaches are schematically compared in Figure 3. 3.3.1

No defined distribution path: normalizing flows

Normalizing flows (NF) [58, 59] can be used as a parametrization of the maps discussed in section 3.2, both in the context of flow-based importance sampling and of (approximate) trivializing maps: here we will discuss applications in lattice field theory of both discrete NFs (for a variety of implementations see refs. [12, 53, 56, 60–80]), and continuous NFs (see refs. [55, 81–84]). A discrete normalizing flow is a diffeomorphism fθ built as a combination of invertible maps fi , usually denoted “coupling layers”. In the following, we collectively indicate with θ the set of parameters defining the transformation fθ . In the case of a lattice gauge theory simulation, it takes a field configuration U0 sampled from some prior distribution q0 and gives as output a new configuration U : U = fθ (U0 ) = (fn ◦ · · · ◦ f1 )(U0 ).

(20)

The training procedure of a normalizing flow requires that the Jacobian of each coupling layer fi be efficiently computable, so that the determinant of the Jacobian of the whole transformation fθ can be calculated as det Jfθ = det Jf1 · . . . · det Jfn . 14

(21)

q Discrete Flow

q0

S1

Target theory p = e −S /Z

≈

p

Continuous Flow

S2

S3

…

m riu b i l i qu ds n-e etho o N M

Diffusion Model

Trivial theory q0 = e −S0 /Z0

Figure 3: Illustrative comparison of the three transport strategies outlined in the main text. To achieve invertibility and a triangular Jacobian, the maps fi usually act separately on two subsets of the field configuration: one is “active”, i.e. the variables are modified by the map, while the other is kept “frozen” and used purely as an input for the map. As a simplified example in lattice gauge theory, we can take ( ′ = Wθi (Ufrozen ) · Uactive Uactive (22) fi : ′ = Ufrozen , Ufrozen where Wθi (·) is a group-valued function. This transformation is invertible regardless of the form of the function Wθi , which can be parametrized by arbitrarily deep and complex neural networks. Simple examples of coupling layers are NICE [85] and RealNVP [86], which can be readily used in a scalar field theory with an even-odd decomposition. For recent implementations of coupling layers enforcing gauge symmetries, please refer to [12, 42, 73] and to [65, 70] for gauge theories with fermions. Crucially, the probability distribution of the output of a discrete normalizing flow can be explicitly written using the change of variables theorem as qθ (U ) = q0 (fθ−1 (U )) |det Jfθ |−1 , −S(U )

(23)

so that the reweighting factors w(U ) = eqθ (U ) can be computed when evaluating the path integral of Eq. (13). 15

An alternative architecture to discrete coupling layers is represented by continuous normalizing flows [87–89]. In this case, the invertible map fθ is defined by integrating the ordinary differential equation (ODE) d U (t) = ∇φθ (U (t); t) (24) dt starting from a given configuration U0 and parameterized by the scalar function φθ . For lattice gauge theories, the gradient ∇ takes values in the Lie algebra, and the resulting field configuration at arbitrary times t may then be written as a path-ordered exponential, Z t  ′ ′ ′ U (t) = P exp dt ∇φθ (U (t ); t ) U0 . (25) 0

The output of the map is conventionally defined as the value at t = 1, i.e., fθ (U0 ) = U (t)|t=1 . Importantly, an analogous ODE can be used to compute the Jacobian determinant as Z 1 log det Jfθ = dt ∇2 φθ (U (t); t), (26) 0

giving qθ (U ) via Eq. (23) as for discrete flows. The choice of the scalar function φθ parametrizes the flow, and may for example be determined analytically or using optimized neural networks. The trivializing map of Ref. [50] is a particular parameter-free choice of a continuous flow, for which the flow ∇φ is evaluated in a perturbative expansion to a given order; unfortunately, the convergence of this expansion limits the range of lattice spacings for which this approach is practical. In order to obtain a flow fθ that can be used in a realistic simulation, a training procedure on the parameters θ must be performed: this is one of the main challenges of flow-based approaches along with the parametrization of fθ itself. For both discrete and continuous normalizing flows, training is done by minimizing certain metrics that measure the level of similarity between qθ (U ) and p(U ). The most used in the context of NF is the so-called (reverse) Kullback–Leibler (KL) divergence, defined as Z qθ (U ) = −⟨log w(U )⟩U ∼qθ + log Z ≥ 0. (27) DKL (qθ ∥p) = DU qθ (U ) log p(U ) This quantity is minimized when qθ = p, so training to minimize the KL divergence drives the model distribution qθ towards the desired Boltzmann distribution. In the last step it is made explicit that this is a case of a self-training procedure: indeed, this metric is computed using samples from the model itself and not from the target p. Note however that this is not the only possibility, as other loss functions can be evaluated in a self-training approach (and can be varied during training as well). Moreover, another typical choice is using the forward KL divergence and training using samples from an MCMC, or a mix of the two approaches. For lattice gauge theories, it has also been found to be helpful to constrain flow models to exactly satisfy certain symmetries, particularly including the high dimensional lattice gauge symmetry. The probability density q(U ) of a flow model is guaranteed to be invariant with respect to any symmetry generator g if [90] 16

1. The prior density is invariant, q0 [g · U ] = q0 [U ], and 2. The flow is equivariant, f ◦ g = g ◦ f 4 . Several classes of discrete and continuous flows have been developed to be equivariant under gauge and other symmetries [12, 55, 63]. Outside of normalizing flows, recent works have also developed and applied new architectures for more general gauge-equivariant transformations [13, 42, 49, 91–95]. An open question in this area is how to achieve a good tradeoff between flexibility of the flow fθ and the implementation of the constraint to be symmetric, and future improvements may emerge from integrating non-flow methods. 3.3.2

Explicitly defined distribution path: non-equilibrium approaches

A similar approach to importance sampling, but purely stochastic in nature, is given by algorithms inspired by Jarzynski’s equality [96, 97] and (equivalently) by Annealed Importance Sampling (AIS) [98]. In these approaches, one builds between the prior q0 ∼ e−S0 and the target p ∼ e−S a sequence of intermediate probability distributions that interpolate between them. The initial configuration U0 is sampled from q0 : then it is updated following a pre-determined sequence of intermediate distributions until reaching the target in a finite number of steps. In particular, we have something like Sλ(1)

Sλ(2)

Sλ(3)

U0 −→ U1 −→ U2 −→ . . .

Sλ(nstep ) ≡S

−→

U

(28)

where each arrow indicates a MCMC update performed using an intermediate distribution, itself defined by the action Sλ(n) , which depends on a (set of) parameter(s) λ(n) varied throughout the process. Since in general we never let the system thermalize at each step, this is a typical example of a non-equilibrium Monte Carlo simulation. Interestingly, such an approach can be used to compute the path integral of Eq. (13) as ! Z nY step ⟨Oe−W ⟩f , (29) ⟨O⟩ = DUi q0 (U0 ) Pf [U0 , U1 , . . . U ] e−(W−∆F ) O(U ) = −W ⟩ ⟨e f i=0 where Pf is the transition probability of the sequence of updates of Eq. (28) and W the work performed on the system when moving from one intermediate action to the next (with Sλ(0) := S0 ), nstep −1 X W= Sλ(i+1) (Ui ) − Sλ(i) (Ui ). (30) i=0

Furthermore, in Eq. (29) we also used Jarzynski’s equality, i.e. e−∆F = ⟨e−W ⟩f , with ∆F being the difference in free energy between q0 and p. Notably, e−W plays the role of the reweighting factor w defined in the previous sections. In the last decade this approach has 4

The term “equivariant” is a property of maps, such as the flow f , whereas “covariant” is typically reserved for fields with well-defined transformation properties under the action of generators g.

17

found applications in lattice gauge theory, both for the computation of partition functions [99, 100] and as a tool to mitigate critical slowing down [101]. This approach can be augmented with the same tools used to define normalizing flows. Namely, one can redefine the sequence of Eq. (28) as Sλ(1)

f1

Sλ(2)

f2

U0 −→ f1 (U0 ) −→ U1 −→ f2 (U1 ) −→ . . .

Sλ(nstep ) ≡S

−→

U

(31)

with the fi being for example the invertible transformations defining the NF in Eq. (20): this defines a stochastic normalizing flow [102] which has been extensively implemented in lattice field theory [103–106]; see also Ref. [107] for related architectures. Note that the Continual Repeated Annealed Flow Transport Monte Carlo (CRAFT) [108] approach follows a very similar philosophy. The work is straightforwardly generalized to W=

nstep −1

X

Sλ(i+1) (fi+1 (Ui )) − Sλ(i) (Ui ) − log det Jfi+1 (Ui ).

(32)

i=0

Unlike NFs, the probability distribution q(U ) of the flow is not analytically accessible and the divergence of Eq. (27) cannot be used for training the parameters. Interestingly, a different divergence can be used to train such architectures, namely DKL (q0 (U0 )Pf ∥p(U )Pr ) = ⟨log

q0 (U0 )Pf [U0 , U1 , . . . U ] ⟩f = ⟨W⟩f − ∆F ≥ 0; p(U )Pr [U0 , U1 , . . . U ]

(33)

here Pr indicates the overall transition probability of a sequence of updates similar to Eq. (28), but starting from p and reaching q0 . So, the KL divergence of Eq. (33) measures the similarity between the stochastic flow moving forward (from q0 → p) and in reverse (from p → q0 ): as for NFs this can be computed using only samples from the forward process. The positivity of the KL divergence also implies the last step of Eq. (33), which has a clear interpretation in terms of the Second Law of thermodynamics. Another promising approach can be obtained starting from AIS and non-equilibrium MCMC and taking the limit in continuous time using annealed Langevin dynamics, i.e., solving the Stochastic Differential Equation (SDE) √ dϕt = −ϵt ∇Sλ(t) (ϕt )dt + 2ϵt dηt , ϕ0 ∼ ρ0 . (34) Here ϕt is the field configuration, Sλ(t) is the action at time t following a certain protocol for the time-dependent parameter(s) λ(t), ϵt is a time-dependent diffusion coefficient, and ηt represents a unit-normalized Wiener noise process. Following Eq. (30), the infinitesimal work can then be naturally defined as dWt = ∂t Sλ(t) (ϕt )dt;

(35)

expectation values and the partition function can then be computed as in Eq. (29). Exactly as the discretized version, such samplers can suffer from high variance when too far from the equilibrium distribution throughout the transformation. To solve or mitigate 18

this problem, the idea proposed in approaches such as NETS [109] (see also [110, 111]) is to introduce an additional transport term. The process now becomes √ dϕt = −ϵt ∇Sλ(t) (ϕt ) dt + bt (ϕt ) dt + 2ϵt dηt , ϕ0 ∼ ρ0 , (36) dWt = −∇ · bt (ϕt ) dt + ∇Sλ(t) (ϕt ) · bt (ϕt ) dt + ∂t Sλ(t) (ϕt ) dt (37) where bt is a velocity field that can be learned. In the perfect transport case, this new field would satisfy: ∇ · (bt ρt ) = −∂t ρt (38) −1 −Sλ(t) with ρt = Zλ(t) e . It is easy to prove that the weights of Eq. (37) would simply become

dWt = ∂t Fλ(t) dt

(39)

with Fλ(t) = − log Zλ(t) . 3.3.3

Implicitly defined distribution path: diffusion models

Diffusion models provide an alternative framework for generating samples from an empirical probability distribution represented by a dataset. Such models can be essentially described as two distinct processes: in the forward process, noise is incrementally added to the data, while the reverse process (the so-called denoising) tries to remove the added noise and recover the (unknown) probability distribution of the starting dataset. The forward process implicitly defines a path through the space of distributions pt that is designed to bridge the target distribution p0 with a simple Gaussian distribution pT . The reverse process is then optimized to approximate this bridge. The forward process is described by a Langevin process in the time t = 0 → T dϕt = v(ϕ, t)dt + g(t)dηt

(40)

with drift term v(ϕ, t) and time-dependent scalar diffusion coefficient given by the square of g(t). The reverse process, on the other hand, can be defined via the reversed time variable t̄ = (T − t) running from 0 → T and it is given by the SDE dϕt̄ = −v(ϕ, t)dt̄ + g(t)2 ∇ log pt dt̄ + g(t)dηt̄

(41)

with the drift term now including pt (ϕ), the probability distribution at time t in the forward diffusion process. Such SDE can be solved after specifying v and g and, most importantly, determining the gradient of the logarithm of pt (ϕ). If this is possible, then Eq. (41) can be used to generate new samples directly from the target distribution the initial data were sampled from. Typically, neither pt (ϕ) nor the required score term ∇ log pt (ϕ) are known. The key for a successful denoising procedure lies then in learning an approximate score function sθ (ϕ, t) ≈ ∇ log pt (ϕ), which can for example be modeled by deep neural networks. This score matching procedure is performed by minimizing a proper loss function [112]. Since the 19

learned score is only approximate, the generated samples are not exact by themselves. We can restore exactness by using Metropolis-adjusted annealed Langevin steps (MAALA) [57], by tracking path weights along the generative process [107,113], or by using the trained model as an energy-based model that defines the intermediate distributions of annealed sampling methods [114,115]. Note that, in contrast to other methods presented in this section, diffusion models require data sampled with the target distribution in order to actually perform score matching. For more details on the applications to lattice gauge theories and the connection with Stochastic Quantization see [57, 116–123]. More general families of interpolation between the target distribution and a simple distribution have also been proposed in works on stochastic interpolants [124,125] and the flow matching procedure [126]. As for diffusions, the implicit paths here are also defined using samples from the target distribution, and the map is designed in a way that provides an efficient training routine for the relevant score function or velocity field.

3.4

Estimating the partition function

The use of trivializing maps for sampling has a pleasant side-effect, as it becomes straightforward to evaluate the partition function itself, which is connected to many interesting observables such as the pressure for the computation of the equation of state [127], the free energy of interfaces [128, 129], ’t Hooft loops [130] and the entanglement entropy [131]. Evaluating Z in an ordinary Monte Carlo is cumbersome and unnatural, precisely because the partition function is not written as the expectation value ⟨O⟩ of any fixed function O; typical determinations of (ratios of) partition functions involve calculating the derivative of log Z with respect to some parameter, and then integrating over the parameter itself. On the other hand, the availability of an (approximate) flow f makes evaluating the partition function straightforward. Recall that an approximate trivializing map for an action S(U ) is defined by the fact that Z −1 e−S(f (U0 )) | det J| ≈ q0 (U0 ), (42) where q0 is a “trivial” distribution whose properties are analytically understood. In particular, the correctly normalized probability distribution q0 (U0 ) can be evaluated efficiently for any U0 . If the trivializing map were exact, we could extract the partition function Z immediately, simply by evaluating both q0 and S ◦ f at the same point U0 : Z=

e−S(f (U0 )) . q0 (U0 )| det J|−1

(43)

For trivializing maps which are only approximate [61, 62], we can still extract the partition function Z, but now as an expectation value:   e−S(U ) Z = ⟨w(U )⟩U ∼q = . (44) q0 (f −1 (U ))| det J|−1 U ∼q 20

Note that in the case where only the unnormalized prior weight e−S0 is known, this formula can still be used to extract the ratio Z/Z0 , which is often the interesting quantity to study. In the case of AIS/SNF approaches, one can use Jarzynski’s equality in the same way, with the weight now given by the exponential of work of eqs. (30) and (32), i.e. Z = ⟨e−W ⟩f Z0

(45)

where the partition function Z0 of the distribution q0 has been made explicit for consistent notation. Flow and flow-related architectures have been applied successfully on a number of phenomenologically interesting studies of partition functions and free energies. An incomplete list includes computations of the equation of state [62, 99, 100], numerical studies of effective string theories [83,104], entanglement entropy and c-functions [132–134] and the critical Casimir effect [135].

4

Observable estimators

Instead of attempting to gather samples more efficiently, we may attempt to reduce the noise in the observables, which reduces the number of samples that are needed. This section concerns such methods: • Control variates learn a total derivative which is correlated with the noise in the observable, and then subtract that noise off explicitly. • Methods which enable the generation of correlated samples, e.g., for different action parameters, reduce the noise in observables related to variations of those parameters. • By evaluating expectation values on a deformed integration contour, we may modify the signal-to-noise ratio without affecting the physics. One particular “observable” has a different form from all others: the partition function (or Helmholtz free energy) itself. Most of the methods discussed in this section assume that the observable is an expectation value, evaluated over a probability distribution given by the Boltzmann factor e−S . This is not the case for the partition function: nevertheless there are ways to evaluate it, and a particularly neat one that uses machine learning methods— specifically normalizing flows—is discussed in section 3.4. Some of the ideas discussed in this section—control variates and deformed contours— operate directly on the level where expectation values are computed. One function of field configurations is replaced by another function which, while related, has more desirable statistical properties. The observables we are most interested in estimating, however, are often not simple expectation values, but are instead to be extracted by performing fits to those expectation values. A common example is the mass, extracted as the asymptotic decay constant of a two-point function: ⟨O(τ )O(0)⟩ ∼ e−M τ . 21

(46)

Methods like control variates will seek to improve the estimate of the mass by replacing the product O(τ )O(0) with an operator which has exactly the same expectation value, but smaller variance. We might also, however, imagine replacing it with an operator that shares the same quantum numbers, but has greater overlap with the state being studied. Recent work has begun to explore this direction, using gauge-equivariant neural networks to parametrize interpolators for static quark-antiquark systems, optimized for overlap with the ground and excited states [136, 137].

4.1

Motivation

Consider computing, in lattice QCD, the nucleon two-point function: CN (τ ) = ⟨N (τ )N̄ (0)⟩ ∼ e−MN τ .

(47)

For the purposes of estimating the nucleon mass, we want to evaluate CN (τ ) for τ as large as possible (to damp away excited state contributions). How large τ we can evaluate at depends (given a sufficiently large lattice) on how the statistical errors behave as a function of τ . The signal —the value of the correlation function itself—is of course known to decay exponentially with decay constant MN . By inspecting the variance of the propagators that form the two-point function N (τ )N̄ (0), we may see that it does as well: ⟨|N (τ )N̄ (0)|2 ⟩ = ⟨(q q̄)3 (τ )(q q̄)3 (0)⟩ ∼ e−3Mπ τ

(48)

The variance of the nucleon two-point function is given by the two-point function of three 3 pions, so the noise (the square root of the variance) decays as e− 2 Mπ τ . The precise value of this noise, like the precise value of the nucleon two-point function, would require a full Monte Carlo to estimate on any given lattice; however we know its asymptotic behavior by inspection. Although both the signal and the noise decay exponentially, the decay constant for the noise is substantially smaller: 32 Mπ < MN . As a result, the relative statistical errors grow (exponentially) at large τ . This situation is shown in Figure 4. This Parisi-Lepage argument [138, 139] shows that calculations of the nucleon two-point function suffer in general from a signal-to-noise problem which is exponential in the separation of the source and sink5 , regardless of which operators are used to extract the nucleon state. Analogous arguments indicate that the signal-to-noise problem occurs for almost all correlation functions considered in lattice studies, including higher n-point functions. A closely related problem, the extensive sign problem [140–145], occurs when performing Monte Carlo on a theory with a complex action, such as systems with a non-zero density of relativistic fermions, or real-time dynamics6 . In such systems the path-integral formula for 5

One should not confuse this exponential scaling in separation with any sort of statement about asymptotic computational complexity. The task at hand is to determine MN to within precision ϵ, and this is not exponentially difficult in any parameter of the theory. 6 In fact sign problems are much more general than this. Almost all quantum Hamiltonians possess sign problems that persist under an arbitrary local change of basis. It is fortunate that many of the systems most interesting in particle and nuclear physics are of the special sort which do not.

22

0.35 0.30

Effective mass

Correlator

100 10−2 10−4 10−6

0.25 0.20 0.15 0.10 0.05 0.00

0

20

40

60

80

0

τ

20

40

60

80

τ

Figure 4: Two views of the signal-to-noise problem, on synthetic data. Left: the correlator on a log-linear plot. The statistical errors grow relative to the signal, although both fall exponentially. Right: an effective mass plot. The smallest mass cannot be effectively determined at early times due to excited-state contamination, and it cannot be determined at late times due to the falling signal-to-noise. the expectation value DU e−S(U ) O(U ) R (49) DU e−S(U ) cannot be interpreted as an expectation value over a probability distribution. To bring Monte Carlo methods to bear on this problem, we most commonly re-express the quantum expectation value as a ratio of classical expectation values taken with respect to the phasequenched system, which is defined by the real part of the action. Concretely, we write Z ⟨Oe−i Im S(U ) ⟩Re S 1 ⟨O⟩ = , ⟨X⟩Re S ≡ DU X(U ) e− Re S(U ) . (50) ⟨e−i Im S(U ) ⟩Re S ZP Q R

⟨O⟩ =

A little consideration of the properties of the partition functions Z and ZP Q is sufficient to convince yourself that their ratio scales exponentially with the (space-time) volume of the system, and therefore that the denominator of the above expression possesses an exponential signal-to-noise problem, similar to those which afflict correlation functions, but far greater in magnitude. In most cases, the signal-to-noise problem for correlation functions is also due to a sign problem. This follows directly from the Parisi-Lepage argument, which indicates that the average magnitude of measurements of a correlation function is exponentially larger than the mean value to be extracted, which requires strong cancellation of positive and negative sign contributions. This theoretical argument has been complemented by empirical evidence for correlation functions in a variety of contexts [146–150]. In this text we will primarily focus on these non-extensive signal-to-noise/sign problems, arising in systems with a real action. There are two related reasons for this. First, the extensive sign problem has exponential scaling in the space-time volume of the lattice7 . For 7

Fermion sign problems scale with the volume measured in physical units (and do not worsen in the

23

QCD-relevant lattice sizes, we are many orders of magnitude away from making meaningful progress on e.g. QCD at T = 0 and large chemical potential. By contrast, signal-to-noise problems have exponential scaling with a length, and are therefore tractable in the near term. Second, the most common computational difficulties encountered in lattice QCD are ordinary signal-to-noise problems. (These are related, of course, because practitioners prefer to work in regimes where we are only a small distance from “success” however defined.)

4.2

Control variates

We begin with what may be the conceptually simplest approach to a signal-to-noise problem; at least, after getting more statistics. When we wish to evaluate ⟨O(U )⟩, we will replace O(U ) with another function Õ(U ), which may be proven to have precisely the same expectation value, but has smaller variance. We accomplish this by subtracting from O a function f , known to have zero expectation. This function is known as a control variate [151]. The key observation is that the integral over a compact8 space, of a total derivative, vanishes. That is, for any function g we have   Z Z
∂g(U ) ∂S(U ) ∂ −S(U ) −S(U ) − g(U ) . (51) g(U )e = DU e 0 = DU ∂U ∂U ∂U This is also true for any measure-preserving derivative operator. In what follows we will simply write ∂ to stand in for any chosen derivative operator, until we come to the point where the precise operator matters. Moreover, P every (sufficiently regular) function f such that ⟨f ⟩ = 0 can be written in the
form e−S f = a ∂a ga e−S , summing over a complete set of such derivative operators ∂a . This is vital: it will mean that nothing is lost by restricting our search for control variates to the space of total derivatives. Note that this decomposition is not unique when there is more than one derivative operator ∂a included in the sum. With this in mind, we can replace our operator O via the equality ⟨O⟩ = ⟨O − ∂g + g∂S⟩,

(52)

for some derivative operator ∂. With this in hand, we are left entirely free to choose g so as to minimize the variance of the improved observable O − ∂g + g∂S. How do we find a favorable function g? An obvious approach is to optimize a neural network [152]. We want to minimize the variance: V [g] = ⟨(O − ∂g + g∂S)2 ⟩ − ⟨O⟩2 ,

(53)

so we can use V [g] as a loss function. continuum limit); real-time sign problems scale with the number of lattice sites (and therefore exhibit exponential scaling in both the continuum and infinite-volume limits). 8 In fact the method does not depend on having a compact configuration space, as long as we’re careful that all functions decay sufficiently quickly at large field values—exponential decay is always sufficient.

24

In practical lattice calculations, the number of configurations available is likely not large relative to the number of parameters in the neural network. If we disregard this and minimize V [g], as defined above, on the ∼ 103 configurations available, we will find that the neural network overfits: it learns a function such that ∂g −g∂S exactly replicates the fluctuations of O on the samples available, but has uncontrolled behavior elsewhere and does not effectively reduce the signal-to-noise problem. Another approach is to write down a basis of candidate control variates, and optimize a linear combination of these candidates. Given a basis of functions G1 , . . . , GK , we define the candidate control variates in the same way as above, according to Fi = ∂Gi − Gi ∂S.

(54)

Then the most general linear combination is defined by K coefficients c1 , . . . , cK : f = ci Fi .

(55)

The chief virtue of this approach is that the minimization of the variance is now a matter of linear algebra. We denote the covariance among the candidates, and the correlations with the observable, as follows: Mij = ⟨Fi Fj ⟩ and bi = ⟨Fi O⟩. (56) With these definitions, the optimal control variate (in this basis) is given by the coefficients c = M −1 b.

(57)

Once again, we have a problem when we consider applying this method to a lattice of a size typical in lattice QCD. Optimistically, we might have 104 configurations, on a lattice of size ∼ 504 . The number of degrees of freedom in the linear ansatz above is far greater than the number of data points on which we can train. Regardless of the family of control variates used, it is valuable to ensure that these control variates respect the same symmetries as does the lattice action. To be clear, this is not necessary for correctness—a control variate constructed as described above will have vanishing expectation regardless of its behavior under the various lattice symmetries. However, ensuring that all appropriate symmetries are respected is a cheap way to alleviate overfitting problems. Control variates have been successfully generalized to both fermionic [153] and discrete [154] systems. The chief obstacle to their application to full-scale lattice QCD is the overfitting problem: it is easy to write down a large family of control variates, but difficult to know which ones are effective without an enormous quantity of training data. Note also that recent advancements [155–157] in the use of flow-based approaches as a variance reduction technique naturally lead to a definition of a control variate.

4.3

Contour deformations

An alternative approach to the signal-to-noise problem starts from complex analyticity9 of the lattice path integrals. Typical lattice path integrals can be written as (high-dimensional) 9

Or, “holomorphy.”

25

iterated integrals over real coordinates that can be analytically continued to complex values of the coordinates. When this is the case, Cauchy’s integral theorem allows the path integral to be deformed from integrating over the original contour M0 to a deformed integration contour M, Z Z 1 1 −S(U ) ⟨O⟩ = DU O(U )e = DŨ O(Ũ )e−S(Ũ ) , (58) Z M0 Z M as illustrated in Figure 5. What properties must M satisfy in order for equation (58) to hold? The answer comes from Stokes’ theorem, of which equation (58) is a special case. The new contour of integration gives an equivalent result whenever the oriented combination of M0 − M bounds a region where the integrand is holomorphic. In practice, it is sufficient (but not necessary) for M0 and M to be homotopic. That is, it is sufficient for there to exist a connected family of manifolds Mt , defined on t ∈ [0, 1], such that the integrand is holomorphic all along each manifold. A practical approach to parameterizing valid integration contours M, is then to define the new contour as the image of a one-to-one map fθ : M0 → M. It is sufficient for the map to be connected continuously to the identity map through variations of the parameters θ, as long as this does not traverse any non-analyticities of the path integral. In the context of lattice QCD, the initial integration contour is usually M0 = SU (3)4V , where SU (3) is viewed as an 8-real-dimensional submanifold of the 8-complex-dimensional space SL(3; C). Deformed contours are then submanifolds of SL(3; C)4V . Deformed integration contours may be parameterized either analytically, using neural networks, or using hybrid approaches. A tradeoff in the parameterization is that the complex Jacobian determinant J(U ) = det(∂ Ũ /∂U ) is required to integrate on the deformed contour. There is a clear analogy with normalizing flow transformations, and complex normalizing flows have been considered to define deformed contours [158]. However, unlike the case of real normalizing flows, composition of discrete transformations only gives a valid deformation with a Jacobian of known factorization if each layer is holomorphic, which strongly restricts the available architectures. Hybrid approaches include restricted classes of deformations with constant or diagonal Jacobians [159–166] which can then be optimized using stochastic gradient descent methods. In some cases, determining the parameters of such deformations as the output of neural networks applied to a constant input significantly improved the stability of optimization and the final results [167]. 4.3.1

Extensive sign problems

This approach has historically been studied as a possible avenue to circumvent the extensive sign problem [158–160, 168–189] (see Ref. [190] for a review). Here, the challenge is a complex-valued action S(U ) which affects both the estimation of the partition function and expectation values of operators. In both cases, Monte Carlo importance sampling according to the absolute value measure |e−S(U ) | requires taking expectation values of severely

26

fluctuating complex phase factors e−i Im S(U ) , for example Z Z 1 = DU |e−S(U ) |e−i Im S(U ) , (59) ZP Q ZP Q R where ZP Q := DU |e−S(U ) | is the normalization of the Monte Carlo probability distribution, often referred to as the “phase quenched” partition function. This integral is difficult to estimate by Monte Carlo precisely due to these phase fluctuations. We are computing the expectation value of e−i Im S over the phase-quenched distribution. Each sample has unit norm; however, the final ratio Z/ZP Q is quite small, and in fact exponentially small in the total space-time volume of the system. We typically refer to 4 Z/ZP Q as the “average phase” ⟨σ⟩ ∼ e−m V , where m is some mass scale characterizing the severity of this sign problem. Integrating on a deformed contour has the potential to mitigate this sign problem. The idea is simple: observables are holomorphic functions of the fields, and are therefore not modified (in expectation) by deforming the contour. The average phase, however, is defined by integrands which are not holomorphic, and therefore can be changed by the choice of contour. We simply wish to find the contour that reduces the sign problem as much as possible. (Although there are still open questions, a great deal is known today about the nature of such contours; see [191] for recent results.) In practice, integrating on a new contour M requires a coordinate system to evaluate the integral. It is practical to simultaneously define the contour and coordinate system by applying a diffeomorphic map that sends each original field configuration U to a complexified configuration Ũ . Leaving Ũ as an implicit function of U , the deformed evaluation of the partition function or an observable can then be written as Z Z 1 −S(Ũ ) DU J(U )O(Ũ ) e−S(Ũ ) , (60) Z = DU J(U ) e , ⟨O⟩ = Z where J(U ) is the (complex) Jacobian of the map from the original integration domain to the deformed contour. How can we select a contour with a reduced sign problem? We can parameterize the map Ũ (U ) by a neural network, but what should be the loss function? We will use as our loss function the negative of the average sign ⟨σ⟩ = Z/ZP Q , which may be defined on any contour. An obstacle immediately presents itself: this average sign is difficult to evaluate precisely, because of the sign problem. A close inspection of the method of gradient descent (and related algorithms) reveals that we do not need to evaluate the loss function at any point during the optimization. Instead, we only want to evaluate the gradient of the loss function with respect to parameters. While the evaluation of the loss function −⟨σ⟩ suffers from precisely the sign problem we are attempting to solve, derivatives of that loss function do not: δ Z [Ũ ] δ δ Z δ Ũ (U ) P Q ⟨σ⟩[Ũ (U )] = = −⟨σ⟩ . δ Ũ (U ) δ Ũ (U ) ZP Q ZP Q [Ũ ]

27

(61)

γ − ≈ Figure 5: A cartoon of a valid contour deformation. The original domain of integration is the real line, which is then deformed to the curve γ. Together, the two form the boundary of a region Ω (symbolically, ∂Ω = R − γ) on which the function being integrated is holomorphic. Here we have used the fact that the partition function Z is contour-independent, and of course we happily disregard an overall factor of ⟨σ⟩, which does not change the direction of gradient descent. Only the phase-quenched partition function ZP Q changes as the contour changes, and derivatives of the phase-quenched partition function are (essentially by definition) sign problem-free. This approach has been used with a good deal of success on low-dimensional fermionic systems [160,161]. More broadly, neural networks have been used to parameterize contours that reduce the sign problem in a range of systems, including finitedensity, Hubbard-model, and real-time problems [91,158,164,165,179,180,182,185,192], and normalizing flows have also been used to compute the density of states for complex actions [193]. 4.3.2

Deformed observables

Deformations can also be applied in the case that the action is real and the observable causes the relevant sign and signal-to-noise problem. In this context, the aim is to integrate on a contour M on which the signal-to-noise problem is milder. It is also possible in this case to evaluate the deformed path integral by reweighting with respect to the original Monte Carlo samples, because the partition function Z itself has no sign problem. This can be done by defining the new observable Q(U ) ≡ J(U ) O(Ũ ) e−S(Ũ )+S(U ) , which has an unchanged expectation value, Z n o 1 DU J(U )O(Ũ )e−S(Ũ )+S(U ) e−S(U ) ⟨Q⟩ = Z Z 1 DU J(U )O(Ũ )e−S(Ũ ) = ⟨O⟩ . = Z 28

(62)

(63)

The variance of Q is in general different from that of O, and the objective is to optimize the manifold parameterization to minimize the variance of the deformed observable. Typically, one is interested in observables that are on average real, in which case the variance Var[Re Q] = ⟨(Re Q)2 ⟩ − (Re ⟨Q⟩)2 determines the effectiveness of this approach. The second term is unaffected by contour deformations, so a useful choice of loss function to minimize is L ≡ (Re Q)2 .

(64)

Both this loss function and its gradient are expressed as an expectation value with respect to the original Monte Carlo distribution, allowing the parameterization to be optimized, and ultimately evaluated, without generating new Monte Carlo samples. This has been evaluated in a variety of contexts, including complex scalar field theory, Abelian gauge theory, and SU(N) gauge theories in 1+1 and 2+1 dimensions [162, 163, 167].

4.4

Correlated samples

For the special case of observables defined in terms of derivatives or differences of primary observables with respect to action parameters, another avenue towards noise reduction lies in the generation of correlated samples [74]. Specifically, if two ensembles of gauge-field configurations are generated at nearby action parameters, α1 and α2 , the derivative or finite difference ⟨O⟩α1 − ⟨O⟩α2 d⟨O⟩ ≈ , (65) dα α1 − α2 benefits from potentially substantial cancellation of statistical fluctuations between the two terms when the ensembles are correlated. For highly correlated ensembles, the variance of the estimator may be dramatically reduced compared to both independent sampling at each parameter value and standard finite-difference reweighting methods. In the limit of a small difference, this has been exactly related to the calculation of control variates of the form of Eq. (52) [155–157]. The generation of correlated ensembles can be achieved using the same flow model architectures as described in Sec. 3 to map between theories at nearby points in parameter space. Physically meaningful applications that involve derivatives with respect to action parameters include chiral and continuum extrapolations, as well as calculations of matrix elements like the chiral condensate, the nucleon sigma term, and other observables using the Feynman–Hellmann theorem. Additional applications arise in evaluating derivatives with respect to the electromagnetic coupling for scale setting or determining isospin-breaking effects in QCD+QED frameworks [194]. There is also relevance for theories exhibiting a sign problem, such as computing derivatives with respect to the baryon chemical potential or the QCD θ-angle. In such cases, the distributions connected by the flow transformations are generally more similar than those encountered in sampling-acceleration contexts, and current flow-based techniques are already viable for typical lattice QCD scales.

29

5

Observable acceleration

This section is concerned with evaluating a fixed observable as quickly as possible. This is most useful in the case where we have a large number of gauge configurations already gathered, but the observable evaluation is sufficiently slow that we will not be able to evaluate the observable on all gauge configurations, without exceeding the available computational resources. There are many ways of optimizing the algorithms used to evaluate various observables, and each method inevitably has some free parameters which can be ‘learned’. To keep this section manageable, we limit the discussion to two particular methods. The first, that of surrogate observables, takes a straightforward ML approach, by searching for a neural network (or other model) which approximates the expensive observable. The second, preconditioners, is already in wide use in its conventional (non-ML) form.

5.1

Surrogate observables

In the case where we have some observable, a function O(U ) of the gauge configurations, which is slow to evaluate, we may accelerate things by finding another function Õ(U ) which is approximately equal, but can be evaluated more quickly. Obviously, if we have exact equality O(U ) = Õ(U ) for all U , we can simply replace the expectation value of the original observable O by its faster equivalent. Nothing is lost in such a replacement. When the equality is not exact, but rather |O − Õ| is merely “typically small”, we can still gain, but more care is required. If we simply replaced O → Õ, we would gain speed but introduce some unknown error. We can handle this error by noticing: ⟨O⟩ = ⟨Õ⟩ − ⟨Õ − O⟩.

(66)

The first term on the right is the fast surrogate for ⟨O⟩; the second term corrects for the bias introduced. If we have N configurations, and both right-hand expectation values are evaluated on each of those N configurations, then there is no speed-up from this method: O still needs to be evaluated N times. The advantage to (66) comes from realizing that there is no need to evaluate both expectation values on all N configurations. If Õ is indeed a good approximation to O, then the square error |O − Õ|2 is small on average, and certainly smaller than the variance of O itself10 : ⟨O2 ⟩ − ⟨O⟩2 ≫ ⟨(Õ − O)2 ⟩ − ⟨Õ − O⟩2 .

(67)

Therefore, if we are evaluating ⟨Õ⟩ on N samples, the number of samples we need to reduce the error in bias-correcting term to be approximately that of the surrogate estimate is Ñ ∼ N

Var (O − Õ) . Var Õ

10

(68)

Note that it is trivial to find a surrogate where the variance in the bias correction is equal to the variance in O: Õ(U ) = 0 will do the job.

30

In the extreme case where the computational resources required to evaluate Õ(U ) are entirely negligible (often a reasonable approximation), the speed-up associated with the use of a surrogate is therefore a factor of Var (O)/Var (O − Õ). This is reminiscent of control variates, but instead of replacing our observable with one of lower variance, we are replacing it with one which may be more efficiently evaluated. Such a substitution is not designed to reduce our sampling costs, but can reduce computational costs associated with post-sampling measurement. The reader will hopefully not be surprised at this point by the suggestion that Õ(U ) may be learned, rather than guessed by a human. The equations above suggest that the relevant figure of merit for a surrogate observable is C[Õ(U )] = Var (Õ − O).

(69)

When the function Õ(U ) is smoothly parameterized as in a neural network, this figure of merit is differentiable and therefore usable as a loss function in gradient descent. This approach has been applied in full-scale lattice QCD calculations only a few times, and with modest success [195, 196]. Typically the improvement corresponds to saving ∼ 20% of the total computational effort—substantial, but in practice insufficient to justify the corresponding human effort in using these techniques.

5.2

Preconditioners

Another promising avenue for accelerating the computation of observables in LQFT involves employing ML techniques to develop improved preconditioners for propagator computations, namely, solutions to the Dirac equation. Preconditioners are transformations applied to linear systems to enhance their numerical properties, facilitating faster convergence of iterative solvers. Formally, given the linear system: Dψ = η,

(70)

a preconditioner is introduced to form a transformed system M Dψ = M η,

(71)

which ideally reduces the condition number and thereby accelerates solver convergence. In lattice field theories, the fields ψ and η represent discretized fermionic fields. Conventionally, iterative solvers such as the conjugate gradient (CG) method or BiCGStab [197] are complemented by analytically or heuristically derived preconditioners. Important examples include even-odd (red-black) preconditioning [198], domain decomposition methods [199], incomplete-factorization approaches such as SSOR preconditioning [200], and adaptive multigrid methods (including the closely related low-mode deflation) [201–205], which are now the standard choice for light-quark propagator computations. These traditional methods exploit the underlying algebraic and geometric structures of the lattice discretization to substantially enhance solver performance and computational efficiency. 31

Recently, there has been a growing interest in leveraging ML to create learned preconditioners, which might surpass traditional heuristics through data-driven optimization. While this field remains nascent, early results demonstrate significant potential. Studies have demonstrated that gauge-equivariant neural networks (GENNs) can successfully function as learned multigrid preconditioners [92, 94, 206]. By explicitly incorporating gauge symmetry—ensuring covariance under transformations Uµ (x) → Ω(x)Uµ (x)Ω(x + µ̂)† —these GENNs respect fundamental symmetries inherent in lattice QCD simulations. Further proofof-principle research has shown that neural-network-based preconditioners trained on small lattice volumes exhibit robust transferability to larger volumes, suggesting promising scalability and potential computational savings compared to classical methods such as even-odd or incomplete Cholesky preconditioners [207, 208]. Despite this encouraging progress, several key challenges remain. Ensuring that MLdriven approaches scale efficiently to larger lattices, maintain strict gauge covariance, and generalize effectively across diverse gauge configurations and lattice parameters are ongoing research priorities. Future research directions will likely continue exploring hybrid methodologies, integrating the strengths of classical numerical techniques with innovative ML-driven insights to push forward computational efficiencies in lattice QCD calculations.

Acknowledgments G.K. was supported by the European Union (ERC, FoundLatt, 101212401). Views and opinions expressed are however those of the authors only and do not necessarily reflect those of the European Union or the European Research Council. Neither the European Union nor the granting authority can be held responsible for them. S.L. was supported by a Richard P. Feynman fellowship from the Los Alamos National Laboratory (LANL) LDRD program; LANL is operated by Triad National Security for the National Nuclear Security Administration of the U.S. Department of Energy (Contract No. 89233218CNA000001). This manuscript has been approved for unlimited release and assigned the identifier LAUR-26-29844. A.N. acknowledges support by the Simons Foundation grant 994300 (Simons Collaboration on Confinement and QCD Strings), by the European Union – Next Generation EU, Mission 4 Component 1, CUP D53D23002970006, under the Italian PRIN “Progetti di Ricerca di Rilevante Interesse Nazionale – Bando 2022” prot. 2022ZTPK4E and by the SFT Scientific Initiative of INFN. PES is supported in part by the U.S. Department of Energy, Office of Science, Office of Nuclear Physics, under grant Contract Number DESC0011090 and by Early Career Award DE-SC0021006, and by Simons Foundation grant 994314 (Simons Collaboration on Confinement and QCD Strings), and has benefited from the QGT Topical Collaboration DE-SC0023646. This work is supported by the U.S. National Science Foundation under Cooperative Agreement PHY-2019786 (The NSF AI Institute for Artificial Intelligence and Fundamental Interactions, http://iaifi.org/).

32

References [1] P. Mehta, M. Bukov, C.-H. Wang, A.G.R. Day, C. Richardson, C.K. Fisher et al., A high-bias, low-variance introduction to Machine Learning for physicists, Phys. Rept. 810 (2019) 1 [1803.08823]. [2] S. Raschka, Y.H. Liu and V. Mirjalili, Machine Learning with PyTorch and Scikit-Learn: Develop machine learning and deep learning models with Python, Packt Publishing, Birmingham, UK, 1st ed. (2022). [3] F. Marquardt, “machine-learning-for-physicists: Code for “machine learning for physicists” lecture series.” https://github.com/FlorianMarquardt/machine-learning-for-physicists, 2025. [4] D. Boyda et al., Applications of Machine Learning to Lattice Quantum Field Theory, in Snowmass 2021, 2, 2022 [2202.05838]. [5] M. Abadi, A. Agarwal, P. Barham, E. Brevdo, Z. Chen, C. Citro et al., TensorFlow: Large-scale machine learning on heterogeneous systems, 2015. [6] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan et al., Pytorch: An imperative style, high-performance deep learning library, 2019. [7] J. Bradbury, R. Frostig, P. Hawkins, M.J. Johnson, C. Leary, D. Maclaurin et al., JAX: composable transformations of Python+NumPy programs, 2018. [8] SciDAC, LHPC, UKQCD collaboration, The Chroma software system for lattice QCD, Nucl. Phys. B Proc. Suppl. 140 (2005) 832 [hep-lat/0409003]. [9] RBC, UKQCD collaboration, Overview of Columbia Physics System, PoS LATTICE2013 (2014) 417. [10] QUDA collaboration, Solving Lattice QCD systems of equations using mixed precision solvers on GPUs, Comput. Phys. Commun. 181 (2010) 1517 [0911.3191]. [11] R.S. Sutton, The bitter lesson, Mar., 2019. [12] G. Kanwar, M.S. Albergo, D. Boyda, K. Cranmer, D.C. Hackett, S. Racanière et al., Equivariant flow-based sampling for lattice gauge theory, Phys. Rev. Lett. 125 (2020) 121601 [2003.06413]. [13] M. Favoni, A. Ipp, D.I. Müller and D. Schuh, Lattice Gauge Equivariant Convolutional Neural Networks, Phys. Rev. Lett. 128 (2022) 032003 [2012.12901]. [14] D.P. Kingma and J. Ba, Adam: A method for stochastic optimization, arXiv preprint arXiv:1412.6980 (2014) . 33

[15] S. Duane, A.D. Kennedy, B.J. Pendleton and D. Roweth, Hybrid Monte Carlo, Phys. Lett. B 195 (1987) 216. [16] ALPHA collaboration, Monte Carlo errors with less errors, Comput. Phys. Commun. 156 (2004) 143 [hep-lat/0306017]. [17] U. Wolff, Critical slowing down, Nucl. Phys. B Proc. Suppl. 17 (1990) 93. [18] L. Del Debbio, G.M. Manca and E. Vicari, Critical slowing down of topological modes, Phys. Lett. B 594 (2004) 315 [hep-lat/0403001]. [19] ALPHA collaboration, Critical slowing down and error analysis in lattice QCD simulations, Nucl. Phys. B 845 (2011) 93 [1009.5228]. [20] M. Hasenbusch, Fighting topological freezing in the two-dimensional CPN-1 model, Phys. Rev. D 96 (2017) 054504 [1706.04443]. [21] C. Bonanno, C. Bonati and M. D’Elia, Large-N SU (N ) Yang-Mills theories with milder topological freezing, JHEP 03 (2021) 111 [2012.14000]. [22] T. Eichhorn, G. Fuwa, C. Hoelbling and L. Varnhorst, Parallel tempered metadynamics: Overcoming potential barriers without surfing or tunneling, Phys. Rev. D 109 (2024) 114504 [2307.04742]. [23] G.G. Batrouni, G.R. Katz, A.S. Kronfeld, G.P. Lepage, B. Svetitsky and K.G. Wilson, Langevin Simulations of Lattice Field Theories, Phys. Rev. D 32 (1985) 2736. [24] C.T.H. Davies, G.G. Batrouni, G.R. Katz, A.S. Kronfeld, G.P. Lepage, K.G. Wilson et al., Fourier Acceleration in Lattice Gauge Theories. 1. Landau Gauge Fixing, Phys. Rev. D 37 (1988) 1581. [25] S. Duane and B.J. Pendleton, Gauge Invariant Fourier Acceleration, Phys. Lett. B 206 (1988) 101. [26] D. Levy, M.D. Hoffman and J. Sohl-Dickstein, Generalizing Hamiltonian Monte Carlo with Neural Networks, 1711.09268. [27] S. Foreman, X.-Y. Jin and J.C. Osborn, Deep Learning Hamiltonian Monte Carlo, in 9th International Conference on Learning Representations, 5, 2021 [2105.03418]. [28] S. Foreman, X.-Y. Jin and J.C. Osborn, LeapfrogLayers: A Trainable Framework for Effective Topological Sampling, PoS LATTICE2021 (2022) 508 [2112.01582]. [29] S. Foreman, X.-Y. Jin and J.C. Osborn, MLMC: Machine Learning Monte Carlo for Lattice Gauge Theory, PoS LATTICE2023 (2023) 036 [2312.08936].

34

[30] U. Wolff, Collective Monte Carlo Updating for Spin Systems, Phys. Rev. Lett. 62 (1989) 361. [31] R.H. Swendsen and J.-S. Wang, Nonuniversal critical dynamics in Monte Carlo simulations, Phys. Rev. Lett. 58 (1987) 86. [32] L. Huang and L. Wang, Accelerated Monte Carlo simulations with restricted Boltzmann machines, Phys. Rev. B 95 (2017) 035105 [1610.02746]. [33] J. Liu, Y. Qi, Z.Y. Meng and L. Fu, Self-learning Monte Carlo method, Phys. Rev. B 95 (2017) 041101 [1610.03137]. [34] J. Liu, H. Shen, Y. Qi, Z.Y. Meng and L. Fu, Self-learning monte carlo method and cumulative update in fermion systems, Phys. Rev. B 95 (2017) 241104 [1611.09364]. [35] X.Y. Xu, Y. Qi, J. Liu, L. Fu and Z.Y. Meng, Self-learning quantum Monte Carlo method in interacting fermion systems, Phys. Rev. B 96 (2017) 041119 [1612.03804]. [36] Y. Nagai, H. Shen, Y. Qi, J. Liu and L. Fu, Self-learning Monte Carlo method: Continuous-time algorithm, Phys. Rev. B 96 (2017) 161102 [1705.06724]. [37] H. Shen, J. Liu and L. Fu, Self-learning monte carlo with deep neural networks, Phys. Rev. B 97 (2018) 205140 [1801.01127]. [38] C. Chen, X.Y. Xu, J. Liu, G. Batrouni, R. Scalettar and Z.Y. Meng, Symmetry-enforced self-learning Monte Carlo method applied to the Holstein model, Phys. Rev. B 98 (2018) 041102 [1802.06177]. [39] Y. Liu, W. Wang, K. Sun and Z.Y. Meng, Designer Monte Carlo simulation for the Gross-Neveu-Yukawa transition, Phys. Rev. B 101 (2020) 064308 [1910.07430]. [40] Y. Nagai, M. Okumura and A. Tanaka, Self-learning Monte Carlo method with Behler-Parrinello neural networks, Phys. Rev. B 101 (2020) 115111 [1807.04955]. [41] M. Medvidovic, J. Carrasquilla, L.E. Hayward and B. Kulchytskyy, Generative models for sampling of lattice field theories, in 34th Conference on Neural Information Processing Systems, 12, 2020 [2012.01442]. [42] Y. Nagai and A. Tomiya, Gauge covariant neural network for quarks and gluons, Phys. Rev. D 111 (2025) 074501 [2103.11965]. [43] Y. Nagai and A. Tomiya, Self-learning Monte Carlo with equivariant Transformer, J. Phys. Soc. Jap. 93 (2024) 114007 [2306.11527]. [44] Y. Nagai, H. Ohno and A. Tomiya, CASK: A Gauge Covariant Transformer for Lattice Gauge Theory, PoS LATTICE2024 (2025) 030 [2501.16955].

35

[45] P. Hasenfratz and F. Niedermayer, Perfect lattice action for asymptotically free theories, Nucl. Phys. B 414 (1994) 785 [hep-lat/9308004]. [46] M. Blatter and F. Niedermayer, New fixed point action for SU(3) lattice gauge theory, Nucl. Phys. B 482 (1996) 286 [hep-lat/9605017]. [47] K. Holland, A. Ipp, D.I. Müller and U. Wenger, Machine learning a fixed point action for SU(3) gauge theory with a gauge equivariant convolutional neural network, Phys. Rev. D 110 (2024) 074502 [2401.06481]. [48] K. Holland, A. Ipp, D.I. Müller and U. Wenger, Machine-Learned Renormalization-Group-Improved Gauge Actions and Classically Perfect Gradient Flows, Phys. Rev. Lett. 136 (2026) 031901 [2504.15870]. [49] J. Aronsson, D.I. Müller and D. Schuh, Geometrical aspects of lattice gauge equivariant convolutional neural networks, 2303.11448. [50] M. Luscher, Trivializing maps, the Wilson flow and the HMC algorithm, Commun. Math. Phys. 293 (2010) 899 [0907.5491]. [51] G.P. Engel and S. Schaefer, Testing trivializing maps in the Hybrid Monte Carlo algorithm, Comput. Phys. Commun. 182 (2011) 2107 [1102.1852]. [52] S.-H. Li and L. Wang, Neural Network Renormalization Group, Phys. Rev. Lett. 121 (2018) 260601. [53] R. Abbott et al., Aspects of scaling and scalability for flow-based sampling of lattice QCD, Eur. Phys. J. A 59 (2023) 257 [2211.07541]. [54] M. Dibak, L. Klein, A. Krämer and F. Noé, Temperature steerable flows and boltzmann generators, Physical Review Research 4 (2022) L042005. [55] M. Gerdes, P. de Haan, C. Rainone, R. Bondesan and M.C.N. Cheng, Learning Lattice Quantum Field Theories with Equivariant Continuous Flows, SciPost Phys. 15 (2023) 238 [2207.00283]. [56] A. Singha, D. Chakrabarti and V. Arora, Conditional normalizing flow for Markov chain Monte Carlo sampling in the critical region of lattice field theory, Phys. Rev. D 107 (2023) 014512 [2207.00980]. [57] Q. Zhu, G. Aarts, W. Wang, K. Zhou and L. Wang, Physics-conditioned diffusion models for lattice gauge theory, JHEP 03 (2026) 111 [2502.05504]. [58] E.G. Tabak and E. Vanden-Eijnden, Density estimation by dual ascent of the log-likelihood, Communications in Mathematical Sciences 8 (2010) 217. [59] E.G. Tabak and C.V. Turner, A family of nonparametric density estimation algorithms, Communications on Pure and Applied Mathematics 66 (2013) 145. 36

[60] M.S. Albergo, G. Kanwar and P.E. Shanahan, Flow-based generative models for Markov chain Monte Carlo in lattice field theory, Phys. Rev. D 100 (2019) 034515 [1904.12072]. [61] K.A. Nicoli, S. Nakajima, N. Strodthoff, W. Samek, K.-R. Müller and P. Kessel, Asymptotically unbiased estimation of physical observables with neural samplers, Phys. Rev. E 101 (2020) 023304 [1910.13496]. [62] K.A. Nicoli, C.J. Anders, L. Funcke, T. Hartung, K. Jansen, P. Kessel et al., Estimation of Thermodynamic Observables in Lattice Field Theories with Deep Generative Models, Phys. Rev. Lett. 126 (2021) 032001 [2007.07115]. [63] D. Boyda, G. Kanwar, S. Racanière, D.J. Rezende, M.S. Albergo, K. Cranmer et al., Sampling using SU(N ) gauge equivariant flows, Phys. Rev. D 103 (2021) 074504 [2008.05456]. [64] L. Del Debbio, J.M. Rossney and M. Wilson, Efficient modeling of trivializing maps for lattice ϕ4 theory using normalizing flows: A first look at scalability, Phys. Rev. D 104 (2021) 094507 [2105.12481]. [65] M.S. Albergo, G. Kanwar, S. Racanière, D.J. Rezende, J.M. Urban, D. Boyda et al., Flow-based sampling for fermionic lattice field theories, Phys. Rev. D 104 (2021) 114507 [2106.05934]. [66] D.C. Hackett, C.-C. Hsieh, M.S. Albergo, D. Boyda, J.-W. Chen, K.-F. Chen et al., Flow-based sampling for multimodal distributions in lattice field theory, 2107.00734. [67] S. Foreman, T. Izubuchi, L. Jin, X.-Y. Jin, J.C. Osborn and A. Tomiya, HMC with Normalizing Flows, PoS LATTICE2021 (2022) 073 [2112.01586]. [68] J. Finkenrath, Tackling critical slowing down using global correction steps with equivariant flows: the case of the Schwinger model, 2201.02216. [69] M.S. Albergo, D. Boyda, K. Cranmer, D.C. Hackett, G. Kanwar, S. Racanière et al., Flow-based sampling in the lattice Schwinger model at criticality, Phys. Rev. D 106 (2022) 014514 [2202.11712]. [70] R. Abbott et al., Gauge-equivariant flow models for sampling in lattice field theories with pseudofermions, Phys. Rev. D 106 (2022) 074506 [2207.08945]. [71] K.A. Nicoli, C.J. Anders, T. Hartung, K. Jansen, P. Kessel and S. Nakajima, Detecting and mitigating mode-collapse for flow-based sampling of lattice field theories, Phys. Rev. D 108 (2023) 114501 [2302.14082]. [72] D. Albandea, L. Del Debbio, P. Hernández, R. Kenway, J. Marsh Rossney and A. Ramos, Learning trivializing flows, Eur. Phys. J. C 83 (2023) 676 [2302.08408]. 37

[73] R. Abbott et al., Normalizing flows for lattice gauge theory in arbitrary space-time dimension, 2305.02402. [74] R. Abbott, A. Botev, D. Boyda, D.C. Hackett, G. Kanwar, S. Racanière et al., Applications of flow models to the generation of correlated lattice QCD ensembles, Phys. Rev. D 109 (2024) 094514 [2401.10874]. [75] A. Bulgarelli, E. Cellini, K. Jansen, S. Kühn, A. Nada, S. Nakajima et al., Flow-Based Sampling for Entanglement Entropy and the Machine Learning of Defects, Phys. Rev. Lett. 134 (2025) 151601 [2410.14466]. [76] D. Schuh, J. Kreit, E. Berkowitz, L. Funcke, T. Luu, K.A. Nicoli et al., Simulating correlated electrons with symmetry-enforced normalizing flows, Phys. Rev. Res. 8 (2026) L012071 [2506.17015]. [77] A. Hasenfratz, E.T. Neil, L. Parato and N. Schwartz, Renormalization-guided cascade upscaling for lattice field generation, 2608.28581. [78] A. Hasenfratz, E.T. Neil, L. Parato and N. Schwartz, Renormalization-guided inverse blocking for lattice field generation: construction and validation, 2608.28580. [79] A. Singha, J. Kauffmann, K. Jansen, J. Finkenrath, V. Arora and S. Nakajima, Sector-Resolved Flow Sampling for Topologically Frozen Lattice Gauge Theories, 2609.07781. [80] A. Singha, J. Finkenrath, K. Jansen, V. Arora and S. Nakajima, Multilevel Plaquette-Space Sampling for Lattice Gauge Theories with Local Constraint Solves, 2609.35722. [81] P. de Haan, C. Rainone, M.C.N. Cheng and R. Bondesan, Scaling Up Machine Learning For Quantum Field Theory with Equivariant Continuous Flows, 2110.02673. [82] S. Bacchio, P. Kessel, S. Schaefer and L. Vaitl, Learning trivializing gradient flows for lattice gauge theories, Phys. Rev. D 107 (2023) L051504 [2212.08469]. [83] M. Caselle, E. Cellini and A. Nada, Sampling the lattice Nambu-Goto string using Continuous Normalizing Flows, JHEP 02 (2024) 048 [2307.01107]. [84] M. Gerdes, P. de Haan, R. Bondesan and M.C.N. Cheng, Nonperturbative trivializing flows for lattice gauge theories, Phys. Rev. D 112 (2025) 094516 [2410.13161]. [85] L. Dinh, D. Krueger and Y. Bengio, NICE: Non-linear independent components estimation, 1410.8516. [86] L. Dinh, J. Sohl-Dickstein and S. Bengio, Density estimation using Real NVP, 1605.08803. 38

[87] R.T.Q. Chen, Y. Rubanova, J. Bettencourt and D. Duvenaud, Neural Ordinary Differential Equations, 1806.07366. [88] L. Zhang, L. Wang et al., Monge-ampère flow for generative modeling, 1809.10188. [89] Q. Zhang and Y. Chen, Diffusion normalizing flow, Advances in Neural Information Processing Systems 34 (2021) 16280 [2110.07579]. [90] J. Köhler, L. Klein and F. Noé, Equivariant flows: exact likelihood generative learning for symmetric densities, in International Conference on Machine Learning, pp. 5361–5370, PMLR, 2020 [2006.02425]. [91] Y. Namekawa, K. Kashiwa, H. Matsuda, A. Ohnishi and H. Takase, Improving efficiency of the path optimization method for a gauge theory, Phys. Rev. D 107 (2023) 034509 [2210.05402]. [92] C. Lehner and T. Wettig, Gauge-equivariant neural networks as preconditioners in lattice QCD, Phys. Rev. D 108 (2023) 034503 [2302.05419]. [93] D. Knüttel, C. Lehner and T. Wettig, Gauge-equivariant multigrid neural networks, PoS LATTICE2023 (2023) 037. [94] C. Lehner and T. Wettig, Gauge-equivariant pooling layers for preconditioners in lattice QCD, Phys. Rev. D 110 (2024) 034517 [2304.10438]. [95] K. Holland, A. Ipp, D.I. Müller and U. Wenger, Fixed point actions from convolutional neural networks, PoS LATTICE2023 (2023) 038 [2311.17816]. [96] C. Jarzynski, Nonequilibrium Equality for Free Energy Differences, Phys. Rev. Lett. 78 (1997) 2690 [cond-mat/9610209]. [97] C. Jarzynski, Equilibrium free-energy differences from nonequilibrium measurements: A master-equation approach, Phys. Rev. E 56 (1997) 5018–5035 [cond-mat/9707325]. [98] R.M. Neal, Annealed importance sampling, Statistics and Computing 11 (2001) 125 [physics/9803008]. [99] M. Caselle, G. Costagliola, A. Nada, M. Panero and A. Toniato, Jarzynski’s theorem for lattice gauge theory, Phys. Rev. D 94 (2016) 034503 [1604.05544]. [100] M. Caselle, A. Nada and M. Panero, QCD thermodynamics from lattice calculations with nonequilibrium methods: The SU(3) equation of state, Phys. Rev. D 98 (2018) 054513 [1801.03110]. [101] C. Bonanno, A. Nada and D. Vadacchino, Mitigating topological freezing using out-of-equilibrium simulations, JHEP 04 (2024) 126 [2402.06561].

39

[102] H. Wu, J. Köhler and F. Noé, Stochastic normalizing flows, Advances in Neural Information Processing Systems 33 (2020) 5933 [2002.06707]. [103] M. Caselle, E. Cellini, A. Nada and M. Panero, Stochastic normalizing flows as non-equilibrium transformations, JHEP 07 (2022) 015 [2201.08862]. [104] M. Caselle, E. Cellini and A. Nada, Numerical determination of the width and shape of the effective string using Stochastic Normalizing Flows, JHEP 02 (2025) 090 [2409.15937]. [105] A. Bulgarelli, E. Cellini and A. Nada, Scaling of stochastic normalizing flows in SU(3) lattice gauge theory, Phys. Rev. D 111 (2025) 074517 [2412.00200]. [106] C. Bonanno, A. Bulgarelli, E. Cellini, A. Nada, D. Panfalone, D. Vadacchino et al., Scaling flow-based approaches for topology sampling in SU(3) gauge theory, JHEP 04 (2026) 051 [2510.25704]. [107] S. Chen, M. Qian, G. Aarts, B. Lucini and K. Zhou, Stochastic Path Sampler For Lattice Field Theory, 2606.13790. [108] A.G.D.G. Matthews, M. Arbel, D.J. Rezende and A. Doucet, Continual Repeated Annealed Flow Transport Monte Carlo, in International Conference on Machine Learning, pp. 15196–15219, PMLR, 2022, https://proceedings.mlr.press/v162/matthews22a.html [2201.13117]. [109] M.S. Albergo and E. Vanden-Eijnden, NETS: A Non-Equilibrium Transport Sampler, 2410.02711. [110] F. Vargas, S. Padhy, D. Blessing and N. Nüsken, Transport meets variational inference: Controlled monte carlo diffusions, in The Twelfth International Conference on Learning Representations, 2024, https://openreview.net/forum?id=PP1rudnxiW [2307.01050]. [111] R. Kapust, F. Ihssen and J.M. Pawlowski, Generative sampling with physics-informed kernels, PoS LATTICE2025 (2026) 055 [2510.26678]. [112] Y. Song, J. Sohl-Dickstein, D.P. Kingma, A. Kumar, S. Ermon and B. Poole, Score-based generative modeling through stochastic differential equations, 2011.13456. [113] L. Wang, Stochastic Quantization as Optimal Control, 2607.21436. [114] D. Sexty, Accurate Sampling from Diffusion Models, 2609.29902. [115] G. Kanwar and O. Vega, Annealed sampling of diffusion models in lattice field theory, 2610.35792. [116] L. Wang, G. Aarts and K. Zhou, Diffusion models as stochastic quantization in lattice field theory, JHEP 05 (2024) 060 [2309.17082]. 40

[117] G. Aarts, D.E. Habibi, L. Wang and K. Zhou, Combining complex Langevin dynamics with score-based and energy-based diffusion models, JHEP 12 (2025) 160 [2510.01328]. [118] O. Vega, J. Komijani, A. El-Khadra and M. Marinkovic, Group-equivariant diffusion models for lattice field theory, JHEP 07 (2026) 213 [2510.26081]. [119] H. Alharazin, J.Y. Panteleeva and B.D. Sun, Diffusion models for SU(2) lattice gauge theory in two dimensions, Phys. Rev. D 114 (2026) 014522 [2602.09045]. [120] G. Aarts, D.E. Habibi, A. Ipp, D.I. Müller, T.R. Ranner, L. Wang et al., Generalizable Equivariant Diffusion Models for Non-Abelian Lattice Gauge Theory, 2601.19552. [121] J. Komijani, M.K. Marinkovic and L. Turgut, Diffusion model for SU(N) gauge theories, 2605.06134. [122] A. Tomiya, Lattice Configuration Generation with a Self-Learning Diffusion Model, 2607.12587. [123] A. Tomiya, Continuous Hasenbusch transport towards gauge diffusion with fermions, arXiv (2026) [2609.15959]. [124] M.S. Albergo and E. Vanden-Eijnden, Building Normalizing Flows with Stochastic Interpolants, 2209.15571. [125] M.S. Albergo, N.M. Boffi and E. Vanden-Eijnden, Stochastic interpolants: A unifying framework for flows and diffusions, 2303.08797. [126] Y. Lipman, R.T.Q. Chen, H. Ben-Hamu, M. Nickel and M. Le, Flow matching for generative modeling, 2210.02747. [127] J. Engels, J. Fingberg, F. Karsch, D. Miller and M. Weber, Nonperturbative thermodynamics of SU(N) gauge theories, Phys. Lett. B 252 (1990) 625. [128] K. Binder, B. Block, S.K. Das, P. Virnau and D. Winter, Monte carlo methods for estimating interfacial free energies and line tensions, J. Stat. Phys. 144 (2011) 690–729 [1103.2241]. [129] M. Caselle, M. Hasenbusch and M. Panero, High precision Monte Carlo simulations of interfaces in the three-dimensional ising model: A Comparison with the Nambu-Goto effective string model, JHEP 03 (2006) 084 [hep-lat/0601023]. [130] P. de Forcrand, M. D’Elia and M. Pepe, A Study of the ’t Hooft loop in SU(2) Yang-Mills theory, Phys. Rev. Lett. 86 (2001) 1438 [hep-lat/0007034]. [131] P. Calabrese and J.L. Cardy, Entanglement entropy and quantum field theory, J. Stat. Mech. 0406 (2004) P06002 [hep-th/0405152]. 41

[132] V. Alba, Out-of-equilibrium protocol for Rényi entropies via the Jarzynski equality, Phys. Rev. E 95 (2017) 062132 [1609.02157]. [133] A. Bulgarelli and M. Panero, Entanglement entropy from non-equilibrium Monte Carlo simulations, JHEP 06 (2023) 030 [2304.03311]. [134] A. Bulgarelli and M. Panero, Duality transformations and the entanglement entropy of gauge theories, JHEP 06 (2024) 041 [2404.01987]. [135] A. Bulgarelli, M. Caselle, A. Nada and M. Panero, Casimir effect in critical O(N) models from nonequilibrium Monte Carlo simulations, Phys. Rev. E 112 (2025) 064126 [2505.20403]. [136] TUMQCD collaboration, Wilson loops with neural networks, Phys. Rev. D 113 (2026) 094501 [2602.02436]. [137] J. Mayer-Steudte, Neural network interpolators for Wilson loops, PoS LATTICE2025 (2026) 186 [2604.07226]. [138] G. Parisi, The strategy for computing the hadronic mass spectrum, Phys. Rept. 103 (1984) 203. [139] G.P. Lepage, The analysis of algorithms for lattice field theory, in Theoretical Advanced Study Institute in Elementary Particle Physics, 6, 1989. [140] P.E. Gibbs, Lattice Monte Carlo Simulations of QCD at Finite Baryonic Density, Phys. Lett. B182 (1986) 369. [141] K. Splittorff and J.J.M. Verbaarschot, Phase of the Fermion Determinant at Nonzero Chemical Potential, Phys. Rev. Lett. 98 (2007) 031601 [hep-lat/0609076]. [142] K. Splittorff and J.J.M. Verbaarschot, The QCD Sign Problem for Small Chemical Potential, Phys. Rev. D75 (2007) 116003 [hep-lat/0702011]. [143] G. Aarts, Can stochastic quantization evade the sign problem? The relativistic Bose gas at finite chemical potential, Phys. Rev. Lett. 102 (2009) 131601 [0810.2089]. [144] P. de Forcrand, Simulating QCD at finite density, PoS LAT2009 (2009) 010 [1005.0539]. [145] A. Alexandru, C. Gattringer, H.P. Schadler, K. Splittorff and J.J.M. Verbaarschot, Distribution of Canonical Determinants in QCD, Phys. Rev. D91 (2015) 074501 [1411.4143]. [146] S.R. Beane, W. Detmold, T.C. Luu, K. Orginos, A. Parreño, M.J. Savage et al., High Statistics Analysis using Anisotropic Clover Lattices: (I) Single Hadron Correlation Functions, Phys. Rev. D79 (2009) 114502 [0903.2990]. 42

[147] S.R. Beane, W. Detmold, K. Orginos and M.J. Savage, Uncertainty Quantification in Lattice QCD Calculations for Nuclear Physics, J. Phys. G 42 (2015) 034022 [1410.2937]. [148] M.L. Wagman and M.J. Savage, Statistics of baryon correlation functions in lattice QCD, Phys. Rev. D96 (2017) 114508 [1611.07643]. [149] M.L. Wagman, Statistical Angles on the Lattice QCD Signal-to-Noise Problem, Ph.D. thesis, U. Washington, Seattle (main), 2017. 1711.00062. [150] Z. Davoudi, W. Detmold, K. Orginos, A. Parreño, M.J. Savage, P. Shanahan et al., Nuclear matrix elements from lattice QCD for electroweak and beyond-Standard-Model processes, 2008.11160. [151] T. Bhattacharya, S. Lawrence and J.-S. Yoo, Control variates for lattice field theory, Phys. Rev. D 109 (2024) L031505 [2307.14950]. [152] P.F. Bedaque and H. Oh, Leveraging neural control variates for enhanced precision in lattice field theory, Phys. Rev. D 109 (2024) 094519 [2312.08228]. [153] S. Lawrence, Schwinger-Dyson control variates for lattice fermions, 2404.10707. [154] S. Lawrence and Y. Yamauchi, Mitigating a discrete sign problem with extreme learning machines, 2312.12636. [155] R. Abbott, D. Boyda, Y. Fu, D.C. Hackett, G. Kanwar, F. Romero-López et al., Variance reduction in lattice QCD observables via normalizing flows, Phys. Rev. D 114 (2026) 014513 [2603.02984]. [156] R. Abbott, Y. Fu, D.C. Hackett, G. Kanwar, F. Romero-López and P.E. Shanahan, Learning the generating functional for variance reduction in lattice QCD, 2606.15986. [157] P. Butti, G. Catumba, A. Nada and L. Spatscheck, A variational framework for variance reduction in lattice field theory, 2609.28142. [158] S. Lawrence and Y. Yamauchi, Normalizing Flows and the Real-Time Sign Problem, Phys. Rev. D 103 (2021) 114509 [2101.05755]. [159] A. Alexandru, P.F. Bedaque and N.C. Warrington, Spin polarized nonrelativistic fermions in 1+1 dimensions, Phys. Rev. D98 (2018) 054514 [1805.00125]. [160] A. Alexandru, P.F. Bedaque, H. Lamm and S. Lawrence, Finite-Density Monte Carlo Calculations on Sign-Optimized Manifolds, Phys. Rev. D97 (2018) 094510 [1804.00697].

43

[161] A. Alexandru, P.F. Bedaque, H. Lamm, S. Lawrence and N.C. Warrington, Fermions at Finite Density in 2+1 Dimensions with Sign-Optimized Manifolds, Phys. Rev. Lett. 121 (2018) 191602 [1808.09799]. [162] W. Detmold, G. Kanwar, M.L. Wagman and N.C. Warrington, Path integral contour deformations for noisy observables, Phys. Rev. D 102 (2020) 014514 [2003.05914]. [163] W. Detmold, G. Kanwar, H. Lamm, M.L. Wagman and N.C. Warrington, Path integral contour deformations for observables in SU(N ) gauge theory, Phys. Rev. D (2021) [2101.12668]. [164] M. Rodekamp, E. Berkowitz, C. Gäntgen, S. Krieg, T. Luu and J. Ostmeyer, Mitigating the Hubbard sign problem with complex-valued neural networks, Phys. Rev. B 106 (2022) 125139 [2203.00390]. [165] M. Rodekamp, E. Berkowitz, M. Dincă, C. Gäntgen, S. Krieg and T. Luu, From Theory to Practice: Applying Neural Networks to Simulate Real Systems with Sign Problems, PoS LATTICE2023 (2024) 031 [2311.18312]. [166] C. Gäntgen, E. Berkowitz, T. Luu, J. Ostmeyer and M. Rodekamp, Reducing the Sign Problem with simple Contour Deformation, PoS LATTICE2023 (2024) 007 [2407.06896]. [167] W. Detmold, G. Kanwar, Y. Lin, P.E. Shanahan and M.L. Wagman, Signal-to-noise improvement through neural network contour deformations for 3D SU (2) lattice gauge theory, PoS LATTICE2023 (2023) 043 [2309.00600]. [168] AuroraScience collaboration, New approach to the sign problem in quantum field theories: High density QCD on a Lefschetz thimble, Phys. Rev. D86 (2012) 074506 [1205.3996]. [169] G. Aarts, Lefschetz thimbles and stochastic quantization: Complex actions in the complex plane, Phys. Rev. D88 (2013) 094501 [1308.4811]. [170] M. Cristoforetti, F. Di Renzo, A. Mukherjee and L. Scorzato, Monte Carlo simulations on the Lefschetz thimble: Taming the sign problem, Phys. Rev. D88 (2013) 051501 [1303.7204]. [171] A. Mukherjee, M. Cristoforetti and L. Scorzato, Metropolis Monte Carlo integration on the Lefschetz thimble: Application to a one-plaquette model, Phys. Rev. D88 (2013) 051502 [1308.0233]. [172] G. Aarts, L. Bongiovanni, E. Seiler and D. Sexty, Some remarks on Lefschetz thimbles and complex Langevin dynamics, JHEP 10 (2014) 159 [1407.2090].

44

[173] M. Cristoforetti, F. Di Renzo, G. Eruzzi, A. Mukherjee, C. Schmidt, L. Scorzato et al., An efficient method to compute the residual phase on a Lefschetz thimble, Phys. Rev. D89 (2014) 114505 [1403.5637]. [174] A. Alexandru, G. Başar and P. Bedaque, Monte Carlo algorithm for simulating fermions on Lefschetz thimbles, Phys. Rev. D93 (2016) 014504 [1510.03258]. [175] A. Alexandru, G. Başar, P.F. Bedaque, G.W. Ridgway and N.C. Warrington, Sign problem and Monte Carlo calculations beyond Lefschetz thimbles, JHEP 05 (2016) 053 [1512.08764]. [176] A. Alexandru, G. Başar, P.F. Bedaque, S. Vartak and N.C. Warrington, Monte Carlo Study of Real Time Dynamics on the Lattice, Phys. Rev. Lett. 117 (2016) 081602 [1605.08040]. [177] H. Fujii, S. Kamata and Y. Kikukawa, Monte Carlo study of Lefschetz thimble structure in one-dimensional Thirring model at finite density, JHEP 12 (2015) 125 [1509.09141]. [178] Y. Tanizaki, Y. Hidaka and T. Hayata, Lefschetz-thimble analysis of the sign problem in one-site fermion model, New J. Phys. 18 (2016) 033002 [1509.07146]. [179] A. Alexandru, P.F. Bedaque, H. Lamm and S. Lawrence, Deep Learning Beyond Lefschetz Thimbles, Phys. Rev. D96 (2017) 094505 [1709.01971]. [180] Y. Mori, K. Kashiwa and A. Ohnishi, Application of a neural network to the sign problem via the path optimization method, PTEP 2018 (2018) 023B04 [1709.03208]. [181] Y. Tanizaki, H. Nishimura and J.J.M. Verbaarschot, Gradient flows without blow-up for Lefschetz thimbles, JHEP 10 (2017) 100 [1706.03822]. [182] K. Kashiwa, Y. Mori and A. Ohnishi, Controlling the model sign problem via the path optimization method: Monte Carlo approach to a QCD effective model with Polyakov loop, Phys. Rev. D99 (2019) 014033 [1805.08940]. [183] M. Fukuma, N. Matsumoto and N. Umeda, Implementation of the HMC algorithm on the tempered Lefschetz thimble method, 1912.13303. [184] M. Fukuma, N. Matsumoto and N. Umeda, Applying the tempered Lefschetz thimble method to the Hubbard model away from half-filling, Phys. Rev. D100 (2019) 114510 [1906.04243]. [185] K. Kashiwa, Y. Mori and A. Ohnishi, Application of the path optimization method to the sign problem in an effective model of QCD with a repulsive vector-type interaction, Phys. Rev. D99 (2019) 114005 [1903.03679].

45

[186] Z.-G. Mou, P.M. Saffin and A. Tranberg, Quantum tunnelling, real-time dynamics and Picard-Lefschetz thimbles, JHEP 11 (2019) 135 [1909.02488]. [187] M. Ulybyshev, C. Winterowd and S. Zafeiropoulos, Lefschetz thimbles decomposition for the Hubbard model on the hexagonal lattice, Phys. Rev. D101 (2020) 014508 [1906.07678]. [188] G. Kanwar and M.L. Wagman, Real-time lattice gauge theory actions: Unitarity, convergence, and path integral contour deformations, Phys. Rev. D 104 (2021) 014513 [2103.02602]. [189] F. Ihssen, R. Kapust and J.M. Pawlowski, Solving sign problems with physics-informed kernels, arXiv (2026) [2603.03159]. [190] A. Alexandru, G. Basar, P.F. Bedaque and N.C. Warrington, Complex paths around the sign problem, Rev. Mod. Phys. 94 (2022) 015006 [2007.05436]. [191] S. Lawrence and Y. Yamauchi, Convex optimization of contour deformations, Phys. Rev. D 110 (2024) 014508 [2311.13002]. [192] J.-L. Wynen, E. Berkowitz, S. Krieg, T. Luu and J. Ostmeyer, Machine learning to alleviate Hubbard-model sign problems, Phys. Rev. B 103 (2021) 125153 [2006.11221]. [193] J.M. Pawlowski and J.M. Urban, Flow-based density of states for complex actions, Phys. Rev. D 108 (2023) 054511 [2203.01243]. [194] N. Hermansson-Truedsson and G. Kanwar, Normalizing flows for all-orders QED corrections in lattice field theory, arXiv (2026) [2605.22444]. [195] B. Yoon, T. Bhattacharya and R. Gupta, Machine Learning Estimators for Lattice QCD Observables, Phys. Rev. D 100 (2019) 014504 [1807.05971]. [196] R. Zhang, Z. Fan, R. Li, H.-W. Lin and B. Yoon, Machine-learning prediction for quasiparton distribution function matrix elements, Phys. Rev. D 101 (2020) 034516 [1909.10990]. [197] A. Frommer, V. Hannemann, B. Nockel, T. Lippert and K. Schilling, Accelerating Wilson fermion matrix inversions by means of the stabilized biconjugate gradient algorithm, Int. J. Mod. Phys. C 5 (1994) 1073 [hep-lat/9404013]. [198] T.A. DeGrand and P. Rossi, Conditioning Techniques for Dynamical Fermions, Comput. Phys. Commun. 60 (1990) 211. [199] M. Luscher, Solution of the Dirac equation in lattice QCD using a domain decomposition method, Comput. Phys. Commun. 156 (2004) 209 [hep-lat/0310048].

46

[200] S. Fischer, A. Frommer, U. Glassner, T. Lippert, G. Ritzenhofer and K. Schilling, A Parallel SSOR preconditioner for lattice QCD, Comput. Phys. Commun. 98 (1996) 20 [hep-lat/9602019]. [201] J. Brannick, R.C. Brower, M.A. Clark, J.C. Osborn and C. Rebbi, Adaptive Multigrid Algorithm for Lattice QCD, Phys. Rev. Lett. 100 (2008) 041601 [0707.4018]. [202] M. Luscher, Local coherence and deflation of the low quark modes in lattice QCD, JHEP 07 (2007) 081 [0706.2298]. [203] R. Babich, J. Brannick, R.C. Brower, M.A. Clark, T.A. Manteuffel, S.F. McCormick et al., Adaptive multigrid algorithm for the lattice Wilson-Dirac operator, Phys. Rev. Lett. 105 (2010) 201602 [1005.3043]. [204] A. Frommer, K. Kahl, S. Krieg, B. Leder and M. Rottmann, Adaptive Aggregation-Based Domain Decomposition Multigrid for the Lattice Wilson–Dirac Operator, SIAM J. Sci. Comput. 36 (2014) A1581 [1303.1377]. [205] QUDA collaboration, Accelerating lattice QCD multigrid on GPUs using fine-grained parallelization, in International Conference for High Performance Computing, Networking, Storage and Analysis, 12, 2016, DOI [1612.07873]. [206] S. Pfahler, D. Knüttel, C. Lehner and T. Wettig, A novel gauge-equivariant neural network architecture for preconditioners in lattice QCD, PoS LATTICE2025 (2026) 027 [2602.23840]. [207] S. Calı̀, D.C. Hackett, Y. Lin, P.E. Shanahan and B. Xiao, Neural-network preconditioners for solving the Dirac equation in lattice gauge theory, Phys. Rev. D 107 (2023) 034508 [2208.02728]. [208] Y. Sun, S. Eswar, Y. Lin, W. Detmold, P. Shanahan, X. Li et al., Matrix-free Neural Preconditioner for the Dirac Operator in Lattice Gauge Theory, arXiv (2025) [2509.10378].

47