\input /home/A10F/paulc/texlib/PRLmac.tex
% \input /net/jack/paulc/texlib/mac.tex
% \input /net/jack/paulc/texlib/moremac.tex
% \input /net/jack/paulc/texlib/fonts.tex
\input /home/A10F/paulc/texlib/peqns.tex

\font\big=cmr10 scaled\magstep1
\font\bigger=cmr10 scaled\magstep2

\line{C3P-945\hfill}
\line{SCCS-12\hfill}
\line{CRPC-TR90159\hfill}
\line{COLO-HEP-229\hfill}
\bigskip
\bigskip
\bigskip
% \centerline{\bf Comparison of Cluster Algorithms for two-dimensional 
% Potts Models}
\centerline{\bf \bigger{COMPARISON OF CLUSTER ALGORITHMS}}
\smallskip
\centerline{\bf \bigger{FOR TWO-DIMENSIONAL POTTS MODELS}}
\bigskip
\bigskip
\bigskip
\centerline{\big{Clive F. Baillie 
\footnote{$^1$}{\tenrm {Current address: Physics Department, University
of Colorado, Boulder CO 80309}}
and Paul D. Coddington}
\footnote{$^2$}{\tenrm {Current address: Physics Department, Syracuse 
University, Syracuse NY 13244}}
}
%
% \baselineskip=24pt
%
\bigskip
\bigskip
\centerline{Caltech Concurrent Computation Program,}
\centerline{California Institute of Technology,}
\centerline{Pasadena, CA 91125, USA}
\bigskip
\bigskip
\centerline{October, 1990}
\bigskip
\bigskip
\bigskip
\centerline{\bf Abstract}
\medskip
\noindent
We have measured the dynamical critical exponent $z$ for the
Swendsen-Wang and the Wolff cluster update algorithms, as well as a 
number of variants of these algorithms, for the $q=2$ and $q=3$ Potts
models in two dimensions. We find that although the autocorrelation
times differ considerably between algorithms, the critical exponents
are the same. 
For $q=2$, we find that although the data are better fitted by a 
logarithmic increase in the autocorrelation time with lattice size,
they are also consistent with a power law with exponent $z \approx 0.25$,
especially if there are non-negligible corrections to scaling.

%
%\bigskip
%\bigskip
%\noindent
%PACS numbers 05.50.+q, 11.15.Ha, 64.60.Ht.
%
\vfill
\eject

\bigskip
{\bf I. \quad INTRODUCTION}
\medskip

Cluster update algorithms can greatly reduce critical slowing down in 
computer simulations of spin models. The original idea of Swendsen
and Wang 
\ref\SW{R.H. Swendsen and J.-S. Wang, Phys. Rev. Lett. {\bf 58}, 86 (1987)}
(SW) was based on the work of Fortuin and Kasteleyn 
\ref\FK{P.W. Kasteleyn and C.M. Fortuin, J. Phys. Soc. Jpn. Suppl. 
{\bf 26s}, 11 (1969); 
C.M. Fortuin and P.W. Kasteleyn, Physica {\bf 57}, 536 (1972)},
which relates the Potts model 
\ref\Potts{R.B. Potts, Proc. Camb. Phil. Soc. {\bf 48}, 106 (1952);
F.Y. Wu, Rev. Mod. Phys. {\bf 54}, 235 (1982)}
to a random bond percolation model
\ref\Perc{D. Stauffer, Phys. Rep. {\bf 54}, 1 (1978); 
J.W. Essam, Rep. Prog. Phys. {\bf 43}, 830 (1980)}.
The $q$-state Potts model consists of a lattice of spins $\sigma_{i}$ 
which can take $q$ different values, and its Hamiltonian is
$$ 
H = K \sum_{<i,j>} \delta_{\sigma_{i}\sigma_{j}},
\eqn\Eqham
$$
where $K$ is a dimensionless coupling constant.
By introducing bonds between neighboring spins, with probability
$$ 
p(\sigma_i,\sigma_j) = \delta_{\sigma_i\sigma_j} (1 - e^{-K}),
\eqn\Eqbond
$$
clusters of bonded spins are created. The SW algorithm consists of generating
all such clusters, then choosing a random spin value for each cluster and 
assigning it to all the sites in that cluster. In the algorithm of Wolff  
\ref\Wolff{U. Wolff, Phys. Rev. Lett. {\bf 62}, 361 (1989)},
a site is chosen at random and a single cluster constructed around it. All 
the spins in this cluster are then flipped, i.e. changed to a single new 
Potts state (different from the current state of the cluster's spins) 
which is chosen at random from the $q-1$ possible values.
Note that this is different from the SW update, which may leave the spins 
unchanged in some clusters. These algorithms, which update whole clusters 
of spins at a time, result in a dramatic decrease in critical slowing down 
over local algorithms which update only one spin at a time.

As shown by Wolff 
\ref\WO{U. Wolff, Nucl. Phys. B {\bf 334}, 581 (1990)}, 
the SW and Wolff algorithms are both special cases of the following general 
cluster update algorithm:

\medskip
\noindent
Step 1. The lattice is decomposed into $N_c$ clusters $c_i$, using the bond
probabilities of equation \Eqbond.

\noindent
Step 2. A subset $C$ of the $c_i$ is chosen using some specified probability 
distribution, and all of the clusters of spins in $C$ are flipped
independently.

\medskip
\noindent
The SW algorithm corresponds to including each cluster in $C$ with probability
$(q-1)/q$, while the Wolff algorithm amounts to picking a single cluster with 
probability proportional to its size. Hence the only differences between the 
Wolff and SW algorithms are that the SW algorithm flips more clusters per 
iteration, while the Wolff algorithm flips relatively larger clusters. In order
to investigate these differences further, we also studied some other methods of
choosing the clusters to be flipped, including the following :
\medskip
\item{(1)} flip all the clusters (hereafter denoted by AC),
\item{(2)} flip one of the clusters, chosen at random with equal probability
for each cluster (denoted by RC), and 
\item{(3)} flip the largest cluster (denoted by LC). 
\medskip
\noindent
Note that flipping all of the clusters will not work in the case of the Ising 
model ($q=2$), since it is equivalent to a global spin flip, and is thus highly 
non-ergodic! However for $q>2$ it is a perfectly acceptable algorithm.

\bigskip
{\bf II. \quad SIMULATIONS}
\medskip

We have made a high statistics study of the above cluster algorithms for the
$q=2$ and $q=3$ Potts models in two dimensions. 
Preliminary results have been reported elsewhere \ref\init{P. D. Coddington 
and C. F. Baillie, Proceedings of the International 
Workshop `Lattice 89', ed. N. Cabibbo {\it et al.}, Nucl. Phys. B (Proc. Suppl.)
{\bf 17}, 305 (1990)}.
We measured the autocorrelations in the energy 
$E = - \sum_{<i,j>} \delta_{\sigma_{i}\sigma_{j}}$
and the absolute magnetization
$|M| = \vert \sum_j e^{2 \pi i \sigma_j / q} \vert$.
%
% DEFINITIONS OF ENERGY AND MAG (with refs for MAG?)
%
%autocorrelations in the energy, magnetization and susceptibility, and also in 
%the improved estimators of these quantities
%\Wolff 
%\ref\Sweeney{M. Sweeny, Phys. Rev. B {\bf 27}, 4445 (1983)}.
%
The normalized autocorrelation function for an operator A is
$$ 
\rho_A(t) = {<A(0)A(t)> - <A(0)>^2 \over <A(0)A(0)> - <A(0)>^2}. 
\eqn\Eqautocorr
$$
At large t this should decay as $e^{-t/\tau_{exp}}$, where $\tau_{exp}$ is the 
exponential autocorrelation time (independent of $A$ unless $A$ is orthogonal 
to the slowest mode).
The quantity which is actually relevant 
to the error in A is the integrated autocorrelation time 
\ref\Sokal{N. Madras and A.D. Sokal, J. Stat. Phys. {\bf 50}, 109 (1988)}
$$ 
\tau_{int,A} = {1 \over 2} + \sum_{t=1}^{\infty} \rho_A(t). 
\eqn\Eqintauto
$$
%If $\rho(t)$ is exactly exponential $\tau_{int} = \tau_{exp}$ to a good 
%approximation. 
The integrated and exponential autocorrelation times have associated exponents
$z_{int,A}$ and $z_{exp}$ governing critical slowing down. These are given by
$\tau(L) \sim L^{z}$ for both the exponential and integrated autocorrelations, 
where $L$x$L$ is the size of the lattice and $\tau(L)$
is measured at the infinite volume critical point.

Autocorrelations are traditionally measured between each update of the
entire lattice, so for the single cluster updates, where only a fraction of
the lattice sites are updated at each iteration, the autocorrelation time
needs to be scaled by a factor of the average cluster size $<|c|>$ divided 
by the lattice size, so that 
%$\tau = \tau' <|c|> L^{-d}$,
$$  
\tau = \tau' <|c|> L^{-d},
\eqn\Eqscale 
$$
where $\tau'$ is the unscaled autocorrelation time and $d$ is the lattice 
dimension (here $d = 2$).
For the Wolff algorithm, the average cluster size scales like the 
susceptibility, hence $<|c|_{W}> \sim L^{\gamma/\nu}$. The average cluster 
size for SW is (apart from finite size effects) independent of $L$, thus
$<|c|_{RC}> \sim L^0$. 
As in the percolation model \Perc,
the maximum cluster size (normalized to the size of the lattice) scales 
% like the magnetization, so that $<|c|_{LC}> \sim L^{2-\beta/\nu}$. 
like the magnetization, so that $<|c|_{LC}> \sim L^{d-\beta/\nu}$. 
From equation \Eqscale\ and the definition of $z$ we have that
$$ \eqalign{
% z_{W} &= z_{W}' - (2 - \gamma / \nu) \cr
% z_{RC} &= z_{RC}' - 2 \cr
z_{W} &= z_{W}' - (d - \gamma / \nu) \cr
z_{RC} &= z_{RC}' - d \cr
z_{LC} &= z_{LC}' - \beta/\nu \cr
}\eqn\Eqzscale
$$
where $z'$ is the unscaled exponent, $\tau'(L) \sim L^{z'}$. In our 
calculations the scaled value of $z$ is actually obtained by scaling the 
autocorrelations by the average cluster size, but we have checked that
using equation \Eqzscale\ instead gives consistent answers for LC and
Wolff. For RC, however, there are large finite size effects on the smaller 
lattices and the average cluster size does not scale well.

Note that the computational requirement for the Wolff algorithm is 
proportional to the average size of the clusters which are flipped \Wolff, 
but for the other single cluster algorithms (RC and LC) it is 
necessary to form all the clusters in order to choose the one to flip.
Hence the computational requirement for the RC and LC algorithms goes like 
$L^d$ (the same as for SW), and in that sense 
the above scaling is unwarranted. However it may be possible to invent
a clever algorithm which constructs single clusters with the required
probability distribution, and in any case the scaling is necessary to
allow a sensible comparison with the Wolff algorithm.

The autocorrelations were measured at the infinite volume critical point
(which is known exactly for 2-$d$ Potts models) over runs of typically 
$5\times10^5$ to $10^6$ sweeps ($10^6$ to $10^7$ sweeps for the RC algorithm),
and 10 or 20 such runs were done for each of the different lattice 
sizes, from $8^2$ up to $128^2$. The results 
from the multiple runs were averaged to give the overall autocorrelation
function $\rho_A(t)$. 
%
% Put number of sweeps in a table?
%
For the SW and Wolff algorithms, we also did runs on a $256^2$ lattice.
For SW we did a total of 7.5 million sweeps for $q=2$ and 3.2 million for $q=3$,
while for Wolff we did 7.5 million sweeps for $q=2$ and 5 million for $q=3$. 
For the SW algorithm with $q=2$ we also did 
a total of 5.6 million sweeps on a $512^2$ lattice.
In the SW case on these larger lattices we used a parallel cluster-finding
algorithm which we have developed 
\ref\parallel{C.F. Baillie and P.D. Coddington,
``Cluster Identification Algorithms for Spin Models --- Sequential and
 Parallel'',
to be published in Concurrency: Practice and Experience}.

The autocorrelation time $\tau_{exp}$ for an operator $A$ was extracted from 
a fit to $\log \rho_A(t)$ in a region where it was linear, typically 
somewhere between $\tau_{exp}$ and $3\tau_{exp}$.
%
%** Say something about Wolff not giving good fits ?? **
%
$\tau_{int,A}$ was calculated by summing $\rho_A(t)$ for $t$ less than some
cutoff value, 
and then using values from the fit to approximate $\rho_A(t)$ for 
the remaining part of 
the sum, which is then just a geometric series and so can be summed 
analytically. 
The cutoff was taken to be the largest point used in the 
fit, which was usually between $2\tau_{exp}$ and $3\tau_{exp}$. 
This method is similar to that used by Wolff
\ref\WIsing{U. Wolff, Phys. Lett. {\bf B228}, 379 (1989)}.
The above analysis was also done for each of the 10 or 20 multiple runs, which 
provided an estimate of the error in the autocorrelation times.

\bigskip
{\bf III. \quad RESULTS}
\medskip

We found as expected that the energy and absolute magnetization have very
similar values of $\tau_{exp}$, however $|M|$ decorrelates faster than $E$, 
giving smaller values of both $\tau_{int}$ and $z_{int}$.
%
%For the SW algorithm the improved estimators gave approximately the same
%autocorrelations as the standard values of the measured quantities. 
%For the Wolff algorithm, the improved estimator of the (unsubtracted)
%susceptibility is just the average cluster size. 
%This quantity seems to have approximately the same 
%$z_{exp}$ as the standard estimators, although $z_{int}$ is very much 
%reduced. In fact for $q=2$ $\tau_{int}$ actually decreases slightly after 
%the scaling of equation \Eqscale.
%Unfortunately this does not mean that the improved 
%estimator gives a greatly reduced error in the susceptibility, in fact at 
%criticality the error is greater than that of the standard estimator due to 
%the very large fluctuations in size of the clusters which are grown.
%
%
Results for $\tau_{int,E}$ for the different algorithms are
shown in Figs. 1 and 2 for $q=2$ and $q=3$ respectively, along with 
straight lines which represent chi-squared fits to a power law. 
The autocorrelation times as a function of $L$ are given for the Wolff and SW 
algorithms in Tables 1 and 2, for $q = 2$ and $q = 3$ respectively.
The values of $z$ from fits to all the data are given in Table 3.
We have not given any values for the RC algorithm, since the data are not 
good enough to reliably extract an exponent.
The exponents for all the other algorithms are approximately the same, and
in particular it appears that $z_{SW} = z_{Wolff}$ 
for $q=2$ and $q=3$ Potts models in two dimensions. 
%
The exponents for the LC algorithm with $q=2$ are rather high, however this
is due to
the fact that the value of $z$ tends to decrease slightly if the fits are done
at larger $L$. The SW and Wolff algorithms also give $z \approx 0.30$ if the 
fits are done only up to $L=128$.
%

Although the exponents are the same, the autocorrelation times differ 
substantially between algorithms. The Wolff algorithm has values
which are less than those of SW by a factor of about 0.45 for $q=2$ and 0.67 
for $q=3$.
It is interesting to note that these ratios are practically the same as the 
probabilities for flipping a cluster in the SW algorithm,
so that if the SW algorithm was scaled (like the single cluster algorithms) by 
the proportion of spins which are actually {\it changed} in the update, then 
the Wolff and SW algorithms would give nearly identical results.
It therefore appears that (at least in two dimensions) the main advantage of a 
single cluster algorithm such as Wolff is that the cluster is always flipped.
Notice however that flipping all the SW clusters (the AC algorithm) is
actually worse than the standard SW algorithm, and in $d > 2$ this argument
does not hold since the exponents (and not just the autocorrelation times) 
are different for the Wolff and SW algorithms
% \WIsing,
$^{ 10,}$
\ref\Brower{P. Tamayo, R.C. Brower and W. Klein,
J. Stat. Phys. {\bf 58}, 1083 (1990)}.
% ``Single Cluster Monte Carlo Dynamics for the Ising Model'', 
% Boston University preprint, (July 1989)}.

Another proposed reason for the superiority of the Wolff algorithm over 
SW is that the average cluster size is larger, and this is supported by the
fact that the LC algorithm, where the largest cluster is flipped,
does slightly better than Wolff. It is therefore quite surprising to see that
the RC algorithm, where the average cluster size is approximately the same as
for SW, does even better than Wolff or LC. This suggests that (at least in
two dimensions) it is the single cluster nature of these updates, rather than 
the size of the clusters, which is most important,
%
and that successive clusters in a single cluster update are somehow less 
correlated than the clusters in a full lattice update such as the SW 
algorithm. The overlap between sites in successive clusters gives some 
indication of how successive
clusters are correlated, and thus should affect the autocorrelations. 
The small RC clusters have a negligible average overlap, which is not the
case for the Wolff clusters in two dimensions. This overlap between 
successive clusters in the Wolff algorithm was studied by Tamayo {\it et al.}
\Brower\ for the Ising model, and found to decrease markedly as $d$ increases. 
It seems
likely that it is a combination of the effect of larger clusters, and a small
overlap between successive clusters, which results in the Wolff algorithm 
having smaller dynamical critical exponents than the SW algorithm for $d > 2$.

% We have also done preliminary investigations of two other possible update
% algorithms, with a single run of $10^5$ to $10^6$ sweeps for each lattice 
% size. The first algorithm is just SW except that the largest cluster 
% is always flipped. It gave autocorrelation times which agreed with those 
% of SW within errors.
We have also done a preliminary investigation, with a single run of $10^5$ 
to $10^6$ sweeps for each lattice size, of another algorithm, for which
all ``large'' clusters (those 
greater than a certain cutoff size) are flipped. In a study of the 
mean-field Ising model, this algorithm was found to give the same 
value of $z$ as SW, but the values of $\tau$ were halved
\ref\MFI{T.S. Ray, P. Tamayo and W. Klein, Phys. Rev. A {\bf 39}, 5949 (1989)}.
For the 2-$d$ Potts model we chose the cutoff
to be $L^{3/2}$, and found that the autocorrelation times for this 
algorithm are decreased only slightly, giving values approximately midway
between those of SW and LC, with the same value of $z$.

A number of other studies have been made of the
dynamical exponents of cluster algorithms for Potts models. 
In their original paper \SW, Swendsen and Wang measured the autocorrelations
on fairly small lattice sizes,
and thus obtained a rather high value of $0.35(1)$ for $z$.
Wolff \WIsing\ 
has compared his algorithm to SW for the 2-$d$ and 3-$d$ Ising models,
and his results agree fairly well with our data.
%
% Do they???
%
A study of the Wolff algorithm for the Ising model has also been carried
out by Tamayo et al. \Brower,
who obtain similar results. 
%
% In particular they show that $z_{W} < z_{SW}$ in three and four dimensions. 
%
Their paper also contains some remarks on 
the RC algorithm, and suggests that $z_{RC} < z_{W}$ for the 2-$d$ model.
From our results it 
appears that this is certainly possible, however the trend at larger lattice 
sizes seems to imply that the exponent is probably the same as for the 
other algorithms.
Much better data would be required to draw a firm conclusion.
Our value for $z_{SW}$ in the $q = 3$ case agrees with that of Swendsen 
and Wang \SW, and the more recent result of Li and Sokal
\ref\Li{X.-J. Li and A.D. Sokal, Phys. Rev. Lett. {\bf 63}, 827 (1989)}.
Kandel et al.
\ref\Kandel{D. Kandel, E. Domany and A. Brandt, Phys. Rev. B {\bf 40}, 
330 (1989)}
obtain a value of 0.4 by measuring the relaxation time from an ordered state 
to equilibrium, claiming to see a crossover to this smaller value of the 
exponent for $L > 32$, however we see no sign of such behavior.  
%
% It is possible that the different results are due to a cluster 
% dynamics equivalent of the difference between the equilibrium and 
% non-equilibrium exponents for local update algorithms
% \ref\Racz{Z. Racz, Phys. Rev. B {\bf 13} 263 (1976); M.E. Fisher and 
% Z. Racz, Phys. Rev. B {\bf 13} 5039 (1976)}.
%

An interesting recent development is the claim by Burkitt and Heermann
\ref\Burkitt{D.W. Heermann and A.N. Burkitt, Physica A {\bf 162}, 210 (1990)}
% ``System Size Dependence of the Autocorrelation Time for the Swendsen-Wang
% Ising Model'', Bergische University preprint, (1989)},
that the autocorrelations in the 2-$d$ Ising model using the SW algorithm 
grow logarithmically rather than as a power law, so that $z_{SW} = 0$. 
It is very difficult to distinguish between a logarithm and a small power.
Since power law or logarithmic behavior is only expected {\it asymptotically},
there is some uncertainty as to which data points should be included in the 
fit (the same is true in fitting the autocorrelation function to obtain $\tau$),
so the fits are somewhat subjective.
If we look at the slopes of lines connecting successive points in Fig. 1
(i.e. the approximations to $z$ using only two successive lattice sizes),
we see a slight trend towards smaller values as $L$ increases.
The difficult problem is to determine whether this trend disappears at large
$L$, resulting in a non-zero value of $z$, or continues, 
signifying logarithmic behavior. 
Note that the $q=3$ data fit well to a power law even down to $L = 16$, and
in this case a logarithmic increase in $\tau$ can be definitely excluded.

In Fig. 3 we plot the exponential and integrated autocorrelations in $|M|$,
along with fits to a power law (solid lines) and a logarithm 
(dashed lines). As seen by Burkitt and Heermann, $\tau_{int,|M|}$ seems to
be better fitted with a logarithm, especially for small values of $L$ (although 
these are not really relevant, since we are looking for asymptotic behavior). 
It is therefore quite possible that $z_{int,|M|}$ is zero. 
% It does not necessarily follow that $z_{exp}$ is zero, or that the 
% the values of $z_{int}$ and $z_{exp}$ for the energy are likewise zero, 
% although it certainly makes it more likely that this is the case.
%
% Does z_int = 0 imply z_exp = 0???
%
The other exponents ($z_{int,E}$ and $z_{exp}$ for $E$ and $|M|$) are much 
harder to ascertain, since the autocorrelations
do not give good fits to either a log or a power law unless smaller values of 
$L$ are excluded. 
For $\tau_{int,E}$, this can be seen in Fig. 1, where the solid 
line is a power law fit, and the dashed line a fit to a logarithm. In the case
of $\tau_{exp}$ for $|M|$, we found substantially 
greater values at large $L$ than did Burkitt and Heermann, and consequently 
were unable to reproduce their good logarithmic fit down to small $L$.
This is possibly due to large errors in their results at large $L$, although
it is difficult to tell since no errors are given.

In Table 4 we give the values of chi-squared per degree of freedom for both 
fits, for the Wolff and SW algorithms. 
The fits include only data with $L \ge 64$ for SW, and $L \ge 50$ for
Wolff, so that we are only fitting to five data points, 
and this may affect the reliability of the chi-squared values. 
The results favor logarithmic behavior, although we would not claim 
that the evidence is conclusive, especially since the differences are so small 
that corrections to scaling could be very important.
% except perhaps for $z_{int}$ of the magnetization for the SW algorithm.
Good data on much larger lattices will probably be necessary to conclusively
differentiate between $z = 0$ and $z \approx 0.25$.

\bigskip
{\bf ACKNOWLEDGEMENTS}
\medskip

The simulations were done using a number of different parallel computers.
The parallel algorithm was run on a 512 node NCUBE-1, a 192 node Symult S2010, 
and a 32 processor Meiko Computing Surface, while the sequential code was
run on a 20 processor Sequent Symmetry, a 20 processor Encore Multimax,
a 96 processor BBN GP1000 Butterfly, and a network of SUN and DEC workstations.
We would like to thank the Caltech Concurrent Supercomputer Facility
and the Advanced Computing Research Facility at Argonne National Laboratory
for the use of these machines.
This work was sponsored in part by Department of Energy grants
DE-FG03-85ER25009 and DE-AC03-81ER40050.
CFB is supported by Air Force Office of Scientific Research Grant AFOSR-89-0422
and by the Department of Energy under contract DE-AC02-86ER40253.

\listrefs\vfill

\eject
\input tab1.tex
\input tab2.tex
\input tab3.tex
\input tab4.tex
\input figcaps.tex

\end
