METROPOLIS ALGORITHM AND EQUIENERGY SAMPLING
FOR TWO MEAN FIELD SPIN SYSTEMS
FEDERICO BASSETTI AND FABRIZIO LEISEN
Abstract. In this paper we study the Metropolis algorithm in connection
with two mean–field spin systems, the so called mean–field Ising model and
the Blume–Emery–Griffiths model. In both this examples the naive choice of
proposal chain gives rise, for some parameters, to a slowly mixing Metropo-
lis chain, that is a chain whose spectral gap decreases exponentially fast (in
the dimension N of the problem). Here we show how a slight variant in the
proposal chain can avoid this problem, keeping the mean computational cost
similar to the cost of the usual Metropolis. More precisely we prove that,
with a suitable variant in the proposal, the Metropolis chain has a spectral
gap which decreases polynomially in 1/N . Using some symmetry structure of
the energy, the method rests on allowing appropriate jumps within the energy
level of the starting state, and it is strictly connected to both the small world
Markov chains of [15, 16] and to the equi-energy sampling of [22] and [26].
1. Introduction.
The Metropolis algorithm, introduced in [29] and later generalized in [18], is
currently (together with other Monte Carlo Markov Chain methods) one of the
most used simulation techniques both in statistics and in physics. See, among
others, [33, 32, 39, 17, 35, 34, 25, 6].
In a finite setting the Metropolis algorithm can be described as follows. Suppose
that, given a probability π(x) on a finite set X , want to approximate
(1.1) µ =
f(x)π(x),
for f : X → R. As a first step, take a reversible Markov chain K(x, y) (the proposal
chain) on X and change its output in order to have a new chain with stationary
distribution π. This can be achieved by constructing a new (π–reversible) chain
(1.2) M(x, y) =
K(x, y)A(x, y) x 6= y
K(x, x) +
z 6=xK(x, z)(1−A(x, z)) x = y
where A(x, y) := min(
π(y)K(y,x)
π(x)K(x,y)
, 1). Then, the metropolis estimate of µ is given by
(1.3) µ̂n =
f(Yi),
where Y0 is generated from some initial distribution π0 and Y1, . . . , Yn fromM(x, y).
It is clear that, from a computational point of view, the speed of convergence to
the stationary distribution and the (asymptotic) variance of the estimate are two
very important features of the Markov chain M .
It is well-known that in some situation a Markov chain can converge very slowly
to its stationary distribution and, moreover, that the asymptotic variance of the es-
timate (1.3) can be much bigger than the variance of f , i.e. V arπ(f) :=
x(f(x)−
Key words and phrases. asymptotic variance, Chain decomposition theorem, fast/slowly mix-
ing chain, mean-field Ising model, Metropolis, spectral gap analysis.
http://arxiv.org/abs/0704.0906v2
2 FEDERICO BASSETTI AND FABRIZIO LEISEN
µ)2π(x), which is equal to the asymptotic variance of the crude Montecarlo estima-
tor. In these cases (1.3) turns out to be a very inefficient estimate of µ.
For the Metropolis chain a classical situation in which the convergence is slow
(and the variance big) is when the target distribution π has many peaks and K is
somehow too “local”.
This is well known in statistical physics, where, typically, a distribution of a
system with energy function h and in thermal equilibrium at temperature T is
described by the Gibbs distribution
πh,T (x) = exp{−h(x)/T }Z−1T
with ZT =
x exp{−h(x)/T }. In point of fact, the Metropolis algorithm has been
proposed in [29] to compute average with respect to such distributions. Indeed,
if h is nice, the Metropolis algorithm is very efficient, but it can perform very
poorly if the energy has many local minima separated by high barriers that cannot
be crossed by the proposal moves K. This problem can be bypassed, for specific
energy, designing appropriate moves that have higher chance to cut across the
energy barrier (see, e.g, [4, 5]), or constructing clever alternative approaches to the
problem, for instance using a reparametrization of the problem (see, e.g., [12, 13])
or using auxiliary variables (see, e.g., [40, 9, 1, 30]). A different kind of solution has
been proposed in [14] and in [28] by introducing the so called simulated tempering,
which essentially means that T is changed (stochastically or not) to flatten h. A
remarkable variant of these methods is the parallel tempering, see, for instance, [19].
More recently new algorithms based on the so called equi–energy levels sampling
have been proposed (see [26] and [22]). In particular, the algorithm proposed in [22]
relies on the so–called equi-energy jump, which enables the chain to reach regions
of the sample space with energy close to the one of the starting state, but that may
be separated by steep energy barriers. In point of fact, even if, according to some
simulations, the method seems to be efficient nothing has been formally proved.
Finally, let us mention a recent algorithm, called small world Markov chains (see
[15, 16]), that combine a local chain with long jumps. In these papers, it has
been shown that a simple modification of the proposal mechanism results in faster
convergence of the chain. That mechanism, which is based on an idea from the
field of small-world networks, amounts to adding occasional wild proposals to any
local proposal scheme.
In the present paper we study two simple examples: the so called mean field Ising
model and the mean field Blume–Emery–Griffiths model. As for the former, it is
well-known that the usual choice of K gives rise, for low temperature, to a slowly
mixing Metropolis chain (see, e.g., [26]). Here we show that a slight variant in the
proposal chain can completely solve this problem, keeping the mean computational
cost similar to the cost of the usual Metropolis. The idea again rests on allowing
appropriate jumps in the same energy level of the starting state. As for the Blume–
Emery–Griffyths mean–field model, we first show that there is a critical region of
the parameters space for which the naive Metropolis chain is slowly mixing. Then
we show how one can modify the proposal chain in order to obtain a better mixing
for the Metropolis chain. The present paper should be intended as a further step in
the direction of a better mathematical understanding of both small world Markov
chains and equi-energy sampling.
The rest of the paper is organized as follows. In Section 2 some general consid-
erations are given. In Section 3 some basic tools concerning Markov chain, which
will be used in the paper, are reviewed. Section 4 contains a warming up example.
In Section 5 the mean field Ising model is treated, while Section 6 deals with the
more complex case of the mean field Blume-Emery-Griffiths model. All the proofs
are deferred to the Appendix.
METROPOLIS ALGORITHM AND EQUIENERGY SAMPLING 3
2. A general strategy
In an abstract setting, what we shall do in the next examples can be summarized
as follows. Let G be a group acting on X for which
(2.1) π(x) = π(g(x)) ∀ x ∈ X , ∀ g ∈ G.
For every x in X let Ox := {y = g(x) : g ∈ G} be the orbit of x (of course if y
belongs to Ox then Ox = Oy).
Assume now that we have a reversible Markov chain KE(x, y) (the proposal) on
X and suppose that the Metropolis chain ME with proposal KE is slowly mixing
(see next section for more details). To speed up the mixing one can try to exploit
(2.1) by taking a proposal of the following form:
(2.2) Kǫ(x, y) = ǫKE(x, y) + (1 − ǫ)KG(x, y)
where
KG(x, y) =
qx(z)Iz(y),
0 < qx(z) < 1 and
qx(z) = 1.
In point of fact, usually KE is “local”; for instance frequently
KE(x, y) = 0
whenever y 6= x belongs to Ox, hence with KG we are adding “long” jumps to the
chain. Moreover, note that if KE is such that KE(x, g(x)) = KE(g(x), x), for every
x in X and g in G, then the Metropolis always accepts the move x→ g(x) and
M(x, g(x)) = ǫKE(x, g(x)) + (1− ǫ)qx(g(x)).
In particular this holds when KE is symmetric.
The heuristics under (2.2) is to combining small world Markov chains and equi-
energy sampling.
Before presenting some examples in which one can actually improve the perfor-
mances of the Metropolis chain using this idea, we collect in the next section some
useful facts concerning Markov chains.
3. Preliminaries
Let P (x, y) be a reversible and ergodic Markov chain on the finite set X with
(unique) stationary distribution p(x). Thus, p(x)P (x, y) = p(y)P (y, x). Let L2(p) =
{f : X → R} with < f, g >p= Ep(fg) =
x f(x)g(x)p(x). Reversibility is equiva-
lent to P : L2 → L2 being self–adjoint. Here Pf(x) =
y f(y)P (x, y). The spec-
tral theorem implies that P has real eigenvalues 1 = λ0(P ) > λ1(P ) ≥ λ2(P ) ≥
· · · ≥ λ|X |−1(P ) > −1 with orthonormal basis of eigen–functions ψi : X → R
(Pψi(x) = λiψi(x), < ψi, ψj >p= δij).
3.1. Spectral gap, variance and speed of convergence. A very important
quantity related to the eigenvalues is the spectral gap, defined by
Gap(P ) = 1−max{λ1, |λ|X |−1|}.
It turns out that the spectral gap is a good index to measure the mixing of a
chain. To better understand this point, assume that f belongs to L2(p) and write
f(x) =
i≥0 aiψi(x) (with ai =< f, ψi >p). Now let Y0 be chosen form some
distribution p0 and Y1, . . . , Yn be a realization of the P (x, y) chain, then
µ̂n =
f(Yi)
4 FEDERICO BASSETTI AND FABRIZIO LEISEN
has asymptotic variance given by
AV ar(f, p, P ) := lim
n · V ar(µ̂n) =
|ak|2
1 + λk
1− λk
See, for instance, Theorem 6.5 in Chapter 6 of [3]. From the last expression, the
classical inequality
(3.1) AV ar(f, p, P ) ≤ 2
1− λ1
V arp(f),
follows easily. The last inequality is the usual way of relating spectral gap to
asymptotic variance and, hence, to the efficency of a chain.
The spectral gap is very important also to give bounds on the speed of con-
vergence to the stationary distribution. For example, if ‖ · ‖TV denotes the total
variation norm, one has
‖δxP k − p‖2TV =
|P k(x,A) − p(x)|
≤ 1− p(x)
4p(x)
(max{λ1, |λ|X |−1|})2k
See, e.g., Proposition 3 in [7]. Another classical bound is
‖p0P k/p− 1‖2,p ≤ Gap(P k)‖p0/p− 1‖2,p
valid for every probability p0. See, for instance, [39].
Roughly speaking one can say that a sequence of Markov chains defined on a
sequence of state space XN is slowly mixing (in the dimension of the problem N)
if the spectral gap decreases exponentially fast in N .
3.2. Cheeger’s inequality. As already recalled, problems of slowly mixing typ-
ically occur when π has two or more peaks and the chain K can only move in
a neighborhood of the starting peak. Usually this phenomenon is called bottle-
neck. A powerful tool to detect the presence of a bottleneck is the conductance and
the related Cheeger’s inequality. Recall that the conductance of a chain P with
stationary distribution p is defined by
h = h(p, P ) := inf
A :p(A)≤ 1
x∈A,y∈Ac
p(x)P (x, y),
and the well-known Cheeger’s inequality is
(3.2) 1− 2h ≤ λ1(P ) ≤ 1−
See, for instance, [3, 37, 7]. Note that, since P is reversible,
(3.3) h ≤ 1
p(x)P (x, y) =
p(y)P (y, x)
for every A such that p(A) ≤ 1/2.
3.3. Chain decomposition theorem. In this subsection we briefly describe a
useful technique to obtain bounds on the spectral gap: the so called chain decom-
position technique. Following [16] assume that A1, . . . , Am is a partition of X .
Moreover, for each i = 1, . . . ,m, define a new Markov chain on Ai by setting
PAi(x, y) := P (x, y) + Ix(y)
P (x, z)
 (x, y ∈ Ai).
PAi is a reversible chain on the state space Ai with respect to the probability
measure
pi(x) := p(x)/p(Ai).
METROPOLIS ALGORITHM AND EQUIENERGY SAMPLING 5
The movement of the original chain among the “pieces” A1, . . . , Am can be de-
scribed by a Markov chain with state space {1, . . . ,m} and transition probabilities
PH(i, j) :=
2p(Ai)
x∈Ai,y∈Aj
P (x, y)p(x)
for i 6= j and
PH(i, i) := 1−
j 6=i
PH(i, j),
which is reversible with stationary distribution
p̄(i) := p(Ai).
A variant of a result of Caracciolo, Pelisetto and Sokal (published in [27]), states
(3.4) Gap(P ) ≥ 1
Gap(PH)
i=1,...,m
Gap(PAi)
holds true, see Theorem 2.2 in [16]. Other results about chain decompositions can
be found, for instance, in [20].
In the next very simple example we shall show how this technique can be used,
starting from a slowly mixing chain, to suggest how to modify the proposal chain
in order to obtain a fast mixing chain.
4. Warming up example
Set X = {−N,−N + 1, . . . , 0, 1, . . . , N} and define a probability measure on X
π(x) =
(θ − 1)θ|x|
2θN+1 + 1− θ
θ being a given parameter bigger than 1. Here we can consider G = {+1,−1} (with
group operation given by the usual product) acting on X by g(x) = gx, hence
Ox = {x,−x}.
Now let KE be a chain defined by
KE(x, x + 1) = 1/2 x 6= N
KE(x, x − 1) = 1/2 x 6= −N
KE(N,N) = KE(−N,−N) = 1/2
KE(x, y) = 0 otherwise
and denote by ME the Metropolis chain with stationary distribution π derived
by KE. It is clear that in this case KE(x, y) = 0 whenever y belongs to Ox.
In this example it is very easy to bound the conductance on ME , indeed, taking
A = {−N, . . . ,−1}, by (3.3), it follows that
h(π,ME) ≤
1− π(0)
Hence,
h(π,ME) ≤ Cθ−N ,
and then (3.2) yields
1− λ1 ≤ 2Cθ−N .
This means that, if f is such that a1 6= 0 and θ > 1, then the asymptotic variance
of f blows up exponentially fast, indeed
AV ar(f, π,ME) ≥ 2Celog(θ)N .
6 FEDERICO BASSETTI AND FABRIZIO LEISEN
Now, instead of KE consider
Kǫ(x, y) = (1− ǫ)KE(x, y) + ǫI{−x}(y)
and let M (ǫ) be the Metropolis chain derived by Kǫ. Decompose X as follows
X = A1 ∪A2 · · · ∪ AN
with A1 = {−1, 0, 1} and Ai = {x ∈ X : |x| = i}, for i > 1. Moreover let
π̄(i) = π(Ai) =
(2θ + 1)/Z for i = 1
2θi/Z for i > 1
where
2θN+1 + 1− θ
(θ − 1)
and set
H (i, j) =
2π(Ai)
l∈Ai,m∈Aj
M (ǫ)(l,m)π(l), M
H (i, i) = 1−
j 6=i
H (i, j).
For i 6= 1, N , one has
H (i, i+ 1) =
2π(Ai)
[M (ǫ)(i, i+ 1)π(i) +M (ǫ)(−i,−i− 1)π(−i)]
and, since π(i) = π(−i) and π(i + 1) ≥ π(i)
H (i, i+ 1) =
In the same way it is easy to see that
H (i, i− 1) =
, i 6= 1, N
H (i, i) = 1−
(1 + θ−1) i 6= 1, N
H (N,N − 1) =
H (N,N) = 1−
H (1, 2) =
4(1 + 1/(2θ))
H (1, 1) = 1−
4(1 + 1/(2θ))
Moreover, for every i 6= 1, M (ǫ)Ai in matrix form is given by
1− ǫ ǫ
ǫ 1− ǫ
and hence
Gap(M
) = 1− |1− 2ǫ|.
While M
is given by
(2θ − 1)(1− ǫ)/(2θ) (1− ǫ)/(2θ) ǫ
(1− ǫ)/2 ǫ (1− ǫ)/2
ǫ (1− ǫ)/(2θ) (2θ − 1)(1− ǫ)/(2θ)
and hence
Gap(M
) = k(θ, ǫ) > −1.
Moreover, since
i6=1,N
H (i, i± 1)),M
H (1, 2),M
H (N,N − 1)
≥ min
(1− ǫ)/(4θ), 1− ǫ
4(1 + 1/(2θ))
=: m(ǫ, θ) > 0
METROPOLIS ALGORITHM AND EQUIENERGY SAMPLING 7
and π̄(i) ≤ 3π̄(j) for every i < j, Lemma A.1 in the appendix yields that
1− λ1(M (ǫ)H ) ≥
m(ǫ, θ)
In the same way, since M
H (i, i+1)+M
H (i, i− 1) ≤ (1− ǫ)M(θ)/4, with M(θ) =
max(1 + θ−1, 2θ/(2θ + 1)) ≤ 2, inequality (A.1) in the Appendix yields that
λN−1(M
H ) ≥ 1−
≥ 1 + ǫ
Hence
Gap(M
H ) ≥
m(ǫ, θ)
and (3.4) yield
Gap(M (ǫ)) ≥
h(θ, ǫ)
for a suitable h. This shows that M (ǫ) is fast mixing for every ǫ > 0 and for every
θ > 1 while ME is slowly mixing for every θ > 1.
5. The mean field Ising model
Let X = {−1, 1}N , N being an even integer. For every β > 0 let π = πβ,N be a
probability on X defined by
π(x) = πβ,N (x) := exp
S2N (x)
Z−1N (β) (x ∈ X )
where
ZN (β) = ZN :=
S2N (x)
is the normalization constant (“partition function”) and
SN (x) :=
xi x = (x1, . . . , xN ).
This is the so called mean field Ising model, or Curie-Weiss model, in which every
particle i, with spin xi, interacts equally with every other particle. It is probably
the most simple but also the most studied example of spin system on a complete
graph. The usual Metropolis algorithm uses as proposal chain
KE(x, y) =
I{x(j)}(y)
where x(j) denotes the vector (x1, . . . ,−xj , . . . , xN ). It has been proved in [26] that,
whenever β > 1,
1− λ1 ≤ Ce−D
where λ1 is the first eigenvalues smaller than 1 of the Metropolis chain ME derived
KE . This yields that the variance of an estimator obtained from this Metropolis
algorithm can blow up exponentially fast in N .
The aim of this section is to show how one can construct a different Metropolis
chain avoiding this problem. In the notation of Section 2, we consider
G = SN × {+1,−1}
(SN being the symmetric group of order N) and we define the action of G on
X = {−1, 1}N by
g(x) = (e · xσ(1), . . . , e · xσ(N)) g = (σ, e).
8 FEDERICO BASSETTI AND FABRIZIO LEISEN
In order to introduce a new proposal, it is useful to write X as the union of its
“energy sets”, that is
X = X0 ∪ X2 ∪ X4 ∪ · · · ∪ XN
where
Xi := {x ∈ X : |SN (x)| = i} (i = 0, 2, . . . , N).
Note that energy takes only even values and that Ox = X|SN (x)|. Moreover, for
i 6= 0, set
X+i := {x ∈ X : SN (x) = i} and X
i := {x ∈ X : SN(x) = −i}.
The new proposal chain will be
K(x, y) = p1KE(x, y) + (1 − p1)K0(x, y) if x ∈ X0
K(x, y) = p1KE(x, y) + p2I{−x}(y) + (1− p1 − p2)Ki(x, y)
if x ∈ Xi, i 6= 0
(5.1)
where p1, p2 belong to (0, 1), p1 + p2 < 1, and
Ki(x, y) = IX+
{x}K+i (x, y) + IX−
{x}K−i (x, y) (i 6= 0).
We shall assume that K±i (K0, respectively) are irreducible, symmetric and aperi-
odic chains on X±i ( X0, respectively).
As a leading example we shall take
K0(x, y) =
) y ∈ X0
K±i (x, y) =
(N−i)/2
) y ∈ X±i ,
(5.2)
that is: a realization of a chain K±i (K0, respectively) is simply a sequence of
independent uniform random sampling from X±i (X0, respectively).
Remark 1. Note that (5.2) is the (n, k)-Bose-Einstein distribution with n = (N +
i)/2 and k = (N − i)/2 + 1 and recall that there is a very easy way to directly
generate Bose-Einstein configurations. One may place n balls sequentially into k
boxes, each time choosing a box with probability proportional to its current content
plus one. Starting from the empty configuration this results in a Bose-Einstein
distribution for every stage.
Now let M be the Metropolis chain defined by the transition kernel (1.2) with
K as in (5.1), i.e. for every x in X±i (i 6= 0)
M(x, y) =
if y = x(j), j = 1...N
p2 if y = −x
(1− p1 − p2)K±i (x, y) if y ∈ X
i , y 6= x
z 6=xM(x, z) if y = x
while for x in X0
M(x, y) =
if y = x(j), j = 1...N
(1− p1)K0(x, y) if y ∈ X0, y 6= x
z 6=xM(x, z) if y = x.
METROPOLIS ALGORITHM AND EQUIENERGY SAMPLING 9
By construction M is an aperiodic, irreducible and reversible chain with stationary
distribution π. Then, when (5.2) holds true,
M(x, y) =
if y = x(j), j = 1...N
p2 if y = −x
(1− p1 − p2) 1( N(N−i)/2)
if y ∈ X±i , y 6= x
z 6=xM(x, z) if y = x
for x in X±i (i 6= 0), while if x belongs to X0
M(x, y) =
if y = x(j), j = 1...N
(1− p1) 1( NN/2)
if y ∈ X0, y 6= x
z 6=xM(x, z) if y = x.
In order to bound the spectral gap ofM we shall use the decomposition theorem
described in Subsection 3.3. To this end, for every i = 0, 2, . . . , N and every j 6= i
P̄ (i, j) :=
2π(Xi)
M(x, y)π(x)
P̄ (i, i) := 1−
j 6=i
P̄ (i, j).
As already noted, P̄ is a reversible chain on {0, 2, . . . , N} with stationary distribu-
π̄(i) := π(Xi).
Moreover define for every i = 0, 2, . . . , N a chain on Xi setting
PXi(x, y) :=M(x, y) + Ix(y)
z∈X c
M(x, z)
where both x and y belong to Xi. In the same way, define chains on X+i and X
for i = 2, . . . , N setting
(x, y) := PXi(x, y) (y 6= x, x, y ∈ X±i )
(x, x) := 1−
y 6=x
PXi(x, y).
These chains are reversible on Xi (X±i , respectively) and have as stationary distri-
butions
πXi(x) :=
π(Xi)
and πX±
(x) :=
πXi(x)
πXi(X±i )
|X±i |
10 FEDERICO BASSETTI AND FABRIZIO LEISEN
respectively. Finally, for every i = 2, 4, . . . , N , define a chain on {+,−} setting
Pi(+,−) :=
2πXi(X+i )
PXi(x, y)πXi(x)
Pi(−,+) :=
2πXi(X−i )
PXi(x, y)πXi(x).
Now the lower bound (3.4), applied two times yields
Gap(M) ≥ 1
Gap(P̄ ) min
i=0,2,...,N
{Gap(PXi)}
Gap(P̄ )min
Gap(PX0),
i=2,...,N
Gap(Pi)min{Gap(PX+
), Gap(PX−
(5.3)
Hence, to get a lower bound on Gap(M) it is enough to obtain bounds on the gaps
of the chains P̄ , PX0 , Pi, PX±
The most important of these bounds is given by the following
Proposition 5.1. P̄ is a birth and death chain on {0, 2, . . . , N}, more precisely
(5.4)
P̄ (0, 2) = p1
P̄ (i, i+ 2) = p1
i 6= N, 0
P̄ (i, i− 2) = p1
exp{2β(1− i)/N} i 6= 0.
Moreover
λ1(P̄ ) ≤ 1−
(N/2 + 1)3
λN/2(P̄ ) ≥ 1− p1.
The proof of the previous proposition is based on a bound for a birth and death
chain, given in the Appendix, which can be of its own interest.
As for the others chains, we have the following
Lemma 5.2. For every i = 2, 4, . . . , N
Gap(PX±
) ≥ (1− p1 − p2)Gap(K±i )
Gap(Pi) = p2,
moreover
Gap(PX0) ≥ (1− p1)Gap(K0).
In this way, using (5.3), we can prove the main result of this section.
Proposition 5.3. Let M be the Metropolis chain derived by the chain K defined
as in (5.1) then
Gap(M) ≥ p1p2
(N/2 + 1)3
[ (1− p1)
Gap(K0),
(1− p1 − p2)
min{Gap(K+i ), Gap(K
i, )}
If K±i and K0 are defined as in (5.2) then
Gap(M) ≥
(N/2 + 1)3
[ (1 − p1 − p2)
(1− p1)
for every β > 0 and N ≥ N0.
METROPOLIS ALGORITHM AND EQUIENERGY SAMPLING 11
Proposition 5.3 shows that the gap is polynomial in 1/N independently of β.
Hence, even when β > 1, the variance of the metropolis estimate obtained with this
proposal can not grow up faster than a polynomial in N .
Note that if in Proposition 5.3 we choose
(5.5) p1 = 1− a/(2N), p2 = a/N
we get
Gap(M) ≥ C
Hence, even with this choice, the Metropolis algorithm is still fast mixing for every
β. It is worth noticing that the mean computational cost of this Metropolis does
not change with respect to the Metropolis which uses the proposal KE. Indeed,
in the case of the usual Metropolis, the computational cost needed to go from Xn
to Xn+1 is O(N), since it is essentially due to a sample of one number among N
numbers (we need to decide which coordinate to flip). In the case of the ”modified”
proposal, things are slight more complex. In this case, at the beginning, we have
an extra “toss”. If with this fist toss we decide to flip at random a coordinate the
cost is still O(N) but if we need to sample from K±i the cost is O(N
2) (in this last
case we need to pick a sample from a Bose-Einstein distribution). Hence, although
our algorithm is ”sometime” more expensive, if we take p1 and p2 as in (5.5), we
get that the mean cost of our algorithm is still O(N).
6. The mean–field Blume-Emery-Griffiths model
The Blume-Emery-Griffiths (BEG) model (see [2]) is an important lattice–spin
model in statistical mechanics, it has been studied extensively as a model of many
diverse systems, including He3 − He4 mixtures as well solid–liquid–gas systems,
microemulsions, semiconductor alloys and electronic conduction models. See, for
instance, [2, 38, 23, 24, 31, 36, 21]. We will focus our attention on a simplified
mean–field version of the BEG model. For a mathematical treatment of this mean–
field model see [10]. In what follows let X := {−1, 0, 1}N , N being an even integer,
and for every β > 0 and K > 0 let πβ,K,N be the probability defined by
π(x) = πβ,K,N(x) = exp{−βRN(x) +
S2N(x)}Z−1N (β,K) (x ∈ X )
where
ZN(β,K) = ZN :=
−βRN (x) +
S2N (x)
is the normalization constant,
SN (x) :=
xi and RN (x) :=
x2i x = (x1, x2, ..., xN ).
A natural Metropolis algorithm can be derived by using the proposal chain
(6.1) KE(x, y) =
[I{x(+j)}(y) + I{x(−j)}(y)]
where x(±j) denotes the vector (x1, . . . , xj ± 1, . . . , xN ), with the convention that
2 = −1 and −2 = 1.
The next proposition shows that there exists a critical region of the parameters
space in which the Metropolis chain is slowly mixing. More precisely, using some
results of [10] it is quite straightforward to proove the following
12 FEDERICO BASSETTI AND FABRIZIO LEISEN
Proposition 6.1. Ler ME be the Metropolis chain (with stationary distribution π)
with proposal chain KE defined in (6.1).Then, there exists a non decreasing function
Γ : (0,+∞) → (0,+∞) with limx→0 Γ(x) = +∞ and limx→∞ Γ(x) = γc ≃ 1.082
such that for every couple of positive parametrs (β,K) with K > Γ(β)
Gap(ME) ≤ Ce−∆N
for suitable constants C = C(γ,K) > 0 and ∆ = ∆(γ,K) > 0.
As in the case of the mean–field Ising model, we intend to by pass the slowly
mixing problem of this Metropolis chain by choosing a different proposal. To un-
derstand which kind of proposal is reasonable, here we choose
G = SN × {+1,−1}
with G acting on X = {−1, 0, 1}N by
g(x) = (e · xσ(1), . . . , e · xσ(N)) g = (σ, e).
At this stage, decompose X as the union of its ”energy sets”, that is
X = X0,0 ∪ X1,1 ∪ X0,2 ∪ X1,3 ∪X3,3 ∪ ... ∪ X0,N ∪ X2,N ∪ ...XN,N
where
Xs,r := {x ∈ X : |SN | = s and RN (x) = r}
r = 0, 1, 2, ..., N and s = 1, 3, ..., r if r is odd and s = 0, 2, ..., N if r is even.
Moreover, for s = 1, 2, ..., N , set
X+s,r := {x ∈ X : SN = s and RN (x) = r}
X−s,r := {x ∈ X : SN = −s and RN (x) = r}.
Note again that Ox = Xs,r with s = SN (x) and r = RN (x). The new proposal
chain will be
K(x, y) = p1KE(x, y) + (1 − p1)K0,r(x, y) if x ∈ X0,r, r = 0, 2, ..., N
K(x, y) = p1KE(x, y) + p2I{−x}(y) + (1− p1 − p2)Ks,r(x, y)
if x ∈ Xs,r, s 6= 0
(6.2)
where p1, p2 belong to (0, 1), p1 + p2 < 1, and
Ks,r(x, y) = IX+s,r{x}K
s,r(x, y) + IX−s,r{x}K
s,r(x, y) (s 6= 0)
K0,r(x, y) =
) y ∈ X0,r
K±s,r(x, y) =
(r−s)/2
) y ∈ X±s,r.
(6.3)
Now let M be the Metropolis chain defined by the transition kernel (1.2) with K
as in (6.2), i.e. for every x in X±s,r (s 6= 0)
M(x, y) =
if y = x(±j), j = 1...N
p2 if y = −x
(1 − p1 − p2) 1(Nr )( r(r−s)/2)
if y ∈ X±s,r, y 6= x
z 6=xM(x, z) if y = x,
METROPOLIS ALGORITHM AND EQUIENERGY SAMPLING 13
while if x belongs to X0,r
M(x, y) =
if y = x(±j), j = 1...N
(1− p1) 1(Nr )( rr/2)
if y ∈ X0,r, y 6= x
z 6=xM(x, z) if y = x.
By construction M is an aperiodic, irreducible and reversible chain with stationary
distribution π.
Also in this case, to bound the spectral gap of M , we shall use the chain decom-
position tools. Let
DN = {(0, 0), (1, 1), (0, 2), (2, 2), (1, 3), (3, 3), (0, 4), (2, 4), (4, 4), ..., (0, N), (2, N), ..., (N,N)}
and, for every couple (s, r), (s̃, r̃) in DN , with (s, r) 6= (s̃, r̃), let
P̄ ((s, r), (s̃, r̃)) :=
2π(Xs,r)
x∈Xs,r
y∈Xs̃,r̃
M(x, y)π(x)
P̄ ((s, r), (s, r)) := 1−
(s̃,r̃) 6=(s,r)
P̄ ((s, r), (s̃, r̃)).
Once again, note that P̄ is a reversible chain on DN with stationary distribution
π̄(s, r) := π(Xs,r).
Moreover, for every (s, r) in DN , define a chain on Xs,r setting
PXs,r (x, y) :=M(x, y) + Ix(y)
z∈X cs,r
M(x, z)
where both x and y belong to Xs,r. In the same way, define chains on X+s,r and X−s,r
for (s, r) in DN , s 6= 0, setting
PX±s,r (x, y) := PXs,r(x, y) (y 6= x, x, y ∈ X
PX±s,r (x, x) := 1−
s,ry 6=x
PXs,r (x, y).
These chains are reversible on Xs,r (X±s,r , respectively) and have as stationary dis-
tributions
πXs,r(x) :=
π(Xs,r)
|Xs,r|
and πX±s,r (x) :=
πXs,r(x)
πXs,r (X±s,r)
|X±s,r|
respectively. Finally, for every (s, r) in DN , s 6= 0, define a chain on {+,−} setting
Ps,r(+,−) :=
2πXs,r(X+s,r)
PXs,r (x, y)πXs,r(x)
Ps,r(−,+) :=
2πXs,r(X−s,r)
PXs,r (x, y)πXs,r(x).
14 FEDERICO BASSETTI AND FABRIZIO LEISEN
At this stage, the lower bound (3.4), applied two times, yields
Gap(M) ≥ 1
Gap(P̄ ) min
(s,r)∈DN
Gap(PXs,r)
Gap(P̄ )min
r=0,2,...,N
Gap(PX0,r )
(s,r)∈DN ,s6=0
Gap(Ps,r)min{Gap(PX+s,r), Gap(PX−s,r)}
(6.4)
To derive from the last bound a more explicit bound we need some preliminary
work. The first result we need is exactly the analogous of Lemma 5.2.
Lemma 6.2. Fore every r = 1, . . . , N
Gap(PX0,r ) ≥ (1− p1)Gap(K0,r) = (1− p1),
moreover, for every (s, r) in DN with s 6= 0,
Gap(P±Xs,r) ≥ (1− p1 − p2)Gap(K
s,r) = (1− p1 − p2).
Finally, for every (s, r) in DN ,
Gap(Ps,r) = p2.
Hence, (6.4) can be rewritten as
(6.5) Gap(M) ≥ Gap(P̄ )p2
min{(1− p1)/2, (1− p1 − p2)/2}.
It remains to bound Gap(P̄ ). Unfortunately the the analogous of Proposition 5.1
is not so simple, hence we shall require an additional hypothesis. In what follows
q|[N ]|(r) :=
(r−2i)2
 if r is even
q|[N ]|(r) : =
(r−2i)2
 if r is odd
r = 0, 1, . . . , N and set
A = {β > 0,K > 0 : ∃N0 such that ∀N ≥ N0, q|[N ]| is unimodal}.
Lemma 6.3. For every (β,K) in A
Gap(P̄ ) ≥ Cp
for a suitable constant C = C(β,K).
Under the same assumptions of the previous Lemma we can state the main
results of this section.
Proposition 6.4. For every (β,K) in A
Gap(M) ≥ C̃p
for a suitable constant C̃ = C̃(β,K).
METROPOLIS ALGORITHM AND EQUIENERGY SAMPLING 15
0 5 10 15
β=2.3, K=0.5,0.6,0.7,0.8
0 5 10 15
β=2.3,2.4,2.5,2.6, K=0.5
0 5 10 15
β=2.3, K=1.2,1.3,1.4,1.5
0 5 10 15
β=2.3,2.4,2.5,2.6, K=1.2
Figure 1. The function q|[N ]| for N = 15 and few values of β and K.
We conjecture that Gap(P̄ ) is polynomial in N for every (β,K) such that β 6=
Γ(K) (where Γ is the function of Proposition 6.1), but we are not able to prove this
conjecture. In point of fact we conjecture that R+×R+ \{(β,K) : Γ(K) = β} ⊂ A.
We plotted q|[N ]| for different N , β andK, and these plotts seem, at least, to confirm
that R+ × R+ \ {(β,K) : |Γ(K)− β| ≤ ǫ} ⊂ A for a suitable small ǫ. In Figure 1
we show the graph of q|[N ]| for few different N , β and K.
Appendix A. The Spectral Gap of a Birth and Death Chain
We derive here some bounds on the eigenvalues of a birth and death chain that we
shall use later. These bounds are obtained using the so called geometric techniques,
see [7]. Let Pn be a birth and death chain on Ωn = {1, . . . , n}. Assume that Pn
is reversible with respect to a probability pn, that is pn(i)Pn(i, j) = pn(j)Pn(j, i).
Moreover let
1 > λ1 ≥ λ2 ≥ . . . λn−1 ≥ −1
the eigenvalues of Pn.
We can now prove the following variant of Proposition 6.3 in [6].
Lemma A.1. If there exist positive constants A, q, B and an integer k such that
Pn(i, i± 1) ≥ An−q (i 6= 1, n)
Pn(1, 2) ≥ An−q
Pn(n, n− 1) ≥ An−q
pn(i) ≤ Bpn(j) i ≤ j ≤ k
pn(j) ≤ Bpn(i) k ≤ i ≤ j
λ1 ≤ 1−
16 FEDERICO BASSETTI AND FABRIZIO LEISEN
Proof. We use the notation and the techniques of [7], see also [3] and [6]. Choose
the set of paths
Γ = {γij = (i, i+ 1, ..., j); i ≤ j; i, j ∈ Ωn}
and for e = (i, i+ 1) (i < n) let
ψ(e) =
pn(i, i+ 1)
γl,m∈Γ
γl,m∋e
|γl,m|
pn(l)pn(m)
pn(i)
where |γ| is the length of the path γ. Setting K := supe ψ(e) one has
λ1 ≤ 1−
(see Proposition 1’ in [7], or Exercise 6.4 page 248 in [3]). So, for our purposes, it
suffices to give an upper bound on K. Assume first that e = (i, i+ 1) with
i < k ≤ n, since |γl,m| ≤ n, it follows that
ψ(e) ≤ n
s≥i+1
pn(r)pn(s)
pn(i)
pn(r)
pn(i)
s≥i+1
pn(s)
pn(s)
≤ nq+2
All the other cases can be treated in the same way. Hence,
ψ(e) ≤ B
and then
λ1 ≤ 1−
As for the smaller eigenvalues, Gershgorin theorem yields that
λn−1 ≥ −1 + 2min
P (i, i).
See, for instance, Corollary 2.1 in the Appendix of [3]. Hence, if there exists a
positive constant D such that
Pn(i, i+ 1) + Pn(i, i− 1) ≤ D/2
for every i, then
(A.1) λn−1 ≥ 1−D.
METROPOLIS ALGORITHM AND EQUIENERGY SAMPLING 17
Appendix B. Proofs
To prove Proposition 5.1 we need first to show that π̄ is essentially unimodal.
Lemma B.1. Let
qN (i) =
i = 0, 2, 4, . . . , N.
For every β < 1 there exists an integer N0 such that for every N ≥ N0
qN (i) ≤ qN (j)
whenever j ≤ i. For every β ≥ 1 there exists an integer N0 such that for every
N ≥ N0
qN (i) ≤ qN (j)
whenever i ≤ j ≤ kN and
qN (i) ≥ qN (j)
whenever kN ≤ i ≤ j, kN being a suitable integer.
Proof. Let ∆N (i) be the ratio
∆N (i) =
qN (i+ 2)
qN (i)
i = 0, 2, 4, ..., N − 2,
so that
∆N (i) =
) exp
(1 + i)
N − i
N + 2 + i
(1 + i)
Setting ∆N (x) =
N+2+x
(1 + x)
, x in [0, N−2], it is enough to prove that
x 7→ ∆N (x) takes the value 1 at most once in [0, N − 2], for sufficiently large N .
To prove this last claim first note that
∆N (0) =
N + 2
1 + 2
= 1− 2
(1 − β) +
− β + 2
Hence, there exists N0 in N such that for N ≥ N0:
β ≥ 1 ⇒ ∆N (0) > 1
β < 1 ⇒ ∆N (0) < 1.
As for the first derivative note that
∆′N (x) =
−2(N + 1) + 2β(N + 2)− 2β
(x2 + 2x)
(N + x+ 2)2
(1 + x)
hence ∆′N (x) = 0 if and only if
−2(N + 1) + 2β(N + 2)− 2β
(x2 + 2x) = 0.
18 FEDERICO BASSETTI AND FABRIZIO LEISEN
Rearranging the last equation as
x2 − 4β
+ 2[(β − 1)N + 2β − 1] = 0
one sees that the roots are
x1,2 = 1±
2β − 1
β − 1
Hence, after setting
r := 1 +
2β − 1
β − 1
N2 and r := 1 +
one has
β < 1 ⇒ ∆′N (x) < 0 ∀x ∈ [0, N − 2]
β > 1 ⇒ ∆′N (x) > 0 for x ∈ [0, r)
∆′N (x) < 0 for x ∈ (r,N − 2]
β = 1 ⇒ ∆′N (x) < 0 for x ∈ [0, r)
∆′N (x) < 0 for x ∈ (r,N − 2]
and this concludes the proof. �
Proof of Proposition 5.1. By direct computations it is easy to prove (5.4). Hence
P (i, i± 2) ≥
4(N + 2)
P (i, i+ 2) + P (i, i− 2) ≤ p1
Now observe that
π̄(0) =
ZN (β)
qN (0)
π̄(i) =
ZN (β)
qN (i) i 6= 0.
Hence, by Lemma B.1, if β < 1
π̄(i) ≤ 2π̄(j)
whenever j ≤ i and N is large enough. While for β > 1
π̄(i) ≤ π̄(j)
whenever i ≤ j ≤ kN and
π̄(i) ≥ π̄(j)
whenever kN ≤ i ≤ j. The thesis follows now by Lemma A.1 and by (A.1). �
In order to prove Lemma 5.2 we recall that by Rayleigh’s theorem
(B.1) 1− λ1(P ) = inf
Ep(f, f)
V arp(f)
: f nonconstant
where
Ep(f, f) :=< (I − P )f, f >p=
(f(x) − f(y))2P (x, y)p(x),
METROPOLIS ALGORITHM AND EQUIENERGY SAMPLING 19
P being a reversible chain w.r.t. p, moreover
(B.2) 1− |λN−1| = inf
x,y(f(x) + f(y))
2P (x, y)p(x)
V arp(f)
: f nonconstant
(see, for instance, Theorem 2.3 in Chapter 6 of [3] and Section 2.1 of [8]). At this
stage set
Pǫ(x, y) := (1 − ǫ)P (x, y) + ǫIx(y).
Hence, (B.1) yields
1− λ1(Pǫ) = inf
f∈L2pf 6=const
x,y(f(x) − f(y))2Pǫ(x, y)p(x)
V arp(f)
= inf
f∈L2pf 6=const
(1− ǫ)
x 6=y(f(x)− f(y))2P (x, y)p(x)
V arp(f)
= (1− ǫ)(1− λ1(P )).
Arguing in the same way and using (B.2) we get
1− |λ|X |−1(Pǫ)| ≥ (1 − ǫ)(1− |λ|X |−1(P )|).
Hence,
(B.3) Gap(Pǫ) ≥ (1− ǫ)Gap(P ).
Proof of Lemma 5.2. Note that
(x, y) = (1 − p1 − p2)K±i (x, y) + (p1 + p2)Ix(y)
and, analogously,
PX0(x, y) = (1 − p1)K0(x, y) + p1Ix(y).
Hence, by (B.3),
Gap(PX±
) ≥ (1− p1 − p2)Gap(K±i )
as well
Gap(PX0) ≥ (1− p1)Gap(K0).
Finally note that Pi is given by
1− p2
1− p2
for every i, hence Gap(Pi) = p2. �
Proof of Proposition 5.3. To prove the first part of the proposition it is enough to
combine Lemma 5.2, Proposition 5.1 and (5.3). To complete the proof observe that
Gap(K±i ) = Gap(K0) = 1, when K
i and K0 are given by (5.2). �
In order to prove Proposition 6.1 we need some results obtained in [10].
Theorem B.2 (Ellis-Otto-Touchette). Let ρN be the distribution of SN (x)/N un-
der πβ,K,N , then ρN satisfies a large deviation principle on [−1, 1] with rate function
Ĩβ,K(z) = Jβ(z)− βKz2 − inf
{Jβ(t)− βKt2}
Jβ(z) = sup
tz − log
1 + e−β(et + e−t
1 + 2e−β
Moreover, if Ẽβ,K := argminĨβ,K, then there exists a non decreasing function Γ :
(0,+∞) → (0,+∞) with limx→0 Γ(x) = +∞ and limx→∞ Γ(x) = γc ≃ 1.082 such
that for every (β,K) with K > Γ(β) then
Ẽβ,K = {±z(β,K) 6= 0}.
20 FEDERICO BASSETTI AND FABRIZIO LEISEN
In particular, for such (β,K) and for every 0 < ǫ < |z(β,K)| there exists a constant
C1 = C1(ǫ, β,K) such that
(B.4) ρ([0, ǫ]) ≤ C1 exp{−
γǫ,β,K}
(B.5) γǫ,β,K = inf
z∈[0,ǫ]
Ĩβ,K(z) > 0.
Proof. For the first part see Theorems 3.3, 3.6 and 3.8 in [10]. As for (B.4)-(B.5),
they are standard consequences of the theory of the large deviations and of the first
part of the proposition, see, e.g., Proposition 6.4 of [11]. �
Proof of Proposition 6.1. We intend to use the Chegeer’s inequality. To do this, let
A := {x : SN (x) < 0}, B := {x : SN (x) > 0}, C := {SN (x) = 0}. First of all note
that, by symmetry, π(A) = π(B) = (1−π(C))/2 ≤ 1/2. The main task is to bound
φ(A) =
π(x)ME(x, y) =
π(y)ME(y, x).
Now, observe that if SN (y) > 1 then ME(y, x) = 0 for every x in A, hence
φ(A) =
y:SN (y)=0
ME(y, x) +
y:SN(y)=1
ME(y, x)
≤ π {y : SN (y) ∈ {0, 1}} .
This yields a bound on the conductance
h = h(π,ME) ≤ φ(A)/π(A) ≤
2π {y : SN (y) ∈ {0, 1}}
1− π{y : SN (y) = 0}
Now by Proposition B.2 we get
h(π,ME) ≤ C2e−∆N
for suitable constants C2 and ∆ > 0. The thesis follows by Cheeger inequality
(3.2). �
Proof of Lemma 6.2. The proof is exactly the same as the proof of Lemma 5.2. �
In order to prove Lemma 6.3 it is convenient to fix some simple properties of the
chain P̄ .
Lemma B.3. P̄ is a random walk on DN . If P̄ ((s, r), (s̃, r̃)) 6= 0,
P̄ ((s, r), (s̃, r̃)) ≥ p1C3
for a suitable constant C3 = C3(β,K), moreover
P̄ ((s, r), (s̃, r̃)) ≤ p1
for every (s, r), (s̃, r̃)) 6= ((0, 0), (1, 1)).
Proof of Lemma B.3. Easy but tedious computations show that
P̄ ((0, 0), (1, 1)) =
1, exp{Kβ
P̄ ((0, N), (1, N − 1)) = p1
P̄ ((0, N), (2, N)) =
METROPOLIS ALGORITHM AND EQUIENERGY SAMPLING 21
P̄ ((0, r), (2, r)) =
r = 0, 2, 4, ..., N − 2
P̄ ((0, r), (1, r − 1)) = p1
r = 0, 2, 4, ..., N − 2
P̄ ((0, r), (1, r + 1)) =
1, exp{Kβ
r = 0, 2, 4, ..., N − 2
P̄ ((s, r), (s + 2, r)) =
(r − s)
(s, r) ∈ DN , 0 < s ≤ N − 2, r ≤ N
P̄ ((s, r), (s − 2, r)) = p1
(r + s) exp{4Kβ
(1 − s)}
(s, r) ∈ DN , 0 < s ≤ N, r ≤ N
P̄ ((s, r), (s + 1, r + 1)) =
(N − r)min
1, exp{Kβ
(2s+ 1)− β}
(s, r) ∈ DN , 0 < s, r ≤ N − 1,
P̄ ((s, r), (s − 1, r + 1)) = p1
(N − r) exp{Kβ
(−2s+ 1)− β}
(s, r) ∈ DN , 0 < s, r ≤ N − 1,
P̄ ((s, r), (s + 1, r − 1)) = p1
(r − s)
(s, r) ∈ DN , 0 < r ≤ N, 0 < s ≤ N − 2
P̄ ((s, r), (s − 1, r − 1)) = p1
(r + s)min
1, exp{Kβ
(2s+ 1)− β}
(s, r) ∈ DN , 0 < r ≤ N, 0 < s ≤ r.
At this stage the statement follows easily. �
Proof of Lemma 6.3. In order to obtain a bound on the gap of P̄ we shall apply
another time the decomposition technique. Write
DN = X̄1 ∪ X̄2 ∪ X̄3 ∪ ... ∪ X̄N ,
where
X̄1 = {(0, 0), (1, 1)} X̄r = {(u1, u2) ∈ Dn : u2 = r}.
On |[N ]| := {1, ..., N} define a chain P|[N ]| setting
P|[N ]|(i, j) :=
2π̄(X̄i)
a∈X̄i
b∈X̄j
P̄ (a, b)π̄(a)
P|[N ]|(i, i) := 1−
j 6=i
P|[N ]|(i, j).
Again P|[N ]| is a reversible chain on |[N ]| with stationary distribution
π̄|[N ]|(i) := π̄(X̄i).
Finally for every r = 1, 2, . . . , N we define a chain on X̄r by setting
PX̄r(a, b) := P̄ (a, b) + Ia(b)
z∈X̄ cr
P̄ (a, z)
where both a and b belong to X̄r. Now note that for every r = 2, 3, . . . , N PX̄r is
a birth and death chain on the state space {(1, r), (3, r), . . . , (r, r)} for r odd and
22 FEDERICO BASSETTI AND FABRIZIO LEISEN
{(0, r), (2, r), . . . , (r, r)} for r even. Let
qr(s) :=
(r − s)/2
and, for r even,
qr(0) := 2
Now observe that PX̄r has stationary distribution
πr(s) ∝ qr(s)
with s = 0, 2, . . . , r if r is even and s = 1, 3, . . . , r if r is odd. First of all let r 6= 1, by
Lemma B.3 and Lemma B.1, it is easy to check that (PX̄r , πr) meets the condition
of Lemma A.1 with
B = 2, n = [(r + 2)/2], A = C3p1[(r + 2)/2]N
([x] being the integer part of x) and then
1− λ1(PX̄r ) ≥
C3p1[(r + 2)/2]
2N [(r + 2)/2]3
Finally, Lemma B.3 with (A.1) yields
λ|X̄r|−1(PX̄r ) ≥ 1− p1.
Hence, for every r 6= 1, we have proved that
(B.6) Gap(PX̄r) ≥ C3/2p1N
For r = 1
PX̄1 =
1− α1/2 α1/2
α2/2 1− α2/2
where
α1 :=
1, exp{
α2 := p1 min
1, exp{Kβ
Gap(PX̄1) ≥ 1− |
2− α1 − α2
| = α1 + α2
where the last equality follows from the fact that α1
and α2
. Hence, for
sufficiently large N , it’s easy to see that
(B.7) Gap(PX̄1) ≥ C4p1N
with C4 = C4(β,K). At this stage (B.6) with (B.7) gives
(B.8) Gap(PX̄r) ≥ C5p1N
for all r ∈ |[N ]|. As for the gap of P|[N ]|, first of all note that P|[N ]| is a birth and
death chain on |[N ]|. From Lemma B.3
P|[N ]|(i, i+1) :=
2π̄(X̄i)
a∈X̄i
b∈X̄i+1
P̄ (a, b)π̄(a) ≥ p1C3
2π̄(X̄i)
a∈X̄i
b∈X̄i+1
π̄(a) ≥ p1C3
and analogously,
P|[N ]|(i, i− 1) ≥
Now, for r 6= 1
π̄|[N ]|(r) = q|[N ]|(r)/(
q|[N ]|(i))
METROPOLIS ALGORITHM AND EQUIENERGY SAMPLING 23
while
π̄|[N ]|(1) = (q|[N ]|(1) + q|[N ]|(0))/(
q|[N ]|(i)).
So, using the unimodality of q|[N ]|, we can apply Lemma B.3 with
which gives
λ1(P|[N ]|) ≤ 1−
≤ 1− p1C3
Using another time Lemma B.3, by (A.1), we get
λN (P|[N ]|) ≥ 1− p1.
Combining this two bounds we have
(B.9) Gap(P|[N ]|) ≥
and so from (3.4)
Gap(P̄ ) ≥
C being a suitable constant that depends by β,K,C3, C4, C5. �
Proof of Proposition 6.4. Combine Lemma 6.3 with 6.5. �
Acknowledgments
We should like to thank Persi Diaconis for useful discussions and for having
encouraged us during this work, Antonietta Mira for suggesting some interesting
references and Claudio Giberti for helping to improve an earlier version of the paper.
References
[1] J. Besag and P. J. Green. Spatial statistics and Bayesian computation. J. Roy. Statist. Soc.
Ser. B, 55(1):25–37, 1993.
[2] M. Blume, V. J. Emery, and R. B. Griffiths. Ising model for the λ transition and phase
separation in he3-he4 mixtures. Phys. Rev. A, 4:1071–1077, 1971.
[3] P. Bremaud. Markov Chains. Springer-Verlag, New York, 1998.
[4] S. Caracciolo, A. Pelissetto, and A. D. Sokal. Nonlocal Monte Carlo algorithm for self-avoiding
walks with fixed endpoints. J. Statist. Phys., 60(1-2):1–53, 1990.
[5] S. Caracciolo, A. Pelissetto, and A. D. Sokal. Dynamic critical exponent of the BFACF
algorithm for self-avoiding walks. J. Statist. Phys., 63(5-6):857–865, 1991.
[6] P. Diaconis and L. Saloff-Coste. What do we know about the Metropolis algorithm? J.
Comput. System Sci., 57(1):20–36, 1998. 27th Annual ACM Symposium on the Theory of
Computing (STOC’95) (Las Vegas, NV).
[7] P. Diaconis and D. Stroock. Geometric bounds for eigenvalues of Markov chains. Ann. Appl.
Probab., 1(1):36–61, 1991.
[8] M. Dyer, L. A. Goldberg, M. Jerrum, and R. Martin. Markov chain comparison. Probab.
Surv., 3:89–111 (electronic), 2006.
[9] R. G. Edwards and A. D. Sokal. Generalization of the Fortuin-Kasteleyn-Swendsen-Wang
representation and Monte Carlo algorithm. Phys. Rev. D (3), 38(6):2009–2012, 1988.
[10] R. S. Ellis, P. T. Otto, and H. Touchette. Analysis of phase transitions in the mean-field
Blume-Emery-Griffiths model. Ann. Appl. Probab., 15(3):2203–2254, 2005.
[11] R.S. Ellis. The Theory of Large Deviation and Applications to Statistical Mechanics.
2006. Lectures for the international seminar on Extreme Events in Complex Dynamics.
http://www.math.umass.edu/ ˜ rsellis/pdf-files/Dresden-lectures.pdf Max-Planck-Institut.
[12] A. E. Gelfand, S. K. Sahu, and B. P. Carlin. Efficient parameterisations for normal linear
mixed models. Biometrika, 82(3):479–488, 1995.
[13] A. Gelman, G. O. Roberts, and W. R. Gilks. Efficient Metropolis jumping rules. In Bayesian
statistics, 5 (Alicante, 1994), Oxford Sci. Publ., pages 599–607. Oxford Univ. Press, New
York, 1996.
http://www.math.umass.edu/
24 FEDERICO BASSETTI AND FABRIZIO LEISEN
[14] C. J. Gey er and E. A. Thompson. Annealing markov chain monte carlo with applications to
ancestral inference. J. Amer. Statist. Assoc., 90:909–920, 1995.
[15] Y. Guan, R. Fleißner, P. Joyce, and S. M. Krone. Markov chain Monte Carlo in small worlds.
Stat. Comput., 16(2):193–202, 2006.
[16] Y. Guan and S. M. Krone. Small-world mcmc and convergence to multi-modal distributions:
From slow mixing to fast mixing. Ann. Appl. Probab., 17(1):284–304, 2007.
[17] J. M. Hammersley and D. C. Handscomb. Monte Carlo methods. Methuen & Co. Ltd., Lon-
don, 1965.
[18] W.K. Hastings. Monte carlo sampling methods using markov chains and their application.
Biometrika, 57:97–109, 1970.
[19] K. Hukushima and K. Nemoto. Exchange monte carlo method and application to spin glass
simulations. J.Phys.Soc.Jpn., 65:1604–1608, 1996.
[20] M. Jerrum, J. Son, P. Tetali, and E. Vigoda. Elementary bounds on Poincaré and log-Sobolev
constants for decomposable Markov chains. Ann. Appl. Probab., 14(4):1741–1765, 2004.
[21] S. A. Kivelson, V. J. Emery, and H. Q. Lin. Doped antiferromagnets in the weak-hopping
limit. Phys. Rev. B, 42:6523–6530, 1990.
[22] S. C. Kou, Qing Zhou, and Wing Hung Wong. Equi-energy sampler with applications in
statistical inference and statistical mechanics. Ann. Statist., 34(4):1581–1619, 2006.
[23] J. Lajzerowicz and J. Sivardiére. Spin 1 lattice gas model. ii. condensation and phase sepa-
ration in a binary fluid. Phys. Rev. A, 11:2090–2100, 1975.
[24] J. Lajzerowicz and J. Sivardiére. Spin 1 lattice gas model. iii. tricritical points in binary and
ternary fluids. Phys. Rev. A, 11:2101–2110, 1975.
[25] J. S. Liu. Monte Carlo strategies in scientific computing. Springer Series in Statistics.
Springer-Verlag, New York, 2001.
[26] N. Madras and M. Piccioni. Importance sampling for families of distributions. Ann. Appl.
Probab., 9(4):1202–1225, 1999.
[27] N. Madras and D. Randall. Markov chain decomposition for convergence rate analysis. Ann.
Appl. Probab., 12(2):581–606, 2002.
[28] E. Marinari and G. Parisi. Simulated tempering: a new monte carlo scheme. Europhy s. Lett,
19:451–458, 1992.
[29] N. Metrpolis, A. Rosenbluth, M. Rosenbluth, A. Teller, and E. Teller. Equation of state
caluclations by fast computing machines. J. Chem. Phys., 21:1087–1092, 1953.
[30] A. Mira and L. Tierney. Efficiency and convergence properties of slice samplers. Scand. J.
Statist., 29(1):1–12, 2002.
[31] K.E. Newman and J.D. Dow. Zinc blende diamond order disorder transition in metastable
crystalline (gaas)1-xge2x alloys. Phys. Rev. B, 27:7495–7508, 1983.
[32] P. H. Peskun. Optimum Monte-Carlo sampling using Markov chains. Biometrika, 60:607–612,
1973.
[33] P. H. Peskun. Guidelines for choosing the transition matrix in Monte Carlo methods using
Markov chains. J. Comput. Phys., 40(2):327–344, 1981.
[34] C. P. Robert and G. Casella. Monte Carlo statistical methods. Springer Texts in Statistics.
Springer-Verlag, New York, second edition, 2004.
[35] R. Y. Rubinstein. Simulation and the Monte Carlo method. John Wiley & Sons Inc., New
York, 1981. Wiley Series in Probability and Mathematical Statistics.
[36] M. Schick and Wei-Heng Shih. Spin 1 model of a microemulsion. Phys. Rev. B, 34:1797–1801,
1986.
[37] A. Sinclair. Algorithms for random generation and counting. Progress in Theoretical Com-
puter Science. Birkhäuser Boston Inc., Boston, MA, 1993. A Markov chain approach.
[38] J. Sivardiére and J. Lajzerowicz. Spin 1 lattice gas model. i. condensation and solidification
of a simple fluid. Phys. Rev. A, 11:2090–2100, 1975.
[39] A. Sokal. Monte Carlo methods in statistical mechanics: foundations and new algorithms. In
Functional integration (Cargèse, 1996), volume 361 of NATO Adv. Sci. Inst. Ser. B Phys.,
pages 131–192. Plenum, New York, 1997.
[40] R.H. Swendsen and J. S. Wang. Non–universal critical dynamics in monte carlo simulations.
Phy. Rev. Lett., pages 86–88, 1987.
Università degli Studi di Pavia, Dipartimento di Matematica, via Ferrata 1, 27100
Pavia, Italy
Università dell’Insubria, Dipartimento di economia, via monte generoso 71, 21100
Varese , Italy
E-mail address: federico.bassetti@unipv.it
E-mail address: leisen.fabrizio@unimore.it
	1. Introduction.
	2. A general strategy
	3. Preliminaries
	3.1. Spectral gap, variance and speed of convergence
	3.2. Cheeger's inequality
	3.3. Chain decomposition theorem
	4. Warming up example
	5. The mean field Ising model
	6. The mean–field Blume-Emery-Griffiths model
	Appendix A. The Spectral Gap of a Birth and Death Chain
	Appendix B. Proofs
	Acknowledgments
	References

  In this paper we study the Metropolis algorithm in connection with two
mean--field spin systems, the so called mean--field Ising model and the
Blume--Emery--Griffiths model. In both this examples the naive choice of
proposal chain gives rise, for some parameters, to a slowly mixing Metropolis
chain, that is a chain whose spectral gap decreases exponentially fast (in the
dimension $N$ of the problem). Here we show how a slight variant in the
proposal chain can avoid this problem, keeping the mean computational cost
similar to the cost of the usual Metropolis. More precisely we prove that, with
a suitable variant in the proposal, the Metropolis chain has a spectral gap
which decreases polynomially in 1/N. Using some symmetry structure of the
energy, the method rests on allowing appropriate jumps within the energy level
of the starting state.

Introduction.
The Metropolis algorithm, introduced in [29] and later generalized in [18], is
currently (together with other Monte Carlo Markov Chain methods) one of the
most used simulation techniques both in statistics and in physics. See, among
others, [33, 32, 39, 17, 35, 34, 25, 6].
In a finite setting the Metropolis algorithm can be described as follows. Suppose
that, given a probability π(x) on a finite set X , want to approximate
(1.1) µ =
f(x)π(x),
for f : X → R. As a first step, take a reversible Markov chain K(x, y) (the proposal
chain) on X and change its output in order to have a new chain with stationary
distribution π. This can be achieved by constructing a new (π–reversible) chain
(1.2) M(x, y) =
K(x, y)A(x, y) x 6= y
K(x, x) +
z 6=xK(x, z)(1−A(x, z)) x = y
where A(x, y) := min(
π(y)K(y,x)
π(x)K(x,y)
, 1). Then, the metropolis estimate of µ is given by
(1.3) µ̂n =
f(Yi),
where Y0 is generated from some initial distribution π0 and Y1, . . . , Yn fromM(x, y).
It is clear that, from a computational point of view, the speed of convergence to
the stationary distribution and the (asymptotic) variance of the estimate are two
very important features of the Markov chain M .
It is well-known that in some situation a Markov chain can converge very slowly
to its stationary distribution and, moreover, that the asymptotic variance of the es-
timate (1.3) can be much bigger than the variance of f , i.e. V arπ(f) :=
x(f(x)−
Key words and phrases. asymptotic variance, Chain decomposition theorem, fast/slowly mix-
ing chain, mean-field Ising model, Metropolis, spectral gap analysis.
http://arxiv.org/abs/0704.0906v2
2 FEDERICO BASSETTI AND FABRIZIO LEISEN
µ)2π(x), which is equal to the asymptotic variance of the crude Montecarlo estima-
tor. In these cases (1.3) turns out to be a very inefficient estimate of µ.
For the Metropolis chain a classical situation in which the convergence is slow
(and the variance big) is when the target distribution π has many peaks and K is
somehow too “local”.
This is well known in statistical physics, where, typically, a distribution of a
system with energy function h and in thermal equilibrium at temperature T is
described by the Gibbs distribution
πh,T (x) = exp{−h(x)/T }Z−1T
with ZT =
x exp{−h(x)/T }. In point of fact, the Metropolis algorithm has been
proposed in [29] to compute average with respect to such distributions. Indeed,
if h is nice, the Metropolis algorithm is very efficient, but it can perform very
poorly if the energy has many local minima separated by high barriers that cannot
be crossed by the proposal moves K. This problem can be bypassed, for specific
energy, designing appropriate moves that have higher chance to cut across the
energy barrier (see, e.g, [4, 5]), or constructing clever alternative approaches to the
problem, for instance using a reparametrization of the problem (see, e.g., [12, 13])
or using auxiliary variables (see, e.g., [40, 9, 1, 30]). A different kind of solution has
been proposed in [14] and in [28] by introducing the so called simulated tempering,
which essentially means that T is changed (stochastically or not) to flatten h. A
remarkable variant of these methods is the parallel tempering, see, for instance, [19].
More recently new algorithms based on the so called equi–energy levels sampling
have been proposed (see [26] and [22]). In particular, the algorithm proposed in [22]
relies on the so–called equi-energy jump, which enables the chain to reach regions
of the sample space with energy close to the one of the starting state, but that may
be separated by steep energy barriers. In point of fact, even if, according to some
simulations, the method seems to be efficient nothing has been formally proved.
Finally, let us mention a recent algorithm, called small world Markov chains (see
[15, 16]), that combine a local chain with long jumps. In these papers, it has
been shown that a simple modification of the proposal mechanism results in faster
convergence of the chain. That mechanism, which is based on an idea from the
field of small-world networks, amounts to adding occasional wild proposals to any
local proposal scheme.
In the present paper we study two simple examples: the so called mean field Ising
model and the mean field Blume–Emery–Griffiths model. As for the former, it is
well-known that the usual choice of K gives rise, for low temperature, to a slowly
mixing Metropolis chain (see, e.g., [26]). Here we show that a slight variant in the
proposal chain can completely solve this problem, keeping the mean computational
cost similar to the cost of the usual Metropolis. The idea again rests on allowing
appropriate jumps in the same energy level of the starting state. As for the Blume–
Emery–Griffyths mean–field model, we first show that there is a critical region of
the parameters space for which the naive Metropolis chain is slowly mixing. Then
we show how one can modify the proposal chain in order to obtain a better mixing
for the Metropolis chain. The present paper should be intended as a further step in
the direction of a better mathematical understanding of both small world Markov
chains and equi-energy sampling.
The rest of the paper is organized as follows. In Section 2 some general consid-
erations are given. In Section 3 some basic tools concerning Markov chain, which
will be used in the paper, are reviewed. Section 4 contains a warming up example.
In Section 5 the mean field Ising model is treated, while Section 6 deals with the
more complex case of the mean field Blume-Emery-Griffiths model. All the proofs
are deferred to the Appendix.
METROPOLIS ALGORITHM AND EQUIENERGY SAMPLING 3
2. A general strategy
In an abstract setting, what we shall do in the next examples can be summarized
as follows. Let G be a group acting on X for which
(2.1) π(x) = π(g(x)) ∀ x ∈ X , ∀ g ∈ G.
For every x in X let Ox := {y = g(x) : g ∈ G} be the orbit of x (of course if y
belongs to Ox then Ox = Oy).
Assume now that we have a reversible Markov chain KE(x, y) (the proposal) on
X and suppose that the Metropolis chain ME with proposal KE is slowly mixing
(see next section for more details). To speed up the mixing one can try to exploit
(2.1) by taking a proposal of the following form:
(2.2) Kǫ(x, y) = ǫKE(x, y) + (1 − ǫ)KG(x, y)
where
KG(x, y) =
qx(z)Iz(y),
0 < qx(z) < 1 and
qx(z) = 1.
In point of fact, usually KE is “local”; for instance frequently
KE(x, y) = 0
whenever y 6= x belongs to Ox, hence with KG we are adding “long” jumps to the
chain. Moreover, note that if KE is such that KE(x, g(x)) = KE(g(x), x), for every
x in X and g in G, then the Metropolis always accepts the move x→ g(x) and
M(x, g(x)) = ǫKE(x, g(x)) + (1− ǫ)qx(g(x)).
In particular this holds when KE is symmetric.
The heuristics under (2.2) is to combining small world Markov chains and equi-
energy sampling.
Before presenting some examples in which one can actually improve the perfor-
mances of the Metropolis chain using this idea, we collect in the next section some
useful facts concerning Markov chains.
3. Preliminaries
Let P (x, y) be a reversible and ergodic Markov chain on the finite set X with
(unique) stationary distribution p(x). Thus, p(x)P (x, y) = p(y)P (y, x). Let L2(p) =
{f : X → R} with < f, g >p= Ep(fg) =
x f(x)g(x)p(x). Reversibility is equiva-
lent to P : L2 → L2 being self–adjoint. Here Pf(x) =
y f(y)P (x, y). The spec-
tral theorem implies that P has real eigenvalues 1 = λ0(P ) > λ1(P ) ≥ λ2(P ) ≥
· · · ≥ λ|X |−1(P ) > −1 with orthonormal basis of eigen–functions ψi : X → R
(Pψi(x) = λiψi(x), < ψi, ψj >p= δij).
3.1. Spectral gap, variance and speed of convergence. A very important
quantity related to the eigenvalues is the spectral gap, defined by
Gap(P ) = 1−max{λ1, |λ|X |−1|}.
It turns out that the spectral gap is a good index to measure the mixing of a
chain. To better understand this point, assume that f belongs to L2(p) and write
f(x) =
i≥0 aiψi(x) (with ai =< f, ψi >p). Now let Y0 be chosen form some
distribution p0 and Y1, . . . , Yn be a realization of the P (x, y) chain, then
µ̂n =
f(Yi)
4 FEDERICO BASSETTI AND FABRIZIO LEISEN
has asymptotic variance given by
AV ar(f, p, P ) := lim
n · V ar(µ̂n) =
|ak|2
1 + λk
1− λk
See, for instance, Theorem 6.5 in Chapter 6 of [3]. From the last expression, the
classical inequality
(3.1) AV ar(f, p, P ) ≤ 2
1− λ1
V arp(f),
follows easily. The last inequality is the usual way of relating spectral gap to
asymptotic variance and, hence, to the efficency of a chain.
The spectral gap is very important also to give bounds on the speed of con-
vergence to the stationary distribution. For example, if ‖ · ‖TV denotes the total
variation norm, one has
‖δxP k − p‖2TV =
|P k(x,A) − p(x)|
≤ 1− p(x)
4p(x)
(max{λ1, |λ|X |−1|})2k
See, e.g., Proposition 3 in [7]. Another classical bound is
‖p0P k/p− 1‖2,p ≤ Gap(P k)‖p0/p− 1‖2,p
valid for every probability p0. See, for instance, [39].
Roughly speaking one can say that a sequence of Markov chains defined on a
sequence of state space XN is slowly mixing (in the dimension of the problem N)
if the spectral gap decreases exponentially fast in N .
3.2. Cheeger’s inequality. As already recalled, problems of slowly mixing typ-
ically occur when π has two or more peaks and the chain K can only move in
a neighborhood of the starting peak. Usually this phenomenon is called bottle-
neck. A powerful tool to detect the presence of a bottleneck is the conductance and
the related Cheeger’s inequality. Recall that the conductance of a chain P with
stationary distribution p is defined by
h = h(p, P ) := inf
A :p(A)≤ 1
x∈A,y∈Ac
p(x)P (x, y),
and the well-known Cheeger’s inequality is
(3.2) 1− 2h ≤ λ1(P ) ≤ 1−
See, for instance, [3, 37, 7]. Note that, since P is reversible,
(3.3) h ≤ 1
p(x)P (x, y) =
p(y)P (y, x)
for every A such that p(A) ≤ 1/2.
3.3. Chain decomposition theorem. In this subsection we briefly describe a
useful technique to obtain bounds on the spectral gap: the so called chain decom-
position technique. Following [16] assume that A1, . . . , Am is a partition of X .
Moreover, for each i = 1, . . . ,m, define a new Markov chain on Ai by setting
PAi(x, y) := P (x, y) + Ix(y)
P (x, z)
 (x, y ∈ Ai).
PAi is a reversible chain on the state space Ai with respect to the probability
measure
pi(x) := p(x)/p(Ai).
METROPOLIS ALGORITHM AND EQUIENERGY SAMPLING 5
The movement of the original chain among the “pieces” A1, . . . , Am can be de-
scribed by a Markov chain with state space {1, . . . ,m} and transition probabilities
PH(i, j) :=
2p(Ai)
x∈Ai,y∈Aj
P (x, y)p(x)
for i 6= j and
PH(i, i) := 1−
j 6=i
PH(i, j),
which is reversible with stationary distribution
p̄(i) := p(Ai).
A variant of a result of Caracciolo, Pelisetto and Sokal (published in [27]), states
(3.4) Gap(P ) ≥ 1
Gap(PH)
i=1,...,m
Gap(PAi)
holds true, see Theorem 2.2 in [16]. Other results about chain decompositions can
be found, for instance, in [20].
In the next very simple example we shall show how this technique can be used,
starting from a slowly mixing chain, to suggest how to modify the proposal chain
in order to obtain a fast mixing chain.
4. Warming up example
Set X = {−N,−N + 1, . . . , 0, 1, . . . , N} and define a probability measure on X
π(x) =
(θ − 1)θ|x|
2θN+1 + 1− θ
θ being a given parameter bigger than 1. Here we can consider G = {+1,−1} (with
group operation given by the usual product) acting on X by g(x) = gx, hence
Ox = {x,−x}.
Now let KE be a chain defined by
KE(x, x + 1) = 1/2 x 6= N
KE(x, x − 1) = 1/2 x 6= −N
KE(N,N) = KE(−N,−N) = 1/2
KE(x, y) = 0 otherwise
and denote by ME the Metropolis chain with stationary distribution π derived
by KE. It is clear that in this case KE(x, y) = 0 whenever y belongs to Ox.
In this example it is very easy to bound the conductance on ME , indeed, taking
A = {−N, . . . ,−1}, by (3.3), it follows that
h(π,ME) ≤
1− π(0)
Hence,
h(π,ME) ≤ Cθ−N ,
and then (3.2) yields
1− λ1 ≤ 2Cθ−N .
This means that, if f is such that a1 6= 0 and θ > 1, then the asymptotic variance
of f blows up exponentially fast, indeed
AV ar(f, π,ME) ≥ 2Celog(θ)N .
6 FEDERICO BASSETTI AND FABRIZIO LEISEN
Now, instead of KE consider
Kǫ(x, y) = (1− ǫ)KE(x, y) + ǫI{−x}(y)
and let M (ǫ) be the Metropolis chain derived by Kǫ. Decompose X as follows
X = A1 ∪A2 · · · ∪ AN
with A1 = {−1, 0, 1} and Ai = {x ∈ X : |x| = i}, for i > 1. Moreover let
π̄(i) = π(Ai) =
(2θ + 1)/Z for i = 1
2θi/Z for i > 1
where
2θN+1 + 1− θ
(θ − 1)
and set
H (i, j) =
2π(Ai)
l∈Ai,m∈Aj
M (ǫ)(l,m)π(l), M
H (i, i) = 1−
j 6=i
H (i, j).
For i 6= 1, N , one has
H (i, i+ 1) =
2π(Ai)
[M (ǫ)(i, i+ 1)π(i) +M (ǫ)(−i,−i− 1)π(−i)]
and, since π(i) = π(−i) and π(i + 1) ≥ π(i)
H (i, i+ 1) =
In the same way it is easy to see that
H (i, i− 1) =
, i 6= 1, N
H (i, i) = 1−
(1 + θ−1) i 6= 1, N
H (N,N − 1) =
H (N,N) = 1−
H (1, 2) =
4(1 + 1/(2θ))
H (1, 1) = 1−
4(1 + 1/(2θ))
Moreover, for every i 6= 1, M (ǫ)Ai in matrix form is given by
1− ǫ ǫ
ǫ 1− ǫ
and hence
Gap(M
) = 1− |1− 2ǫ|.
While M
is given by
(2θ − 1)(1− ǫ)/(2θ) (1− ǫ)/(2θ) ǫ
(1− ǫ)/2 ǫ (1− ǫ)/2
ǫ (1− ǫ)/(2θ) (2θ − 1)(1− ǫ)/(2θ)
and hence
Gap(M
) = k(θ, ǫ) > −1.
Moreover, since
i6=1,N
H (i, i± 1)),M
H (1, 2),M
H (N,N − 1)
≥ min
(1− ǫ)/(4θ), 1− ǫ
4(1 + 1/(2θ))
=: m(ǫ, θ) > 0
METROPOLIS ALGORITHM AND EQUIENERGY SAMPLING 7
and π̄(i) ≤ 3π̄(j) for every i < j, Lemma A.1 in the appendix yields that
1− λ1(M (ǫ)H ) ≥
m(ǫ, θ)
In the same way, since M
H (i, i+1)+M
H (i, i− 1) ≤ (1− ǫ)M(θ)/4, with M(θ) =
max(1 + θ−1, 2θ/(2θ + 1)) ≤ 2, inequality (A.1) in the Appendix yields that
λN−1(M
H ) ≥ 1−
≥ 1 + ǫ
Hence
Gap(M
H ) ≥
m(ǫ, θ)
and (3.4) yield
Gap(M (ǫ)) ≥
h(θ, ǫ)
for a suitable h. This shows that M (ǫ) is fast mixing for every ǫ > 0 and for every
θ > 1 while ME is slowly mixing for every θ > 1.
5. The mean field Ising model
Let X = {−1, 1}N , N being an even integer. For every β > 0 let π = πβ,N be a
probability on X defined by
π(x) = πβ,N (x) := exp
S2N (x)
Z−1N (β) (x ∈ X )
where
ZN (β) = ZN :=
S2N (x)
is the normalization constant (“partition function”) and
SN (x) :=
xi x = (x1, . . . , xN ).
This is the so called mean field Ising model, or Curie-Weiss model, in which every
particle i, with spin xi, interacts equally with every other particle. It is probably
the most simple but also the most studied example of spin system on a complete
graph. The usual Metropolis algorithm uses as proposal chain
KE(x, y) =
I{x(j)}(y)
where x(j) denotes the vector (x1, . . . ,−xj , . . . , xN ). It has been proved in [26] that,
whenever β > 1,
1− λ1 ≤ Ce−D
where λ1 is the first eigenvalues smaller than 1 of the Metropolis chain ME derived
KE . This yields that the variance of an estimator obtained from this Metropolis
algorithm can blow up exponentially fast in N .
The aim of this section is to show how one can construct a different Metropolis
chain avoiding this problem. In the notation of Section 2, we consider
G = SN × {+1,−1}
(SN being the symmetric group of order N) and we define the action of G on
X = {−1, 1}N by
g(x) = (e · xσ(1), . . . , e · xσ(N)) g = (σ, e).
8 FEDERICO BASSETTI AND FABRIZIO LEISEN
In order to introduce a new proposal, it is useful to write X as the union of its
“energy sets”, that is
X = X0 ∪ X2 ∪ X4 ∪ · · · ∪ XN
where
Xi := {x ∈ X : |SN (x)| = i} (i = 0, 2, . . . , N).
Note that energy takes only even values and that Ox = X|SN (x)|. Moreover, for
i 6= 0, set
X+i := {x ∈ X : SN (x) = i} and X
i := {x ∈ X : SN(x) = −i}.
The new proposal chain will be
K(x, y) = p1KE(x, y) + (1 − p1)K0(x, y) if x ∈ X0
K(x, y) = p1KE(x, y) + p2I{−x}(y) + (1− p1 − p2)Ki(x, y)
if x ∈ Xi, i 6= 0
(5.1)
where p1, p2 belong to (0, 1), p1 + p2 < 1, and
Ki(x, y) = IX+
{x}K+i (x, y) + IX−
{x}K−i (x, y) (i 6= 0).
We shall assume that K±i (K0, respectively) are irreducible, symmetric and aperi-
odic chains on X±i ( X0, respectively).
As a leading example we shall take
K0(x, y) =
) y ∈ X0
K±i (x, y) =
(N−i)/2
) y ∈ X±i ,
(5.2)
that is: a realization of a chain K±i (K0, respectively) is simply a sequence of
independent uniform random sampling from X±i (X0, respectively).
Remark 1. Note that (5.2) is the (n, k)-Bose-Einstein distribution with n = (N +
i)/2 and k = (N − i)/2 + 1 and recall that there is a very easy way to directly
generate Bose-Einstein configurations. One may place n balls sequentially into k
boxes, each time choosing a box with probability proportional to its current content
plus one. Starting from the empty configuration this results in a Bose-Einstein
distribution for every stage.
Now let M be the Metropolis chain defined by the transition kernel (1.2) with
K as in (5.1), i.e. for every x in X±i (i 6= 0)
M(x, y) =
if y = x(j), j = 1...N
p2 if y = −x
(1− p1 − p2)K±i (x, y) if y ∈ X
i , y 6= x
z 6=xM(x, z) if y = x
while for x in X0
M(x, y) =
if y = x(j), j = 1...N
(1− p1)K0(x, y) if y ∈ X0, y 6= x
z 6=xM(x, z) if y = x.
METROPOLIS ALGORITHM AND EQUIENERGY SAMPLING 9
By construction M is an aperiodic, irreducible and reversible chain with stationary
distribution π. Then, when (5.2) holds true,
M(x, y) =
if y = x(j), j = 1...N
p2 if y = −x
(1− p1 − p2) 1( N(N−i)/2)
if y ∈ X±i , y 6= x
z 6=xM(x, z) if y = x
for x in X±i (i 6= 0), while if x belongs to X0
M(x, y) =
if y = x(j), j = 1...N
(1− p1) 1( NN/2)
if y ∈ X0, y 6= x
z 6=xM(x, z) if y = x.
In order to bound the spectral gap ofM we shall use the decomposition theorem
described in Subsection 3.3. To this end, for every i = 0, 2, . . . , N and every j 6= i
P̄ (i, j) :=
2π(Xi)
M(x, y)π(x)
P̄ (i, i) := 1−
j 6=i
P̄ (i, j).
As already noted, P̄ is a reversible chain on {0, 2, . . . , N} with stationary distribu-
π̄(i) := π(Xi).
Moreover define for every i = 0, 2, . . . , N a chain on Xi setting
PXi(x, y) :=M(x, y) + Ix(y)
z∈X c
M(x, z)
where both x and y belong to Xi. In the same way, define chains on X+i and X
for i = 2, . . . , N setting
(x, y) := PXi(x, y) (y 6= x, x, y ∈ X±i )
(x, x) := 1−
y 6=x
PXi(x, y).
These chains are reversible on Xi (X±i , respectively) and have as stationary distri-
butions
πXi(x) :=
π(Xi)
and πX±
(x) :=
πXi(x)
πXi(X±i )
|X±i |
10 FEDERICO BASSETTI AND FABRIZIO LEISEN
respectively. Finally, for every i = 2, 4, . . . , N , define a chain on {+,−} setting
Pi(+,−) :=
2πXi(X+i )
PXi(x, y)πXi(x)
Pi(−,+) :=
2πXi(X−i )
PXi(x, y)πXi(x).
Now the lower bound (3.4), applied two times yields
Gap(M) ≥ 1
Gap(P̄ ) min
i=0,2,...,N
{Gap(PXi)}
Gap(P̄ )min
Gap(PX0),
i=2,...,N
Gap(Pi)min{Gap(PX+
), Gap(PX−
(5.3)
Hence, to get a lower bound on Gap(M) it is enough to obtain bounds on the gaps
of the chains P̄ , PX0 , Pi, PX±
The most important of these bounds is given by the following
Proposition 5.1. P̄ is a birth and death chain on {0, 2, . . . , N}, more precisely
(5.4)
P̄ (0, 2) = p1
P̄ (i, i+ 2) = p1
i 6= N, 0
P̄ (i, i− 2) = p1
exp{2β(1− i)/N} i 6= 0.
Moreover
λ1(P̄ ) ≤ 1−
(N/2 + 1)3
λN/2(P̄ ) ≥ 1− p1.
The proof of the previous proposition is based on a bound for a birth and death
chain, given in the Appendix, which can be of its own interest.
As for the others chains, we have the following
Lemma 5.2. For every i = 2, 4, . . . , N
Gap(PX±
) ≥ (1− p1 − p2)Gap(K±i )
Gap(Pi) = p2,
moreover
Gap(PX0) ≥ (1− p1)Gap(K0).
In this way, using (5.3), we can prove the main result of this section.
Proposition 5.3. Let M be the Metropolis chain derived by the chain K defined
as in (5.1) then
Gap(M) ≥ p1p2
(N/2 + 1)3
[ (1− p1)
Gap(K0),
(1− p1 − p2)
min{Gap(K+i ), Gap(K
i, )}
If K±i and K0 are defined as in (5.2) then
Gap(M) ≥
(N/2 + 1)3
[ (1 − p1 − p2)
(1− p1)
for every β > 0 and N ≥ N0.
METROPOLIS ALGORITHM AND EQUIENERGY SAMPLING 11
Proposition 5.3 shows that the gap is polynomial in 1/N independently of β.
Hence, even when β > 1, the variance of the metropolis estimate obtained with this
proposal can not grow up faster than a polynomial in N .
Note that if in Proposition 5.3 we choose
(5.5) p1 = 1− a/(2N), p2 = a/N
we get
Gap(M) ≥ C
Hence, even with this choice, the Metropolis algorithm is still fast mixing for every
β. It is worth noticing that the mean computational cost of this Metropolis does
not change with respect to the Metropolis which uses the proposal KE. Indeed,
in the case of the usual Metropolis, the computational cost needed to go from Xn
to Xn+1 is O(N), since it is essentially due to a sample of one number among N
numbers (we need to decide which coordinate to flip). In the case of the ”modified”
proposal, things are slight more complex. In this case, at the beginning, we have
an extra “toss”. If with this fist toss we decide to flip at random a coordinate the
cost is still O(N) but if we need to sample from K±i the cost is O(N
2) (in this last
case we need to pick a sample from a Bose-Einstein distribution). Hence, although
our algorithm is ”sometime” more expensive, if we take p1 and p2 as in (5.5), we
get that the mean cost of our algorithm is still O(N).
6. The mean–field Blume-Emery-Griffiths model
The Blume-Emery-Griffiths (BEG) model (see [2]) is an important lattice–spin
model in statistical mechanics, it has been studied extensively as a model of many
diverse systems, including He3 − He4 mixtures as well solid–liquid–gas systems,
microemulsions, semiconductor alloys and electronic conduction models. See, for
instance, [2, 38, 23, 24, 31, 36, 21]. We will focus our attention on a simplified
mean–field version of the BEG model. For a mathematical treatment of this mean–
field model see [10]. In what follows let X := {−1, 0, 1}N , N being an even integer,
and for every β > 0 and K > 0 let πβ,K,N be the probability defined by
π(x) = πβ,K,N(x) = exp{−βRN(x) +
S2N(x)}Z−1N (β,K) (x ∈ X )
where
ZN(β,K) = ZN :=
−βRN (x) +
S2N (x)
is the normalization constant,
SN (x) :=
xi and RN (x) :=
x2i x = (x1, x2, ..., xN ).
A natural Metropolis algorithm can be derived by using the proposal chain
(6.1) KE(x, y) =
[I{x(+j)}(y) + I{x(−j)}(y)]
where x(±j) denotes the vector (x1, . . . , xj ± 1, . . . , xN ), with the convention that
2 = −1 and −2 = 1.
The next proposition shows that there exists a critical region of the parameters
space in which the Metropolis chain is slowly mixing. More precisely, using some
results of [10] it is quite straightforward to proove the following
12 FEDERICO BASSETTI AND FABRIZIO LEISEN
Proposition 6.1. Ler ME be the Metropolis chain (with stationary distribution π)
with proposal chain KE defined in (6.1).Then, there exists a non decreasing function
Γ : (0,+∞) → (0,+∞) with limx→0 Γ(x) = +∞ and limx→∞ Γ(x) = γc ≃ 1.082
such that for every couple of positive parametrs (β,K) with K > Γ(β)
Gap(ME) ≤ Ce−∆N
for suitable constants C = C(γ,K) > 0 and ∆ = ∆(γ,K) > 0.
As in the case of the mean–field Ising model, we intend to by pass the slowly
mixing problem of this Metropolis chain by choosing a different proposal. To un-
derstand which kind of proposal is reasonable, here we choose
G = SN × {+1,−1}
with G acting on X = {−1, 0, 1}N by
g(x) = (e · xσ(1), . . . , e · xσ(N)) g = (σ, e).
At this stage, decompose X as the union of its ”energy sets”, that is
X = X0,0 ∪ X1,1 ∪ X0,2 ∪ X1,3 ∪X3,3 ∪ ... ∪ X0,N ∪ X2,N ∪ ...XN,N
where
Xs,r := {x ∈ X : |SN | = s and RN (x) = r}
r = 0, 1, 2, ..., N and s = 1, 3, ..., r if r is odd and s = 0, 2, ..., N if r is even.
Moreover, for s = 1, 2, ..., N , set
X+s,r := {x ∈ X : SN = s and RN (x) = r}
X−s,r := {x ∈ X : SN = −s and RN (x) = r}.
Note again that Ox = Xs,r with s = SN (x) and r = RN (x). The new proposal
chain will be
K(x, y) = p1KE(x, y) + (1 − p1)K0,r(x, y) if x ∈ X0,r, r = 0, 2, ..., N
K(x, y) = p1KE(x, y) + p2I{−x}(y) + (1− p1 − p2)Ks,r(x, y)
if x ∈ Xs,r, s 6= 0
(6.2)
where p1, p2 belong to (0, 1), p1 + p2 < 1, and
Ks,r(x, y) = IX+s,r{x}K
s,r(x, y) + IX−s,r{x}K
s,r(x, y) (s 6= 0)
K0,r(x, y) =
) y ∈ X0,r
K±s,r(x, y) =
(r−s)/2
) y ∈ X±s,r.
(6.3)
Now let M be the Metropolis chain defined by the transition kernel (1.2) with K
as in (6.2), i.e. for every x in X±s,r (s 6= 0)
M(x, y) =
if y = x(±j), j = 1...N
p2 if y = −x
(1 − p1 − p2) 1(Nr )( r(r−s)/2)
if y ∈ X±s,r, y 6= x
z 6=xM(x, z) if y = x,
METROPOLIS ALGORITHM AND EQUIENERGY SAMPLING 13
while if x belongs to X0,r
M(x, y) =
if y = x(±j), j = 1...N
(1− p1) 1(Nr )( rr/2)
if y ∈ X0,r, y 6= x
z 6=xM(x, z) if y = x.
By construction M is an aperiodic, irreducible and reversible chain with stationary
distribution π.
Also in this case, to bound the spectral gap of M , we shall use the chain decom-
position tools. Let
DN = {(0, 0), (1, 1), (0, 2), (2, 2), (1, 3), (3, 3), (0, 4), (2, 4), (4, 4), ..., (0, N), (2, N), ..., (N,N)}
and, for every couple (s, r), (s̃, r̃) in DN , with (s, r) 6= (s̃, r̃), let
P̄ ((s, r), (s̃, r̃)) :=
2π(Xs,r)
x∈Xs,r
y∈Xs̃,r̃
M(x, y)π(x)
P̄ ((s, r), (s, r)) := 1−
(s̃,r̃) 6=(s,r)
P̄ ((s, r), (s̃, r̃)).
Once again, note that P̄ is a reversible chain on DN with stationary distribution
π̄(s, r) := π(Xs,r).
Moreover, for every (s, r) in DN , define a chain on Xs,r setting
PXs,r (x, y) :=M(x, y) + Ix(y)
z∈X cs,r
M(x, z)
where both x and y belong to Xs,r. In the same way, define chains on X+s,r and X−s,r
for (s, r) in DN , s 6= 0, setting
PX±s,r (x, y) := PXs,r(x, y) (y 6= x, x, y ∈ X
PX±s,r (x, x) := 1−
s,ry 6=x
PXs,r (x, y).
These chains are reversible on Xs,r (X±s,r , respectively) and have as stationary dis-
tributions
πXs,r(x) :=
π(Xs,r)
|Xs,r|
and πX±s,r (x) :=
πXs,r(x)
πXs,r (X±s,r)
|X±s,r|
respectively. Finally, for every (s, r) in DN , s 6= 0, define a chain on {+,−} setting
Ps,r(+,−) :=
2πXs,r(X+s,r)
PXs,r (x, y)πXs,r(x)
Ps,r(−,+) :=
2πXs,r(X−s,r)
PXs,r (x, y)πXs,r(x).
14 FEDERICO BASSETTI AND FABRIZIO LEISEN
At this stage, the lower bound (3.4), applied two times, yields
Gap(M) ≥ 1
Gap(P̄ ) min
(s,r)∈DN
Gap(PXs,r)
Gap(P̄ )min
r=0,2,...,N
Gap(PX0,r )
(s,r)∈DN ,s6=0
Gap(Ps,r)min{Gap(PX+s,r), Gap(PX−s,r)}
(6.4)
To derive from the last bound a more explicit bound we need some preliminary
work. The first result we need is exactly the analogous of Lemma 5.2.
Lemma 6.2. Fore every r = 1, . . . , N
Gap(PX0,r ) ≥ (1− p1)Gap(K0,r) = (1− p1),
moreover, for every (s, r) in DN with s 6= 0,
Gap(P±Xs,r) ≥ (1− p1 − p2)Gap(K
s,r) = (1− p1 − p2).
Finally, for every (s, r) in DN ,
Gap(Ps,r) = p2.
Hence, (6.4) can be rewritten as
(6.5) Gap(M) ≥ Gap(P̄ )p2
min{(1− p1)/2, (1− p1 − p2)/2}.
It remains to bound Gap(P̄ ). Unfortunately the the analogous of Proposition 5.1
is not so simple, hence we shall require an additional hypothesis. In what follows
q|[N ]|(r) :=
(r−2i)2
 if r is even
q|[N ]|(r) : =
(r−2i)2
 if r is odd
r = 0, 1, . . . , N and set
A = {β > 0,K > 0 : ∃N0 such that ∀N ≥ N0, q|[N ]| is unimodal}.
Lemma 6.3. For every (β,K) in A
Gap(P̄ ) ≥ Cp
for a suitable constant C = C(β,K).
Under the same assumptions of the previous Lemma we can state the main
results of this section.
Proposition 6.4. For every (β,K) in A
Gap(M) ≥ C̃p
for a suitable constant C̃ = C̃(β,K).
METROPOLIS ALGORITHM AND EQUIENERGY SAMPLING 15
0 5 10 15
β=2.3, K=0.5,0.6,0.7,0.8
0 5 10 15
β=2.3,2.4,2.5,2.6, K=0.5
0 5 10 15
β=2.3, K=1.2,1.3,1.4,1.5
0 5 10 15
β=2.3,2.4,2.5,2.6, K=1.2
Figure 1. The function q|[N ]| for N = 15 and few values of β and K.
We conjecture that Gap(P̄ ) is polynomial in N for every (β,K) such that β 6=
Γ(K) (where Γ is the function of Proposition 6.1), but we are not able to prove this
conjecture. In point of fact we conjecture that R+×R+ \{(β,K) : Γ(K) = β} ⊂ A.
We plotted q|[N ]| for different N , β andK, and these plotts seem, at least, to confirm
that R+ × R+ \ {(β,K) : |Γ(K)− β| ≤ ǫ} ⊂ A for a suitable small ǫ. In Figure 1
we show the graph of q|[N ]| for few different N , β and K.
Appendix A. The Spectral Gap of a Birth and Death Chain
We derive here some bounds on the eigenvalues of a birth and death chain that we
shall use later. These bounds are obtained using the so called geometric techniques,
see [7]. Let Pn be a birth and death chain on Ωn = {1, . . . , n}. Assume that Pn
is reversible with respect to a probability pn, that is pn(i)Pn(i, j) = pn(j)Pn(j, i).
Moreover let
1 > λ1 ≥ λ2 ≥ . . . λn−1 ≥ −1
the eigenvalues of Pn.
We can now prove the following variant of Proposition 6.3 in [6].
Lemma A.1. If there exist positive constants A, q, B and an integer k such that
Pn(i, i± 1) ≥ An−q (i 6= 1, n)
Pn(1, 2) ≥ An−q
Pn(n, n− 1) ≥ An−q
pn(i) ≤ Bpn(j) i ≤ j ≤ k
pn(j) ≤ Bpn(i) k ≤ i ≤ j
λ1 ≤ 1−
16 FEDERICO BASSETTI AND FABRIZIO LEISEN
Proof. We use the notation and the techniques of [7], see also [3] and [6]. Choose
the set of paths
Γ = {γij = (i, i+ 1, ..., j); i ≤ j; i, j ∈ Ωn}
and for e = (i, i+ 1) (i < n) let
ψ(e) =
pn(i, i+ 1)
γl,m∈Γ
γl,m∋e
|γl,m|
pn(l)pn(m)
pn(i)
where |γ| is the length of the path γ. Setting K := supe ψ(e) one has
λ1 ≤ 1−
(see Proposition 1’ in [7], or Exercise 6.4 page 248 in [3]). So, for our purposes, it
suffices to give an upper bound on K. Assume first that e = (i, i+ 1) with
i < k ≤ n, since |γl,m| ≤ n, it follows that
ψ(e) ≤ n
s≥i+1
pn(r)pn(s)
pn(i)
pn(r)
pn(i)
s≥i+1
pn(s)
pn(s)
≤ nq+2
All the other cases can be treated in the same way. Hence,
ψ(e) ≤ B
and then
λ1 ≤ 1−
As for the smaller eigenvalues, Gershgorin theorem yields that
λn−1 ≥ −1 + 2min
P (i, i).
See, for instance, Corollary 2.1 in the Appendix of [3]. Hence, if there exists a
positive constant D such that
Pn(i, i+ 1) + Pn(i, i− 1) ≤ D/2
for every i, then
(A.1) λn−1 ≥ 1−D.
METROPOLIS ALGORITHM AND EQUIENERGY SAMPLING 17
Appendix B. Proofs
To prove Proposition 5.1 we need first to show that π̄ is essentially unimodal.
Lemma B.1. Let
qN (i) =
i = 0, 2, 4, . . . , N.
For every β < 1 there exists an integer N0 such that for every N ≥ N0
qN (i) ≤ qN (j)
whenever j ≤ i. For every β ≥ 1 there exists an integer N0 such that for every
N ≥ N0
qN (i) ≤ qN (j)
whenever i ≤ j ≤ kN and
qN (i) ≥ qN (j)
whenever kN ≤ i ≤ j, kN being a suitable integer.
Proof. Let ∆N (i) be the ratio
∆N (i) =
qN (i+ 2)
qN (i)
i = 0, 2, 4, ..., N − 2,
so that
∆N (i) =
) exp
(1 + i)
N − i
N + 2 + i
(1 + i)
Setting ∆N (x) =
N+2+x
(1 + x)
, x in [0, N−2], it is enough to prove that
x 7→ ∆N (x) takes the value 1 at most once in [0, N − 2], for sufficiently large N .
To prove this last claim first note that
∆N (0) =
N + 2
1 + 2
= 1− 2
(1 − β) +
− β + 2
Hence, there exists N0 in N such that for N ≥ N0:
β ≥ 1 ⇒ ∆N (0) > 1
β < 1 ⇒ ∆N (0) < 1.
As for the first derivative note that
∆′N (x) =
−2(N + 1) + 2β(N + 2)− 2β
(x2 + 2x)
(N + x+ 2)2
(1 + x)
hence ∆′N (x) = 0 if and only if
−2(N + 1) + 2β(N + 2)− 2β
(x2 + 2x) = 0.
18 FEDERICO BASSETTI AND FABRIZIO LEISEN
Rearranging the last equation as
x2 − 4β
+ 2[(β − 1)N + 2β − 1] = 0
one sees that the roots are
x1,2 = 1±
2β − 1
β − 1
Hence, after setting
r := 1 +
2β − 1
β − 1
N2 and r := 1 +
one has
β < 1 ⇒ ∆′N (x) < 0 ∀x ∈ [0, N − 2]
β > 1 ⇒ ∆′N (x) > 0 for x ∈ [0, r)
∆′N (x) < 0 for x ∈ (r,N − 2]
β = 1 ⇒ ∆′N (x) < 0 for x ∈ [0, r)
∆′N (x) < 0 for x ∈ (r,N − 2]
and this concludes the proof. �
Proof of Proposition 5.1. By direct computations it is easy to prove (5.4). Hence
P (i, i± 2) ≥
4(N + 2)
P (i, i+ 2) + P (i, i− 2) ≤ p1
Now observe that
π̄(0) =
ZN (β)
qN (0)
π̄(i) =
ZN (β)
qN (i) i 6= 0.
Hence, by Lemma B.1, if β < 1
π̄(i) ≤ 2π̄(j)
whenever j ≤ i and N is large enough. While for β > 1
π̄(i) ≤ π̄(j)
whenever i ≤ j ≤ kN and
π̄(i) ≥ π̄(j)
whenever kN ≤ i ≤ j. The thesis follows now by Lemma A.1 and by (A.1). �
In order to prove Lemma 5.2 we recall that by Rayleigh’s theorem
(B.1) 1− λ1(P ) = inf
Ep(f, f)
V arp(f)
: f nonconstant
where
Ep(f, f) :=< (I − P )f, f >p=
(f(x) − f(y))2P (x, y)p(x),
METROPOLIS ALGORITHM AND EQUIENERGY SAMPLING 19
P being a reversible chain w.r.t. p, moreover
(B.2) 1− |λN−1| = inf
x,y(f(x) + f(y))
2P (x, y)p(x)
V arp(f)
: f nonconstant
(see, for instance, Theorem 2.3 in Chapter 6 of [3] and Section 2.1 of [8]). At this
stage set
Pǫ(x, y) := (1 − ǫ)P (x, y) + ǫIx(y).
Hence, (B.1) yields
1− λ1(Pǫ) = inf
f∈L2pf 6=const
x,y(f(x) − f(y))2Pǫ(x, y)p(x)
V arp(f)
= inf
f∈L2pf 6=const
(1− ǫ)
x 6=y(f(x)− f(y))2P (x, y)p(x)
V arp(f)
= (1− ǫ)(1− λ1(P )).
Arguing in the same way and using (B.2) we get
1− |λ|X |−1(Pǫ)| ≥ (1 − ǫ)(1− |λ|X |−1(P )|).
Hence,
(B.3) Gap(Pǫ) ≥ (1− ǫ)Gap(P ).
Proof of Lemma 5.2. Note that
(x, y) = (1 − p1 − p2)K±i (x, y) + (p1 + p2)Ix(y)
and, analogously,
PX0(x, y) = (1 − p1)K0(x, y) + p1Ix(y).
Hence, by (B.3),
Gap(PX±
) ≥ (1− p1 − p2)Gap(K±i )
as well
Gap(PX0) ≥ (1− p1)Gap(K0).
Finally note that Pi is given by
1− p2
1− p2
for every i, hence Gap(Pi) = p2. �
Proof of Proposition 5.3. To prove the first part of the proposition it is enough to
combine Lemma 5.2, Proposition 5.1 and (5.3). To complete the proof observe that
Gap(K±i ) = Gap(K0) = 1, when K
i and K0 are given by (5.2). �
In order to prove Proposition 6.1 we need some results obtained in [10].
Theorem B.2 (Ellis-Otto-Touchette). Let ρN be the distribution of SN (x)/N un-
der πβ,K,N , then ρN satisfies a large deviation principle on [−1, 1] with rate function
Ĩβ,K(z) = Jβ(z)− βKz2 − inf
{Jβ(t)− βKt2}
Jβ(z) = sup
tz − log
1 + e−β(et + e−t
1 + 2e−β
Moreover, if Ẽβ,K := argminĨβ,K, then there exists a non decreasing function Γ :
(0,+∞) → (0,+∞) with limx→0 Γ(x) = +∞ and limx→∞ Γ(x) = γc ≃ 1.082 such
that for every (β,K) with K > Γ(β) then
Ẽβ,K = {±z(β,K) 6= 0}.
20 FEDERICO BASSETTI AND FABRIZIO LEISEN
In particular, for such (β,K) and for every 0 < ǫ < |z(β,K)| there exists a constant
C1 = C1(ǫ, β,K) such that
(B.4) ρ([0, ǫ]) ≤ C1 exp{−
γǫ,β,K}
(B.5) γǫ,β,K = inf
z∈[0,ǫ]
Ĩβ,K(z) > 0.
Proof. For the first part see Theorems 3.3, 3.6 and 3.8 in [10]. As for (B.4)-(B.5),
they are standard consequences of the theory of the large deviations and of the first
part of the proposition, see, e.g., Proposition 6.4 of [11]. �
Proof of Proposition 6.1. We intend to use the Chegeer’s inequality. To do this, let
A := {x : SN (x) < 0}, B := {x : SN (x) > 0}, C := {SN (x) = 0}. First of all note
that, by symmetry, π(A) = π(B) = (1−π(C))/2 ≤ 1/2. The main task is to bound
φ(A) =
π(x)ME(x, y) =
π(y)ME(y, x).
Now, observe that if SN (y) > 1 then ME(y, x) = 0 for every x in A, hence
φ(A) =
y:SN (y)=0
ME(y, x) +
y:SN(y)=1
ME(y, x)
≤ π {y : SN (y) ∈ {0, 1}} .
This yields a bound on the conductance
h = h(π,ME) ≤ φ(A)/π(A) ≤
2π {y : SN (y) ∈ {0, 1}}
1− π{y : SN (y) = 0}
Now by Proposition B.2 we get
h(π,ME) ≤ C2e−∆N
for suitable constants C2 and ∆ > 0. The thesis follows by Cheeger inequality
(3.2). �
Proof of Lemma 6.2. The proof is exactly the same as the proof of Lemma 5.2. �
In order to prove Lemma 6.3 it is convenient to fix some simple properties of the
chain P̄ .
Lemma B.3. P̄ is a random walk on DN . If P̄ ((s, r), (s̃, r̃)) 6= 0,
P̄ ((s, r), (s̃, r̃)) ≥ p1C3
for a suitable constant C3 = C3(β,K), moreover
P̄ ((s, r), (s̃, r̃)) ≤ p1
for every (s, r), (s̃, r̃)) 6= ((0, 0), (1, 1)).
Proof of Lemma B.3. Easy but tedious computations show that
P̄ ((0, 0), (1, 1)) =
1, exp{Kβ
P̄ ((0, N), (1, N − 1)) = p1
P̄ ((0, N), (2, N)) =
METROPOLIS ALGORITHM AND EQUIENERGY SAMPLING 21
P̄ ((0, r), (2, r)) =
r = 0, 2, 4, ..., N − 2
P̄ ((0, r), (1, r − 1)) = p1
r = 0, 2, 4, ..., N − 2
P̄ ((0, r), (1, r + 1)) =
1, exp{Kβ
r = 0, 2, 4, ..., N − 2
P̄ ((s, r), (s + 2, r)) =
(r − s)
(s, r) ∈ DN , 0 < s ≤ N − 2, r ≤ N
P̄ ((s, r), (s − 2, r)) = p1
(r + s) exp{4Kβ
(1 − s)}
(s, r) ∈ DN , 0 < s ≤ N, r ≤ N
P̄ ((s, r), (s + 1, r + 1)) =
(N − r)min
1, exp{Kβ
(2s+ 1)− β}
(s, r) ∈ DN , 0 < s, r ≤ N − 1,
P̄ ((s, r), (s − 1, r + 1)) = p1
(N − r) exp{Kβ
(−2s+ 1)− β}
(s, r) ∈ DN , 0 < s, r ≤ N − 1,
P̄ ((s, r), (s + 1, r − 1)) = p1
(r − s)
(s, r) ∈ DN , 0 < r ≤ N, 0 < s ≤ N − 2
P̄ ((s, r), (s − 1, r − 1)) = p1
(r + s)min
1, exp{Kβ
(2s+ 1)− β}
(s, r) ∈ DN , 0 < r ≤ N, 0 < s ≤ r.
At this stage the statement follows easily. �
Proof of Lemma 6.3. In order to obtain a bound on the gap of P̄ we shall apply
another time the decomposition technique. Write
DN = X̄1 ∪ X̄2 ∪ X̄3 ∪ ... ∪ X̄N ,
where
X̄1 = {(0, 0), (1, 1)} X̄r = {(u1, u2) ∈ Dn : u2 = r}.
On |[N ]| := {1, ..., N} define a chain P|[N ]| setting
P|[N ]|(i, j) :=
2π̄(X̄i)
a∈X̄i
b∈X̄j
P̄ (a, b)π̄(a)
P|[N ]|(i, i) := 1−
j 6=i
P|[N ]|(i, j).
Again P|[N ]| is a reversible chain on |[N ]| with stationary distribution
π̄|[N ]|(i) := π̄(X̄i).
Finally for every r = 1, 2, . . . , N we define a chain on X̄r by setting
PX̄r(a, b) := P̄ (a, b) + Ia(b)
z∈X̄ cr
P̄ (a, z)
where both a and b belong to X̄r. Now note that for every r = 2, 3, . . . , N PX̄r is
a birth and death chain on the state space {(1, r), (3, r), . . . , (r, r)} for r odd and
22 FEDERICO BASSETTI AND FABRIZIO LEISEN
{(0, r), (2, r), . . . , (r, r)} for r even. Let
qr(s) :=
(r − s)/2
and, for r even,
qr(0) := 2
Now observe that PX̄r has stationary distribution
πr(s) ∝ qr(s)
with s = 0, 2, . . . , r if r is even and s = 1, 3, . . . , r if r is odd. First of all let r 6= 1, by
Lemma B.3 and Lemma B.1, it is easy to check that (PX̄r , πr) meets the condition
of Lemma A.1 with
B = 2, n = [(r + 2)/2], A = C3p1[(r + 2)/2]N
([x] being the integer part of x) and then
1− λ1(PX̄r ) ≥
C3p1[(r + 2)/2]
2N [(r + 2)/2]3
Finally, Lemma B.3 with (A.1) yields
λ|X̄r|−1(PX̄r ) ≥ 1− p1.
Hence, for every r 6= 1, we have proved that
(B.6) Gap(PX̄r) ≥ C3/2p1N
For r = 1
PX̄1 =
1− α1/2 α1/2
α2/2 1− α2/2
where
α1 :=
1, exp{
α2 := p1 min
1, exp{Kβ
Gap(PX̄1) ≥ 1− |
2− α1 − α2
| = α1 + α2
where the last equality follows from the fact that α1
and α2
. Hence, for
sufficiently large N , it’s easy to see that
(B.7) Gap(PX̄1) ≥ C4p1N
with C4 = C4(β,K). At this stage (B.6) with (B.7) gives
(B.8) Gap(PX̄r) ≥ C5p1N
for all r ∈ |[N ]|. As for the gap of P|[N ]|, first of all note that P|[N ]| is a birth and
death chain on |[N ]|. From Lemma B.3
P|[N ]|(i, i+1) :=
2π̄(X̄i)
a∈X̄i
b∈X̄i+1
P̄ (a, b)π̄(a) ≥ p1C3
2π̄(X̄i)
a∈X̄i
b∈X̄i+1
π̄(a) ≥ p1C3
and analogously,
P|[N ]|(i, i− 1) ≥
Now, for r 6= 1
π̄|[N ]|(r) = q|[N ]|(r)/(
q|[N ]|(i))
METROPOLIS ALGORITHM AND EQUIENERGY SAMPLING 23
while
π̄|[N ]|(1) = (q|[N ]|(1) + q|[N ]|(0))/(
q|[N ]|(i)).
So, using the unimodality of q|[N ]|, we can apply Lemma B.3 with
which gives
λ1(P|[N ]|) ≤ 1−
≤ 1− p1C3
Using another time Lemma B.3, by (A.1), we get
λN (P|[N ]|) ≥ 1− p1.
Combining this two bounds we have
(B.9) Gap(P|[N ]|) ≥
and so from (3.4)
Gap(P̄ ) ≥
C being a suitable constant that depends by β,K,C3, C4, C5. �
Proof of Proposition 6.4. Combine Lemma 6.3 with 6.5. �
Acknowledgments
We should like to thank Persi Diaconis for useful discussions and for having
encouraged us during this work, Antonietta Mira for suggesting some interesting
references and Claudio Giberti for helping to improve an earlier version of the paper.
References
[1] J. Besag and P. J. Green. Spatial statistics and Bayesian computation. J. Roy. Statist. Soc.
Ser. B, 55(1):25–37, 1993.
[2] M. Blume, V. J. Emery, and R. B. Griffiths. Ising model for the λ transition and phase
separation in he3-he4 mixtures. Phys. Rev. A, 4:1071–1077, 1971.
[3] P. Bremaud. Markov Chains. Springer-Verlag, New York, 1998.
[4] S. Caracciolo, A. Pelissetto, and A. D. Sokal. Nonlocal Monte Carlo algorithm for self-avoiding
walks with fixed endpoints. J. Statist. Phys., 60(1-2):1–53, 1990.
[5] S. Caracciolo, A. Pelissetto, and A. D. Sokal. Dynamic critical exponent of the BFACF
algorithm for self-avoiding walks. J. Statist. Phys., 63(5-6):857–865, 1991.
[6] P. Diaconis and L. Saloff-Coste. What do we know about the Metropolis algorithm? J.
Comput. System Sci., 57(1):20–36, 1998. 27th Annual ACM Symposium on the Theory of
Computing (STOC’95) (Las Vegas, NV).
[7] P. Diaconis and D. Stroock. Geometric bounds for eigenvalues of Markov chains. Ann. Appl.
Probab., 1(1):36–61, 1991.
[8] M. Dyer, L. A. Goldberg, M. Jerrum, and R. Martin. Markov chain comparison. Probab.
Surv., 3:89–111 (electronic), 2006.
[9] R. G. Edwards and A. D. Sokal. Generalization of the Fortuin-Kasteleyn-Swendsen-Wang
representation and Monte Carlo algorithm. Phys. Rev. D (3), 38(6):2009–2012, 1988.
[10] R. S. Ellis, P. T. Otto, and H. Touchette. Analysis of phase transitions in the mean-field
Blume-Emery-Griffiths model. Ann. Appl. Probab., 15(3):2203–2254, 2005.
[11] R.S. Ellis. The Theory of Large Deviation and Applications to Statistical Mechanics.
2006. Lectures for the international seminar on Extreme Events in Complex Dynamics.
http://www.math.umass.edu/ ˜ rsellis/pdf-files/Dresden-lectures.pdf Max-Planck-Institut.
[12] A. E. Gelfand, S. K. Sahu, and B. P. Carlin. Efficient parameterisations for normal linear
mixed models. Biometrika, 82(3):479–488, 1995.
[13] A. Gelman, G. O. Roberts, and W. R. Gilks. Efficient Metropolis jumping rules. In Bayesian
statistics, 5 (Alicante, 1994), Oxford Sci. Publ., pages 599–607. Oxford Univ. Press, New
York, 1996.
http://www.math.umass.edu/
24 FEDERICO BASSETTI AND FABRIZIO LEISEN
[14] C. J. Gey er and E. A. Thompson. Annealing markov chain monte carlo with applications to
ancestral inference. J. Amer. Statist. Assoc., 90:909–920, 1995.
[15] Y. Guan, R. Fleißner, P. Joyce, and S. M. Krone. Markov chain Monte Carlo in small worlds.
Stat. Comput., 16(2):193–202, 2006.
[16] Y. Guan and S. M. Krone. Small-world mcmc and convergence to multi-modal distributions:
From slow mixing to fast mixing. Ann. Appl. Probab., 17(1):284–304, 2007.
[17] J. M. Hammersley and D. C. Handscomb. Monte Carlo methods. Methuen & Co. Ltd., Lon-
don, 1965.
[18] W.K. Hastings. Monte carlo sampling methods using markov chains and their application.
Biometrika, 57:97–109, 1970.
[19] K. Hukushima and K. Nemoto. Exchange monte carlo method and application to spin glass
simulations. J.Phys.Soc.Jpn., 65:1604–1608, 1996.
[20] M. Jerrum, J. Son, P. Tetali, and E. Vigoda. Elementary bounds on Poincaré and log-Sobolev
constants for decomposable Markov chains. Ann. Appl. Probab., 14(4):1741–1765, 2004.
[21] S. A. Kivelson, V. J. Emery, and H. Q. Lin. Doped antiferromagnets in the weak-hopping
limit. Phys. Rev. B, 42:6523–6530, 1990.
[22] S. C. Kou, Qing Zhou, and Wing Hung Wong. Equi-energy sampler with applications in
statistical inference and statistical mechanics. Ann. Statist., 34(4):1581–1619, 2006.
[23] J. Lajzerowicz and J. Sivardiére. Spin 1 lattice gas model. ii. condensation and phase sepa-
ration in a binary fluid. Phys. Rev. A, 11:2090–2100, 1975.
[24] J. Lajzerowicz and J. Sivardiére. Spin 1 lattice gas model. iii. tricritical points in binary and
ternary fluids. Phys. Rev. A, 11:2101–2110, 1975.
[25] J. S. Liu. Monte Carlo strategies in scientific computing. Springer Series in Statistics.
Springer-Verlag, New York, 2001.
[26] N. Madras and M. Piccioni. Importance sampling for families of distributions. Ann. Appl.
Probab., 9(4):1202–1225, 1999.
[27] N. Madras and D. Randall. Markov chain decomposition for convergence rate analysis. Ann.
Appl. Probab., 12(2):581–606, 2002.
[28] E. Marinari and G. Parisi. Simulated tempering: a new monte carlo scheme. Europhy s. Lett,
19:451–458, 1992.
[29] N. Metrpolis, A. Rosenbluth, M. Rosenbluth, A. Teller, and E. Teller. Equation of state
caluclations by fast computing machines. J. Chem. Phys., 21:1087–1092, 1953.
[30] A. Mira and L. Tierney. Efficiency and convergence properties of slice samplers. Scand. J.
Statist., 29(1):1–12, 2002.
[31] K.E. Newman and J.D. Dow. Zinc blende diamond order disorder transition in metastable
crystalline (gaas)1-xge2x alloys. Phys. Rev. B, 27:7495–7508, 1983.
[32] P. H. Peskun. Optimum Monte-Carlo sampling using Markov chains. Biometrika, 60:607–612,
1973.
[33] P. H. Peskun. Guidelines for choosing the transition matrix in Monte Carlo methods using
Markov chains. J. Comput. Phys., 40(2):327–344, 1981.
[34] C. P. Robert and G. Casella. Monte Carlo statistical methods. Springer Texts in Statistics.
Springer-Verlag, New York, second edition, 2004.
[35] R. Y. Rubinstein. Simulation and the Monte Carlo method. John Wiley & Sons Inc., New
York, 1981. Wiley Series in Probability and Mathematical Statistics.
[36] M. Schick and Wei-Heng Shih. Spin 1 model of a microemulsion. Phys. Rev. B, 34:1797–1801,
1986.
[37] A. Sinclair. Algorithms for random generation and counting. Progress in Theoretical Com-
puter Science. Birkhäuser Boston Inc., Boston, MA, 1993. A Markov chain approach.
[38] J. Sivardiére and J. Lajzerowicz. Spin 1 lattice gas model. i. condensation and solidification
of a simple fluid. Phys. Rev. A, 11:2090–2100, 1975.
[39] A. Sokal. Monte Carlo methods in statistical mechanics: foundations and new algorithms. In
Functional integration (Cargèse, 1996), volume 361 of NATO Adv. Sci. Inst. Ser. B Phys.,
pages 131–192. Plenum, New York, 1997.
[40] R.H. Swendsen and J. S. Wang. Non–universal critical dynamics in monte carlo simulations.
Phy. Rev. Lett., pages 86–88, 1987.
Università degli Studi di Pavia, Dipartimento di Matematica, via Ferrata 1, 27100
Pavia, Italy
Università dell’Insubria, Dipartimento di economia, via monte generoso 71, 21100
Varese , Italy
E-mail address: federico.bassetti@unipv.it
E-mail address: leisen.fabrizio@unimore.it
	1. Introduction.
	2. A general strategy
	3. Preliminaries
	3.1. Spectral gap, variance and speed of convergence
	3.2. Cheeger's inequality
	3.3. Chain decomposition theorem
	4. Warming up example
	5. The mean field Ising model
	6. The mean–field Blume-Emery-Griffiths model
	Appendix A. The Spectral Gap of a Birth and Death Chain
	Appendix B. Proofs
	Acknowledgments
	References
