Swendsen–Wang algorithm

The Swendsen–Wang algorithm is the first non-local or cluster algorithm for Monte Carlo simulation for large systems near criticality. It has been introduced by Robert Swendsen and Jian-Sheng Wang in 1987 at Carnegie Mellon. The original algorithm was designed for the Ising and Potts models, and it was later generalized to other systems as well, such as the XY model by Wolff algorithm and particles of fluids. The key ingredient was the random cluster model, a representation of the Ising or Potts model through percolation models of connecting bonds, due to Fortuin and Kasteleyn. It has been generalized by Barbu and Zhu to arbitrary sampling probabilities by viewing it as a Metropolis–Hastings algorithm and computing the acceptance probability of the proposed Monte Carlo move.

The Swendsen–Wang algorithm is the first non-local or cluster algorithm for Monte Carlo simulation for large systems near criticality. It has been introduced by Robert Swendsen and Jian-Sheng Wang in 1987 at Carnegie Mellon. The original algorithm was designed for the Ising and Potts models, and it was later generalized to other systems as well, such as the XY model by Wolff algorithm and particles of fluids. The key ingredient was the random cluster model, a representation of the Ising or Potts model through percolation models of connecting bonds, due to Fortuin and Kasteleyn. It has been generalized by Barbu and Zhu to arbitrary sampling probabilities by viewing it as a Metropolis–Hastings algorithm and computing the acceptance probability of the proposed Monte Carlo move. 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 \tau usually increases as L^z with z\simeq 2 or greater; since, to be accurate, the simulation time must be 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: z=0.35 for the 2D Ising model ( z=2.125 for standard simulations); z=0.75 for the 3D Ising model, as opposed to z=2.0 for standard simulations. Starting from a given configuration of spins, we associate to each pair of nearest neighbours on sites n,m a random variable b_{n,m}\in \lbrace 0,1\rbrace which is interpreted in the following way: if b_{n,m}=0 then there is no link between the sites n and m (the bond is closed); if b_{n,m}=1 then there is a link connecting the spins \sigma_n \text{ and } \sigma_m (the bond is open). These values are assigned according to the following (conditional) probability distribution: P\left[b_{n,m}=0|\sigma_n\neq\sigma_m\right]=1 ; P\left[b_{n,m}=1|\sigma_n\neq\sigma_m\right]=0 ; P\left[b_{n,m}=0|\sigma_n=\sigma_m\right]=e^{-2\beta J_{nm}} ; P\left[b_{n,m}=1|\sigma_n=\sigma_m\right]=1-e^{-2\beta J_{nm}} ; where 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 H[\sigma]=\sum\limits_{\langle i,j\rangle }-J_{i,j}\sigma_i\sigma_j , and the partition function is Z=\sum\limits_{\lbrace\sigma\rbrace}e^{-\beta H[\sigma]} . Consider the interaction between a pair of selected sites n and m and eliminate it from the total Hamiltonian, defining H_{nm}[\sigma]=\sum\limits_{\langle i,j \rangle \neq \langle n,m \rangle}-J_{i,j}\sigma_i\sigma_j. Define also the restricted sums: Z_{n,m}^{same}=\sum\limits_{\lbrace\sigma\rbrace}e^{-\beta H_{nm}[\sigma]}\delta_{\sigma_n,\sigma_m} ; Z_{n,m}^{diff}=\sum\limits_{\lbrace\sigma\rbrace}e^{-\beta H_{nm}[\sigma]}\left(1-\delta_{\sigma_n,\sigma_m}\right). Z=e^{\beta J_{nm}}Z_{n,m}^{same}+e^{-\beta J_{nm}}Z_{n,m}^{diff}. Introduce the quantity Z_{nm}^{ind}=Z_{n,m}^{same}+Z_{n,m}^{diff} ; the partition function can be rewritten as 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: P_{\langle n,m \rangle\;link}=1-e^{-2\beta J_{nm}}. The process can be easily adapted to antiferromagnetic spin systems, as it is sufficient to eliminate Z_{n,m}^{same} in favor of 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. 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 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 p ). So the ratio of the transition probabilities of going from one state to another is \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_{\langle l,m \rangle}\delta_{\sigma_l,\sigma_m}J_{lm}\right]}{p\cdot \exp\left[-2\beta\sum\limits_{\langle l,m \rangle}\delta_{\sigma'_l,\sigma'_m}J_{lm}\right]} =e^{-\beta\Delta E} since \Delta E=-\sum\limits_{\langle l,m \rangle}J_{lm}\left(\sigma'_l \sigma'_m - \sigma_l \sigma_m\right)=-\sum\limits_{\langle l,m \rangle}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_{\langle l,m \rangle}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.

References

Kasteleyn P. W. and Fortuin (1969) J. Phys. Soc. Jpn. Suppl. 26s:11

Text derived from https://en.wikipedia.org/w/index.php?oldid=1375977191 under CC BY-SA 4.0; modifications have been made.