Ergodicity of the LLR method for the Density of States

The LLR method is a novel algorithm that enables us to evaluate the density of states in lattice gauge theory. We present our study of the ergodicity properties of the LLR algorithm for the model of Yang Mills SU(3). We show that the use of the replica exchange method alleviates significantly the topological freeze-out that severely affects other algorithms.


Introduction
Monte Carlo (MC) importance sampling simulations are at the basis of lattice investigations and are one of the most powerful methods to study non-perturbative phenomena in quantum field theory. Despite their success, there is a number of known cases where they perform very poorly, for example in presence of an overlap problem, i.e. when the measure distribution induces on the update algorithm a poor sample of the portion of phase space relevant for the evaluation of the observables, or in case of quantities that cannot be expressed as expectation values over a probability measure like partition functions and free energies. Over the years many algorithms has been proposed that try to overcome some of these problems. In this work we investigate the ergodicity properties of a recently proposed method, the LLR method.

The method
Our method is based on the evaluation of the density of states, and has been recently proposed in [1,2]. The density of states ρ(S) is defined via the functional integral where φ is some generic dynamical field variable, [Dφ] is the path integral over the field configurations and S [φ] the action of the Euclidean quantum field theory. From the density of states it is possible to Speaker, e-mail: antonio.rago@plymouth.ac.uk compute the partition function and expectation values of energy dependent observables via a simple one dimensional integration: while a bit more of care needs to be payed to evaluate generic observables (see sect.2.1). The conceptual difference of this formulation with respect to classical Montecarlo techniques is that all the path integrals are evaluated on a micro-canonical ensemble, i.e. at fixed action.

The LLR method in a nutshell
The LLR algorithm gives access to a controllable approximation of the ρ(S) featuring strong convergence properties. The procedure can be summarised in the steps below 1 .
For a given action range, partition the interval between S min and S max in N subintervals and call S i the center action value and δ S = (S max − S min )/N, for evenly spread actions intervals 2 . For each of the intervals evaluate and from the knowledge of the coefficients a i reconstruct the log-linear approximation of the density of states as It can be shown that our approximation converges to the correct function: ρ(S) = ρ LLR (S) e c δ 2 S "almost everywhere" for δ S → 0 (the ρ(S) is supposed to be almost everywhere C 2 ).
To evaluate the a i we first define: where the double brackets . . . represent the insertion of a strongly localised support function centred around the S i , in the present case a gaussian function. The a i are the zeros of the stochastic equation The stochastic nature of the latter is due to the determination of the double bracket expectation value, which can be accessed only through Monte Carlo estimates. The numerical solution of stochastic equations can be found using the iterative Robbins-Monro procedure [3]. This begins with a guessed value a (0) i , which is updated using the iteration and converges in L 2 , and hence in probability, to root of the equation: for more details see [3].
Once the a's are computed the density of states is reconstructed using eq.5 with error-bars determined using bootstrap resampling.
A few remarks are in order: the ρ LLR is an approximation of the correct density of state O(δ 2 S ) and it shows exponential error suppression: the relative approximation error does not depend on the magnitude, as a consequence this method works over several orders of magnitude. Once the a i have been evaluated, energy observables can be evaluated through direct integration of the density of states while generic observables can be evaluated through and also in this case the convergence in δ S is O(δ 2 S ). Note also that being the support function analytic this method is amenable to HMC simulation at practically no additional coding and computational cost, see [4].
To summarise for each step of our iteration we have to evaluate our observables over configurations distributed according to the weight Once the Robins-Monro procedure has converged we can reconstruct the observables of the full theory by an appropriate unidimensional integration. It is however clear that the strong localisation in action values imposed by the presence of the gaussian could slow down the update dynamic. Indeed the probability of visiting states with action far from the peak of the gaussian will be parametrically small in δ S and this will lead to a slow dynamic of the Markov Chain. The proposed solution is to simultaneously simulate multiple overlapping intervals with fixed central action and periodically propose a swap of the configurations belonging to pairs of them with probability: Such step preserves the detailed balance of action of the entire system and hence the resulting algorithm is still ergodic. Subsequent exchanges allow any configuration sequence to travel through all the action intervals, hence overcoming any potential action barrier. In the following we will refer to this method as replica 3 exchange, however it has to be noted that in literature it has also been referred to as umbrella sampling or parallel tempering. This method has been applied to the study of systems with strong metastabilities like the q-state Potts model at large q [5]. In this work we want to investigate the behaviour of the algorithm for theories with topological sectors by studying autocorrelation time of observables sensitive to these sectors.
We studied a pure SU(3) Yang Mills model in d = 4 with Wilson action. To simulate the model we used an HMC update, with Molecular Dynamics trajectory length τ = 1, we used a 2 nd order Omelyan integrator tuned to obtain an acceptance ∼ 98% for all our runs. For the choice of the MC parameter and a direct comparison of our data we refer to [6]. Our simulations were performed for a range of action values corresponding to 5.789 ≤ β ≤ 6.2 or 0.140 f m ≥ a ≥ 0.068 f m for fixed number of lattice points V = 16 4 hence with a extension ranging from 2.2 f m to 1.1 f m.
The observables that we measured are: • Action observables: observables that can be written as a function of the only action (these are a byproduct of the LLR method).
• The flowed action density E defined through the average of the plaquettes evaluated over the link evolved with the Wilson Flow at flow time t and its more symmetric clover definition E clover [7].
• The topological charge, using a bosonic estimator computed along the Wilson flow: where G t,µν (x) is again the clover term for the field strength tensor on the lattice constructed from the gauge links at flow time t [8]. For all the above observables we use the Madras-Sokal definition of the autocorrelation time [9], implemented according to Refs. [10,11].Γ(ξ) is the autocorrelation as function of the time lag ξ between measurements: where N is total number of measurements, B i denotes the i-th measurement of the observable B, and B its average. The integrated autocorrelation time τ int is defined as: where W is the size of the Madras-Sokal window. The swapping probability of two replicas will depend on the size of their overlapping probability distributions, hence typically it will depend on S i and on the fluctuation size of the constrained evolution σ LLR = S 2 i − S 2 i . In order to tune the swapping probability among different replicas we noted that it is possible to relate the σ LLR to the width of the support function (δ S ) and to the variance of the unconstrained evolution σ for identical action values. This study is reported in figure  2 where the asymptotic behaviour both for small δ S (σ δ S ) and large δ S (σ LLR σ) can be clearly highlighted. An additional feature of this study is the independence within our numerical precision of the scaling function with respect to the central action value. Thanks to this feature, it becomes possible, with the only knowledge of the variance of the unconstrained action as function of β, to generate a sequence of S i and (δ S ) i for which the overlapping spread of action among neighbouring replicas is equal. If we define the swap probability as the probability for fixed energy value to swap configurations with any other replica, it is possible to show that the tuning just described gives rise to a flat swap probability among all the replicas. By changing the size of the overlapping region it is possible to tune the swapping probability to the desired value, in figure 3 the case of probability ∼ 50%.

Tuning of the replica exchange method
A final note is in order: the Markov Chain process defined in this way is still distributed according the detailed balance of energy when considering the whole set of replicas, as such the autocorrelation time can be measured only on the observables reweighed over all the replicas together and not on a single replica at a time.

Results
We have performed an extensive investigation of the properties of our algorithm and compared our findings against the results of traditional HMC simulations for all the observables described in sec.3, together with their autocorrelation times. For a rough cost budget we have generated 0.7 · 10 5 configurations per replica and divided the action interval over 26 replicas, while for each traditional simulation we have generated 1.1 · 10 5 configurations. We produced estimates of our observables only in an internal region defined such that the reweighed observables could have a contribution from the boundary 10 −4 times smaller than the central contribution. All our estimates agree within the statistical error with their reference comparisons, and depending on the observables the statistical error associated with the observables measured with the LLR approach is at least a factor 2 smaller then traditional simulations. More importantly the autocorrelation time of all the investigated observables, once the replica exchange is taken into account, is smaller than the equivalent traditional simulation, signalling a good ergodicity property of the algorithm.

Conclusions
In this work we have investigated the ergodicity properties of the LLR method to evaluate the density of state. We focused in our study on the SU(3) pure Yang Mills theory, to test the behaviour of our method with theories that possess topological vacua. We measured observables sensitive to the topological structure alongside other less sensitive observables, and compared our findings with equivalent investigation done with traditional techniques. We have shown that the use of a replica exchange technique alleviates the problem of slow update dynamic evolution also in presence of the topological vacua. We have shown that the LLR method is capable of comparable precision in all the estimates at a comparable resource cost once the autocorrelation is taken into account. Simulation are still ongoing to investigate the scaling properties of this algorithm towards the continuum limit.