Measuring the Decorrelation Times of Fourier Modes in Simulations

R. Ben-Av^{1}

Physics Department, Princeton University, Princeton, N.J. 08540,USA and G. Bhanot^{2} Thinking Machines Corporation 245, First Street, Cambridge, MA 02142, USA and

Institute for Advanced Study, Princeton, N.J. 08540, USA

February 8, 1993

Abstract

We describe a method to study the rate at which modes decorrelate in numerical simulations. We study the XY model updated with the Metropolis and Wolff dynamics respectively and compute the rate at which each eigenvector of the dynamics decorrelates. Our method allows us to identify the decorrelation time for each mode separately. We find that the autocorrelation function of the various modes is markedly different for the ‘local’ Metropolis compared to the ‘non-local’ Wolff dynamics. Equipped with this new insight, it may be possible to devise highly efficient algorithms.

IAS preprint: IASSNS-HEP-93/2

E-mail: 1 - radi@puhep1.princeton.edu ; 2 - gyan@think.com

1 Introduction

In the last few years there has been a lot of progress in finding new Monte-Carlo dynamics with little or no critical slowing down. This progress was initiated by the paper of Swendsen and Wang [1] who described a novel cluster dynamics which was not in the universality class of local dynamics such as the Metropolis or heat-bath algorithms and had z < 1 for Potts models. Their results were later extended to the ϕ^{4} model [2] and to O(n) models [3]. For a recent survey see [4].

However there has not been much progress in understanding what makes these improved algorithms work. The features that would contribute to improvement are easy to state: the dynamics must work on ‘relevant degrees of freedom’. Unfortunately, these are difficult to isolate for most models. The Potts models are in a special, simple class where we expect clusters of spins to constitute the important degrees of freedom. However, even here, one does not know why the Swendsen-Wang algorithm has the dynamical critical exponent that it does. In our opinion, such an understanding would lead to the invention of better algorithms. One should note the work of Lee and Sokal [5] which established a lower bound on the dynamical critical exponent z for Ising models. However, many fundamental questions remain unanswered, such as: Is there a universality class of ‘non-local’ dynamics? What is the mechanism behind the success of acceleration algorithms? Does the new dynamics apply to general gauge models, models with frustrations, and if so, how?

In this paper we take a first step in this direction by describing a method to study the rate at which modes at different length scales decorrelate for any given dynamics. We use the XY model in 2-dimensions as a platform to describe our method. We introduce a large set of operators and measure their time-time autocorrelation function for Metropolis and Wolff dynamics. Translation invariance allows us to identify Fourier modes as eigenvectors of the autocorrelation matrix. We compute the eigenvalues of the Fourier modes as a function of time. These directly measure the rate at which the various modes decorrelate.

Our results show that there is a marked difference between the two methods. In the Metropolis case, high frequency mode converge much faster than low frequency modes and it is the low frequency modes that govern the slow decorrelation in this dynamics. For the Wolff update, all modes decorrelate much faster than the Metropolis update. However, some high frequency modes of the Wolff dynamics are slower than some low frequency modes. This latter phenomena, to our knowledge, has never been demonstrated before.

Because our method provides a direct way to study the rate at which modes equilibrate, it could be used to both find the slow degrees of freedom and also to check the effectiveness of different algorithms in accelerating them. We expect that our method will help in fine tuning the search for effective algorithms.

2 The Model

The action of the two dimensional ferromagnetic XY model is:

E = ∑_{<i,j>} cos(θ_{i} − θ_{j})  (1)

where i = (n_{x}, n_{y}), 1 ≤ n_{x} ≤ L_{x} and 1 ≤ n_{y} ≤ L_{y}. The lattice is square with periodic boundary conditions, β is the inverse temperature. This model has the famous K-T phase transition [6] and has been extensively investigated numerically [7, 8, 9].

The ‘local’ dynamics we use is a parallel Metropolis algorithm. The dynamics is defined as follows: The lattice is decomposed into a red-black checkerboard of even and odd sites. All the even (odd) spins are updated simultaneously. A sweep is an update of the even spins followed by an update of the odd spins. The Wolff dynamics is implemented in its single cluster version [3]. The origin of the cluster is chosen randomly and this fact makes the dynamics translation invariant. For the moment, we will define one sweep to be one cluster update.

3 Autocorrelation Function

Let us define the elements of the autocorrelation matrix M(t) to be:

M_{i,j}(t) =< cos(θ_{i}(0))cos(θ_{j}(t)) > − < cos(θ_{i}) >^{2}  (2)

where i = (n_{x}, n_{y}) as above. In various limits, M_{i,j}(t) is related to static quantities. M_{i,j}(0) is the static correlation function of the observables cos(θ_{i}) and hence it is independent of the dynamics. ∑_{i,j} M_{i,j}(0) is the magnetic susceptibility times the volume. M_{i,j}(∞) is zero. For t between 0 and ∞, M_{i,j}(t) depends strongly on the dynamics. It is always easier to analyze the case where the dynamics is translational invariant. When this is true, M_{i,j}(t) = M_{i+d,j+d}(t) where d⃗ is any lattice vector. In this case it is obvious that the eigenvectors of M_{i,j} are the Fourier modes for any t.

It is clear that ‘typewrite order’ Metropolis dynamics in two dimensions is not translationaly invariant [10]. The parallel checkerboard Metropolis algorithm we use is translation invariant for even k⃗; that is, for k⃗ = (k_{1}, k_{2}) such that both k_{1} and k_{2} are even. However, it is easy to define a set of variables for which the dynamics is translation invariant. These are block spin variables which are the sums of cos(θ_{i}) on the four sites of a 2 × 2 lattice square. We will compute M_{i,j} for these variables. The Wolff dynamics we use is obviously translation invariant. We also use the same block variables for this dynamics to make the comparison between the two straightforward.

Let us denote the eigenvalue of M(t) by λk⃗(t) where k⃗ is a vector in the first Brillouin zone. For each k⃗ one can assign a decorrelation time τk⃗ via:

λk⃗(t) ∼ e^{−t/τk⃗} , t → ∞  (3)

4 The Results

We have measured the autocorrelation eigenvalues λk⃗(t) for both the Metropolis and Wolff dynamics for lattices with L=16 and β = 0.81, 0.83, 0.85, 0.87, 0.89, 0.91. Note that for ξ(β = 0.91) = 20 > L [8, 9] and so this temperature is already in the ‘spinwave’ phase for this lattice size. We did five sets of 120, 000 sweeps each, discarding the first 20,000. We took 10, 000 measurements of the autocorrelation function of Eq.2 with ten sweeps between each measurement. The reason we were limited to small lattices is that it is quite time consuming to measure the correlation function for all pairs i, j in Eq. 2 (which is what we did). If one knows the eigenvectors it is of course possible to measure only the autocorrelation matrix for them directly. This results in a much smaller matrix (in our case, it reduces the matrix M from 64 × 64 to 8 × 8). This is the preferred approach for future analysis.

In Fig. 1a-f we show λ_{k} for both dynamics as a function of update time T for each β value. The time T is the number of updates of both colors of the checkerboard in the case of the Metropolis algorithm. For the Wolff algorithm, the appropriate definition for T is the number of cluster updates times the average cluster size divided by the number of lattice sites. We have only plotted the data for T ≤ 40. In reality, our measurements extended out to T = 100. We plot only the first three modes corresponding to k⃗ = (0, 0), k⃗ = (0, 1) and k⃗ = (0, 2) where we measure k⃗ in units of 2π/L. The error bars on the correlation functions were estimated by repeating each measurement up to five times and computing the standard deviation.

We fit the data for λ_{k} to Eq.3 in the region where T was large enough for this relationship to be valid. From these fits, we extract the autocorrelation time τ which is shown as a function of β in Figs 2, 3a, and 3b.

5 Discussion

First it is clear from Fig.1 that the Wolff algorithm is much superior to the Metropolis algorithm. All Fourier modes decorrelate faster at all temperatures for it. It is also evident from Figs.1 and 2 that for the Wolff algorithm, the zero mode is much faster than the other modes. This unusual fact can be seen most clearly in Fig.2 where the τ for the zero mode is the smallest for all β values. One also sees in Fig. 1 that the modes exhibit different short time behavior. In Fig.2 we clearly see another interesting property of the Wolff dynamics. The autocorrelation time is generally larger for shorter distance modes. This should be contrasted with the opposite behavior of the ‘local-Metropolis’ where the shorter distance modes (modes with higher k) generally decorrelate faster than longer distance modes (modes with lower k). This pattern is what one would expect given the fact that the Metropolis algorithm is local and so should accelerate short distance modes and the Wolff algorithm is more global and so should do a better job for longer distance modes. However, the rich pattern that is evident in Figs.2,3 is not completely intuitive. For example, in Fig3b, one would expect the k = (0, 3) mode to be faster than the k = (0, 1) mode but the data clearly shows that the reverse is true. On the other hand, the relative order of the modes (0, 2) and (0, 1) is what one would expect for a local algorithm.

One might be able to use the information in Figs. 1-3 to arrive at an algorithm that can take advantage of the best features of the Metropolis and Wolff dynamics. This would involve doing some number n_{1} of Metropolis sweeps followed by some number n_{2} of Wolff sweeps. By monitoring the dynamics of a few modes as we have done in this paper and changing n_{1} and n_{2} one would arrive at the optimum algorithm. One can obviously also add in other methods such as the multigrid method [9] and the over-relaxation method [11] to get a truly hybrid algorithm. In fact, it may be that one need not actually have to adjust the n_{i}’s to find the optimum combination if the function λk⃗ is known for all t and for all modes. One might then be able to do an analytic calculation to find the optimum n_{i}’s.

It is well known that there are other degrees of freedom which are important in the XY model; viz., the non-local vortex excitations [6]. We believe that for a complete treatement of this model one should include the vortices as well as spins in the definition of the auto-correlation matrix M(t). In fact, we have tried to do this. However, we found that the contribution of the vortex-vortex and spin-vortex part of M(t) was insignificant when compared to the spin-spin part in Eq. 2. This comes from the fact that the vortex density is almost everywhere zero and any non-trivial vortex contribution to M(t) is swamped by the spin contribution. One might study the vortices by themselves and see what can be learned from them before attempting to get information about the correlation between the dynamics of vortices and spins. However, the issue of the relative normalization of spin and vortex densities remains.

We hope to take up the issues discussed in the last two paragraphs in a future publication. Acknowledgement

The research of GB was partly supported by U.S. DOE Grant DE-FG02-90ER40542 and the Ambrose Monell Foundation. We thank Thinking Machines Corporation for lots of computer time on their CM-2 and CM-200 supercomputers. We thank Herbert Neuberger for discussions and Steve Adler for a critical reading of the manuscript.

References

[1] R.H. Swendsen and J.S. Wang, Phys. Rev. Lett. 58, 86 (1987)

[2] R. C. Brower and P. Tamayo, Phys. Rev. Letts 62(1989) 62.

[3] U. Wolf, Phys. Rev. Lett. 62, 361 (1989); Nucl. Phys. B 322, 759 (1989); Phys. Lett. B 228, 379 (1989).

[4] A. D. Sokal, Nucl. Phys. B20 (Proc. Supp.) (1991) 55.

[5] X.J. Li and A.D. Sokal, Phys. Rev. Lett. 63, 827 (1989); X.J. Li and A.D. Sokal, Phys. Rev. Lett. 67, 1482 (1991).

[6] J. M. Kosterlitz and D. J. Thouless, J. Phys. C 6, 1182 (1973).

[7] J. Tobochnic and G. V. Chester, Phys. Rev. B 20, 3761 (1979).

[8] R. Gupta and C. Baillie, Phys. Rev. B45 (1992) 2883.

[9] A. Hulsebos, J. Smit, J. Vink, Nucl. Phys. B 356, 775 (1991); R.G. Edwards, J. Goodman, A.D. Sokal, Nucl. Phys. B 354, 289 (1991)

[10] H. Neuberger, private communication.

[11] S. L. Adler, Phys. Rev. 23 2901 (1981); C. Whitmer, Phys. Rev. D29 (1984) 306.

Figure Captions

• Figure 1a-f: λk⃗ against the time T to update L_{x} × Ł_{y} spins for different temperatures. For each k⃗, the slower dynamics is Metropolis and the faster is Wolff.

• Fig2: Decorrelation time τ as a function of β for various k⃗ for the Wolff algorithm. Only the error bars on the least accurate data are shown. Note that the zero mode decorrelates the fastest (has the smallest τ).

• Fig3a: τ versus β for the zero mode for the Metropolis algorithm. Note the relative scale of this figure and Fig.2.

Fig3b: τ versus β for some other modes for the Metropolis algorithm. Again, only one set of error bars are shown to give a feeling for the accuracy of the data.

8 9 10 11 12 13 14 15 16