A novel density of state method for complex action systems
Recently, a new and efficient algorithm (the LLR method) has been proposed for computing densities of states in statistical systems and gauge theories. In this talk, we explore whether this novel density of states method can be applied to numerical computations of observables in systems for which the action is complex. To this purpose, we introduce a generalised density of states, in terms of which integrals of oscillating observables can be determined semi-analytically, and we define a strategy to compute it with the LLR method. As a case study, we apply these ideas to the spin model at finite density, finding a remarkable agreement of our results for the phase twist with those obtained with the worm algorithm for all explored chemical potentials, including values for which there are cancellations over sixteen orders of magnitude. These findings open new perspectives for dealing with the sign problem on physically more relevant systems.
A novel density of state method for complex action systems
\FullConferenceThe 32nd International Symposium on Lattice Field Theory
23-28 June, 2014
Columbia University New York, NY
1 Introduction and motivations
Monte-Carlo simulations of the theory discretised on a spacetime lattice provide a first principle method to compute observables in QCD at zero baryon density. At the heart of this approach is the interpretation of the Euclidean path integral measure as a Boltzmann weight. This relies on the measure itself being positive. However, as soon as a chemical potential coupled to the baryon number is switched on, the measure becomes complex and importance sample methods are no longer of use. In fact, because of the action being complex, the path integral has contributions alternating in sign that give rise to severe numerical cancellations. This phenomenon is called the sign problem and characterises not only QCD at finite density, but dense quantum systems in general. Other cases in which the sign problem hinder numerical simulations with importance sampling methods include interactions with external electromagnetic fields, rotating frames and real-time dynamics.
For numerical simulations of those systems, radically different approaches are needed. Significant progress has been achieved recently, with the introduction of several methods based on a wide variety of ideas (see [Aarts:2013bla, Gattringer:2014cpa] for recent reviews). Most of these algorithms are still in their infancy. In particular, both their numerical and mathematical properties need to be better understood. For this reason, it is particularly useful from an empirical point of view to compare them on the same models, in order to understand and assess their strengths. In this work, we give an overview of a recently proposed method [Langfeld:2014nta] based on a novel algorithm [Langfeld:2012ah] for the computation of the density of states. The underlying idea is that if a suitable generalised density of states is defined, the complex integral can be reduced to a one-dimensional oscillatory integral. If sufficient precision is available on the computation of the density of states, these types of integrals can be done with the desired accuracy. We remark that the general strategy is not new [Gocksch:1988iz] and has been revisited several times (see e.g. [Azcoiti:2011ei] and references therein). Our main original contribution is the new method for an accurate determination of the density of states over the orders of magnitude requested by the problem for a statistically meaningful evaluation of the integral.
In order to test our proposal, we have studied the spin model at finite density [Karsch:1985cb, Kim:2005ck], in which the sign problem disappears in a dual reformulation. The existence of a dual representation that is free from the sign problem enables one to devise an exact algorithm [Mercado:2012yf] that can be used to cross-check results obtained with our method. The results presented below extend our original work [Langfeld:2014nta], on which this contribution is mostly based.
2 The spin model
At strong coupling and for large fermion mass, for finite temperature and non-zero chemical potential QCD is described by the three-dimensional spin model
with the spin variable. The spin interaction is nearest-neighbour and is weighted by the coupling . The couplings and are related to the chemical potential and to the fermion hopping parameter . In the expression defining the partition function, indicates the complex conjugate of the spin-spin interaction. and are respectively the spin-spin interaction and the contribution weighted by .
At , the model reduces to the three-state Potts model, which is known to have a first order phase transition at . For the system is in the symmetric phase, while for the system is in a broken symmetry phase. At , the transition persists for a small non-zero , to terminate in a tricritical point. More details on the phase structure are provided in [Mercado:2011ua] and references therein.
At non-zero , the action is complex. In fact, while , for any . However, because of the symmetry , the partition function is real. This can be seen explicitly if we reformulate the problem in terms of the number of spins that are equal to each of the three cubic roots of unity . Introducing the variables
which fulfil the constraint ( being the number of sites), can be written as
and takes the form
with being the difference of the spins aligned along the roots with positive and negative imaginary part. Although real, contains an oscillating term, which is at the origin of the sign problem in this approach to the model.
The model admits a reformulation in terms of which there is no sign problem. This dual model has been used to construct an algorithm that allows to simulate the system without incurring in cancellations. The existence of an algorithm that can be trusted in the high- regime makes the model an ideal testbed for alternative approaches.
3 Computing the density of states
If we define a generalised density of state as
the partition function takes the form
being always positive, with this definition we have isolated the oscillating contribution from the non-oscillating one. Hence, if we are able to determine with sufficiently high accuracy, we can in principle try to do the oscillating sum. The needed accuracy is set by the simulation parameters. However, the stronger the oscillation, the higher is the accuracy requested to overcome the noise coming from cancellations of positive and negative contributions.
The Local Linear Relaxation (LLR) algorithm [Langfeld:2012ah] (see [Pellegrini:2014dva] for recent developments) has been proven to give an accurate determination for the density of states in gauge theories and in spin systems [Guagnelli:2012dk]. The method is based on a linear approximation of the logarithm of the density of states in a sufficiently small interval of variation of the independent variable, with the angular coefficient determined with a recursive relation. The full density of states can be reconstructed by imposing continuity of the piecewise approximations at the edges of the intervals.
Following [Langfeld:2014nta], we determine the generalised density of states by using a modification of this algorithm. Since implies , we only need to determine for . For obtaining this quantity, we formulate the ansatz
and use the LLR algorithm to determine the . The procedure goes as follows. We define -restricted expectation values of a function as
where for and otherwise. is a normalisation factor such that . The double bracket expectation values can be computed using standard Monte-Carlo methods. In particular, we can use and to obtain the in each interval starting from a trial . This can be achieved using the Newton-Raphson recursion
that has been adapted from [Langfeld:2012ah]. In this way, is determined up to the free parameter , which can be fixed by imposing the normalisation condition .
Since a Monte-Carlo averaging is involved in the recursion, the will be determined up to a statistical error. In order to account for this error in a realistic way, a bootstrap procedure has been used employing around 100 independent determinations of for each .
Fig. 1 (a) provides an example of a typical determination of the density of states. In order to check the correctness of the calculation, we have performed a simulation with the worm algorithm. The two methods agree in the whole range for which the histogram of the density of states can be reliably extracted from the latter simulation (up to ). The LLR method allows us to go well beyond this value (we stopped our determination at ), obtaining a density of states that spans well over 60 orders of magnitude (see Fig. 1 (b)).
4 Performing oscillating sums
In this section, we provide evidence that the density of states determined with the LLR method has sufficient accuracy for performing directly oscillating sums even in regions in which the sign problem is severe. The severity of the sign problem is measured by the expectation value of the phase factor . This is given by
Values of close to one mean that the sign problem is mild; conversely, means that the system is afflicted by a severe sign problem.
can be computed using a worm update at various and a snake algorithm [deForcrand:2000fi] to solve the overlap problem that arises when taking the ratio of partition functions at significantly different values of . This calculation is not afflicted by the sign problem. Within the LLR method, can be computed directly using the numerical determination of . We have performed a simulation at for . We have found that the approach of reconstructing directly from the numerical data for does not provide the required precision on the final result to make it statistically different from zero when the sign problem is significant. A better technique uses a polynomial interpolation of the logarithm of the density of states111A similar method has been used in [Azcoiti:2002vk].. More in details, we can write , where we have imposed the constraint that the logarithm of the density of states is even for . We have used interpolations up to , finding that the result is very stable and only and are significantly different from zero (within errors). All fits provide acceptable values of /dof. Using this semi-analytical procedure, we were able to obtain an agreement of the phase factor up to the maximum simulated value , where , indicating a strong sign problem. We refer to [Langfeld:2014nta] for details.
Another oscillating quantity is the phase twist
which in our formalism can be expressed as
Fig. 2 (a) shows a plot of the numerator and of the denominator in the definition of , both normalised by dividing them by . These two quantities vary over several orders of magnitude and following strong cancellations become of order at . Hence, even if their final ratio is of order , this is still a non-trivial test of how cancellations are resolved by our method. Fig 2 (b) shows good agreement between a semi-analytical determination with the fitted density of states and a direct Monte-Carlo simulation using a sign problem-free algorithm.
5 Discussion and conclusions
In this contribution, we have proposed a general method for simulating systems with the sign problem. The method consists in a precise determination of a generalised density of states using the LLR algorithm and a semi-analytical calculation of oscillating quantities. The method has been tested on the spin model, where it has been shown to reproduce the numerical results obtained with a dual algorithm, which does not suffer from the sign problem.
The same model has been studied with similar techniques in another contribution [Mercado:2014dva] on lattices for two different sets of parameters corresponding to the deep symmetric phase (which is the regime of our investigation) and to a situation in which the system is close to the phase transition. Ref. [Mercado:2014dva] confirms our conclusions deep in the symmetric phase. Near the phase transition, the authors observe good agreement for the phase twist obtained with the density of states and with the dual algorithm in a wide range of , with some deviations appearing at very high values of . These deviations could be due to the fact that the proposed ansatz for the density of states might be not appropriate in that regime or to a loss of efficiency of the dual algorithm. In any case, the discrepancy needs to be understood (and resolved) by performing dedicated simulations.
We thank V. Azcoiti, Ph. de Forcrand, C. Gattringer, J. Greensite, Y. Mercado, R. Pellegrini, A. Rago and P. Törek for discussions. This work is supported by STFC under the DiRAC framework. We are grateful for the support from the HPCC Plymouth, where the numerical computations have been carried out. KL is supported by the Leverhulme Trust (grant RPG-2014-118) and STFC (grant ST/L000350/1). BL is supported by STFC (grant ST/G000506/1).
Aarts:2013blaG. Aarts, Complex Langevin dynamics and other approaches at finite chemical potential, PoS LATTICE2012 (2012) 017, [arXiv:1302.3028].
Gattringer:2014cpaC. Gattringer, New developments for dual methods in lattice field theory at non-zero density, PoS LATTICE2013 (2013) 002.
Langfeld:2014ntaK. Langfeld and B. Lucini, The density of states approach to dense quantum systems, Phys. Rev. D to appear (2014) [arXiv:1404.7187].
Langfeld:2012ahK. Langfeld, B. Lucini, and A. Rago, The density of states in gauge theories, Phys.Rev.Lett. 109 (2012) 111601, [arXiv:1204.3243].
Gocksch:1988izA. Gocksch, Simulating lattice QCD at finite density, Phys.Rev.Lett. 61 (1988) 2054.
Azcoiti:2011eiV. Azcoiti, E. Follana, and A. Vaquero, Progress in numerical simulations of systems with a vacuum like term: The two and three-dimensional Ising model within an imaginary magnetic field, Nucl.Phys. B851 (2011) 420–442, [arXiv:1105.1020].
Karsch:1985cbF. Karsch and H. Wyld, Complex Langevin Simulation of the SU(3) Spin Model With Nonzero Chemical Potential, Phys.Rev.Lett. 55 (1985) 2242.
Kim:2005ckS. Kim, P. de Forcrand, S. Kratochvila, and T. Takaishi, The 3-state Potts model as a heavy quark finite density laboratory, PoS LAT2005 (2006) 166, [hep-lat/0510069].
Mercado:2012yfY. D. Mercado, H. G. Evertz, and C. Gattringer, Worm algorithms for the 3-state Potts model with magnetic field and chemical potential, Comput.Phys.Commun. 183 (2012) 1920–1927, [arXiv:1202.4293].
Mercado:2011uaY. D. Mercado, H. G. Evertz, and C. Gattringer, The QCD phase diagram according to the center group, Phys.Rev.Lett. 106 (2011) 222001, [arXiv:1102.3096].
Pellegrini:2014dvaR. Pellegrini, K. Langfeld, B. Lucini, and A. Rago, The density of states from first principles, PoS LATTICE2014 (2014) 229.
Guagnelli:2012dkM. Guagnelli, Sampling the density of states, arXiv:1209.4443.
deForcrand:2000fiP. 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].
Azcoiti:2002vkV. Azcoiti, G. Di Carlo, A. Galante, and V. Laliena, New proposal for numerical simulations of theta vacuum - like systems, Phys.Rev.Lett. 89 (2002) 141601, [hep-lat/0203017].
Mercado:2014dvaY. D. Mercado, P. Torek, and C. Gattringer, The Z3 model with the density of states method, arXiv:1410.1645.