Comparison among HB-inspired algorithms for continuous-spin systems and gauge fields

Comparison among HB-inspired algorithms for continuous-spin systems and gauge fields

A. Cucchieri; R. B. Frigori; T. Mendes; A. Mihara

Instituto de Física de São Carlos – Universidade de São Paulo C.P. 369, 13560-970, São Carlos, SP, Brazil

ABSTRACT

We propose a new local algorithm for the thermalization of n-vector spin models, which can also be used in the numerical simulation of SU(N) lattice gauge theories. The algorithm combines heat-bath (HB) and micro-canonical updates in a single step - as opposed to the hybrid overrelaxation method, which alternates between the two kinds of update steps - while preserving ergodicity. We test our proposed algorithm in the case of the one-dimensional 4-vector spin model and compare its performance with the standard HB algorithm and with other HB-inspired algorithms.

Keywords: Monte Carlo algorithms, Heat-bath; Overrelaxation; Critical slowing-down

I. INTRODUCTION

The lattice formulation of quantum field theories provides a first principles approach to investigate non-perturbative aspects of high-energy physics. In this formulation the theory becomes equivalent to a statistical mechanical model, which can be studied numerically using Monte Carlo simulations (see for example [1,2] and references therein). As a consequence, the system considered evolves according to a Markov process in the so-called Monte Carlo time and the action-weighted configuration-space average of the observables is substituted by a time average over successive (independent) field configurations of the system. The possibility to use Monte Carlo simulations to study the theory nonperturbatively is especially important in the case of quantum chromodynamics, the nonabelian SU(3) gauge theory describing the strong interactions between hadrons. These simulations are computationally very demanding and must be done using local thermalization algorithms, since global methods (such as cluster algorithms) do not work well in this case.

We should notice that, based on the above-mentioned equivalence, when the continuum limit of the lattice quantum field theory is taken, the corresponding statistical mechanical model approaches its critical point. More precisely, the distance from the critical point is given by the lattice parameter b (directly related to the lattice bare coupling constant), which corresponds to an inverse temperature. For nonabelian gauge theories (as a consequence of asymptotic freedom) the continuum limit is given by b ® ¥, i.e. the critical point corresponds to temperature zero. The expected critical behavior is therefore similar to the one of a two-dimensional continuous-spin model (e.g. the classical Heisenberg or n-vector spin model with n > 2) or of a one-dimensional spin model, which show a critical point only at temperature zero, or b ® ¥.

The process of obtaining independent field configurations is called thermalization and is usually carried out by applying at each link of the lattice a local algorithm, such as Metropolis or heat bath (HB). When a critical point is approached, this process is afflicted by the well-known phenomenon of critical slowing-down (CSD) [3], which increases the correlation among successive field configurations. This implies that the integrated auto-correlation time t_int increases as a power of the lattice side N. In particular, for the Metropolis or HB algorithms one has t_int ~ N², i.e. the dynamic critical exponent z is equal to 2. Since statistical Monte-Carlo errors are proportional to , numerical simulations become increasingly inefficient close to a critical point. In order to reduce the problem of CSD one can combine the standard Metroplis and HB algorithms with so-called micro-canonical updates, allowing larger jumps in the configuration space and therefore improving the generation of independent samples. This is the idea behind the so-called hybrid overrelaxation method [4]. In general, adding a few micro-canonical sweeps greatly reduces CSD and, correspondly, the computational work.

A modification of the heat-bath algorithm, called overheat-bath (OHB), was introduced some years ago in Ref. [5]. The basic idea was to incorporate a micro-canonical move directly into the heat-bath step, thus reducing the computational cost while preserving the large moves in the configuration space. As it turns out, combining the two moves (heat-bath and micro-canonical) in a single step leads to a significant improvement in performance when compared to the hybrid overrelaxation method described above, which is based on alternating the two kinds of update moves. The resulting algorithm was indeed able to speed up the thermalization process, but, as already stressed in Ref. [5], it is not clear if it preserves ergodicity (especially when working at small temperatures). The OHB algorithm is used today in numerical simulations [6], usually combined with other algorithms in order to ensure ergodicity. In this work we propose a modification of the overheat-bath algorithm, which we call the modified heat-bath algorithm (MHB), incorporating a micro-canonical move into the heat-bath step without compromising ergodicity. In order to test the MHB algorithm — and compare its performance with the standard HB algorithm and the OHB algorithm — we consider the 4-vector spin model on a 1-dimensional lattice. Note that the model is exactly solvable, which makes it especially suited for testing the algorithm and comparing results to the exact solution.

Let us also note that, due to the isometry between the groups O(4) and SU(2), it is possible to study the SU(2) case with the same algorithm used for the 4-vector case. More precisely, the local update — i.e. the update of a single spin or gauge field variable while holding the rest of the lattice fixed — of the SU(2) lattice gauge theory is identical to the one for the 4-vector model. This does not mean, however, that the global update of an SU(2) gauge field configuration can be obtained from the corresponding update of a 4-vector model. Indeed, the latter update can be done very efficiently using the Swendsen-Wang-Wolff algorithm, while no such class of algorithms works well for lattice gauge theories. The method proposed here is therefore intended mostly for use in simulations of lattice gauge theories and not of n-vector models, although hybrid overrelaxation methods have also been applied in large-scale simulations of n-vector models [7]. Of course, an efficient thermalization in the SU(2) case is of great physical interest, since the quenched SU(N_c) case (for N_c > 3) is usually studied applying the SU(2)-embedding technique introduced in Ref. [8]. The SU(2) embedding is also used for simulating other Lie groups, such as the Sp(2) and Sp(3) groups, in studies of the deconfinement phase transition [9]. Preliminary results for the 2d SU(2) case have been presented in [10].

II. THE 4-VECTOR SPIN MODEL AND THE ALGORITHMS

The 4-vector spin model (on a 1-d lattice) is defined by the Hamiltonian

where the spins S_x are four-dimensional unit vectors, b ~ 1/Temperature , N is the lattice side and · indicates a scalar product.

In the case of a local algorithm, one has to consider the contribution to the Hamiltonian due to a single spin S_x. This gives _ss = -b S_x · H_x + constant , where the "effective magnetic field" H_x is given by H_x = S_x-1 + S_x+1 . (Here we consider periodic boundary conditions.) The above equation can also be written as

where S_x and _x are now interpreted as SU(2) matrices in the fundamental representation and N_x = . Note that Eq. (2) is clearly analogous to the expression of the single-link action obtained by considering the contribution of a single link variable to the SU(2) Wilson action [10].

Using the single-side action (2) and the invariance of the group measure under group multiplication, one obtains the HB update [11]

Here, U = u₀ I + i å_j u_j s_j is an SU(2) matrix, I is the 2 ×2 identity matrix, s_j are the three Pauli matrices and u₀ is randomly generated according to the distribution

At the same time, the vector - with components u_j normalized to - points along a uniformly chosen random direction in three-dimensional space [2].

In the hybrid over-relaxed algorithm one does m micro-canonical (or energy-conserving) update sweeps over the lattice, followed by one HB sweep. The micro-canonical update is a local deterministic transformation applied to the SU(2) matrix S_x, which does not change the value of the Hamiltonian. This is achieved by considering the update

As stressed in the Introduction, this update represents a large move in the configuration space, allowing the hybrid over-relaxed algorithm to reduce CSD at the price of a greater computational cost, due to the micro-canonical sweeps. Usually, by setting m = 2,3 one obtains a strong reduction in the value of t_int while increasing the computational cost by a factor smaller than 2.

In the OHB [5] one tries to include the micro-canonical step (5) directly into the heat-bath algorithm. To this end one generates u₀ according to the distribution (4) while the components of are not randomly chosen but are set using the relation = -, where W = w₀I + i· = . As a final step, the vector is normalized to . Clearly, the idea here is the same one applied in the standard hybrid over-relaxed algorithm: one tries to maximize the move in the configuration space by changing the sign of the component of W that is orthogonal to the effective magnetic field. (Note that the action S = -b N / 2 Tr W can be viewed as the action of a matrix W in an effective magnetic field given by the identity matrix I.) Clearly, this step does not obey the uniform distribution for but is a microcanonical move. Indeed we can think of the algorithm as a two-step process: a HB move followed by a microcanonical step. Thus, the algorithm is exact but may not be ergodic. We analyze the conditions for applicability of the OHB algorithm in a separate work [12], but it is clear that for some initial configurations the algorithm can get "trapped" inside a subset of the space of configurations, compromising the ergodicity condition. One such configuration is, for example, a "cold" start for an n-vector model, in which all spins are aligned along a fixed direction. In this case, the OHB is not able to change the lattice configuration at all. There is a clear problem also if the initial configuration in the SU(2)-lattice-gauge-theory case is given by variables in an Abelian subgroup. As can be easily seen, the update will change the initial configuration in this case, but the resulting Markov chain will remain restricted to the Abelian sector, without exploring the full space of configurations.

In this work we propose a modified HB algorithm (MHB) in which the generation of the SU(2) matrix U is done as in the HB case, but followed by one final step: if the scalar product · is positive then one does ® —, where the matrix W has been defined above. This may also be thought of as a modification of the overheat-bath (OHB). As said above, the basic idea, in both cases, is to incorporate a micro-canonical move into the heat-bath step. The difference is that in our case the direction of is randomly chosen and only its sign is modified, according to the above rule, while in the OHB algorithm one fixes µ - directly. Our modification should ensure ergodicity while keeping the advantage of having a micro-canonical step included into the HB step. We have indeed verified that the MHB algorithm has no problem with cold starts and shows a better t_int than that of the OHB algorithm for energy-related observables. On the other hand, in the OHB case, the move in configuration space is more optimized and the iteration cost is lower, since the algorithm is simpler and it needs fewer random numbers. Our results are presented in the next section.

III. RESULTS AND CONCLUSIONS

In order to study the CSD of the new algorithm we have to investigate if, and with what exponent, its relaxation time t_int (for certain quantities we are interested in) diverges as the lattice size N increases. To this end, we have to evaluate t_int for different lattice sides N at "constant physics", namely as the lattice size is increased, the physical size of the lattice should remain fixed. This is done by introducing a correlation length x and by keeping the ratio x/ N fixed, with N x in order to keep finite-size effects under control. For the 1-d 4-vector spin model one has x µ b [13], thus one should tune N µ b in order to keep the ratio x/ N constant. In our simulations we considered the lattice sizes N = 32, 64, 96, 128, 160, 224, 256 and b = 2.5 N / 32 and b = 5.0 N / 32 . In all cases we did 1.1 ×10⁶ thermalization sweeps and for the analysis we discarded the first 10⁵ sweeps.

In order to test the new algorithm we have evaluated the energy density e = N^-1

S_x· S_x+1 , the magnetization ² = ( S_x )² and the magnetic susceptibility c = N^-1á²ñ . Let us recall that there exists an exact solution [13] for the 1-d 4-vector spin model, even at finite lattice side N and with periodic boundary conditions. Thus, we can check the numerical results for the various quantities against the exact solution (see also [14,Table 1]).

In all cases we have evaluated the integrated auto-correlation time t_int using an automatic windowing procedure [3] with two different window factors (c = 4 and 8). We also employ a method based on a comparison between the naive statistical error with a binning error [15]. We checked that these three estimates are always in agreement. The error on the integrated auto-correlation time t_int has been evaluated using the Madras-Sokal formula [3].

Results are reported in Figs. 1 and 2. The data obtained for t_int have been fitted to the Ansatz t_int = aN^z in order to estimate the dynamic critical exponent z. Our study shows that the MHB algorithm has essentially the same critical exponent z of the standard HB algorithm but the value of the integrated auto-correlation time t_int is about 30 % smaller compared to the standard HB algorithm. This implies a reduction in the statistical error and in the computational cost by a factor of about 20%. Similar results have been obtained in the SU(2) case [10]. Thus, a good decorrelation can be reached without the use of multiple micro-canonical sweeps. In particular, from our simulations it appears that the HB algorithm with m micro-canonical steps is essentially equivalent to the MHB algorithm with m-1 micro-canonical steps. This feature may be useful in parallel computing [2] due to the reduction of the inter-node communication.

At the same time, from our data one sees that the MHB algorithm has a critical exponent z larger than that of the OHB algorithm. In order to reduce CSD for the MHB algorithm, while keeping its ergodicity properties, we are now testing an algorithm that interpolates between MHB and OHB.

A more extensive study of these algorithms, both in the O(4) and in the SU(2) cases, will be presented in a forthcoming work [12].

Acknowledgments

Research supported by FAPESP (Projects No. 00/05047-5 and 03/00928-1). Partial support from CNPq is also acknowledged (AC, TM).

Received on 30 September, 2005

• [1] T. Mendes, "Computational aspects of lattice QCD", arXiv:hep-lat/0309127.
• [2] A. Cucchieri, T. Mendes, G. Travieso, and A. R. Taurines, arXiv:hep-lat/0308005.
• [3] A. D. Sokal, "Monte Carlo methods in statistical mechanics: foundations and new algorithms", http://citeseer.nj.nec.com/sokal96monte.html, (1996).
• [4] U. Wolff, Phys. Lett. B 288, 166 (1992).
• [5] R. Petronzio and E. Vicari, Phys. Lett. B 254, 444 (1991).
• [6] R. Fiore, A. Papa, and P. Provero, Phys. Rev. D 67, 114508 (2003).
• [7] K. Chen and D.P. Landau, Phys. Rev. B 49, 3266 (1994).
• [8] N. Cabibbo and E. Marinari, Phys. Lett. B 119, 387 (1982).
• [9] K. Holland, M. Pepe, and U. J. Wiese, Nucl. Phys. B 694, 35 (2004).
• [10] R. B. Frigori, A. Cucchieri, T. Mendes, and A. Mihara, AIP Conf. Proc. 739, 593 (2005).
• [11] M. Creutz, Phys. Rev. D 21, 2308 (1980).
• [12] A. Cucchieri, R. B. Frigori, T. Mendes, and A. Mihara, "Improved heat-bath dynamics for continuous-spin systems and gauge theories", in preparation.
• [13] A. Cucchieri, T. Mendes, A. Pelissetto, and A. D. Sokal, J. Statist. Phys. 86, 581 (1997).
• [14] T. Mendes and A. D. Sokal, Phys. Rev. D 53, 3438 (1996).
• [15] G. de Divitiis et al., Nucl. Phys. B 437, 447 (1995).

Publication Dates

• Publication in this collection
  16 Oct 2006
• Date of issue
  Sept 2006

History

• Received
  30 Sept 2005