% -- LATEX --
\documentstyle{article}
\pagestyle{empty}

\title{\bf IV.3. Cluster Algorithms for Spin Models}
\author{\it Clive F. Baillie and Paul D. Coddington}

\begin{document}
\maketitle

\bigskip
\noindent
{\bf Monte Carlo Calculations of Spin Models}\hfill

The goal of computer simulations of spin models is to generate configurations 
of spins typical of statistical equilibrium and measure physical observables on
this ensemble of configurations. The generation of configurations is 
traditionally performed by Monte Carlo methods such as the Metropolis algorithm
\cite{Metropolis:53a}, which produce configurations with a probability given by the Boltzmann 
distribution $e^{- \beta S(\phi)}$, where $S(\phi)$ is the action, or energy,
of the system in configuration $\phi$, and $\beta$ is the inverse temperature.
One of the main problems with these methods in practice is that successive 
configurations are not statistically independent, but rather are correlated, 
with some autocorrelation time $\tau$ between effectively independent
configurations. 

A key feature about traditional (Metropolis-like) Monte Carlo algorithms is
that the updates are {\it local}, that is, one spin at a time is updated, and 
its new value depends only on the values of spins which affect its contribution
to the action, i.e. only on local (usually nearest neighbor) spins.
Thus in a single step of the algorithm, information about the state of a 
spin is transmitted only to its near neighbors. In order for the system
to reach a new effectively independent configuration, this information
must travel a distance of order the (static or spatial) correlation length
$\xi$. As the information executes a random walk around the lattice, one would
suppose that the autocorrelation time $\tau \sim \xi^2$.
However, in general $\tau \sim \xi^z$, where $z$ is called the dynamical 
critical exponent. Almost all numerical simulations of spin models
have measured $z \approx 2$ for local update algorithms.

REFER TO OTHER MONTE CARLOS IN BOOK (IV.5,12,26,27,28,30,32).

For a spin model with a phase transition, as the inverse temperature $\beta$ 
approaches the critical value, $\xi$ diverges to 
infinity, so that the computational efficiency rapidly goes to zero! This 
behavior is called critical slowing down (CSD), and until very recently it has 
plagued Monte Carlo simulations of statistical mechanical systems,
in particular spin models, at or near their phase transitions. 
Recently, however, several new ``cluster algorithms'' have been introduced which 
decrease $z$ dramatically by performing {\it non-local} spin updates, thus
reducing (or even eliminating) CSD and facilitating much more efficient
computer simulations.

\bigskip
\noindent
{\bf Cluster algorithms}\hfill

The aim of the cluster update algorithms is to find a suitable collection of 
spins which can be flipped with relatively little cost in energy.
We could obtain non-local updating very simply by using the standard 
Metropolis Monte Carlo algorithm to flip randomly selected bunches of spins,
but then the acceptance would be tiny. Therefore we need a method which picks 
sensible bunches or clusters of spins to be updated. 
The first such algorithm was proposed by Swendsen and Wang \cite{Swendsen:87a},
and was based on an equivalence between a Potts spin model \cite{Potts:52a}
\cite{Wu:82a} and percolation models \cite{Stauffer:78a} \cite{Essam:80a}, 
for which cluster properties play a fundamental role.

The Potts model is a very simple spin model of a ferromagnet, in which the 
spins can take $q$ different values.
The case $q = 2$ is just the well-known Ising model.
In the Swendsen and Wang algorithm, clusters of spins are created by 
introducing bonds between neighboring spins with probability $1-e^{-\beta}$
if the two spins are the same, and zero if they are not.
All such clusters are created, and then updated by choosing a random new 
spin value for each cluster and assigning it to all the spins in that cluster.

A variant of this algorithm, for which only a single cluster is constructed
and updated at each sweep, has been proposed by Wolff \cite{Wolff:89a}. The implementation 
of this algorithm is shown in Figs. 1-3, which shows a q = 3 Potts model
at its critical temperature, with different colors representing the three 
different spin values. From the starting configuration (Fig. 1), we choose a 
site at random, and construct a cluster around it by bonding together 
neighboring sites with the appropriate probabilities (Fig. 2).
All sites in this cluster are then given the same new spin value, producing
the new configuration shown in Fig. 3, which is obviously far less correlated
with the initial configuration than the result of a single Metropolis update (Fig. 4).
Although Wolff's method is probably the best {\it sequential} cluster algorithm,
the Swendsen and Wang algorithm seems to be better suited for 
parallelization, since it involves the entire lattice rather than just
a single cluster. We have therefore concentrated our attention on the
method of Swendsen and Wang, where {\it all} the clusters must be identified 
and labeled.

First we outline a sequential method for labeling clusters, the so-called 
``ants in the labyrinth'' algorithm.
The reason for its name is that we can visualize the algorithm as follows \cite{Dewar:87a}.
An ant is put somewhere on the lattice, and notes which of the neighboring 
sites are connected to the site it is on. At the next time-step this ant places
children on each of these connected sites which are not already occupied.
The children then proceed to reproduce likewise until the entire cluster is 
populated. 
In order to label all the clusters, we start by giving every site a negative
label, set the initial cluster label to be zero, and then loop through all the 
sites in turn. If a site's label is negative then the site has not already been 
assigned to a cluster, so we place an ant on this site, give it the current
cluster label, and let it reproduce, passing the label on to all its offspring.
When this cluster is identified we increment the cluster label and carry on,
repeating the ant-colony birth, growth and death cycle until all the clusters
have been identified.

\bigskip
\noindent
{\bf Parallel cluster algorithms}\hfill

As with the percolation models upon which the cluster algorithms are based, 
the phase transition in a spin model occurs when the clusters of bonded 
spins become large enough to span the entire lattice.
Thus near criticality (which in most cases is where we want to perform
the simulation) clusters come in all sizes, from order N (where N is the 
number of sites in the lattice) right down to a single site.
The highly irregular and non-local nature of the clusters means that
cluster update algorithms do not vectorize, and hence give poor 
performance on vector machines. On this problem a CRAY X-MP is only 
about ten times faster than a Sun4 workstation.
The irregularity of the problem also
means that SIMD machines are not well suited to this problem, whereas
for the Metropolis type algorithms they are perhaps the best machines 
available. It therefore appears that the optimum performance for this type of 
problem will come from MIMD parallel computers.

A parallel cluster algorithm involves distributing the lattice onto
an array of processors using the usual domain decomposition.
Clearly a sequential algorithm can be used to label the clusters on each 
processor, but we need a procedure for converting these labels to their
correct {\it global} values. We need to be able to tell many processors,
which may be any distance apart, that some of their clusters are actually
the same, to agree on which of the many different 
local labels for a given cluster should be assigned to be the global 
cluster label, and to pass this label to all the processors containing a 
part of that cluster. 
We have implemented two such algorithms, ``self-labeling'' and
``global equivalencing''.

\bigskip
\noindent
{\bf Self-labeling}\hfill

We shall refer to this algorithm as ``self-labeling'', since each site
figures out which cluster it is in by itself from local information.
We begin by assigning each site $i$ a unique cluster label $S_i$.
In practice this is simply chosen as the position of that site in the lattice.
At each step of the algorithm, in parallel, every site looks in turn at each 
of its neighbors in the positive directions. 
If it is bonded to a neighboring site $n$ which
has a different cluster label $S_n$, then both $S_i$ and $S_n$ are set to 
the minimum of the two. This is continued until nothing changes, by which
time all the clusters will have been labeled with the minimum initial label of 
all the sites in the cluster. Note that to check termination of the algorithm
involves each processor sending a termination flag (finished or not finished)
to every other processor after each step, which can become very costly for
a large processor array.
This is an SIMD algorithm and can therefore be run on machines like the
AMT DAP and TMC Connection Machine.
However, the SIMD nature of these computers leads to very poor
load balancing. Most processors end up waiting for the few in the largest
cluster which are the last to finish. We implemented this on the AMT DAP
and obtained only about 20\% efficiency.

We can improve this method, on a MIMD machine,
by using a faster sequential algorithm, such as 
``ants in the labyrinth'', to label the clusters in the sublattice on each
processor, and then just use self-labeling on the sites at the edges of each
processor to eventually arrive at the global cluster labels. 
The number of steps required to do the self-labeling will depend on the
largest cluster, which at the phase transition will generally span the 
entire lattice. The number of self-labeling steps will therefore be of the 
order of the maximum distance between processors, which for a square array
of $P$ processors is just $2\sqrt P$. 
Hence the amount of communication (and calculation) involved in doing the 
self-labeling, which is proportional to the number of iterations times the 
perimeter of the sublattice, goes like $L$ for an $L$x$L$ lattice, 
whereas the time taken on each processor to do the local cluster labeling goes 
like the area of the sublattice, which is $L^2/P$. Therefore as long as $L$ is 
substantially greater than the number of processors we can expect to obtain a 
reasonable speedup. 
Of course, this algorithm suffers from the same type of load imbalance as does the
SIMD version, however in this case it is much less severe since
most of the work is done with ``ants in the labyrinth'' which is well load balanced.
The speedups obtained on the Symult 2010 for a variety of lattice
sizes are shown in Fig. 5. The dashed line indicates perfect speedup (i.e.
100\% efficiency).  The lattice sizes for
which we actually need large numbers of processors are of the order of $512^2$
or greater, and we can see that running on 64 nodes (or running multiple
simulations of 64 nodes each) gives us quite acceptable efficiencies of
about 70\% for $512^2$ and 80\% for $1024^2$.

\bigskip
\noindent
{\bf Global equivalencing}\hfill

In this method we again use the fastest sequential algorithm to identify the 
clusters in the sublattice on every node. Each node then looks at
the labels of sites along the edges of the neighboring nodes in the positive
directions, and works out which ones are connected and should be matched up.
These lists of ``equivalences'' are all passed to one of the nodes,
which uses an algorithm for finding equivalence classes \cite{Knuth:68a}
\cite{Press:86a}
(which in this case
are the global cluster labels) to match up the connected clusters. 
This node then broadcasts the results back to all the other nodes. 

This part of the algorithm is purely sequential, and is thus a potentially
disastrous bottleneck for large numbers of processors. The amount of work 
involved in doing the 
global matchup goes like $P$ times the perimeter of the sublattice on each 
node, so that the efficiency should be less than for self-labeling,
although we might still expect reasonable speedups if the number of nodes
is not extremely large.
The speedups obtained on the Symult 2010 for a variety of lattice
sizes are shown in Fig. 6.
Global equivalencing gives about the same speedups as self-labeling for
small numbers of processors, but as expected self-labeling
does much better as the number of nodes increases .


\bigskip
\noindent
{\bf Other algorithms}\hfill

Currently the only other parallel cluster algorithm proposed for spin models 
is due to Burkitt and Heermann \cite{Burkitt:89a}, but it is much more complicated, and
less efficient, than the self-labeling algorithm. However
the problem of labeling clusters of spins is very similar to an important
problem in image analysis, that of identifying and labeling the connected
components in a binary or multi-colored image composed of an array of pixels, 
and there have been a number of parallel algorithms implemented for this 
problem \cite{Embrechts:89a} \cite{Lim:87a} \cite{Cypher:89a}.
The most promising of these parallel algorithms for spin models
has a hierarchical divide-and-conquer approach \cite{Embrechts:89a}.
The processor array is divided up into smaller subarrays of, for 
example, $2$x$2$ processors. In each subarray, the processors look at the edges
of their neighbors for clusters which are connected across processor boundaries.
As in global equivalencing, these equivalences are all passed to one of the nodes of the sub-array,
which places them in equivalence classes.
The results of these partial matchings are similarly combined
on each $4$x$4$ subarray, and this process is continued until finally all the
partial results are merged together on a single processor to give the global 
cluster values.

Finally, we should mention the trivial
parallelization technique of running independent Monte Carlo simulations 
on different processors.
This method works well until the lattice size gets too big to fit into the 
memory of each node.
In the case of the Potts model, for example, only lattices of size 
less than about $300^2$ or $50^3$ will fit into 1 Mbyte, though most other spin 
models are more complicated and more memory intensive.
The smaller lattices which are seen to give poor speedups in
Figs. 5 and 6 can be run with 100\% efficiency in this way. Note, of course,
that this requires an MIMD computer.
In fact, we have used this method to calculate the
dynamical critical exponents of various cluster algorithms for Potts models 
\cite{Baillie:89l} on lattices up to $128^2$.
In order to distinguish whether the dynamical critical exponent is small
or is in fact zero, we need to go to larger lattices.
Thus, we are now simulating lattices of size $256^2$ and $512^2$ using
the self-labeling parallel cluster algorithm, on a number of MIMD computers.

\newpage
\bibliography{/home/swan/kea/bib/doc.bib,/home/swan/kea/bib/outside.bib,/home/swan/kea/bib/spcp.bib}
\bibliographystyle{/usr/lib/publisher/inputs/ccsf_sort}

\end{document}
