Swendsen–Wang algorithm

The problem of the critical slowing-down affecting local processes is of fundamental importance in the study of second-order phase transitions (like ferromagnetic transition in the Ising model), as increasing the size of the system in order to reduce finite-size effects has the disadvantage of requiring a far larger number of moves to reach thermal equilibrium. Indeed the correlation time ${\displaystyle \tau }$ usually increases as ${\displaystyle L^{z}}$ with ${\displaystyle z\simeq 2}$ or greater; since, to be accurate, the simulation time must be ${\displaystyle t\gg \tau }$, this is a major limitation in the size of the systems that can be studied through local algorithms. SW algorithm was the first to produce unusually small values for the dynamical critical exponents: ${\displaystyle z=0.35}$ for the 2D Ising model (${\displaystyle z=2.125}$ for standard simulations); ${\displaystyle z=0.75}$ for the 3D Ising model, as opposed to ${\displaystyle z=2.0}$ for standard simulations.

The algorithm is non-local in the sense that a single sweep updates a collection of spin variables based on the Fortuin–Kasteleyn representation. The update is done on a "cluster" of spin variables connected by open bond variables that are generated through a percolation process, based on the interaction states of the spins.

Consider a typical ferromagnetic Ising model with only nearest-neighbor interaction.

• Starting from a given configuration of spins, we associate to each pair of nearest neighbours on sites ${\displaystyle n,m}$ a random variable ${\displaystyle b_{n,m}\in \lbrace 0,1\rbrace }$ which is interpreted in the following way: if ${\displaystyle b_{n,m}=0}$ then there is no link between the sites ${\displaystyle n}$ and ${\displaystyle m}$ (the bond is closed); if ${\displaystyle b_{n,m}=1}$ then there is a link connecting the spins ${\displaystyle \sigma _{n}{\text{ and }}\sigma _{m}}$(the bond is open). These values are assigned according to the following (conditional) probability distribution:

 ${\displaystyle P\left[b_{n,m}=0|\sigma _{n}\neq \sigma _{m}\right]=1}$;
 ${\displaystyle P\left[b_{n,m}=1|\sigma _{n}\neq \sigma _{m}\right]=0}$;
 ${\displaystyle P\left[b_{n,m}=0|\sigma _{n}=\sigma _{m}\right]=e^{-2\beta J_{nm}}}$;
 ${\displaystyle P\left[b_{n,m}=1|\sigma _{n}=\sigma _{m}\right]=1-e^{-2\beta J_{nm}}}$;

where ${\displaystyle J_{nm}>0}$ is the ferromagnetic coupling strength.

This probability distribution has been derived in the following way: the Hamiltonian of the Ising model is

${\displaystyle H[\sigma ]=\sum \limits _{<i,j>}-J_{i,j}\sigma _{i}\sigma _{j}}$,

and the partition function is

${\displaystyle Z=\sum \limits _{\lbrace \sigma \rbrace }e^{-\beta H[\sigma ]}}$.

Consider the interaction between a pair of selected sites ${\displaystyle n}$ and ${\displaystyle m}$ and eliminate it from the total Hamiltonian, defining ${\displaystyle H_{nm}[\sigma ]=\sum \limits _{<i,j>\neq <n,m>}-J_{i,j}\sigma _{i}\sigma _{j}.}$

Define also the restricted sums:

${\displaystyle Z_{n,m}^{same}=\sum \limits _{\lbrace \sigma \rbrace }e^{-\beta H_{nm}[\sigma ]}\delta _{\sigma _{n},\sigma _{m}}}$;

${\displaystyle Z_{n,m}^{diff}=\sum \limits _{\lbrace \sigma \rbrace }e^{-\beta H_{nm}[\sigma ]}\left(1-\delta _{\sigma _{n},\sigma _{m}}\right).}$

${\displaystyle Z=e^{\beta J_{nm}}Z_{n,m}^{same}+e^{-\beta J_{nm}}Z_{n,m}^{diff}.}$

Introduce the quantity

${\displaystyle Z_{nm}^{ind}=Z_{n,m}^{same}+Z_{n,m}^{diff}}$;

the partition function can be rewritten as

${\displaystyle Z=\left(e^{\beta J_{nm}}-e^{-\beta J_{nm}}\right)Z_{n,m}^{same}+e^{-\beta J_{nm}}Z_{n,m}^{ind}.}$

Since the first term contains a restriction on the spin values whereas there is no restriction in the second term, the weighting factors (properly normalized) can be interpreted as probabilities of forming/not forming a link between the sites: ${\displaystyle P_{<n,m>\;link}=1-e^{-2\beta J_{nm}}.}$ The process can be easily adapted to antiferromagnetic spin systems, as it is sufficient to eliminate ${\displaystyle Z_{n,m}^{same}}$ in favor of ${\displaystyle Z_{n,m}^{diff}}$ (as suggested by the change of sign in the interaction constant).

• After assigning the bond variables, we identify the same-spin clusters formed by connected sites and make an inversion of all the variables in the cluster with probability 1/2. At the following time step we have a new starting Ising configuration, which will produce a new clustering and a new collective spin-flip.

It can be shown that this algorithm leads to equilibrium configurations. To show this, we interpret the algorithm as a Markov chain, and show that the chain is both ergodic (when used together with other algorithms) and satisfies detailed balance, such that the equilibrium Boltzmann distribution is equal to the stationary distribution of the chain.

Ergodicity means that it is possible to transit from any initial state to any final state with a finite number of updates. It has been shown that the SW algorithm is not ergodic in general (in the thermodynamic limit).[2] Thus in practice, the SW algorithm is usually used in conjunction with single spin-flip algorithms such as the Metropolis–Hastings algorithm to achieve ergodicity.

The SW algorithm does however satisfy detailed-balance. To show this, we note that every transition between two Ising spin states must pass through some bond configuration in the percolation representation. Let's fix a particular bond configuration: what matters in comparing the probabilities related to it is the number of factors ${\displaystyle q=e^{-2\beta J}}$ for each missing bond between neighboring spins with the same value; the probability of going to a certain Ising configuration compatible with a given bond configuration is uniform (say ${\displaystyle p}$). So the ratio of the transition probabilities of going from one state to another is

${\displaystyle {\frac {P_{\lbrace \sigma \rbrace \rightarrow \lbrace \sigma '\rbrace }}{P_{\lbrace \sigma '\rbrace \rightarrow \lbrace \sigma \rbrace }}}={\frac {Pr\left(\lbrace \sigma '\rbrace |B.C.\right)Pr\left(B.C.|\lbrace \sigma \rbrace \right)}{Pr\left(\lbrace \sigma \rbrace |B.C.\right)Pr\left(B.C.|\lbrace \sigma '\rbrace \right)}}={\frac {p\cdot \exp \left[-2\beta \sum \limits _{<l,m>}\delta _{\sigma _{l},\sigma _{m}}J_{lm}\right]}{p\cdot \exp \left[-2\beta \sum \limits _{<l,m>}\delta _{\sigma '_{l},\sigma '_{m}}J_{lm}\right]}}=e^{-\beta \Delta E}}$

since ${\displaystyle \Delta E=-\sum \limits _{<l,m>}J_{lm}\left(\sigma '_{l}\sigma '_{m}-\sigma _{l}\sigma _{m}\right)=-\sum \limits _{<l,m>}J_{lm}\left[\delta _{\sigma '_{l},\sigma '_{m}}-\left(1-\delta _{\sigma '_{l},\sigma '_{m}}\right)-\delta _{\sigma _{l},\sigma _{m}}+\left(1-\delta _{\sigma _{l},\sigma _{m}}\right)\right]=-2\sum \limits _{<l,m>}J_{lm}\left(\delta _{\sigma '_{l},\sigma '_{m}}-\delta _{\sigma _{l},\sigma _{m}}\right)}$.

This is valid for every bond configuration the system can pass through during its evolution, so detailed balance is satisfied for the total transition probability. This proves that the algorithm is correct.

Although not analytically clear from the original paper, the reason why all the values of z obtained with the SW algorithm are much lower than the exact lower bound for single-spin-flip algorithms (${\displaystyle z\geq \gamma /\nu }$) is that the correlation length divergence is strictly related to the formation of percolation clusters, which are flipped together. In this way the relaxation time is significantly reduced. Another way to view this is through the correspondence between the spin statistics and cluster statistics in the Edwards-Sokal representation.[3]

The algorithm is not efficient in simulating frustrated systems, because the correlation length of the clusters is larger than the correlation length of the spin model in the presence of frustrated interactions.[4] Currently, there are two main approaches to addressing this problem, such that the efficiency of cluster algorithms is extended to frustrated systems.

The first approach is to extend the bond-formation rules to more non-local cells, and the second approach is to generate clusters based on more relevant order parameters. In the first case, we have the KBD algorithm for the fully-frustrated Ising model, where the decision of opening bonds are made on each plaquette, arranged in a checkerboard pattern on the square lattice.[5] In the second case, we have replica cluster move for low-dimensional spin glasses, where the clusters are generated based on spin overlaps, which is believed to be the relevant order parameter.

• Random cluster model
• Monte Carlo method
• Wolff algorithm
• http://www.hpjava.org/theses/shko/thesis_paper/node69.html
• http://www-fcs.acs.i.kyoto-u.ac.jp/~harada/monte-en.html

• Swendsen, Robert H.; Wang, Jian-Sheng (1987-01-12). "Nonuniversal critical dynamics in Monte Carlo simulations". Physical Review Letters. American Physical Society (APS). 58 (2): 86–88. Bibcode:1987PhRvL..58...86S. doi:10.1103/physrevlett.58.86. ISSN 0031-9007. PMID 10034599.
• Kasteleyn P. W. and Fortuin (1969) J. Phys. Soc. Jpn. Suppl. 26s:11
• Fortuin, C.M.; Kasteleyn, P.W. (1972). "On the random-cluster model". Physica. Elsevier BV. 57 (4): 536–564. Bibcode:1972Phy....57..536F. doi:10.1016/0031-8914(72)90045-6. ISSN 0031-8914.
• Wang, Jian-Sheng; Swendsen, Robert H. (1990). "Cluster Monte Carlo algorithms". Physica A: Statistical Mechanics and Its Applications. Elsevier BV. 167 (3): 565–579. Bibcode:1990PhyA..167..565W. doi:10.1016/0378-4371(90)90275-w. ISSN 0378-4371.
• Barbu, A. (2005). "Generalizing Swendsen-Wang to sampling arbitrary posterior probabilities". IEEE Transactions on Pattern Analysis and Machine Intelligence. Institute of Electrical and Electronics Engineers (IEEE). 27 (8): 1239–1253. doi:10.1109/tpami.2005.161. ISSN 0162-8828. PMID 16119263. S2CID 410716.