Framework and properties of cluster algorithms

I~l[|llM||l|la[IN'!

PROCEEDINGS SUPPLEMENTS

Nuclear PhysicsB (Proc. Suppl.) 30 (1993) 289-292 North-Holland

Framework and properties of cluster algorithms Werner Kerler Fachbereich Physik, Universit~it Marburg W-3550 Marburg, Germany A general framework for cluster algorithms is given, considering the available freedom and implications for the transition spectra. A numerical study of Ising systems, allowing to get hold of the eigenvalues and of their weights, leads to further rules governing cluster mechanisms and reveals various dependences on the lattice size.

1. G E N E R A L F O R M U L A T I O N In the transition probability

W(o., o.') = Z A(o', n)B(n, a, o.')

(1.1)

{n}

one gets detailed balance. Using the joint distribution u j ( n , cr) = U(o.)A(o., n) and the marginal distribution iS(n) = ~{o} U~,(n, o.) one may introduce the special conditioned probability

for A, describing the generation of bonds between variables, as well as for B, responsible for the change of the values of the variables, optimal choices are to be found. With probabilities Pa,~, where na = 0 denotes deleting and na > 0 freezing, one gets

A(n, o.) - us(n, o.)

A(o., n) -- H p a , o ( o . )

for which detailed balance follows immediately without further assumptions. If a projection function 0(a, n) exists such that ~(o., n)A(a, n) = A(o., n), for the choice

(1.2)

The index a can be determined by any subdivision

H(o.) -- Z H a ( o . )

(1.3)

ot

in U(o-) = Z - l e x p ( - H ( o - ) ) . The bonds occur within the set of variables in Ha (if this are more than two also na > 1 becomes possible). Imposing in A the condition Pa0(o.) = exp(Ha(o.) - g a - An)

(1.4)

(1.8)

~(-)

and consider the particular transition probability

ITV(o.,o.') --" Z A(o., n)fi(n, o.')

B(n, o., o.') = tg(o., n)ft(n, o.') ,

e -H`'(°) = Z~an,,Dan,,(o.)

U(o.')A(o.', n) = U(o.)A(o., n)

into p(o.) to obtain the joint distribution

provided that Pa..(o.') = p.,,°(o.) for na > 0. This can be guaranteed in (1.1) by the condition

(1.10)

provided it satisfies (1.6) and (1.7), W reduces to Izd. This situation can be verified to occur in the original approach of Swendsen and Wang [1]. It should be realized that the probabilities Pan, for freezing and deleting are distinct from the percolation probabilities/~a,,o. The latter ones enter if one inserts a decomposition of type [2, 3]

where Ka = max{o} (Ha(o.)) and where An _> 0 are arbitrary constants, one obtains (1.5)

(1.9)

{,,}

(1.11)

n~

1

uj(n,o.) =-~HpanoDano(O.)

(1.12)

a

B(n, o., o.') = 0 if Pano,(o') ¢ Pano(o.')

(1.6)

for na > 0. Then by requiring

B(,, o", o.) = B(n, o3 o")

From this one gets the marginal distribution 1

(1.7)

~(n)=gH~an~7(n) a

0920-5632/93/$06.00 © 1993 - Elsevier Science Publishers B.V. All rights reserved.

,

(1.13)

290

W. Kerler~Frameworkand properties of cluster algorithnu

where 7(n) = E ~ H,~ D~no(cr), which describes a generalized percolation model. On the other hand, in terms of the p ~ , . the joint distribution reads ~

p.n.(a)

(1.14)

where ~ = I-Ia e x p ( - K , - Aa). The generation of joint distributions allows to consider more general observables. Ones of practical interest are related to the decomposition

B(n,a,a') = E F(n,z)C(n,z,o',a')

(1.15)

where {z} describes partitions of the clusters in two sets, with values of variables to be changed or to be kept. For observables of type f(n, z, a) one has to use the transition probabilities

Wc(n,z,a,n',z',~') = G(n, z, cr, a')A( a', n')F(n', z')

(1.18)

(1.19)

while for f(n, z, a) = f(n, z) in general no simplification occurs. To discuss spectral properties a Hilbert space with (f, g) = (fig) is appropriate [4]. For (1.18) in the case of f(n, z, a) or f(n, z)

(f' g) = E #c(n, z, o')f*(n, z, a)g(n, z, o') (1.20) is to be used while for (1.19)in the case o f f ( a ) (f, g) =

Z#(a)f*(a)g(o')

2. S T U D Y O F I S I N G S Y S T E M S In the special case with Aij -- 0 and

(1.17)

for f(n, z, a) = f(~) reduces to

u( gf(a)W(a,

with projectors Fv , eigenvalues Av , and ,~1 = 1. In the indicated metric detailed balance means self-adjointnes which implies real eigenvalues and positive weights. In addition, e.g., in the case of (1.9) only positive eigenvalues occur. In more general cases (described by a biorthogonal system) also conjugate complex pairs of eigenvalues are possible and the weights need no longer to be positive. The time dependences in contributions to (1.22) have the general forms e -t/r , e - q ~ ( - 1 ) t , and e-t/~cos(¢t + ¢).

one gets pij,o = 1 - 6o,,oj(1 - e -s~j) = 1 - Pij,1. Then A is proportional to the projection function O(a, n) = 1-I(ij):.m,=l 6o,,0, • This special feature allows to consiaer

f(n, z, a)Wc(n, z, a, n', z', a') f(n', z', a')

=

(1.22)

H(a) : ~(-Jij6~,,~j (#)

invariant, however, in general do not satisfy detailed balance. It is to be noted that

(ftf,+l) = E pc(n, z, o')

(f'Fvf)

(1.16)

which leave the joint distribution

pc(n, z, o') = .(o')A(o', n)F(n, z)

R(t) = (fufu+t)- (fu)~ = Z At

(1.21)

is more convenient. Spectral properties manifest themselves in autocorrelation functions

B(n, a, ~') =

+ Cij)

,

(2.1)

(2.2)

8j8 t

with 0c(a, s, n) = I-Ir=l N~ 1-1iel. 6,,,,,,, where s~ denotes the spin and L(n) the set of sites of the rth cluster. Thus it remains to investigate the freedom in the choices ofb with b(n, s, s') = b(n, s', s). Monte Carlo simulations with high statistics have been performed for Ising systems at the critical point in 2, 3, and 4 dimensions and with observables E, M, X, and C, where C is the sum of relative sizes of the clusters with spin flips. Also autocorrelations have been determined in each case. All simulations have been done for four choices of b, corresponding to the flipping rule of Swendsen and Wang [1] (SW), to the probability of Wolff [5] (WO), to flipping the spin of the largest cluster [6] (LC), and to flipping that of the largest one and of those not touching it [7] (SC). The vectorization of the program developed also extends to cluster labeling, (1) to parallel

W. Kerler ~Framework and properties of cluster algorithms

propagation of the label front in each cluster and (2) to a separate identification of small clusters using the bond numbers of sites, which introduces parallel features also with respect to different clusters. The analysis of the data is based on X2-fits of the normalized autocorrelation functions p(t) = R ( t ) / R ( O ) . The covariance matrix needed for this is obtained from p using that higher cumulants vanish. Stability against variations of the fit intervals is required and X2-values are observed to check reliability. The fit function rl(t) is used to extrapolate the estimator of the integrated autocorrelation time giving

=

1

M

+

p(t) + t=l

(2.3) t>M

It turns out that ri,~t and var(rin~) are practically independent of the window M. The fit functions used are ce - t / ~ for SC and SW, and ce - t l ~ q-c2e -t/ra for LC and WO. These main modes are seen for all observables though with different weights (except for C in the case of SW). Fast modes with small weights occur close to t = 0. They are studied considering p(t) - rl(t ) below the fit interval. It turns out that 1 - 7/(0) is smallest for E and smaller for X than for M. 3. D I S C U S S I O N

OF RESULTS

3.1. Spectral properties With respect to the form of the spectrum it turns out that only a few eigenvalues play a role. A broad spectrum with a threshold (as related to time dependences of type t - a e - t / r ) can be excluded. The rule is found that the weight of the leading eigenvalue e -1]r dominates if the pattern of clusters subject to spin flips is extended. This is reflected by the fact that for LC and WO a second main mode is observed the decay of which is not fast as compared to that of the leading one , while no such mode occurs for SC and SW. This rule can be understood by noting that the spectral representation relates changes of configurations to the spectrum in an analogous way as

291

is familiar for the conjugate variables in Fourier transforms. Similar trends are observed for the observables. The weight of the leading mode relative to that of the second one is larger for E than for the other observables. Further, for E the smallest fraction of fast modes occurs. Considering fast contributions for the observables E, X, and M, which are of type f ( a ) , due to detailed balance only real eigenvalues and positive weights are observed. For SW, being based on W, in addition only positive eigenvalues occur. On the other hand, for SC and LC also negative eigenvalues are possible and related oscillations are in fact seen for the observables M and X. In the case of fast contributions for the observable C, which is of type f ( n , z), one has to note that there is no detailed balance. For SW and WO contributions proportional to 6t,0 occur. For SC in d = 2 a relatively large negative contribution is found. For SC and LC superpositions of oscillating and decaying fast modes are seen. 3.2. L a t t i c e size d e p e n d e n c e s The dependences on the lattice size L have been analyzed for the quantities and fits to k L z calculated where possible. For the autocorrelation times in d = 3 and 4 within errors power laws are found. In d = 2 in general different behaviors are observed (logarithmic fits only possible for Vin × t and r~n Mt of SC and SW as well as a few cases with power laws appear to be exceptions). It turns out that also the coefficients depend on L. Various behaviors are observed (depending on observables and dimensions). Frequently power laws are seen. In general simultaneous power behavior for r, c, and tint cannot be expected (unless L is very large) because the term ce - t / r of p(t) contributes c(e I/~ - 1) -1 to tint. If the leading mode dominates such that ri~t ~ cv the power laws combine. Empirically one has primarily power behavior of r and c. For LC and WO (C) decreases with L according to laws k L z . For SC it first decreases and then increases. This indicates that the sum of the rel-

292

W. Kerler/ Framework and properties of cluster algorithms

ative sizes of the clusters not in contact with the largest one increases. Figure 1 shows the related behavior of the difference ( C ) s c - ( C ) L c and also

(c)Lc. 1

I

LC, d=2

I

I

I

I

I 20

I 50

I

LC, d=3 ,,,~,,,~+.

0.5

0.2 0.1 0.05 DC, d=3

0.02 0.01 -

0.005

0,002

DC, d=2

0.001

I 2

I 5

I 10

L

I 100

Figure 1. ( C ) s c - (C)Lc (denoted by DC) and (C)Lc in d = 2, 3, 4 as functions of L.

For the separate labeling of small clusters used here, with parallel features not possible in a onecluster approach, the gain is such that with SC one improves also with respect to the effective WO times. Further improvement is possible by extending the parallel labeling part. Thus the measured autocorrelation times are compared. By SC one gains a factor of about 2 with respect to SW (the gain is in k and the z-values are about the same). LC is between S C a n d SW, in higher dimension closer to SW. For larger L and d the autocorrelation times measured for WO are much larger than for SW. The rule is found that faster decay of autocorrelations is related to larger total size of the pattern of disconnected clusters with spin flips. For SC, LC, and WO this is immediately seen. For SW the effect comes into play that flipping the spin of the largest cluster simultaneously with spins of clusters touching it makes the decay slower. It becomes obvious that improvement is possible by avoiding simultaneous flips in touching clusters as much as possible. Clearly going from LC to SC then is only a very first step. One has also to account for connectivity among the clusters not in contact with the largest one. From the behavior of ( C ) s c - (C)Lc in Figure 1 it is seen that by this considerable improvement can be expected in particular for larger L and d, promising also a more favorable universality class. REFERENCES

1 3.3.

Decay

of autocorrelations

The time scale is such that each time step corresponds to a new cluster configuration. Thus all choices of B are treated on equal footing. For WO it is possible to grow only one cluster. To account for this effective autocorrelation times have been introduced [8] by multiplying the measured autocorrelation times by (C). The comparison with these effective times is appropriate if the cluster labeling takes most of the time and if it is done sequentially with respect to different clusters. Otherwise the measured autocorrelation times are to be compared.

2 3 4 5 6 7

R.H. Swendsen and J.-S. Wang, Phys. Rev. Lett. 58 (1987) 86. R.G. Edwards and A.D. Sokal, Phys. Rev. D 38 (1988) 2009. R.C. Brower and P. Tamayo, Boston University preprint, hep-lat/9206007 (1992). N. Madras and A.D. Sokal, J. Stat. Phys.50 (1988) 109. U. Wolff, Phys. Rev. Lett. 62 (1989) 361. C.F. Baillie and P.D. Coddington, Phys. Rev. B 43 (1991) 10617. W. Kerler, Nucl. Phys. B (Proc. Suppl.) 26

(1992) 628. 8

U. Wolff, Phys. Lett. B 228 (1989) 379.