Swendsen–Wang algorithm: Difference between revisions

Revision as of 14:49, 18 December 2016

Latest revision as of 20:10, 22 August 2025

(41 intermediate revisions by 19 users not shown)

Line 1:

The '''Swendsen–Wang algorithm''' is the first non-local or cluster [[algorithm]] for [[Monte Carlo simulation]] for large systems near [[Critical point (thermodynamics)|criticality]]. It has been introduced by [[Robert Swendsen]] and [[Jian-Sheng Wang]] in 1987 at [[Carnegie Mellon University|Carnegie Mellon]].

The original algorithm was designed for the Ising and Potts models, and later it was generalized to other systems as well, such as the XY model by [[Wolff algorithm]] and particles of fluids. A key ingredient of the method is based on the 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 (2005) to sampling arbitrary probabilities by viewing it as a [[Metropolis–Hastings algorithm]] and computing the acceptance probability of the proposed Monte Carlo move.

The original algorithm was designed for the [[Ising model|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|Potts]] model through percolation models of connecting bonds, due to Fortuin and Kasteleyn. It has been generalized by Barbu and Zhu<ref>{{Cite journal|last1=Barbu|first1=Adrian|last2=Zhu|first2=Song-Chun|date=August 2005|title=Generalizing Swendsen-Wang to sampling arbitrary posterior probabilities|journal=IEEE Transactions on Pattern Analysis and Machine Intelligence|volume=27|issue=8|pages=1239–1253|doi=10.1109/TPAMI.2005.161|issn=0162-8828|pmid=16119263|bibcode=2005ITPAM..27.1239B |s2cid=410716}}</ref> to arbitrary sampling probabilities by viewing it as a [[Metropolis–Hastings algorithm]] and computing the acceptance probability of the proposed Monte Carlo move.

== Motivation ==

The problem of the critical slowing-down affecting local processes is of fundamental importance in the study of second-order [[phase transitionstransition]]s (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 <math>\tau</math> usually increases as <math>L^z</math> with <math>z\simeq 2</math> or greater; since, to be accurate, the simulation time must be <math>t\gg\tau</math>, this is a major limitation in the size of the systems that can be studied through [[Local algorithm|local algorithms]]. SW algorithm was the first to produce unusually small values for the dynamical critical exponents: <math>z=0.35</math> for the 2D Ising model (<math>z=2.125</math> for standard simulations); <math>z=0.75</math> for the 3D Ising model, as opposed to <math>z=2.0</math> for standard simulations.

Indeed the correlation time <math>\tau</math> usually increases as <math>L^z</math> with <math>z\simeq 2</math> or greater; since, to be accurate, the simulation time must be <math>t\gg\tau</math>, 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: <math>z=0.35</math> for the 2D Ising model (<math>z=2.125</math> for standard simulations); <math>z=0.75</math> for the 3D Ising model, as opposed to <math>z=2.0</math> for standard simulations.

== Description ==

{{Main|Random cluster model}}

The algorithm is non-local in the sense that in a single sweep of moves a collective update of the spin variables of the system is done. The key idea is to take an additional number of 'bond' variables, as suggested by Fortuin and Kasteleyn, who mapped the Potts model onto a percolation model.

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

Let's considerConsider a typical ferromagnetic Ising model with only nearest-neighbourneighbor interaction.

* Starting from a given configuration of spins, we associate to each couplepair of nearest neighbours on sites <math>n,m</math> a random variable <math>b_{n,m}\in \lbrace 0,1\rbrace</math> which is interpreted in the following way: if <math>b_{n,m}=0</math> then there is no link between the sites <math>n</math> and <math>m</math> (the bond is '''closed'''); if <math>b_{n,m}=1</math>, then there is a link connecting the spins <math>\sigma_n</math> \text{ and <math>} \sigma_m</math>(the arebond is connected'''open'''). These values are assigned according to the following (conditional) probability distribution:

<math>P\left[b_{n,m}=0|\sigma_n\neq\sigma_m\right]=1</math>;

Line 21 ⟶ 18:

<math>P\left[b_{n,m}=0|\sigma_n=\sigma_m\right]=e^{-2\beta J_{nm}}</math>;

<math>P\left[b_{n,m}=1|\sigma_n=\sigma_m\right]=1-e^{-2\beta J_{nm}}</math>;

where <math>J_{nm}>0</math> is the ferromagnetic interactioncoupling intensitystrength.

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

<math>H[\sigma]=\sum\limits_{<i,j>}-J_{i,j}\sigma_i\sigma_j</math>,

and the [[Partition function (statistical mechanics)|partition function]] is

<math>Z=\sum\limits_{\lbrace\sigma\rbrace}e^{-\beta H[\sigma]}</math>.

Consider the interaction between a couplepair of selected sites $<math>n$</math> and $<math>m$</math> and eliminate it from the total Hamiltonian, defining

<math>H_{nm}[\sigma]=\sum\limits_{<i,j>\neq<n,m>}-J_{i,j}\sigma_i\sigma_j.</math>

Define also the restricted sums:

<math>Z_{n,m}^{same}=\sum\limits_{\lbrace\sigma\rbrace}e^{-\beta H_{nm}[\sigma]}\delta_{\sigma_n,\sigma_m}</math>;

<math>Z_{n,m}^{diff}=\sum\limits_{\lbrace\sigma\rbrace}e^{-\beta H_{nm}[\sigma]}\left(1-\delta_{\sigma_n,\sigma_m}\right).</math>

<math>Z=e^{\beta J_{nm}}Z_{n,m}^{same}+e^{-\beta J_{nm}}Z_{n,m}^{diff}.</math>

Introduce the quantity

<math>Z_{nm}^{ind}=Z_{n,m}^{same}+Z_{n,m}^{diff}</math>;

the partition function can be rewritten as

<math>Z=\left(e^{\beta J_{nm}}-e^{-\beta J_{nm}}\right)Z_{n,m}^{same}+e^{-\beta J_{nm}}Z_{n,m}^{ind}.</math>

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: <math>P_{<n,m>\;link}=1-e^{-2\beta J_{nm}}.</math>

The process can be easily adapted to antiferromagnetic spin systems, as it is sufficient to eliminate <math>Z_{n,m}^{same}</math> in favor of <math>Z_{n,m}^{diff}</math>, in(as accordancesuggested toby 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 $\frac{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.

* 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 $\frac{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.

== Correctness ==

It can be shown that this algorithm leads to equilibrium configurations. TheTo firstshow waythis, towe proveinterpret itthe isalgorithm usingas the theory ofa [[Markov chainschain]], eitherand notingshow that the equilibriumchain (describedis byboth Boltzmann-Gibbs[[Ergodicity|ergodic]] distribution)(when mapsused intotogether itself,with orother showingalgorithms) thatand insatisfies a[[detailed singlebalance]], sweepsuch ofthat the latticeequilibrium there[[Boltzmann distribution]] is aequal non-zeroto probabilitythe of[[stationary going from any statedistribution]] of the Markov chain to any other; thus the corresponding irreducible ergodic Markov chain has an asymptotic probability distribution satisfying detailed balance.

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).<ref>{{Cite journal|last1=Gore|first1=Vivek K.|last2=Jerrum|first2=Mark R.|date=1999-10-01|title=The Swendsen–Wang Process Does Not Always Mix Rapidly|url=https://doi.org/10.1023/A:1004610900745|journal=Journal of Statistical Physics|language=en|volume=97|issue=1|pages=67–86|doi=10.1023/A:1004610900745|bibcode=1999JSP....97...67G|s2cid=189821827|issn=1572-9613|url-access=subscription}}</ref> 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.

Alternatively, we can show explicitily that detailed balance is satisfied. Every transition between two Ising configurations 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 <math>q=e^{-2\beta J}</math> 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 <math>p</math>). So the ratio of the transition probabilities of going from one state to another is

Alternatively,The weSW canalgorithm showdoes explicitilyhowever thatsatisfy detailed -balance. isTo satisfied.show this, we note that Everyevery transition between two Ising configurationsspin 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 <math>q=e^{-2\beta J}</math> 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 <math>p</math>). So the ratio of the transition probabilities of going from one state to another is

<math>\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]}

={-\beta\Delta E}</math>

since <math>\Delta E=-\sum\limits_frac{<l,m>}J_P_{lm}\left(lbrace\sigma\rbrace\rightarrow\lbrace\sigma'_l \rbrace}}{P_{\lbrace\sigma'_m - \sigma_l rbrace\sigma_mrightarrow\right)=-lbrace\sumsigma\limits_{<l,m>rbrace}J_{lm}=\frac{Pr\left[(\delta_{lbrace\sigma'_l,\rbrace|B.C.\right)Pr\left(B.C.|\lbrace\sigma'_m\rbrace\right)}-{Pr\left(1-\delta_{lbrace\sigma'_l,\rbrace|B.C.\right)Pr\left(B.C.|\lbrace\sigma'_m}\rbrace\right)-}=\delta_frac{p\sigma_l,cdot \sigma_m}+exp\left(1[-2\beta\sum\limits_{<l,m>}\delta_{\sigma'_lsigma_l,\sigma'_msigma_m}J_{lm}\right)]}{p\right]=cdot \exp\left[-2\beta\sum\limits_{<l,m>}J_{lm}\left(\delta_{\sigma'_l,\sigma'_m}-\delta_J_{\sigma_l,\sigma_mlm}\right)</math>.]}

=e^{-\beta\Delta E}</math>

since <math>\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)</math>.

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 works.

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 worksis correct.

== Efficiency ==

Although not analiticallyanalytically 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 algortithmsalgorithms (<math>z\geq\gamma/\nu</math>) 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 [[Random cluster model|Edwards-Sokal representation]].<ref>{{Cite journal|last1=Edwards|first1=Robert G.|last2=Sokal|first2=Alan D.|date=1988-09-15|title=Generalization of the Fortuin-Kasteleyn-Swendsen-Wang representation and Monte Carlo algorithm|url=https://link.aps.org/doi/10.1103/PhysRevD.38.2009|journal=Physical Review D|volume=38|issue=6|pages=2009–2012|doi=10.1103/PhysRevD.38.2009|pmid=9959355|bibcode=1988PhRvD..38.2009E|url-access=subscription}}</ref> Some mathematically rigorous results on the mixing time of this process have been obtained by Guo and Jerrum [https://projecteuclid.org/journals/annals-of-applied-probability/volume-28/issue-2/Random-cluster-dynamics-for-the-Ising-model-is-rapidly-mixing/10.1214/17-AAP1335.full].

== Generalizations ==

The algorithm is not particularly good for simulating frustrated systems.

The algorithm is not efficient in simulating [[Geometrical frustration|frustrated systems]], because the [[Percolation critical exponents|correlation length of the clusters]] is larger than the [[Correlation function (statistical mechanics)|correlation length of the spin model]] in the presence of frustrated interactions.<ref>{{Cite journal|last1=Cataudella|first1=V.|last2=Franzese|first2=G.|last3=Nicodemi|first3=M.|last4=Scala|first4=A.|last5=Coniglio|first5=A.|date=1994-03-07|title=Critical clusters and efficient dynamics for frustrated spin models|url=https://link.aps.org/doi/10.1103/PhysRevLett.72.1541|journal=Physical Review Letters|volume=72|issue=10|pages=1541–1544|doi=10.1103/PhysRevLett.72.1541|pmid=10055635|bibcode=1994PhRvL..72.1541C|hdl=2445/13250|hdl-access=free}}</ref> 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 [[Domino tiling|fully-frustrated Ising model]], where the decision of opening bonds are made on each plaquette, arranged in a checkerboard pattern on the square lattice.<ref>{{Cite journal|last1=Kandel|first1=Daniel|last2=Ben-Av|first2=Radel|last3=Domany|first3=Eytan|date=1990-08-20|title=Cluster dynamics for fully frustrated systems|url=https://link.aps.org/doi/10.1103/PhysRevLett.65.941|journal=Physical Review Letters|volume=65|issue=8|pages=941–944|doi=10.1103/PhysRevLett.65.941|pmid=10043065|bibcode=1990PhRvL..65..941K|url-access=subscription}}</ref> In the second case, we have [[replica cluster move]] for low-dimensional [[Spin glass|spin glasses]], where the clusters are generated based on spin overlaps, which is believed to be the relevant order parameter.

==References==

*Swendsen, R. H., and Wang, J.-S. (1987), ''Nonuniversal critical dynamics in Monte Carlo simulations'', Phys. Rev. Lett., 58(2):86&ndash;88.

*Kasteleyn P. W. and Fortuin (1969) J. Phys. Soc. Jpn. Suppl. 26s:11; Fortuin C. M. and Kasteleyn P.W. (1972), Physica (Utrecht) 57:536.

*Wang J.-S. and Swendsen, R. H. (1990),''Cluster Monte Carlo algorithms,'' Physica A 167:565.

*Barbu, A., Zhu, S. C. (2005), ''Generalizing Swendsen-Wang to sampling arbitrary posterior probabilities'', IEEE Trans Patt. Anal. Mach. Intell., 27(8):1239-1253.

==External linksSee also ==

*[http://www.hpjava.org/theses/shko/thesis_paper/node69.html]

* [[Random cluster model]]

*[http://www-fcs.acs.i.kyoto-u.ac.jp/~harada/monte-en.html]

* [[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]

==References==

{{Reflist}}

*{{cite journal | last1=Swendsen | first1=Robert H. | last2=Wang | first2=Jian-Sheng | title=Nonuniversal critical dynamics in Monte Carlo simulations | journal=Physical Review Letters | publisher=American Physical Society (APS) | volume=58 | issue=2 | date=1987-01-12 | issn=0031-9007 | doi=10.1103/physrevlett.58.86 | pages=86–88| pmid=10034599 | bibcode=1987PhRvL..58...86S }}

*Kasteleyn P. W. and Fortuin (1969) J. Phys. Soc. Jpn. Suppl. 26s:11; Fortuin C. M. and Kasteleyn P.W. (1972), Physica (Utrecht) 57:536.

*{{cite journal | last1=Fortuin | first1=C.M. | last2=Kasteleyn | first2=P.W. | title=On the random-cluster model | journal=Physica | publisher=Elsevier BV | volume=57 | issue=4 | year=1972 | issn=0031-8914 | doi=10.1016/0031-8914(72)90045-6 | pages=536–564| bibcode=1972Phy....57..536F }}

*{{cite journal | last1=Wang | first1=Jian-Sheng | last2=Swendsen | first2=Robert H. | title=Cluster Monte Carlo algorithms | journal=Physica A: Statistical Mechanics and Its Applications | publisher=Elsevier BV | volume=167 | issue=3 | year=1990 | issn=0378-4371 | doi=10.1016/0378-4371(90)90275-w | pages=565–579| bibcode=1990PhyA..167..565W }}

*{{cite journal | last=Barbu | first=A. | title=Generalizing Swendsen-Wang to sampling arbitrary posterior probabilities | journal=IEEE Transactions on Pattern Analysis and Machine Intelligence | publisher=Institute of Electrical and Electronics Engineers (IEEE) | volume=27 | issue=8 | year=2005 | issn=0162-8828 | doi=10.1109/tpami.2005.161 | pages=1239–1253| pmid=16119263 | bibcode=2005ITPAM..27.1239B | s2cid=410716 }}

{{DEFAULTSORT:Swendsen-Wang algorithm}}

[[Category:Monte Carlo methods]]

[[Category:Statistical mechanics]]

[[Category:Critical phenomena]]

[[Category:Phase transitions]]

{{physics-stub}}