The Swendsen-Wang process does not always mix rapidly

Vivek K. Gore and Mark R. Jerrum^{∗} Department of Computer Science University of Edinburgh

21st October, 1996

Abstract

The Swendsen-Wang process provides one possible dynamics for the Qstate Potts model in statistical physics. Computer simulations of this process are widely used to estimate the expectations of various observables (random variables) of a Potts system in the equilibrium (or Gibbs) distribution. The legitimacy of such simulations depends on the rate of convergence of the process to equilibrium, often known as the mixing rate. Empirical observations suggest that the Swendsen-Wang process mixes rapidly in many instances of practical interest. In spite of this, we show that there are occasions on which the Swendsen-Wang process requires exponential time (in the size of the system) to approach equilibrium.

∗Address for correspondence: Department of Computer Science, University of Edinburgh, The King’s Buildings, Edinburgh EH9 3JZ, United Kingdom.

1 Introduction

The Potts model is a natural generalisation of the Ising model to an arbitrary number Q ≥ 2 of states or “spins.” Configurations in the Potts model can be viewed as Q-colourings (not in general proper) of the vertices of an undirected n-vertex graph. One is interested in sampling configurations from a certain distribution, known as the Gibbs distribution, with the aim of obtaining estimates for certain random variables on configurations. For the ferromagnetic Potts model — the focus of this article — the Gibbs distribution assigns greater probability to configurations in which a larger number of pairs of adjacent spins are alike. A precise definition of the Potts model is provided in Section 2.

In the absence of effective direct methods, the usual approach to sampling configurations is via the “Markov chain Monte Carlo” method [11]. The idea is to provide the model with a dynamics by defining an ergodic random walk on configurations whose stationary distribution is the required Gibbs distribution. Provided the walk is rapidly mixing, i.e., converges rapidly to equilibrium, configurations may be efficiently sampled by simulating the walk for a sufficient, but not excessive, number of steps.

A number of different dynamics are possible. The simplest is to move between configurations by changing one spin at a time, with transition probabilities determined by the “Metropolis rule” [17]. It is fairly easy to demonstrate situations in which this random walk takes exponential time (in n, the size of the graph) to approach equilibrium, even in the ferromagnetic case.^{1} A more complicated dynamics, which allows many spins to change in one step, was proposed by Swendsen and Wang [22] and is now widely used in computer simulations.

The Swendsen-Wang process (as we shall call it) appears to converge rapidly to equilibrium in many instances of practical interest. This empirical observation might encourage us to attempt to prove that the mixing time of the process grows not too quickly as a function of n, specifically that it is bounded by a fixed polynomial in n, independent of the other parameters of the system. Such a result would establish the existence of an efficient approximation algorithm — more precisely, a “fully polynomial randomised approximation scheme” or fpras [13, 11] — for computing the partition function of a Q-state ferromagnetic Potts system. Such an algorithm is only known to exist in the case Q = 2 [9].

Our main result (see Proposition 7 for a precise statement) demonstrates that this is a vain hope. For a certain particularly simple family of Potts systems based on the complete graph K_{n} on n vertices (the so-called “Curie-Weiss model”) the Swendsen-Wang process is still far from equilibrium after exponentially many steps. This counterexample is valid for all Q ≥ 3 and for a suitably chosen “coupling constant.” It is an open question whether rapid mixing obtains when Q = 2 (the Ising model), or if the negative result can be extended to more physically realistic instances of the Potts model, for example, 2- or 3-dimensional lattices.

1The antiferromagnetic model, in which adjacent spins tend to be unlike, includes graph colouring as a limit, so rapid convergence cannot be expected for any reasonable dynamics.

2 The Potts model

The Potts model was introduced by R. B. Potts [19] in 1952, has been a focus of much attention in the physics and mathematics communities ever since. Rather than present a detailed historical account of the model here, we refer the interested reader to Baxter [2, Chap. 12].

The problem is easily stated. Consider a collection of sites {1, 2, …, n}, denoted by [n], each pair i, j of which has an associated interaction energy, which, for simplicity, we assume takes on one of just two values, either 0 or J. In most cases of physical interest, the set E of unordered pairs of sites with nonzero interaction energy forms a regular lattice graph ([n], E). A configuration σ = (σ_{1}, …, σ_{n}) is an assignment of spins to sites, where σ_{i} denotes the spin at site i. The number of spins is denoted by Q, where Q ≥ 2; individual spins may simply be denoted by the numbers in [Q]. The energy of a configuration σ is given by the Hamiltonian

H(σ) = ∑_{(i,j)∈E} J(1 − δ(σ_{i}, σ_{j}))

where δ is the Kronecker-δ function which is 1 if its arguments are equal, and 0 otherwise.

The central problem is to compute the partition function

Z = ∑_{σ}exp(−βH(σ))

where β > 0 is what is called the inverse temperature (to be precise, β = 1/kT, where T is the temperature and k is Boltzmann’s constant) and the sum is over all possible configurations σ. Many of the physical properties of the system can be computed from the knowledge of Z. Essentially, Z is the normalising factor in the calculation of probabilities: according to the fundamental theory of statistical mechanics, the probability that the system in equilibrium is found in state σ (the steady state probability) is π(σ) = Z^{−1} exp(−βH(σ)). Moreover, certain logarithmic derivatives of Z correspond to quantities such as mean energy and mean magnetic moment. Singularities in these derivatives (in the limit, as n → ∞) generally correspond to phase transitions, when a small change in a parameter has an observable effect on the macroscopic properties of the system. If a small change in temperature causes a phase transition, then that temperature is called the critical temperature.

Consider the effect of the parameter J in the Hamiltonian. The high probability configurations are those for which H(σ) is low. Let d(σ) denote the number of edges in E that connect sites with different spins in the configuration σ. It is easy to see that H(σ) = Jd(σ), and that

Z = ∑_{σ}exp(−βJd(σ))

We let K = βJ; K is usually known as the coupling constant in the statistical physics community. By the definition of H(σ), if J (or K, since β > 0) is positive, then configurations in which neighbouring spins (spins associated with a pair of sites with non-zero interaction energy) are the same are preferred: this is the ferromagnetic (attractive) case. On the other hand, if J is negative, the neighbouring spins will tend to be different, and this is the antiferromagnetic (repulsive) case.

The search for efficient computational solutions to these problems has proved extremely hard and has generated a vast body of literature. (The reader is referred to Chapter 12 of Baxter [2], where some special cases have been considered.) A huge amount of computational effort has also been poured into numerical solutions, especially for regular lattices.

Although this problem arises in statistical physics, it has a very interesting connection with theoretical computer science. It is another example of a significant combinatorial enumeration problem which is #P-complete, and is hence apparently intractable in exact form. This is an intriguing class of problems, and includes the problems of computing the volume of a convex body and the permanent of a 0-1 matrix. The Potts model also turns out to be one of the many specialisations of the famous Tutte polynomial in graph theory. The reader is referred to Welsh [23, p. 62] for more on this interesting connection.

A lot of research effort has been devoted to finding efficient approximation algorithms for #P-complete problems, where by efficient we mean that the algorithm is guaranteed to run in time polynomial in the number of sites n. Randomness has played a major role in this area, and efficient randomised approximation algorithms have been given for computing the volume of a convex body and estimating the permanent of a dense 0-1 matrix, as well as for many other problems. Each of these algorithms is a fully polynomial randomised approximation scheme (fpras), i.e., one that produces solutions which, with very high probability, fall within arbitrarily small error bounds specified by the user, the price of greater accuracy being a modest increase in runtime.

Most of these algorithms use Markov chain simulation. This approach has been used extensively over the years in the field of physics, and in the last ten years or so, it has been used by researchers in computer science to provide fully polynomial randomised approximation schemes for many problems. For further information on this approach and its applications, refer to the surveys by Kannan [12] and Jerrum and Sinclair [11].

2.1 Single spin flip process

It is, in fact, very easy to use the Markov chain approach to approximate the partition function of a Potts system. The state space Σ is simply the set of all possible configurations. Let π(σ) denote the steady state probability of the configuration σ as described earlier. Transition probabilities from the current state σ are modelled by the following procedure:

(1) choose a site i ∈ [n] and a spin s ∈ [Q], uniformly at random;

(2) assign spin s to site i to get a new configuration σ′, and let the probability of accepting the move, p_{acc} = min{1, π(σ′)/π(σ)};

(3) with probability p_{acc} let the next state be σ′, and with probability 1 − p_{acc} let the next state be σ itself.

The procedure described above is an example of what is known in the statistical physics community as the “Glauber dynamics” or the “single spin flip dynamics” of the Potts model. The acceptance condition used here is called the Metropolis rule, which is familiar in the computer simulation of models in statistical physics [17]. Unfortunately, the Markov chain described above has not been proved to be rapidly mixing. In fact, when Q = 2, i.e., when we have an Ising system, it can be shown that the chain is not rapidly mixing in the ferromagnetic case, when K is sufficiently large. It is well known that ferromagnetic Ising systems typically exhibit a phase transition at a certain value of the parameter β; for values of β above this critical value, the system settles into a state in which there is a preponderance of one or the other of the two spins.

A significant theoretical advance came in 1990 when Jerrum and Sinclair [9] described the first provably efficient approximation algorithm for the partition function of an arbitrary ferromagnetic Ising system. They used a completely different Markov chain, one in which the state space consists of all the subgraphs of the interaction graph G. They proved that this Markov chain is rapidly mixing at all temperatures, and this fact can be used to estimate the partition function of the system. Unfortunately, their approach doesn’t seem to generalise to Q > 2.

The antiferromagnetic case seems to be even harder: in fact, for the Ising model, Jerrum and Sinclair [9] (following Barahona [1]) proved that the existence of even an fpras is highly unlikely. Similar results for the Potts model have been proved by Welsh [23, p. 138]. However, it is worth mentioning that if Q ≥ 2Δ, where Δ is the maximum degree of a vertex in the interaction graph G, then the Markov chain described above can be shown to be rapidly mixing and yields an fpras. This was shown for the zero temperature case (β → ∞) by Jerrum [10], and then extended to arbitrary temperature (at least for Q > 2Δ), by Sokal [21]. The latter result follows directly from the standard proof of the Dobrushin uniqueness theorem [4, 5], combined with Salas and Sokal’s verification [20] of the hypotheses of that theorem for the case of the antiferromagnetic Potts model with Q > 2Δ.

2.2 The random cluster model

Before introducing the Swendsen-Wang approach for approximating the partition function for the Q-state Potts model, it would be instructive to look at a related model, called the Random Cluster model, which was introduced by Fortuin and Kasteleyn [8] in 1972. We have the interaction graph G = ([n], E) just as before, however, there are no spins on the sites. The behaviour of the system depends on the formation of bonds between pairs of interacting sites, and the clusters (connected components in graph terminology) formed by these bonds. There are two parameters associated with the model, a probability p of the formation of a bond between two interacting sites, and a weighting factor Q. (When this Q is a positive integer it corresponds to the Q in the Potts model; however, in the random cluster model, Q may be an arbitrary non-negative real number.) The partition function for this model is given by

Z = ∑_{A⊆E} p^{|A|}(1 − p)^{|E|−|A|}Q^{C(A)}

where A denotes the set of interactions that form a bond, and C(A) is the number of clusters (connected components) in the bond graph A. The sum is over all possible subgraphs of G.

It turns out that the Q-state Potts model is equivalent to the random cluster model with p = 1 − exp(−K), where K is the coupling constant as described earlier, and the parameter Q is common to both models. The equivalence is close, even at the microstate level: to obtain a Potts configuration from a random cluster configuration, simply assign a spin from [Q] independently and u.a.r. to each cluster. Note that the random cluster model effectively generalises the Potts model to an arbitrary (possibly nonintegral) positive number of spins.

2.3 The Swendsen-Wang approach

This approach is based on the “bond graph and cluster” view of the Potts model described in the previous section. Unlike the single spin flip dynamics, this approach is not local in that a single transition can affect a large number of sites. Each transition can be described as a two-step process, and the connection with the random cluster model is very clear. Let the current Potts configuration be denoted by σ. The next configuration σ′ is obtained as follows:

(1) Let A ⊆ E be the subset of edges that form a bond, i.e., ones with the same spin on both incident sites. Each of the edges in A is retained independently with a probability p = 1 − exp(−K); this gives a subset A′ of A.

(2) Using A′ as the set of edges, connected components (clusters) are formed. For each cluster, a spin is chosen uniformly at random from [Q], and all sites within the cluster are assigned that spin.

That the Markov chain with transitions defined by this experiment is ergodic is immediate; that it has the correct distribution is not too difficult to show. (See, for example, Edwards and Sokal [6].)

This Markov chain seems to work very well in practice, however, it has not been proved to be rapidly mixing for any class of graphs. The non-local nature of the transitions seems to allow the chain to move more freely within the state space, thus avoiding the possible constrictions that might result at low temperatures. Our result deals with the case where the interaction graph is the complete graph K_{n}, and we prove that the chain is not rapidly mixing in this case for Q ≥ 3.

3 A first-order phase transition

The slow mixing rate of the Swendsen-Wang process is connected with what is known as a “first-order phase transition,” which we now investigate in the context of the Potts model on the complete graph (Curie-Weiss model). We exhibit two distinct kinds of configurations that account for all but an exponentially small fraction of the partition function Z. In fact, by tuning the coupling constant, we arrange that the two kinds of configurations make a roughly equal contribution to Z. Such a system is said to be in a mixed phase and the two kinds of configurations are called coexisting phases [7]. This value of the coupling constant is usually referred to as its critical value (denoted here by K_{cr}). The phase transition reflects a crucial instability in the model in the following sense. When K > K_{cr}, the system prefers the so called ordered phase (one of the spins dominates). As K is decreased (i.e., temperature is increased), the system goes into a mixed phase at K = K_{cr} and then makes an abrupt transition to the disordered phase (each of the Q spins appears roughly the same number of times) when K < K_{cr}.

We consider the case where K = K_{cr} and in Section 4 we show that the Swendsen-Wang process only very infrequently makes a transition between the coexisting phases, which results in a slow mixing rate.

Consider an arbitrary configuration σ ∈ Σ of the Q-state Potts system. Recall that the equilibrium probability of σ is given by

π(σ) = Z^{−1} exp(−Kd(σ)),

where d(σ) is the number of pairs of sites in E with different spins. We choose K = c/n, where 1 < c < Q is a constant depending on Q. (See (4) for an explicit expression for c.) For the complete graph K_{n}, since all interactions are present, the only relevant observable quantities are the number of different spins and the sizes of these spin classes. Let n = (n_{1}, …, n_{Q}) be the vector whose ith component is the size of the ith spin class of σ; we say that n is the type of σ.

Note that

d(σ) = \frac{1}{2}(n^{2} − ∑_{i=1}^{Q} n_{i}^{2}).  (1)

Since π(σ) is a function only of type, we may write π(n) to denote the value of π(σ) for any configuration σ of type n.

In equilibrium, the probability of being in a configuration of type n is N(n)× π(n), where

N(n) = \binom{n}{n_{1}, n_{2}, …, n_{Q}}

denotes the number of configurations of type n. Using Stirling’s approximation,

N(n) = n^{−(Q−1)/2} exp {( − ∑_{i=1}^{Q} a_{i} ln a_{i})n + Δ(a)}  (2)

where a = (a_{1}, …, a_{n}) = n/n, and Δ(a) is an error term; in general, |Δ(a)| = O(log n), but the tighter estimate |Δ(a)| = O(1) holds if it is known that a ≥ (ε, …, ε) for some constant ε > 0. (The implicit constants depend of Q and ε only.)

From (1), recalling K = c/n, we have

π(n) = Z^{−1} exp \{ − \frac{c}{2}(1 − ∑_{i=1}^{Q} a_{i}^{2})n\}.

Therefore,

Pr(σ has type n) = Z^{−1}n^{−(Q−1)/2} exp{f(a)n + Δ(a)},  (3)

where

f(a) = ∑_{i=1}^{Q} g(a_{i}) − \frac{c}{2}

and g(x) = \frac{1}{2}cx^{2} − x ln x.

In order to identify the configurations that have the largest weights, we need to maximise f, which in turn means that we need to maximise ∑_{i=1}^{Q} g(a_{i}), in the region defined by a_{i} ≥ 0 for all i, and ∑_{i=1}^{Q} a_{i} = 1. This is clearly a closed region (viewed as a set in (Q − 1)-dimensional Euclidean space), and we use R to denote it. We now proceed to look at the behaviour of g in the interval [0, 1]. The following are easy observations:

• if we define g(0) = 0, then g is continuous in [0, 1].

• g′(x) = cx − ln x − 1 is defined in (0, 1], and tends to ∞ as x → 0^{+}.

• g′(x) has a unique minimum in (0, 1] at x = c^{−1}, such that g′(c^{−1}) = ln c.

• g′′(x) > 0 in (c^{−1}, 1] and g′′(x) < 0 in (0, c^{−1}).

Proposition 1 Let a = (a_{1}, ⋯, a_{Q}) be a local maximum point of f. Then a satisfies the following properties:

(i) a lies in the interior of R.

(ii) Either a_{i} = Q^{−1} for all i, or there are α and β such that 0 < α < c^{−1} < β < 1, and a_{i} ∈ {α, β}, for all i.

(iii) If a is such that the a_{i} are not all equal, then there is a unique component a_{j} such that a_{j} = β; the other components a_{i} with i ≠ j satisfy a_{i} = α. Furthermore, g′(α) = g′(β).

Proof

(i) Suppose, on the contrary, that a is such that a_{i} = 0 and a_{j} > 0. Since g′(x) → ∞ as x → a_{i}^{+} and g′(a_{j}) is finite, we can increase f by setting a_{i} = ε and a_{j} = a_{j} − ε, where ε > 0 is sufficiently small.

(ii) At any local maximum, it must be the case that g′(a_{i}) = g′(a_{j}), for all i and j. For suppose g′(a_{i}) ≠ g′(a_{j}), for some i ≠ j. Then a small perturbation of ε to a_{i} and a_{j} (either a_{i} ← a_{i} + ε and a_{j} ← a_{j} − ε or the other way round, depending on the values of g′(a_{i}) and g′(a_{j})) would cause f(a) to increase. Since g′(x) is unimodal in (0, 1], a_{i} ∈ {α, β} for all i where α and β are on either side of the minimum of g′.

(iii) Suppose on the contrary, that a_{j} = a_{k} = β, where j ≠ k. Since g′′(β) > 0, setting a_{j} ← a_{j} − ε and a_{k} ← a_{k} + ε would cause f(a) to increase.

If we now set^{2}

c = \frac{2(Q − 1) ln(Q − 1)}{Q − 2}  (4)

it is routine to verify that the following three choices for a satisfy properties (i)–(iii) in the statement of Proposition 1:

(S1) a_{i} = Q^{−1} for all i = 1, …, Q;

(S2) a_{i} = (Q(Q − 1))^{−1} for all i = 1, …, Q − 1, and a_{Q} = (Q − 1)/Q;

(S3) a_{i} = (2(Q − 1))^{−1} for all i = 1, …, Q − 1, and a_{Q} = 1/2.

2Beyond a certain value c_{0} > 1 of c, the function f can be shown to have two local maximum points and a local minimum point. Our chosen value of c is the unique value at which the two local maxima become equal to yield a global maximum. Much of the analysis has been done in [15], albeit in a somewhat different context.

We shall shortly see that they are the only solutions satisfying properties (i)– (iii). Note that the ordering of subscripts on the spin classes is not significant, as explained earlier.

Claim 2 The first two solutions (S1)–(S2) above are the only local maximum points of f in R, and they both correspond to the global maximum of f. (The final solution (S3) is a local minimum.)

Proof It is clear from Proposition 1 that any maximum point of f should have the form a_{i} = α for all i = 1, …, Q − 1, and a_{Q} = β = 1 − (Q − 1)α, for some α in (0, Q^{−1}] satisfying

h(α) = g′(α) − g′(1 − (Q − 1)α) = 0.

Now h(α) = cQα − c + ln(1 − (Q − 1)α) − ln α, and

h′(α) = cQ − \frac{Q − 1}{1 − (Q − 1)α} − \frac{1}{α} = cQ − \frac{1}{α(1 − (Q − 1)α)}

Setting h′(α) = 0, we get the quadratic equation α(1 − (Q − 1)α) = 1/cQ, which implies that h(α) has at most two turning points (and hence, at most three zeros) in the interval (0, Q^{−1}]. Since α = (Q(Q − 1))^{−1}, α = (2(Q − 1))^{−1} and α = Q^{−1} all satisfy h(α) = 0, they are in fact the only solutions to that equation. We conclude that (S1)–(S3) are the only choices for a consistent with the conditions of Proposition 1, and hence they must cover all the local maximum points of f(a).

We now proceed to show that solutions (S1) and (S2) correspond to the global maximum of f, and that (S3) does not. (In fact it is a local minimum point.) Since, by Proposition 1, we are only interested in solutions of the form a_{i} = α for all i = 1, …, Q − 1, and a_{Q} = β = 1 − (Q − 1)α. we may view f as a function of the single variable α. Accordingly, define f̂(α) = f(α, …, α, β). For the value of c chosen above (4), we have, by direct calculation,

f̂(Q^{−1}) = f̂((Q(Q − 1))^{−1}) = ln Q − \frac{(Q − 1)^{2} ln(Q − 1)}{Q(Q − 2)}  (5)

and

f̂((2(Q − 1))^{−1}) = ln 2 − \frac{Q ln(Q − 1)}{4(Q − 2)}.  (6)

Denote by d(Q) the difference between (5) and (6):

d(Q) = f̂(Q^{−1}) − f̂((2(Q − 1))^{−1}) = ln Q − ln 2 − \frac{(3Q − 2) ln(Q − 1)}{4Q}.

Then d(Q) > 0 for all Q ≥ 3. To verify this, observe that for Q ≥ 16,

d(Q) > ln Q − ln 2 − \frac{3Q ln Q}{4Q} = ln Q − ln 2 − \frac{3}{4}ln Q ≥ 0.

The cases 3 ≤ Q ≤ 15 can be checked separately. Thus (S1) and (S2) are global maximum points of f and are the only such.

Note that d(Q) = f(Q^{−1}) − f((2(Q − 1))^{−1}) is very small for smaller values of Q (e.g., of the order of 10^{−3} for Q = 3); however, since f appears in the exponent and is multiplied by n, the actual difference in weights of these two types of configurations is quite large even for fairly small values of n.

Denote by B_{=}(ε) the set of all points in R that are within (Euclidean) distance ε of the balanced maximum point (Q^{−1}, …, Q^{−1}), and by B_{≠}(ε) the points that are within distance ε of any one of the other maximum points. (B_{=}(ε) is a (Q − 1)-dimensional ball and B_{≠}(ε) a union of Q such balls.) Let Σ_{=}(ε) ⊂ Σ (respectively Σ_{≠}(ε)) be the set of configurations whose type n lies in nB_{=}(ε) (respectively nB_{≠}(ε)). The following result summarises what we have discovered.

Proposition 3 For any ε > 0:

(i) Pr(σ ∈ Σ_{=}(ε)) = Ω(n^{−(Q−1)});

(ii) Pr(σ ∈ Σ_{≠}(ε)) = Ω(n^{−(Q−1)}); and

(iii) Pr(σ ∉ Σ_{=}(ε) ∪ Σ_{≠}(ε)) = e^{−Ω(n)}.

The implicit constants depend only on Q and ε.

Proof Let a ∈ R be chosen so that every component of a has value either i/n or (i + 1)/n for some integer i, and let n_{0} = na; observe that all configurations σ of type n_{0} are in Σ_{=}(ε), provided n is sufficiently large. Then |f(a) − f(Q^{−1}, …, Q^{−1})| = O(n^{−1}), and hence, by (3), n_{0} comes within a constant factor of maximising Pr(σ is of type n) over all types n. (Note that we are operating within the |Δ(a)| = O(1) regime.) Since the total number of distinct types is O(n^{(Q−1)}) we have part (i). Part (ii) is proved in a similar manner.

Finally note that the supremum of f(a) over the region a ∈ R\(Σ_{=}(ε)∪Σ_{≠}(ε)) is strictly less than the supremum over the whole of R. Part (iii) follows by combining this observation with (3).

4 Dynamics

It is clear from Proposition 3 that the single spin flip process described in Section 2.1 will converge only very slowly to equilibrium, since it is difficult to escape from either of the neighbourhoods Σ_{=}(ε) or Σ_{≠}(ε) using small steps. However, the Swendsen-Wang process is able to change large blocks of spins in one step, which at first sight seems to give it a significant advantage. Our main result (Proposition 7 below) suggests that this advantage may on occasion be illusory.

Before we present a formal proof, it would be useful get an intuitive feel as to why we expect Proposition 7 to be true. Let σ denote the current Potts configuration. Note that the Swendsen-Wang process only considers edges that form a bond, so that the configuration may be viewed as a collection of smaller complete graphs, one for each spin in [Q]. Let n_{1}, n_{2}, …, n_{Q} denote the sizes of these graphs. Let G_{ν,p} denote the standard random graph model, in which an undirected ν-vertex graph is formed by adding, independently for each unordered pair of vertices u, v, an edge connecting u and v with probability p. Step (1) of the process essentially creates Q random graphs, G_{n_{i},p}, one of each size n_{i}, 1 ≤ i ≤ Q, where the probability of retaining an edge is p = 1 − exp(−K). Recall that K = c/n, where 1 < c < Q, so that for large n, p ≈ c/n. For a carefully chosen value of c, we expect two possibilities:

• Just prior to Step (1), if all the spin classes in σ have roughly the same size ν (≈ n/Q), then, for any such class, the probability p of retaining an edge can be written as d/ν, where d ≈ c/Q, so that d < 1. A well-known result in the theory of random graphs [3] tells us that with very high probability, the spin class will break up into very small components (of size O(log ν)) so that at the end of Step (2), after assigning random spin values to these very small components, with high probability, we again end up with all spin classes having roughly equal size.

• Just prior to Step (1), if there is one very large spin class in σ and all the other spin classes are very small, then the value of d (as above) would be > 1 for the large class whereas it would be ≪ 1 for the other classes. We can now appeal to results about “the giant component” in random graphs [3] to say that with high probability, at the end of Step (1) there would be one large component and all the other components would be very small. This means that at the end of Step (2), we expect, with high probability, a configuration similar to the one before Step (1). The choice of c determines how close the new configuration would be to the previous one.

Proposition 7 shows that for a carefully chosen value of c (depending on Q), a Q-state Potts system tends to settle in one of the two kinds of configurations mentioned above, and the probability of making a transition from one kind to the other is very small. Our proof utilises some standard bounds on the tails of distributions of sums of independent r.v’s that we state here for convenient reference.

Lemma 4 [Chernoff.] Let the random variable Z have distribution Bin(ν, p), where Bin(ν, p) is the binomial distribution with parameters ν and p. Then for any real γ > 0,

Pr(Z > γνp) < (\frac{e^{γ−1}}{γ^{γ}} )^{νp}.

Proof See [18, Theorem 4.1].

Lemma 5 [Hoeffding.] Let Z_{1}, …, Z_{k} be independent r.v’s with a_{i} ≤ Z_{i} ≤ b_{i}, for suitable constants a_{i}, b_{i}, and all 1 ≤ i ≤ k. Also let Ẑ = ∑_{i=1}^{k} Z_{i}. Then for any λ > 0,

Pr (|Ẑ − Exp Ẑ| ≥ λn) ≤ exp(−2λ^{2}n^{2} / ∑_{i=1}^{k}(b_{i} − a_{i})^{2})

Proof See [16, Theorem 5.7].

Let G_{ν,p} denote the standard random graph model as before. Suppose that p < d/ν, with d < 1 a constant, and G is selected according to the model G_{ν,p}. It is a classical result that, with probability tending to 1 as ν → ∞, the connected components of G all have size O(log ν). We require a (fairly crude) large deviation version of this result.

Lemma 6 Let G be selected according to the model G_{ν,p}, where p < d/ν and 0 < d < 1 is a constant. Then the probability that G contains a component of size exceeding \sqrt{ν} is exp(−Ω(\sqrt{ν})).

Proof Following Karp [14], we consider a simple stochastic procedure for growing a connected component of G from specified vertex s. Let D_{0} = {s} and P_{0} = ∅. At step t, D_{t} will be the set of “discovered” vertices (those that have been shown to be connected to s), and P_{t} ⊆ D_{t} the set of “processed” vertices. If D_{t} = P_{t} we are done: D_{t} is the connected component of G containing s. Otherwise, we select v ∈ D_{t} \ P_{t} and let D_{t+1} = D_{t} ∪ G(v) and P_{t+1} = P_{t} ∪ {v}, where G(v) denotes the set of neighbours of v in G. Note that the termination condition is equivalent to |D_{t}| = t.

We must show that this process terminates within \sqrt{ν} steps with very high probability. We do this by comparing the evolution of D_{i} against another sequence of random variables (X_{i}) defined by X_{0} = 1 and X_{t+1} = X_{t} + Bin(ν, d/ν), where Bin(ν, p) is the binomial distribution with parameters ν and p. If Y and Z are random variables taking non-negative integer values, then Y is said to stochastically dominate Z if Pr(Y ≤ i) ≤ Pr(Z ≤ i) for all i ≥ 0. It is an elementary fact that stochastic domination is preserved by addition. Now it is easily checked that X_{t+1} − X_{t} stochastically dominates |D_{t+1}| − |D_{t}|, and so X_{t} stochastically dominates |D_{t}|. But X_{t} − 1 is clearly distributed as Bin(tν, d/ν) so we can estimate the probability that X_{t} is large using a Chernoff bound. In particular, letting γ = (t − 1)/td in Lemma 4, the probability that X_{t} exceeds t is bounded as follows:

Pr(X_{t} > t) = Pr(X_{t} − 1 > t − 1)
= Pr(X_{t} − 1 > γtd)
= Pr(X_{t} − 1 > γ Exp(X_{t} − 1))
≤ (\frac{e^{γ−1}}{γ^{γ}} )^{td}.

Setting t = ⌊\sqrt{ν}⌋, and noting that γ → d^{−1} > 1 as ν → ∞, we obtain

Pr(|D_{t}| > t) ≤ Pr(X_{t} > t) = exp(−Ω(\sqrt{ν})),

where we have used the fact that X_{t} stochastically dominates |D_{t}|. But |D_{t}| > t is the event that the component building procedure has not terminated at or before time t, i.e., that the connected component of G containing s has size greater than t = ⌊\sqrt{ν}⌋. Multiplying by ν we obtain a bound on the probability that any connected component in G has size exceeding \sqrt{ν}; this small extra factor may be absorbed by the Ω-notation.

Proposition 7 Suppose Q ≥ 3 is an integer, c = 2(Q − 1)(Q − 2)^{−1} ln(Q − 1), and consider a Potts system on K_{n} with coupling constant K = c/n. Let ε > 0 be sufficiently small, and let Σ_{=}(ε) and Σ_{≠}(ε) be as in Proposition 3. Starting at any configuration σ(0) ∈ Σ_{=}(ε), the expected time T for the Swendsen-Wang process to reach a configuration σ(T) ∈ Σ_{≠}(ε) is exp(Ω(\sqrt{n})).

Proof Suppose the configuration σ(t) at time t is an arbitrary member of Σ_{=}(ε). By definition of Σ_{=}(ε), the size ν of any spin-class of σ(0) is bounded above by ν ≤ (Q^{−1} + ε)n. Focussing attention on a particular spin-class of size ν, the set A constructed in step (1) of the Swendsen-Wang process is the edge set of a complete graph K_{ν} on ν vertices, and the set A′ is the edge set of a random graph G selected according to the model G_{ν,p}, where

p = 1 − e^{−K} ≤ \frac{c}{n} ≤ (\frac{1}{Q} + ε)\frac{c}{ν}

Since cQ^{−1} < 1 we have p ≤ d/ν where d < 1, provided ε is sufficiently small. By Lemma 6, with probability 1 − exp(−Ω(\sqrt{ν})), all connected components of G have size at most \sqrt{ν}. Since (Q^{−1} − ε)n ≤ ν ≤ n, the same statement holds with n replacing ν. Similar arguments apply to the other spin-classes, so, with probability 1 − exp(−Ω(\sqrt{n})), all the connected components formed in step (1) of the Swendsen-Wang process have size at most \sqrt{n}.

Let κ be the number of such components, and s_{1}, …, s_{κ} be their respective sizes. The expected size of a spin-class constructed in step (2) of the Swendsen-Wang process is n/Q, and because there are many components (at least \sqrt{n}) we expect the actual size of each spin-class to be close to the expectation. We quantify this intuition by appealing to a Hoeffding bound. Fix a spin γ, and define the random variables Y_{1}, …, Y_{κ} and Ŷ by

Y_{i} = \{ \begin{matrix}s_{i}, & if the ith component receives spin γ in step (2); \\ 0, & otherwise, \end{matrix}

and Ŷ = ∑_{i=1}^{κ} Y_{i}. Then Exp Ŷ = n/Q and, by Lemma 5, for any λ > 0

Pr (|Ŷ − Exp Ŷ | ≥ λn) ≤ exp(−2λ^{2}n^{2} / ∑_{i=1}^{κ} s_{i}^{2})
≤ exp(−2λ^{2}\sqrt{n}),

since

∑_{i=1}^{κ} s_{i}^{2} ≤ ∑_{i=1}^{κ} s_{i}\sqrt{n} = n^{3/2}.

Similar bounds apply, of course, to the other spins. Choosing λ = ε/\sqrt{Q} we see that, with probability 1 − exp(−Ω(\sqrt{n})), the size of every spin-class in σ(t + 1) lies in the range ((Q^{−1} − ε/\sqrt{Q})n, (Q^{−1} + ε/\sqrt{Q})n); but this condition implies σ(t + 1) ∈ Σ_{=}(ε). The claimed result follows easily.

Note that Proposition 7 only proves that starting from any configuration in Σ_{=}(ε), the Swendsen-Wang process requires exponential time to reach a configuration in Σ_{≠}(ε). However, since the Swendsen-Wang process is reversible, it is clear that the same statement also holds if we interchange the roles of Σ_{=}(ε) and Σ_{≠}(ε) in the Proposition.

Acknowledgements

We thank Alan Frieze and Boris Pittel for invaluable pointers to relevant literature, Peter Winkler for an enlightening discussion on Gibbs measures and phase transitions, and Alan Sokal for his careful reading and constructive criticism of an early version of this paper.

References

[1] F. Barahona, On the computational complexity of Ising spin glass models Journal of Physics A: Mathematical and General 15 (1982), pp. 3241–3253.

[2] R. J. Baxter, Exactly Solved Models in Statistical Mechanics, Academic Press, London, 1982.

[3] B. Bollobás, Random Graphs, Academic Press, London, 1985.

[4] R. L. Dobrushin, The description of a random field by means of conditional probability and conditions of its regularity, Theory of Probability and its Applications 13 (1968), pp. 197–224.

[5] R. L. Dobrushin, Prescribing a system of random variables by conditional distributions, Theory of Probability and its Applications 15 (1970), pp. 458– 486.

[6] R. G. Edwards and A. D. Sokal, Generalizations of the Fortuin-Kasteleyn-Swendsen-Wang representation and Monte Carlo algorithm, Physical Review D 38 (1988), pp. 2009–2012.

[7] R. S. Ellis, Entropy, large deviations, and statistical mechanics, Grundlehren der mathematischen Wissenschaften 271, Springer-Verlag, 1985.

[8] C. M. Fortuin and P. W. Kasteleyn, On the random cluster model I: introduction and relation to other models, Physica 57 (1972) pp. 536–564.

[9] M. Jerrum and A. Sinclair, Polynomial-time approximation algorithms for the Ising model, SIAM Journal on Computing 22 (1993), pp. 1087–1116.

[10] M. Jerrum, A very simple algorithm for estimating the number of kcolourings of a low-degree graph, Random Structures and Algorithms 7 (1995), pp. 157–165.

[11] M. Jerrum and A. Sinclair, The Markov chain Monte Carlo method: an approach to approximate counting and integration. In Approximation Algorithms for NP-hard Problems (Dorit Hochbaum, ed.), PWS, 1996, pp. 482– 520.

[12] R. Kannan, Markov chains and polynomial time algorithms. Proceedings of the 35th Annual IEEE Symposium on Foundations of Computer Science, Computer Society Press, 1994, pp. 656–671.

[13] R. M. Karp and M. Luby, Monte-Carlo algorithms for enumeration and reliability problems, Proceedings of the 24th IEEE Symposium on Foundations of Computer Science, Computer Society Press, 1983, pp. 56–64.

[14] R. M. Karp, The transitive closure of a random digraph, Random Structures and Algorithms 1 (1990), pp. 73–93.

[15] H. Kesten and R. H. Schonmann, Behavior in large dimensions of the Potts and Heisenberg models, Reviews in Mathematical Physics 1 (1990), pp. 147– 182.

[16] C. McDiarmid, On the method of bounded differences, London Mathematical Society Lecture Note Series 141, Cambridge University Press, 1989, pp. 148– 188.

[17] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller and E. Teller, Equation of state calculations by fast computing machines, Journal of Chemical Physics 21 (1953), pp. 1087–1092.

[18] R. Motwani and P. Raghavan, Randomized Algorithms, Cambridge University Press, 1995.

[19] R. B. Potts, Some generalised order-disorder transformations, Proceedings of the Cambridge Philosophical Society 48 (1952), pp. 106–109.

[20] J. Salas and A. D. Sokal, Absence of phase transition for antiferromagnetic Potts models via the Dobrushin uniqueness theorem, Journal of Statistical Physics. (To appear.)

[21] A. Sokal, Personal communication, January 1995.

[22] R. H. Swendsen and J-S. Wang, Non-universal critical dynamics in Monte-Carlo simulations, Physical Review Letters 58 (1987), pp. 86–88.

[23] D. J. A. Welsh, Complexity: knots, colourings and counting, London Mathematical Society Lecture Notes 186, Cambridge University Press, 1993.