On Estimation of Thermodynamic Observables in Lattice Field Theories with Deep Generative Models

References

Busza et al. (2018) W. Busza, K. Rajagopal, and W. Van Der Schee, Heavy ion collisions: the big picture and the big questions, Annual Review of Nuclear and Particle Science 68, 339 (2018).
Philipsen (2013) O. Philipsen, The qcd equation of state from the lattice, Progress in Particle and Nuclear Physics 70, 55 (2013).
De Forcrand et al. (2001) P. De Forcrand, M. D’Elia, and M. Pepe, ’t hooft loop in su(2) yang-mills theory, Physical Review Letters 86, 1438 (2001).
Jarzynski (1997a) C. Jarzynski, Equilibrium free-energy differences from nonequilibrium measurements: A master-equation approach, Physical Review E 56, 5018 (1997a).
Jarzynski (1997b) C. Jarzynski, Nonequilibrium equality for free energy differences, Physical Review Letters 78, 2690 (1997b).
Caselle et al. (2018) M. Caselle, A. Nada, and M. Panero, Qcd thermodynamics from lattice calculations with nonequilibrium methods: The su(3) equation of state, Physical Review D 98, 054513 (2018).
Goodfellow et al. (2016) I. Goodfellow, Y. Bengio, and A. Courville, Deep learning (MIT press, 2016).
Albergo et al. (2019) 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, 034515 (2019).
Kanwar et al. (2020) G. Kanwar, M. S. Albergo, D. Boyda, K. Cranmer, D. C. Hackett, S. Racanière, D. J. Rezende, and P. E. Shanahan, Equivariant flow-based sampling for lattice gauge theory, arXiv preprint arXiv:2003.06413 (2020).
Urban and Pawlowski (2018) J. M. Urban and J. M. Pawlowski, Reducing autocorrelation times in lattice simulations with generative adversarial networks, arXiv preprint arXiv:1811.03533 (2018).
Wirnsberger et al. (2020) P. Wirnsberger, A. J. Ballard, G. Papamakarios, S. Abercrombie, S. Racanière, A. Pritzel, D. J. Rezende, and C. Blundell, Targeted free energy estimation via learned mappings, arXiv preprint arXiv:2002.04913 (2020).
Zwanzig (1955) R. W. Zwanzig, High-temperature equation of state by a perturbation method. ii. polar gases, The Journal of Chemical Physics 23, 1915 (1955).
Akiyama et al. (2020) S. Akiyama, D. Kadoh, Y. Kuramashi, T. Yamashita, and Y. Yoshimura, Tensor renormalization group approach to four-dimensional complex $ϕ^{4}$ theory at finite density, arXiv preprint arXiv:2005.04645 (2020).
Asakawa et al. (2014) M. Asakawa, T. Hatsuda, E. Itou, M. Kitazawa, H. Suzuki, F. Collaboration, et al., Thermodynamics of s u (3) gauge theory from gradient flow on the lattice, Physical Review D 90, 011501 (2014).
Giusti and Pepe (2014) L. Giusti and M. Pepe, Equation of state of a relativistic theory from a moving frame, Physical review letters 113, 031601 (2014).
Papamakarios et al. (2019) G. Papamakarios, E. Nalisnick, D. J. Rezende, S. Mohamed, and B. Lakshminarayanan, Normalizing flows for probabilistic modeling and inference, arXiv preprint arXiv:1912.02762 (2019).
Dinh et al. (2014) L. Dinh, D. Krueger, and Y. Bengio, Nice: Non-linear independent components estimation, arXiv preprint arXiv:1410.8516 (2014).
MacKay and Mac Kay (2003) D. J. MacKay and D. J. Mac Kay, Information theory, inference and learning algorithms (Cambridge university press, 2003).
(19) More precisely, the distributions $p$ and $q$ have to be identically only almost everywhere, i.e. up to a set of measure zero.
Kingma and Welling (2013) D. P. Kingma and M. Welling, Auto-encoding variational bayes, arXiv preprint arXiv:1312.6114 (2013).
Rezende et al. (2014) D. J. Rezende, S. Mohamed, and D. Wierstra, Stochastic backpropagation and approximate inference in deep generative models, arXiv preprint arXiv:1401.4082 (2014).
Rezende and Mohamed (2015) D. J. Rezende and S. Mohamed, Variational inference with normalizing flows, arXiv preprint arXiv:1505.05770 (2015).
Müller et al. (2019) T. Müller, B. Mcwilliams, F. Rousselle, M. Gross, and J. Novák, Neural importance sampling, ACM Transactions on Graphics (TOG) 38, 1 (2019).
Noé et al. (2019) F. Noé, S. Olsson, J. Köhler, and H. Wu, Boltzmann generators: Sampling equilibrium states of many-body systems with deep learning, Science 365 (2019).
Wu et al. (2019) D. Wu, L. Wang, and P. Zhang, Solving statistical mechanics using variational autoregressive networks, Physical Review Letters 122, 080602 (2019).
Nicoli et al. (2020)
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, Physical Review E
101, 023304 (2020).
Nicoli et al. (2019) K. Nicoli, P. Kessel, N. Strodthoff, W. Samek, K.-R. Müller, and S. Nakajima, Comment on ”solving statistical mechanics using vans”: Introducing savant-vans enhanced by importance and mcmc sampling, arXiv preprint arXiv:1903.11048 (2019).
Adler (1981) S. L. Adler, Over-relaxation method for the monte carlo evaluation of the partition function for multiquadratic actions, Physical Review D 23, 2901 (1981).
Whitmer (1984) C. Whitmer, Over-relaxation methods for monte carlo simulations of quadratic and multiquadratic actions, Physical Review D 29, 306 (1984).
Callaway and Rahman (1983) D. J. Callaway and A. Rahman, Lattice gauge theory in the microcanonical ensemble, Physical Review D 28, 1506 (1983).
Fodor and Jansen (1994) Z. Fodor and K. Jansen, Overrelaxation algorithm for coupled gauge-higgs systems, Physics Letters B 331, 119 (1994).
Wolff et al. (2004) U. Wolff, A. Collaboration, et al., Monte carlo errors with less errors, Computer Physics Communications 156, 143 (2004).
Burda et al. (2015) Y. Burda, R. Grosse, and R. Salakhutdinov, Importance weighted autoencoders, arXiv preprint arXiv:1509.00519 (2015).
Nowozin (2018)
S. Nowozin, Debiasing evidence approximations: On importance-weighted autoencoders and jackknife variational inference, ICLR 2018 (2018).
Teh et al. (2007)
Y. W. Teh, D. Newman, and M. Welling, A collapsed variational bayesian inference algorithm for latent dirichlet allocation, in
Advances in neural information processing systems (2007) pp. 1353–1360.
Gattringer and Lang (2009) C. Gattringer and C. Lang, Quantum chromodynamics on the lattice: an introductory presentation, Vol. 788 (Springer Science & Business Media, 2009).
Bickel and Doksum (2015) P. J. Bickel and K. A. Doksum, Mathematical statistics: basic ideas and selected topics, volume I, Vol. 117 (CRC Press, 2015).

Appendix A Conventions for Action

The form of the action $S$ used in the main text is

	$S (ϕ) = \sum x \in Λ - 2 κ 2 \sum^μ = 1 φ (x)$	$φ (x +^μ) + (1 - 2 λ) φ (x)^{2}$
		$+ λ φ (x)^{4} .$		(11)

It can be obtained by starting from the (more standard) action

	$S (φ) = \sum x \in Λ a^{2} \frac{1}{2} 2 \sum^μ = 1$	$\frac{φ (x + a^μ) - φ (x)}{a^{2}}$
		$+ \frac{m_{0}^{2}}{2} φ^{2} (x) + \frac{g_{0}}{4!} φ^{4} (x)$		(12)

and performing the following re-definitions

$φ$	$= (2 κ)^{\frac{1}{2}} ϕ,$	(13)
$(a m_{0})^{2}$	$= \frac{1 - 2 λ}{κ} - 4,$	(14)
$a^{2} g_{0}$	$= \frac{6 λ}{κ^{2}} .$	(15)

Appendix B A brief overview of Deep Learning

Neural Networks

Neural networks are a machine learning algorithm which has proven to be particularly powerful. A neural network is build of layers which are defined by

y^{(l)} (x) = σ (W^{l} x + b^{l}),

where $x \in R^{n}$ , $y^{l} \in R^{m}$ are input and output of the layer. The output of the layer is also often called the activation of the layer. The weights $W^{l} \in R^{m, n}$ and the bias $b^{l} \in R^{n}$ are the learnable parameters of the neural network. The non-linearity $σ : R \to R$ is a non-linear function which is applied element-wise to the components of $W^{l} x + b^{l}$

. Widely-used activation function are

σ (x) = max (x, 0)

σ (x) = tanh (x)

A neural network consists of $L$ such layers, i.e.

g (x) = (y^{(L)} \circ \dots \circ y^{(1)}) (x) .

It is important to note that the weights $W^{l}$ and biases $b^{l}$ do not have to be of the same dimensionality for each layer (although we did not make this explicit in our notation). It is also important to note that we merely described the most simple type of neural network, namely a fully-connected neural network. There is a zoo of other neural networks but we will refrain from a more detailed discussion as it is not needed for our purposes (see Goodfellow et al. (2016) for an overview).

Learning Parameters with Backpropagation

The parameters of the neural networks,

W = {(W^{l}, b^{l}), i = 1, \dots, L},

are determined by minimizing a certain loss function

L

by gradient descent (see (7) for the particular loss function used in this work). It is important to emphasize that the number of parameters are typically large (of order

10^{3} - 10^{6}

for typical modern neural networks). It is therefore clear that one cannot determine the gradient by finite-difference (as we would need calculate the finite difference ratio for each of these parameters which is prohibitively expensive).

The basic idea for calculating the gradient $\nabla_{W} L$ is to use the fact that we know the functional form of the neural network: the gradient of the loss is given by

\nabla_{W} L = \frac{\partial L}{\partial y^{L}} \frac{\partial y^{L}}{\partial y^{L - 1}} \dots \frac{\partial y^{1}}{\partial x} .

Each term in this expression is known, for example

\frac{\partial y^{l + 1}}{\partial y^{l}} = σ^{'} (y^{l}) W^{l} .

For a fixed non-linearity, we know the analytical form of the derivative $σ^{'}$ . This observation leads to the following algorithm: we first perform a forward pass of the neural network, i.e. starting from the input $x$ , we calculate the activations $y^{l}$ for each layer and store them in memory. This process ends with the final activation $y^{L}$ which is, by definition, the output of the neural network. The gradient of each layer $\frac{\partial y^{l + 1}}{\partial y^{l}} = σ^{'} (y^{l}) W^{l}$ can then directly be calculated (as we have stored the activation $y^{l}$ ). Crucially, we only need the matrix product of these Jacobians and it is efficient to start by calculating the gradient with respect to the output layer $L$ , then the layer $L - 1$ and so forth. This is because the loss function has a scalar output value and therefore the matrix product of the Jacobians

\frac{\partial L}{\partial y^{L}} \frac{\partial y^{L}}{\partial y^{L - 1}} \dots \frac{\partial y^{l + 1}}{\partial y^{l}}

is a vector with the same number of components as

y^{l}

. We can therefore save memory by simply overwriting the stored activation

y^{l}

. This algorithm is called backpropagation and allows us to calculate the gradient

\nabla_{W} L

for roughly the same cost as a forward pass of the neural network.

Appendix C Relation between $Var (C)$ and KL divergence

Theorem.

Let $C (ϕ) = S (ϕ) + ln q (ϕ)$ . The following relation between the KL divergence and the variance of $C$ holds:

KL (q_{θ} | | p) = \frac{1}{2} {Var}_{q} (C) + O (E_{q} [| w - 1 |^{3}]),

where $w (ϕ) = \frac{p (ϕ)}{q (ϕ)}$ is the normalized importance weight.

Proof.

The expectation value of the normalized importance weight is

E_{q} [\frac{p}{q}] = \int D [ϕ] p (ϕ) = 1 .

The Kullback-Leibler divergence can be rewritten in terms of the normalized importance weight

KL (q | p) = E_{q} [ln \frac{q}{p}] = - E_{q} [ln w] .

We now expand the KL divergence around the expectation value of the normalized importance weight

	$KL (q \| p)$	$= - E_{q} [ln (1 - (w - 1))]$
		$= E_{q} [\infty \sum j = 0 \frac{(- 1)^{j}}{j} (w - 1)^{j}]$
		$= E_{q} [w - 1]      = 0 + \frac{1}{2} E_{q} [(w - 1)^{2}] + O (E_{q} [\| w - 1 \|^{3}]) .$

We now relate this expression to the variance of $C$ . To this end, we first observe that $C = - ln ~ w$ , where $~ w = exp (- S) / q$ is the unnormalized importance weight. We then rewrite the expectation value of $C$ as

	$E_{q} [C]$	$= - E_{q} [ln ~ w]$
		$= - E_{q} [ln w + ln Z]$
		$= - ln Z - E_{q} [ln w]$
		$= - ln Z + KL (q \| p)$
		$= - ln Z + O (E_{q} [(w - 1)^{2}]),$

where the last step uses the expansion for the KL divergence derived above. It then follows its variance is given by

	${Var}_{q} (C)$	$= E_{q} [{(C - E_{q} [C])}_{q}^{2}]$
		$= E_{q} [{(- ln ~ w + E_{q} [ln ~ w])}_{q}^{2}]$
		$= E_{q} [(- ln ~ w + ln Z      = - ln w + O (E_{q} [(w - 1)^{2}]))^{2}]$

Expanding the logarithm around $E_{q} [w] = 1$ again, we obtain

{Var}_{q} [C] = E_{q} [(w - 1)^{2}] + O (E_{q} [| w - 1 |^{3}]) .

Combining this expression with the expansion derived for the KL divergence, we obtain the claim of the theorem. ∎

The higher-order moments will be small towards the end of the training process for which

q \approx p

and thus

w \approx 1

. Thus, the variance of

C

will become small. We indeed observe this behaviour in our numerical experiments, see Figure 3 for an example.

Figure 3: The variance ${Var}_{q} (C)$ decreases during training. We use hopping parameter $κ = 0.3$ and bare coupling $λ = 0.022$ for a $16 \times 8$ lattice.

Appendix D Analytic Solution for Partition Function

Figure 4: Integrated autocorrelation time of the free energy during refinement of the step size. The shaded red areas refers to the interval in hopping parameter $κ \in [0.265, 0.28]$ for which the refinement is applied. Darker shading indicates narrower step size. The experiments were performed using the overrelaxed HCM algorithm.

The action of the scalar field theory is given by

S = \sum x \in Λ - 2 κ 2 \sum^μ = 1 φ (x) φ (x +^μ) + (1 - 2 λ) φ (x)^{2} + λ φ (x)^{4} .

We want to calculate the partion function

Z = \int D [ϕ] exp (- S (ϕ)),

which for vanishing hopping parameter $κ$ decouples in independent integrals of each lattice site of the lattice $Λ$ :

Z = \prod x \in Λ (\int d ϕ (x) exp (- λ ϕ (x)^{4} - (1 - 2 λ) ϕ (x)^{2}))

The partition function can then be calculated analytically using the integral

\int exp (- a x^{4} - b x^{2}) d x = \sqrt{\frac{b}{4 a}} exp (\frac{b^{2}}{8 a}) K_{\frac{1}{4}} (\frac{b^{2}}{8 a}),

where $K_{n}$ is the modified Bessel function of the second kind. Using this formula, we obtain the following analytic form of the free energy

F = - T | Λ | ln (z),

where we have defined

z (λ) = \sqrt{\frac{1 - 2 λ}{4 λ}} exp (\frac{(1 - 2 λ)^{2}}{8 λ}) K_{\frac{1}{4}} (\frac{(1 - 2 λ)^{2}}{8 λ}) .

This result corresponds to the zeroth order of the hopping expansion Gattringer and Lang (2009) of the partition function $Z$ and one may, in principle, also calculate higher-order corrections. However, they are not needed for our purposes.

Appendix E Error Analysis for Free Energy Estimator

Our discussion is based on Nowozin (2018) which discussed the same results in the context of variational inference. We provide a review here since these results may be hard to extract for physicist not familiar with variational inference. We also point out a subtlety that was not discussed in the previous work.

The estimator for the free energy is given by

^F = - T ln \frac{1}{N} N \sum i = 0 ~ w (ϕ_{i}),

ϕ_{i} \sim q_{θ} .

(16)

Theorems for the variance and bias of this estimator are discussed in the following. For this, we use the delta method of moments which is summarized in the following theorem.

Theorem.

Let ${^X}_{N} = \frac{1}{N} \sum_{i = 1}^{N} X_{i}$ be the sample mean of independent and identically distributed random variables $X_{i}$ with $E [X_{i}^{2 k + 2}] < \infty$ for $k \in {0, 1}$ . Let $h$ be a real-valued function with uniformly bounded derivatives. It then holds that

E [h ({^X}_{N})] = c_{0} + \frac{c_{1}}{N} + O (\frac{1}{N^{2}}),

where

c_{0} = h (μ),

c_{1} = h^{''} (μ) \frac{σ^{2}}{2},

with $σ^{2} = E [(X - E X)^{2}]$ and $μ = E [X]$ .

We refer to Chapter 5.3 of Bickel and Doksum (2015) for a proof and more details.

The application of the delta method to the free energy estimator $^F$ is, in practise, subject to a subtlety regarding the bounded differentiablity of the function $h$ . We will ignore this subtlety in the following and return to it at the end of the section.

Theorem.

The bias of $^F$ is given by

B [- β^F] = - \frac{1}{2 N} \frac{E_{q} [(~ w - E_{q} [~ w])^{2}]}{E_{q} [~ w]^{2}} + O (N^{- 2}),

assuming that $E_{q} [{~ w}^{2 k + 2}] < \infty$ for $k = 0, 1$ .

Proof.

The bias of $- β^F = ln^Z$ is given by

B [- β^F] = E_{q} [ln^Z] - ln Z .

Using the delta method for moments, we derive that

E_{q} [ln^Z] = ln Z - \frac{1}{2 N Z^{2}} E_{q} [(~ w - E_{q} [~ w])^{2}] + O (N^{- 2}),

(17)

where we have used that $h (x) = ln (x)$ has second derivative $h^{''} (x) = - \frac{1}{x^{2}}$ . The proof then concludes by observing that

E_{q} [~ w] = Z .

∎

Theorem.

The variance of $^F$ is given by

assuming that $E_{q} [{~ w}^{2 k + 2}] < \infty$ for $k = 0, 1$ .

Proof.

The variance can be written as

Var [- β^F] = E_{q} [(ln^Z)^{2}] - E_{q} [ln^Z]^{2} .

We now evaluate both terms on the right-hand-side individually using the delta method. For the first term, we use the delta method with $h (x) = (ln x)^{2}$ which has second derivative

h^{''} (x) = \frac{2}{x^{2}} - 2 \frac{ln (x)}{x^{2}} .

Using this expression, we then obtain that $E_{q} [(ln^Z)^{2}]$ is equal to

(ln Z)^{2} + (\frac{1}{Z^{2}} - \frac{log Z}{Z^{2}}) \frac{E_{q} [(~ w - E_{q} [~ w])^{2}]}{N} + O (N^{- 2})

For the squared expectation value, we use the expansion (17) derived in the proof for the bias. This gives that $(E_{q} ln^Z)^{2}$ is equal to

	${(ln Z - \frac{1}{2 N Z^{2}} E_{q} [(~ w - E_{q} [~ w])^{2}] + O (N^{- 2}))}_{q}^{2}$
	$= (ln Z)^{2} - \frac{1}{N Z^{2}} E_{q} [(~ w - E_{q} [~ w])^{2}] + O (N^{- 2}) .$

Subtracting these two expressions, it then follows

Var [- β^F] = \frac{1}{Z^{2}} \frac{E_{q} [(~ w - E_{q} [~ w])^{2}]}{N} + O (N^{- 2}),

and the proof concludes by observing that $Z = E_{q} [~ w]$ . ∎

A few remarks are in order: from the theorems, it follows that the standard deviation of the estimator

^F

is of order

O (1 / \sqrt{N})

. In the large

N

limit, we can therefore neglect the bias correction as it is of order

O (N^{- 1})

. Furthermore, we can replace the expectation values in the theorems by the sample mean up to (negligible) higher-order corrections. In practise, we therefore use these results to estimate the variance and bias of

^F

. Alternatively, one can use a standard Jackknife analysis to estimate variance and bias (see for example Gattringer and Lang (2009)). In our experiments, we use both methods to estimate the errors and check that they lead to consistent results. Lastly, we remark that error estimators for general observables involving the partition function can be derived, see Nicoli et al. (2020).

As mentioned above, the delta method requires that the derivatives of the function $h$ are (uniformly) bounded. For a generic LQFT, this will not be the case for $h (x) = ln (x)$ since its derivatives diverge for $x \to 0^{+}$ . To the best of our knowledge, the same problem will generically arise in the context of variational inference but seems to have not been discussed in the literature.

To address this subtlety, one could require that the action of the lattice quantum field theory is bounded. For example, this can be ensured by putting the field theory in a box potential. Since only very high energy configurations are affected by this (for suitably large choice of the box potential) and since these configurations are extremely unlikely to be sampled, this modification will have no practical effect on the numerical experiments. After this modification, $^Z$ is bounded from below and $h^{(n)} (^Z)$ is also bounded as a result.

More rigorously, the result for the variance can be derived without assumptions on a bound for the derivatives by using the delta method for in law approximation which takes the following form

Theorem.

Let ${^X}_{N} = \frac{1}{N} \sum_{i = 1}^{N} X_{i}$ be the sample mean of independent and identically distributed random variables $X_{i}$ with $E [X_{i}^{k}] < \infty$ for $k \in {1, 2}$ . Let $h$ be a differential function at $μ = E_{q} [X]$ . Then

\sqrt{n} (h ({^X}_{N}) - h (μ)) D \to N (0, σ^{2} (h)),

where $σ^{2} (h) = h^{'} (μ) Var (X)$ .

For a proof, we again refer to Bickel and Doksum (2015), see Theorem 5.3.3. Applying this theorem to the free energy estimator $- β^F = ln^Z$ , we obtain the same expression for its variance as derived above. However, the theorem does not require any bound on the derivatives of $h (x) = ln (x)$ .

Appendix F Step Size Analysis

Figure 5: Free energy at $κ = 0.3$ obtained by both MCMC-based and flow-based method. For the MCMC-based method, different step size $δ κ$ of the hopping parameter were used. Details on the step sizes $δ κ$ are summarized in Table 1. For the flow-based method, we use the same number of samples as for the corresponding refined MCMC method. As a result, also the error of the flow’s estimate decreases.

Figure 6: Free energy at $κ = 0.3$ and $λ = 0.022$ using both $Z_{i} / Z_{i + 1}$ (down) and $Z_{i + 1} / Z_{i}$ (up). We use the same setup (i.e. number of steps, hopping parameter change $δ κ$ , etc) as for Figure 2.

As explained in the main text, the free energy difference $Δ F_{e b}$ is calculated in steps

Δ F_{e b} = Δ F_{e, i_{K}} + Δ F_{i_{K} i_{K - 1}} + \dots + Δ F_{i_{1} b} .

In this appendix, we will analyze the dependency of our results on the chosen steps.

We start from an initial step size corresponding to a change in hopping parameter $κ$ of $δ κ = 0.05$ . Between $κ = 0.2$ and $κ = 0.3$ , we however take a finer step size of $δ κ = 0.01$ . Since we are interested in the free energy difference $Δ F_{e b}$ between $κ_{b} = 0.0$ and $κ_{e} = 0.3$ , this corresponds to running $K = 14$ Markov chains. We focus on the $16 \times 8$ lattice and use an overelaxed HMC algorithm to sample $400 k$ configurations for each chain. The overrelaxation is performed every $10$ steps.

We then repeatedly refine the step size in a certain subregion around the critical $κ$ value. The details can be found in Table 1.

The results of this analysis are shown in Figure 5. We observe that the error of the estimator does not significantly decrease. We note that the error of the flow decreases during refinement because its free energy estimation uses the same number of samples as all Markov chains combined (and this number increases by the additional refinement steps).

In order to ensure that the distributions $p_{i}$ and $p_{i + 1}$ in

have sufficient overlap, we also estimate $\frac{Z_{i}}{Z_{i + 1}}$ by exchanging $p_{i}$ with $p_{i + 1}$ in the relation above. We then check that this leads to compatible results, see Figure 6. We note that this consistency check is relatively cheap as it requires running one additional Markov chain.

We also study the dependence of the integrated autocorrelation of the free energy on this refinement procedure, see Figure 4.

$δ κ$	$#$ chains	$#$ samples	refined $κ$ region
$0.01$	14	5.6 M	0.20-0.30
$0.005$	24	9.6 M	0.20-0.30
$0.0025$	40	16 M	0.22-0.30
$0.001$	76	30.4 M	0.24-0.30
$0.0005$	88	35.2 M	0.267-0.279

Table 1: Details on the refinement analysis. In each refinement stage, we take smaller steps

δ κ

in a certain subregion of the hopping parameter

κ

trajectory (see last column). The step size taken in this region is shown in the first column. Outside of the most refined region, the same step sizes as in the previous refinement stage are taken. Since each chain is taken to be of the same length, the total number of samples (third column) grows proportional to the number of chains (second column).

Appendix G Details on Numerical Experiments

Figure 7: Free energy estimation with error analysis by both Jackknife and delta method. Both lead to compatible results. We use the same data as in Figure 2.

Hmc:

We use a HMC algorithm with overrelaxation. Each Markov chain has 5k thermalization steps followed by 400k estimation steps. The sign of the field configuration is flipped every ten steps.

Training of flow:

For every lattice, we use a normalizing flow with six coupling layers. Each coupling layer (5) has neural network $m$ with five fully-connected layers with no bias and Tanh non-linearities. The hidden layers of $m$

consist of 1000 neurons each. We train the flow for 1M steps using an 8k mini-batch. We use ReduceLROnPlateau learning rate scheduler of PyTorch with an initial learning rate of

5 \times 10^{- 4}

and patience of 3k steps. The minimum learning rate was set to

1 \times 10^{- 7}

Estimation:

As described in the main text, for HMC-based estimation we use a step size of $δ κ = 0.01$ for $κ \in [0.2, 0.3]$ and a step size of $δ κ = 0.05$ for all other values of the hopping parameter. As a result 14 Markov chains are run. In total, the HMC-based method therefore uses $14 \times 400 k = 5.6 M$ configurations. We use the same number of samples for the flow-based estimation. For efficiency, we sample these configurations in mini-batches of 3k samples.

Error estimation:

we use both the uwerr Wolff et al. (2004) and jackknife method to estimate uncertainties for HMC. In order to deal with autocorrelation for jackknife, we perform binning with a 1k bin size. Error estimation for flow is performed by the delta method and also by jackknife, see Figure 7.

From the free energy estimates, one can then derive other thermodynamic observables such as the entropy. We refer to the main text for a discussion of this. Figure 8 shows estimation of entropy. Errors were estimated using both the Jackknife and uwerr method. Both error analysis methods lead to consistent results.

On Estimation of Thermodynamic Observables in Lattice Field Theories with Deep Generative Models

Stochastic normalizing flows as non-equilibrium transformations

Machine Learning of Thermodynamic Observables in the Presence of Mode Collapse

Computing Absolute Free Energy with Deep Generative Models

Deep Involutive Generative Models for Neural MCMC

Machine Learning and Variational Algorithms for Lattice Field Theory

Markov-Chain Monte Carlo Approximation of the Ideal Observer using Generative Adversarial Networks

Hierarchical autoregressive neural networks for statistical systems

References

Appendix A Conventions for Action