12.6 Cluster Algorithms for Spin Models

12.6.1 Monte Carlo Calculations of Spin Models

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 [?], which produce configurations with a probability given by the Boltzmann distribution e^{−βS(ϕ)}, where S(ϕ) is the action, or energy, of the system in configuration ϕ, and β 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, τ, between effectively independent configurations.

A key feature about traditional (Metropolis-like) Monte Carlo algorithms is that the updates are 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 nearest 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 ξ. As the information executes a random walk around the lattice, one would suppose that the autocorrelation time, τ ∼ ξ^{2}. However, in general τ ∼ ξ^{z}, where z is called the dynamical critical exponent. Almost all numerical simulations of spin models have measured z ≈ 2 for local update algorithms. (See also Sections ??, ??, ??, 12.6, and ??).

For a spin model with a phase transition, as the inverse temperature β approaches the critical value, ξ 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 non-local spin updates, thus reducing (or even eliminating) CSD and facilitating much more efficient computer simulations.

12.6.2 Cluster Algorithms

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 [?], and was based on an equivalence between a Potts spin model [?], [?] and percolation models [?], [?], 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^{−β} 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 [?]. The implementation of this algorithm is shown in Figure 12.1 to Figure 12.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 (Figure 12.1), we choose a site at random, and construct a cluster around it by bonding together neighboring sites with the appropriate probabilities (Figure 12.2). All sites in this cluster are then given the same new spin value, producing the new configuration shown in Figure 12.3, which is obviously far less correlated with the initial configuration than the result of a single Metropolis update (Figure 12.4). Although Wolff’s method is probably the best 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 parallelizing the method of Swendsen and Wang, where 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 [?]. 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 Figure 12.1: Initial Configuration of 3-state Potts Spins on which Wolff Algorithm is to be Applied

Figure 12.2: Bonds of Cluster Constructed by Wolff Algorithm Indicated in Yellow Figure 12.3: Result of Wolff Algorithm Applied to Spin Configuration in Figure 12.1—all spins in cluster flipped to same new value (in this case from blue to red)

Figure 12.4: Result of Metropolis Algorithm Applied to Spin Configuration in Figure 12.1—only a few single spins flipped the ant-colony birth, growth and death cycle until all the clusters have been identified.

12.6.3 Parallel Cluster Algorithms

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 well, and hence give poor performance on vector machines. On this problem, a CRAY X-MP is only about ten times faster than a Sun 4 workstation. The irregularity of the clusters 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 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” [?], [?].

12.6.4 Self-labeling

We shall refer to this algorithm as “self-labeling”, since each site figures out which cluster it is in by itself from local information. This method has also been referred to as “local label propagation” [?], [?]. 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 sub-lattice 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 sub-lattice, goes like L for an L × L lattice, whereas the time taken on each processor to do the local cluster labeling goes like the area of the sub-lattice, 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 Figure 12.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}.

12.6.5 Global Equivalencing

In this method, we again use the fastest sequential algorithm to identify the clusters in the sub-lattice on every processor. Each processor then looks at the labels of sites along the edges of the neighboring processors 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 processors, which uses an algorithm for finding equivalence classes [?], [?]

Figure 12.5: Speedups for Self-Labeling Algorithm

(which, in this case, are the global cluster labels) to match up the connected clusters. This processor then broadcasts the results back to all the other processors.

This part of the algorithm is purely sequential, and is thus a potentially disastrous bottleneck for large numbers of processors. It also requires this processor to have a large amount of memory in which to store all the labels from every other processor. The amount of work involved in doing the global matchup goes like P times the perimeter of the sub-lattice on each processor, i.e., like L\sqrt{P} so that the efficiency should be less than for selflabeling, although we might still expect reasonable speedups if the number of processors is not extremely large. The speedups obtained on the Symult 2010 for a variety of lattice sizes are shown in Figure 12.6. The figure for 512^{2} on 128 processors is missing due to memory constraints. 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 processors increases.

12.6.6 Other Algorithms

The problem of labeling clusters of spins is an example of a standard graph problem known as connected component labeling [?]. Another important instance occurs in image analysis, in identifying and labeling the connected components of a binary or multi-colored image composed of an array of pixels [?]. There have been a number of parallel algorithms implemented for this problem [?], [?], [?], [?]. The most promising of these parallel algorithms for spin models has a hierarchical divide-and-conquer approach

Figure 12.6: Speedups for Global Equivalencing Algorithm

[?], [?]. The processor array is divided up into smaller sub-arrays of, for example, 2 × 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 × 4 sub-array, 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 Figure 12.5 and Figure 12.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 [?], [?], [?] (see Section ??).

12.6.7 Credits

This research was performed by C. F. Baillie, P. D. Coddington, J. Apostolakis and E. Marinari.