Synthesis of MCMC and Belief Propagation
Sungsoo Ahn

SUNGSOO . AHN @ KAIST. AC . KR

Michael Chertkov

CHERTKOV @ LANL . GOV

arXiv:1605.09042v5 [stat.ML] 21 Nov 2016

Jinwoo Shin

JINWOOS @ KAIST. AC . KR

Abstract
Markov Chain Monte Carlo (MCMC) and Belief Propagation (BP) are the most popular algorithms for
computational inference in Graphical Models (GM). In principle, MCMC is an exact probabilistic method
which, however, often suffers from exponentially slow mixing. In contrast, BP is a deterministic method, which
is typically fast, empirically very successful, however in general lacking control of accuracy over loopy graphs.
In this paper, we introduce MCMC algorithms correcting the approximation error of BP, i.e., we provide a
way to compensate for BP errors via a consecutive BP-aware MCMC. Our framework is based on the Loop
Calculus approach which allows to express the BP error as a sum of weighted generalized loops. Although
the full series is computationally intractable, it is known that a truncated series, summing up all 2-regular
loops, is computable in polynomial-time for planar pair-wise binary GMs and it also provides a highly accurate
approximation empirically. Motivated by this, we first propose a polynomial-time approximation MCMC
scheme for the truncated series of general (non-planar) pair-wise binary models. Our main idea here is to use
the Worm algorithm, known to provide fast mixing in other (related) problems, and then design an appropriate
rejection scheme to sample 2-regular loops. Furthermore, we also design an efficient rejection-free MCMC
scheme for approximating the full series. The main novelty underlying our design is in utilizing the concept
of cycle basis, which provides an efficient decomposition of the generalized loops. In essence, the proposed
MCMC schemes run on transformed GM built upon the non-trivial BP solution, and our experiments show that
this synthesis of BP and MCMC outperforms both direct MCMC and bare BP schemes.

1. Introduction
GMs express factorization of the joint multivariate probability distributions in statistics via graph of relations
between variables. The concept of GM has been used successfully in information theory, physics, artificial
intelligence and machine learning [1, 2, 3, 4, 5, 6]. Of many inference problems one can set with a GM,
computing partition function (normalization), or equivalently marginalizing the joint distribution, is the most
general problem of interest. However, this paradigmatic inference problem is known to be computationally
intractable in general, i.e., formally it is #P-hard even to approximate [7, 8].
To address this obstacle, extensive efforts have been made to develop practical approximation methods,
among which MCMC- [9] based and BP- [10] based algorithms are, arguably, the most popular and practically
successful ones. MCMC is exact, i.e., it converges to the correct answer, but its convergence/mixing is, in
general, exponential in the system size. On the other hand, message passing implementations of BP typically
demonstrate fast convergence, however in general lacking approximation guarantees for GM containing loops.
Motivated by this complementarity of the MCMC and BP approaches, we aim here to synthesize a hybrid
approach benefiting from a joint use of MCMC and BP.
At a high level, our proposed scheme uses BP as the first step and then runs MCMC to correct for the
approximation error of BP. To design such an “error-correcting" MCMC, we utilize the Loop Calculus approach
[11] which allows, in a nutshell, to express the BP error as a sum (i.e., series) of weights of the so-called
generalized loops (sub-graphs of a special structure). There are several challenges one needs to overcome.
First of all, to design an efficient Markov Chain (MC) sampler, one needs to design a scheme which allows
efficient transitions between the generalized loops. Second, even if one designs such a MC which is capable of
accessing all the generalized loops, it may mix slowly. Finally, weights of generalized loops can be positive or
negative, while an individual MCMC can only generate non-negative contributions.
1

Since approximating the full loop series (LS) is intractable in general, we first explore whether we can
deal with the challenges at least in the case of the truncated LS corresponding to 2-regular loops. In fact, this
problem has been analyzed in the case of the planar pairwise binary GMs [12, 13] where it was shown that
the 2-regular LS is computable exactly in polynomial-time through a reduction to a Pfaffian (or determinant)
computation [14]. In particular, the partition function of the Ising model without external field (i.e., where
only pair-wise factors present) is computable exactly via the 2-regular LS. Furthermore, the authors show that
in the case of general planar pairwise binary GMs, the 2-regular LS provides a highly accurate approximation
empirically. Motivated by these results, we address the same question in the general (i.e., non-planar) case of
pairwise binary GMs via MCMC. For the choice of MC, we adopt the Worm algorithm [15]. We prove that
with some modification including rejections, the algorithm allows to sample (with probabilities proportional to
respective weights) 2-regular loops in polynomial-time. Then, we design a novel simulated annealing strategy
using the sampler to estimate separately positive and negative parts of the 2-regular LS. Given any ε > 0, this
leads to a ε-approximation polynomial-time scheme for the 2-regular LS under a mild assumption.
We next turn to estimating the full LS. In this part, we ignore the theoretical question of establishing the
polynomial mixing time of a MC, and instead focus on designing an empirically efficient MCMC scheme. We
design an MC using a cycle basis of the graph [16] to sample generalized loops directly, without rejections.
It transits from one generalized loop to another by adding or deleting a random element of the cycle basis.
Using the MC sampler, we design a simulated annealing strategy for estimating the full LS, which is similar
to what was used earlier to estimate the 2-regular LS. Notice that even though the prime focus of this paper
is on pairwise binary GMs, the proposed MCMC scheme allows straightforward generalization to general
non-binary GMs.
In summary, we propose novel MCMC schemes to estimate the LS correction to the BP contribution to the
partition function. Since already the bare BP provides a highly non-trivial estimation for the partition function,
it is naturally expected and confirmed in our experimental results that the proposed algorithm outperforms
other standard (not related to BP) MCMC schemes applied to the original GM. We believe that our approach
provides a new angle for approximate inference on GM and is of broader interest to various applications
involving GMs.

2. Preliminaries
2.1 Graphical models and belief propagation
Given undirected graph G = (V, E) with |V | = n, |E| = m, a pairwise binary Markov Random Fields (MRF)
defines the following joint probability distribution on x = [xv ∈ {0, 1} : v ∈ V ]:
p(x) =

Y
1 Y
ψv (xv )
ψu,v (xu , xv ),
Z
v∈V

Z :=

X

Y

ψv (xv )

x∈{0,1}n v∈V

(u,v)∈E

Y

ψu,v, (xu , xv )

(u,v)∈E

where ψv , ψu,v are some non-negative functions, called compatibility or factor functions, and the normalization
constant Z is called the partition function. Without loss of generality, we assume G is connected. It is known
that approximating the partition function is #P-hard in general [8]. Belief Propagation (BP) is a popular
message-passing heuristic for approximating marginal distributions of MRF. The BP algorithm iterates the
following message updates for all (u, v) ∈ E:
X
Y
mt+1
ψu,v (xu , xv )ψu (xu )
mtw→u (xu ),
u→v (xv ) ∝
xu ∈{0,1}

w∈N (u)\v

where N (v) denotes the set of neighbors of v. In general BP may fail to converge, however in this case one
may substitute it with a somehow more involved algorithm provably convergent to its fixed point [22, 23, 24].
Estimates for the marginal probabilities are expressed via the fixed-point messages {mu→v : (u, v) ∈ E} as

2

follows: τv (xv ) ∝ ψv (xv )

Q

u∈N (v)

mu→v (xv ) and


τu,v (xu , xv ) ∝ ψu (xu )ψv (xv )ψu,v (xu , xv ) 


Y

mw→v (xu ) 

w∈N (u)


Y

mw→v (xv ) .

w∈N (v)

2.2 Bethe approximation and loop calculus
BP marginals also results in the following Bethe approximation for the partition function Z:
XX
X X
log ZBethe =
τv (xv ) log ψv (xv ) +
τu,v (xu , xv ) log ψu,v (xu , xv )
v∈V xv

−

(u,v)∈E xu ,xv

XX

X

τv (xv ) log τv (xv ) −

v∈V xv

X

τu,v (xu , xv ) log

(u,v)∈E xu ,xv

τu,v (xu , xv )
τu (xu )τv (xv )

If graph G is a tree, the Bethe approximation is exact, i.e., ZBethe = Z. However, in general, i.e. for the graph
with cycles, BP algorithm provides often rather accurate but still an approximation.
Loop Series (LS) [11] expresses, Z/ZBethe , as the following sum/series:
Z
ZBethe
w(F ) :=



Y
(u,v)∈EF

= ZLoop :=

X

w(F ),

w(∅) = 1,

F ∈L

 Y
τu,v (1, 1)
−1
τu (1)τv (1)

τv (1) + (−1)

dF (v)

v∈VF



τv (1)
1 − τv (1)

dF (v)−1

!
τv (1)

where each term/weight is associated with the so-called generalized loop F and L denotes the set of all
generalized loops in graph G (including the empty subgraph ∅). Here, a subgraph F of G is called generalized
loop if all vertices v ∈ F have degree dF (v) (in the subgraph) no smaller than 2.
Since the number of generalized loops is exponentially large, computing ZLoop is intractable in general.
However, the following truncated sum of ZLoop , called 2-regular loop series, is known to be computable in
polynomial-time if G is planar [12]:1
X
Z2-Loop :=
w(F ),
F ∈L2-Loop

where L2-Loop denotes the set of all 2-regular generalized loops, i.e., F ∈ L2-Loop if dF (v) = 2 for every vertex
v of F . One can check that ZLoop = Z2-Loop for the Ising model without the external fields. Furthermore, as
stated in [12, 13] for the general case, Z2-Loop provides a good empirical estimation for ZLoop .

3. Estimating 2-regular loop series via MCMC
In this section, we aim to describe how the 2-regular loop series Z2-Loop can be estimated in polynomial-time.
To this end, we first assume that the maximum degree ∆ of the graph G is at most 3. This degree constrained
assumption is not really restrictive since any pairwise binary model can be easily expressed as an equivalent
one with ∆ ≤ 3, e.g., see the supplementary material. The rest of this section consists of two parts. We first
propose an algorithm generating a 2-regular loop sample with the probability proportional to the absolute
value of its weight, i.e.,
π2-Loop (F ) :=

|w(F )|
†
Z2-Loop

,

†
where Z2-Loop
=

X
F ∈L2-Loop

1. Note that the number of 2-regular loops is exponentially large in general.

3

|w(F )|.

Note that this 2-regular loop contribution allows the following factorization: for any F ∈ L2-Loop ,
|w(F )| =

Y

w(e),

where w(e) := p

e∈F

τu,v (1, 1) − τu (1)τv (1)
τu (1)τv (1)(1 − τu (1))(1 − τv (1))

.

(1)

In the second part we use the sampler constructed in the first part to design a simulated annealing scheme to
estimate Z2-Loop .
3.1 Sampling 2-regular loops
We suggest to sample the 2-regular loops distributed according to π2-Loop through a version of the Worm
S algorithm proposed by Prokofiev and Svistunov [15]. It can be viewed as a MC exploring the set, L2-Loop L2-Odd ,
where L2-OddSis the set of all subgraphs of G with exactly two odd-degree vertices. Given current state
F ∈ L2-Loop L2-Odd , it chooses the next state F 0 as follows:
1. If F ∈ L2-Odd , pick a random vertex v (uniformly) from V . Otherwise, pick a random odd-degree vertex
v (uniformly) from F .
2. Choose a random neighbor u of v (uniformly) within G, and set F 0 ← F initially.
3. Update F 0 ← F ⊕ {u, v} with the probability



1 |w(F ⊕{u,v})|

,
1
min

|w(F )|


n

n |w(F ⊕{u,v})|
min 4
,1
|w(F
)|





min d(v) |w(F ⊕{u,v})| , 1
2d(u)
|w(F )|

if F ∈ L2-Loop
else if F ⊕ {u, v} ∈ L2-Loop
else if F, F ⊕ {u, v} ∈ L2-Odd

Here,
⊕ denotes the symmetric difference and for F ∈ L2-Odd , its weight is defined according to w(F ) =
Q
e∈F w(e). In essence, the Worm algorithm consists in either deleting or adding an edge to the current
subgraph F . From the Worm algorithm, we transition to the following algorithm which samples 2-regular
loops with probability π2-Loop simply by adding rejection of F if F ∈ L2-Odd .
Algorithm 1 Sampling 2-regular loops
1:
2:
3:
4:
5:
6:
7:
8:
9:

Input: Number of trials N ; number of iterations T of the Worm algorithm
Output: 2-regular loop F .
for i = 1 → N do
Set F ← ∅ and update it T times by running the Worm algorithm
if F is a 2-regular loop then
BREAK and output F .
end if
end for
Output F = ∅.
The following theorem states that Algorithm 1 can generate a desired random sample in polynomial-time.

Theorem 1. Given δ > 0, choose inputs of Algorithm 1 as
N ≥ 1.2 n log(3δ −1 ),

and

T ≥ (m − n + 1) log 2 + 4∆mn4 log(3nδ −1 ).

Then, it follows that
1
2

X



P Algorithm 1 outputs F − π2-Loop (F ) ≤ δ.

F ∈L2-Loop

namely, the total variation distance between π2-Loop and the output distribution of Algorithm 1 is at most δ.
4

The proof of the above theorem is presented in the supplementary material due to the space constraint. In
the proof, we first show that MC induced by the Worm algorithm mixes in polynomial time, and then prove
that acceptance of a 2-regular loop, i.e., line 6 of Algorithm 1, occurs with high probability. Notice that the
uniform-weight version of the former proof, i.e., fast mixing, was recently proven in [18]. For completeness
of the material exposition, we present the general case proof of interest for us. The latter proof, i.e., high
acceptance, requires to bound |L2-Loop | and |L2-Odd | to show that the probability of sampling 2-regular loops
under the Worm algorithm is 1/poly(n) for some polynomial function poly(n).
3.2 Simulated annealing for approximating 2-regular loop series
Here we utilize Theorem 1 to describe an algorithm approximating the 2-regular LS Z2-Loop in polynomial time.
To achieve this goal, we rely on the simulated annealing strategy [19] which requires to decide a monotone
cooling schedule β0 , β1 , . . . , β`−1 , β` , where β` corresponds to the target counting problem and β0 does to its
relaxed easy version. Thus, designing an appropriate cooling strategy is the first challenge to address. We
will also describe how to deal with the issue that Z2-Loop is a sum of positive and negative terms, while most
simulated annealing strategies in the literature mainly studied on sums of non-negative terms. This second
challenge is related to the so-called ‘fermion sign problem’ common in statistical mechanics of quantum
systems [25]. Before we describe the proposed algorithm in details, let us provide its intuitive sketch.
†
The proposed algorithm consists of two parts: a) estimating Z2-Loop
via a simulated annealing strategy
†
and b) estimating Z2-Loop /Z2-Loop
via counting samples corresponding to negative terms in the 2-regular loop
series. First consider the following β-parametrized, auxiliary distribution over 2-regular loops:

π2-Loop (F : β) =

1
†
Z2-Loop
(β)

|w(F )|β ,

for 0 ≤ β ≤ 1.

(2)

Note that one can generate samples approximately with probability (2) in polynomial-time using Algorithm 1
by setting w ← wβ . Indeed, it follows that for β 0 > β,
†
Z2-Loop
(β 0 )
†
(β)
Z2-Loop

=

X

|w(F )|β

F ∈L2-Loop

0

−β

|w(F )|β
†
(β)
Z2-Loop

h
i
0
= Eπ2-Loop (β) |w(F )|β −β ,

where the expectation can be estimated using O(1) samples if it is Θ(1), i.e., β 0 is sufficiently close to β.
Then, for any increasing sequence β0 = 0, β1 , . . . , βn−1 , βn = 1, we derive
†
Z2-Loop

=

†
Z2-Loop
(βn )
†
Z2-Loop
(βn−1 )

·

†
Z2-Loop
(βn−1 )
†
Z2-Loop
(βn−2 )

···

†
†
Z2-Loop
(β2 ) Z2-Loop
(β1 )
†
Z2-Loop
(β1 )

†
Z2-Loop
(β0 )

†
Z2-Loop
(0),

†
where it is know that Z2-Loop
(0), i.e., the total number of 2-regular loops, is exactly 2m−n+1 [16]. This allows


†
us to estimate Z2-Loop simply by estimating Eπ2-Loop (βi ) |w(F )|βi+1 −βi for all i.
†
Our next step is to estimate the ratio Z2-Loop /Z2-Loop
. Let L−
2-Loop denote the set of negative 2-regular loops,

i.e.,
L−
2-Loop := {F : F ∈ L2-Loop , w(F ) < 0}.
Then, the 2-regular loop series can be expressed as
P
!



|w(F )|
F ∈L−
2-Loop
†
†
Z2-Loop = 1 − 2
Z2-Loop = 1 − 2Pπ2–Loop w(F ) < 0 Z2-Loop
,
†
Z2-Loop


where we estimate Pπ2–Loop w(F ) < 0 again using samples generated by Algorithm 1.
We provide the formal description of the proposed algorithm and its error bound as follows.

5

Algorithm 2 Approximation for Z2-Loop
Input: Increasing sequence β0 = 0 < β1 < · · · < βn−1 < βn = 1; number of samples s1 , s2 ; number
of trials N1 ; number of iterations T1 for Algorithm 1.
2: for i = 0 → n − 1 do
3:
Generate 2-regular loops F1 , . . . , Fs1 for π2-Loop (βi ) using Algorithm 1 with input N1 and T1 , and set

1:

Hi ←
4:
5:

end for
Generate 2-regular loops F1 , . . . , Fs2 for π2-Loop using Algorithm 1 with input N2 and T2 , and set
κ←

6:

1 X
w(Fj )βi+1 −βi .
s1 j

|{Fj : w(Fj ) < 0}|
.
s2

b2-Loop ← (1 − 2κ)2m−n+1 Q Hi .
Output: Z
i

Theorem 2. Given ε, ν > 0, choose inputs of Algorithm 2 as βi = i/n for i = 1, 2, . . . , n − 1,
−1
s1 ≥ 18144n2 ε−2 wmin
dlog(6nν −1 )e,

−1
N1 ≥ 1.2n log(144nε−1 wmin
),

−1
T1 ≥ (m − n + 1) log 2 + 4∆mn4 log(48nε−1 wmin
),

s2 ≥ 18144ζ(1 − 2ζ)−2 ε−2 dlog(3ν −1 )e,

N2 ≥ 1.2n log(144ε−1 (1 − 2ζ)−1 ),

T2 ≥ (m − n + 1) log 2 + 4∆mn4 log(48ε−1 (1 − 2ζ)−1 )
where wmin = mine∈E w(e) and ζ = Pπ2–Loop [w(F ) < 0]. Then, the following statement holds
#
"
b2-Loop − Z2-Loop |
|Z
≤ ε ≤ 1 − ν,
P
Z2-Loop
which means Algorithm 2 estimates Z2-Loop within approximation ratio 1 ± ε with high probability.
The proof of the above theorem is presented in the supplementary material due to the space constraint.
We note that all constants entering in Theorem 2 were not optimized. Theorem 2 implies that complexity of
−1
Algorithm 2 is polynomial with respect to n, 1/ε, 1/ν under assumption that wmin
and 1 − 2Pπ2–Loop [w(F ) <
−1
0] are polynomially small. Both wmin and 1 − 2Pπ2–Loop [w(F ) < 0] depend on the choice of BP fixed
point, however it is unlikely (unless a degeneracy) that these characteristics become large. In particular,
Pπ2–Loop [w(F ) < 0] = 0 in the case of attractive models [20].

4. Estimating full loop series via MCMC
In this section, we aim for estimating the full loop series ZLoop . To this end, we design a novel MC sampler for
generalized loops, which adds (or removes) a cycle basis or a path to (or from) the current generalized loop.
Therefore, we naturally start this section introducing necessary backgrounds on cycle basis. Then, we turn to
describe the design of MC sampler for generalized loops. Finally, we describe a simulated annealing scheme
similar to the one described in the preceding section. We also report its experimental performance comparing
with other methods.
4.1 Sampling generalized loops with cycle basis
The cycle basis C of the graph G is a minimal set of cycles which allows to represent every Eulerian subgraph
of G (i.e., subgraphs containing no odd-degree vertex) as a symmetric difference of cycles in the set [16]. Let
6

us characterize the combinatorial structure of the generalized loop using the cycle basis. To this end, consider
a set of paths between any pair of vertices:
P = {Pu,v : u 6= v, u, v ∈ V, Pu,v is a path from u to v},

i.e., |P| = n2 . Then the following theorem allows to decompose any generalized loop with respect to any
selected C and P.
Theorem 3. Consider any cycle basis C and path set P. Then, for any generalized loop F , there exists a
decomposition, B ⊂ C ∪ P, such that F can be expressed as a symmetric difference of the elements of B, i.e.,
F = B1 ⊕ B2 ⊕ · · · Bk−1 ⊕ Bk for some Bi ∈ B.
The proof of the above theorem is given in the supplementary material due to the space constraint. Now
given any choice of C, P, consider the following transition from F ∈ L, to the next state F 0 :
1. Choose, uniformly at random, an element B ∈ C ∪ P, and set F 0 ← F initially.
n
o
(
⊕B|)|
F ⊕ B with probability min 1, |w(F
|w(F )|
0
2. If F ⊕ B ∈ L, update F ←
.
F
otherwise
Due to Theorem 3, it is easy to check that the proposed MC is irreducible and aperiodic, i.e., ergodic, and the
distribution of its t-th state converges to the following stationary distribution as t → ∞:
πLoop (F ) =

|w(F )|
†
ZLoop

,

†
where ZLoop
=

X

|w(F )|.

F ∈LLoop

One also has a freedom in choosing C, P. To accelerate mixing of MC, we suggest to choose the minimum
weighted cycle basis C and the shortest paths P with respect to the edge weights {log w(e)} defined in (1),
which are computable using the algorithm in [16] and the Bellman-Ford algorithm [21], respectively. This
encourages transitions between generalized loops with similar weights.
4.2 Simulated annealing for approximating full loop series
Algorithm 3 Approximation for ZLoop
1:
2:

3:
4:
5:
6:

7:

Input: Decreasing sequence β0 > β1 > · · · > β`−1 > β` = 1; number of samples s0 , s1 , s2 ; number of
iterations T0 , T1 , T2 for the MC described in Section 4.1
Generate generalized loops F1 , · · · , Fs0 by running T0 iterations of the MC described in Section 4.1 for
πLoop (β0 ), and set
s0
U ← ∗ |w(F ∗ )|β0 ,
s
where F ∗ = arg maxF ∈{F1 ,··· ,Fs0 } |w(F )| and s∗ is the number of F ∗ sampled.
for i = 0 → ` − 1 do
Generate generalized loops P
F1 , · · · , Fs1 by running T1 iterations of the MC described in Section 4.1
for πLoop (βi ), and set Hi ← s11 j |w(Fj )|βi+1 −βi .
end for
Generate generalized loops F1 , · · · Fs2 by running T2 iterations of the MC described in Section 4.1 for
πLoop , and set
|{Fj : w(Fj ) < 0}|
κ←
.
s2
bLoop ← (1 − 2κ) Q Hi U .
Output: Z
i

7

(a)

(b)

(c)

Figure 1: Plots of the log-partition function approximation error with respect to (average) interaction strength.
Each point is averaged over 20 (random) models.
Now we are ready to describe a simulated annealing scheme for estimating ZLoop . It is similar, in
principle, with that in Section 3.2. First, we again introduce the following β-parametrized, auxiliary probability
†
distribution πLoop (F : β) = |w(F )|β /ZLoop
(β). For any decreasing sequence of annealing parameters,
β0 , β1 , · · · , β`−1 , β` = 1, we derive
†
ZLoop
=

†
ZLoop
(β` )

·

†
ZLoop
(β`−1 )

†
†
ZLoop
(β`−1 ) ZLoop
(β`−2 )

···

†
†
ZLoop
(β2 ) ZLoop
(β1 ) †
·
ZLoop (β0 ).
†
†
ZLoop (β1 ) ZLoop (β0 )
0

†
†
Following similar procedures in Section 3.2, one can estimate ZLoop
(β 0 )/ZLoop
(β) = EπLoop (β) [|w(F )|β −β ]
†
using the sampler described in Section 4.1. Moreover, ZLoop
(β0 ) = |w(F ∗ )|/PπLoop (β0 ) (F ∗ ) is estimated by
∗
sampling generalized loop F with the highest probability PπLoop (β0 ) (F ∗ ). For large enough β0 , the approximation error becomes relatively small since PπLoop (β0 ) (F ∗ ) ∝ |w(F ∗ )|β0 dominates over the distribution. In
combination, this provides a desired approximation for ZLoop . The result is stated formally in Algorithm 3.

4.3 Experimental results
In this section, we report experimental results for computing partition function of the Ising model and the
hard-core model. We compare Algorithm 2 in Section 3 (coined MCMC-BP-2reg) and Algorithm 3 in Section
4.2 (coined MCMC-BP-whole), with the bare Bethe approximation (coined BP) and the popular Gibbs-sampler
(coined MCMC-Gibbs). To make the comparison fair, we use the same annealing scheme for all MCMC
schemes, thus making their running times comparable. More specifically, we generate each sample after
running T1 = 1, 000 iterations of an MC and take s1 = 100 samples to compute each estimation (e.g., Hi ) at
intermediate steps. For performance measure, we use the log-partition function approximation error defined
as | log Z − log Zapprox |/| log Z|, where Zapprox is the output of the respective algorithm. We conducted 3
experiments on the 4 × 4 grid graph. In our first experimental setting, we consider the Ising model with varying
interaction strength and no external (magnetic) field. To prepare the model of interest, we start from the Ising
model with uniform (ferromagnetic/attractive and anti-ferromagnetic/repulsive) interaction strength and then
add ‘glassy’ variability in the interaction strength modeled via i.i.d Gaussian random variables with mean 0
and variance 0.52 , i.e. N (0, 0.52 ). In other words, given average interaction strength 0.3, each interaction
strength in the model is independently chosen as N (0.3, 0.52 ). The second experiment was conducted by
adding N (0, 0.52 ) corrections to the external fields under the same condition as in the first experiment. In
this case we observe that BP often fails to converge, and use the Concave Convex Procedure (CCCP) [23]
for finding BP fixed points. Finally, we experiment with the hard-core model on the 4 × 4 grid graph with
varying a positive parameter λ > 0, called ‘fugacity’ [26]. As seen clearly in Figure 1, BP and MCMC-Gibbs
are outperformed by MCMC-BP-2reg or MCMC-BP-whole at most tested regimes in the first experiment with
no external field, where in this case, the 2-regular loop series (LS) is equal to the full one. Even in the regimes
8

where MCMC-Gibbs outperforms BP, our schemes correct the error of BP and performs at least as good as
MCMC-Gibbs. In the experiments, we observe that advantage of our schemes over BP is more pronounced
when the error of BP is large. A theoretical reasoning behind this observation is as follows. If the performance
of BP is good, i.e. the loop series (LS) is close to 1, the contribution of empty generalized loop, i.e., w(∅), in
LS is significant, and it becomes harder to sample other generalized loops accurately.

5. Conclusion
In this paper, we propose new MCMC schemes for approximate inference in GMs. The main novelty of our
approach is in designing BP-aware MCs utilizing the non-trivial BP solutions. In experiments, our BP based
MCMC scheme also outperforms other alternatives. We anticipate that this new technique will be of interest to
many applications where GMs are used for statistical reasoning.

Acknowledgements.
The work of MC was carried out under the auspices of the National Nuclear Security Administration of the
U.S. Department of Energy at Los Alamos National Laboratory under Contract No. DE-AC52-06NA25396,
and was partially supported by the Advanced Grid Modeling Program in the U.S. Department of Energy Office
of Electricity.

References
[1] J. Pearl, “Probabilistic reasoning in intelligent systems: networks of plausible inference,” Morgan Kaufmann, 2014.
[2] R. G. Gallager, “Low-density parity-check codes,” Information Theory, IRE Transactions 8(1): 21-28,
1962.
[3] R. F. Kschischang, and J. F. Brendan, “Iterative decoding of compound codes by probability propagation
in graphical models,” Selected Areas in Communications, IEEE Journal 16(2): 219-230, 1998.
[4] M. I. Jordan, ed. “Learning in graphical models,” Springer Science & Business Media 89, 1998.
[5] R.J. Baxter, “Exactly solved models in statistical mechanics,” Courier Corporation, 2007.
[6] W.T. Freeman, C.P. Egon, and T.C. Owen, “Learning low-level vision.” International journal of computer
vision 40(1): 25-47, 2000.
[7] V. Chandrasekaran, S. Nathan, and H. Prahladh, “Complexity of Inference in Graphical Models,” Association for Uncertainty and Artificial Intelligence, 2008
[8] M. Jerrum, and A. Sinclair, “Polynomial-time approximation algorithms for the Ising model,” SIAM
Journal on computing 22(5): 1087-1116, 1993.
[9] C. Andrieu, N. Freitas, A. Doucet, and M. I. Jordan, “An introduction to MCMC for machine learning,”
Machine learning 50(1-2), 5-43, 2003.
[10] J. Pearl, “Reverend Bayes on inference engines: A distributed hierarchical approach,” Association for the
Advancement of Artificial Intelligence, 1982.
[11] M. Chertkov, and V. Y. Chernyak, “Loop series for discrete statistical models on graphs,” Journal of
Statistical Mechanics: Theory and Experiment 2006(6): P06009, 2006.
[12] M. Chertkov, V. Y. Chernyak, and R. Teodorescu, “Belief propagation and loop series on planar graphs,”
Journal of Statistical Mechanics: Theory and Experiment 2008(5): P05003, 2008.
9

[13] V. Gomez, J. K. Hilbert, and M. Chertkov, “Approximate inference on planar graphs using Loop Calculus
and Belief Propagation,” The Journal of Machine Learning Research, 11: 1273-1296, 2010.
[14] P. W. Kasteleyn, “The statistics of dimers on a lattice,” Classic Papers in Combinatorics. Birkhäuser
Boston, 281-298, 2009.
[15] N. Prokof’ev, and B. Svistunov, “Worm algorithms for classical statistical models,” Physical review
letters 87(16): 160601, 2001.
[16] J.D. Horton, “A polynomial-time algorithm to find the shortest cycle basis of a graph.” SIAM Journal on
Computing 16(2): 358-366, 1987. APA
[17] H. A. Kramers, and G. H. Wannier, “Statistics of the two-dimensional ferromagnet. Part II,” Physical
Review 60(3): 263, 1941.
[18] A. Collevecchio, T. M. Garoni, T.Hyndman, and D. Tokarev, “The worm process for the Ising model is
rapidly mixing,” arXiv preprint arXiv:1509.03201, 2015.
[19] S. Kirkpatrick, “Optimization by simulated annealing: Quantitative studies.” Journal of statistical physics
34(5-6): 975-986, 1984.
[20] R. Nicholas, “The Bethe partition function of log-supermodular graphical models,” Advances in Neural
Information Processing Systems. 2012.
[21] J. Bang, J., and G. Z. Gutin. “Digraphs: theory, algorithms and applications.” Springer Science &
Business Media, 2008.
[22] Y. W. Teh and M. Welling, “Belief optimization for binary networks: a stable alternative to loopy belief
propagation,” Proceedings of the Eighteenth conference on Uncertainty in artificial intelligence, 493-500,
2001.
[23] A. L. Yuille, “CCCP algorithms to minimize the Bethe and Kikuchi free energies: Convergent alternatives
to belief propagation,” Neural Computation, 14(7): 1691-1722, 2002.
[24] J. Shin, “The complexity of approximating a Bethe equilibrium,” Information Theory, IEEE Transactions
on, 60(7): 3959-3969, 2014.
[25] https://www.quora.com/Statistical-Mechanics-What-is-the-fermion-sign-problem
[26] Dyer, M., Frieze, A., and Jerrum, M. “On counting independent sets in sparse graphs,” SIAM Journal on
Computing 31(5): 1527-1541, 2002.
[27] J. Schweinsberg, “An O(n2 ) bound for the relaxation time of a Markov chain on cladograms.” Random
Structures & Algorithms 20(1): 59-70, 2002.

10

Appendix A. Transformation to an equivalent binary pairwise model with maximum
degree at most 3

Figure 2: Demonstration of building an equivalent model with maximum degree ∆ ≤ 3 via ‘expanding’
vertices (in grey). In the new model, one can introduce edge factor ψu,v between the duplicated vertices u, v
(in bold) such that ψu,v (xu , xv ) = 1 if xu = xv and ψu,v (xu , xv ) = 0 otherwise.

Appendix B. Proof of Theorem 1
First, note that the MC induced by the worm algorithm converges to the following stationary distribution
Y
πWA (F ) ∝ Ψ(F )
w(e),
e∈F

where
(
Ψ(F ) =

n,
2,

∀F ∈ L2-Loop ,
∀F ∈ L2-Odd .

We first prove its polynomial mixing, i.e. it produces a sample from a distribution with the desired total
variation distance from πWA in a polynomial number of iterations.
Lemma 1. Given any δ > 0 and any F0 ∈ L2-Loop ∪ L2-Loop , choose
Tmix ≥ w(F0 )−1 + (m − n + 1) log 2 + 12∆mn4 log δ −1 ,
t
and let πWA
(·) denote the resulting distribution of after updating t times by the worm algorithm with initial
state F0 . Then, it follows that

1
2

X

Tmix
πWA
(F ) − πWA (F ) ≤ δ,

F ∈L2-Loop ∪L2-Loop

namely, the mixing time of the MC is bounded above by Tmix .
The proof of the above lemma is given in Section B.1. Collevecchio et al. [18] recently proved that the
worm algorithm mixes in polynomial time when the weights are uniform, i.e., equal. We extend the result to
our case of non-uniform weights. The proof is based on the method of canonical path, which views the state
space as a graph and constructs a path between every pair of states having certain amount of flow defined by
πWA . From Lemma 1 with parameters
N ≤ 1.2n log(3δ −1 ),

T ≤ (m − n + 1) log 2 + 4∆mn4 log(3nδ −1 ),

and

F0 ← ∅,

we obtain that the total variation distance between πWA and the distribution of updated states in line 4 of
δ
Algorithm 1 is at most 3n
. Next, we prove that the probability of acceptance in line 6 of Algorithm 1 is
sufficiently large.

11

Lemma 2. The probability of sampling a 2-regular loop from distribution πW A is bounded below by n−1 , i.e.
πWA (L2-Loop ) ≥ n1 .
The proof of the above lemma is given in Section B.2. The proof relies on the fact that the size of L2-Loop
is bounded by a polynomial of the size of L2-Odd .
Now we are ready to complete the proof of Theorem 1. Let π
b2-Loop denote the distribution of 2-regular
loops from line 6 of Algorithm 1 under parameters as in Theorem 1. We say Algorithm 1 fails if it outputs
F = ∅ from line 9. Choose a set of 2-regular loops Lb2-Loop := {F ∈ L2-Loop : π
b2-Loop (F ) > π2-Loop (F )}.
Then the total variation distance between π2-Loop and π
b2-Loop can be expressed as:
1
2

X

|b
π2-Loop (F ) − π2-Loop (F )| = π
b2-Loop (Lb2-Loop ) − π2-Loop (Lb2-Loop ).

F ∈L2-Loop

By applying Lemma 1 and Lemma 2, we obtain the following under parameters as in Theorem 1:
π
b2-Loop (Lb2-Loop ) − π2-Loop (Lb2-Loop )
(a) π
bWA (Lb2-Loop )
− (1 − π
bWA (L2-Loop ))N − π2-Loop (Lb2-Loop )
≥
π
bWA (L2-Loop )
b2-Loop ) + δ
(b) πWA (L
δ N
3n
− (1 − πWA (L2-Loop ) −
≥
) − π2-Loop (Lb2-Loop )
δ
3n
πWA (L2-Loop ) − 3n
(c)

≥ −
(d)

≥ −

δ
2δ
− e−(πWA (L2-Loop )+ 3n )N
3n πWA (L2-Loop )

δ
2δ
− = −δ.
3
3

In the above, (a) comes from the fact that a sample from line 6 of Algorithm 1 follows the distribution
b2-Loop )
π
bWA (L
π
bWA (L2-Loop )

and the failure probability of Algorithm 1 is (1 − π
bWA (L2-Loop ))N . For (b), we use the variation
distance between π
bWA and πWA due to Lemma 1 and parameters as in Theorem 1, i.e.,
|b
πWA (S) − πWA (S)| ≤

δ
3n

∀ S ⊆ L2-Loop ∪ L2-Odd .

For (c), we use (1 − x) ≤ e−x for any x ≥ 0 and (d) follows from Lemma 2 and N ≤ n ln(3δ −1 ). The
converse π
b2-Loop (Lb2-Loop ) − π2-Loop (Lb2-Loop ) ≤ δ can be done similarly by considering the complementary set
b
L2-Loop \L2-Loop . This completes the proof of Theorem 1.
B.1 Proof of Lemma 1
First, let PWA denote the transition matrix of MC induced by the worm algorithm in Section 3.1. Then we are
able to define the corresponding transition graph GWA = (L2-Loop ∪ L2-Odd , EWA ), where each vertex is a state
of the MC, and edges are defined on state pairs with nonzero transition probability, i.e.
EWA = {(A, A0 ) : (A, A0 ) ∈ (L2-Loop ∪ L2-Odd ) × (L2-Loop ∪ L2-Odd ), PπWA (A, A0 ) > 0}.
Our proof makes use of the following result proved in [18].
Theorem 4 (Schweinsberg 2002 [18]). Consider an irreducible and lazy MC, with finite state space Ω,
transition matrix P and transition graph GP , which is reversible with respect to the distribution π. Let O ⊆ Ω
be nonempty, and for each pair (I, J) ∈ Ω × O, specify a path γI,J in GP from I to J. Let
Γ = {γI,J : (I, J) ∈ Ω × O}
12

denote the collection of all such paths, and let L(Γ) be the length of longest path in Γ. For any transition
T ∈ EP , let
HT = {(I, F ) ∈ Ω × O : T ∈ γI,J }.
Then

"
τA (δ) ≤ log

where
Φ(Γ) =

max

(A,A0 )∈EP

1
π(A) + log



X



I,J∈H(A,A0 )

!#
1
δ



4L(T )Φ(Γ)



π(I)π(J)
.
π(O)π(A)P (A, A0 ) 

To this end, we choose O = L2-Loop and we show that there exists a choice of paths Γ = {γI,J : (I, J) ∈
(L2-Loop ∪ L2-Loop ) × L2-Odd } such that
Φ(Γ) ≤ ∆n4 ,

L(Γ) ≤ m.

Then we obtain the statement in Lemma 1 immediately.
We begin by specifying Γ, and then proceed to the bound of Φ(Γ). To this end, we fix an [n]-valued vertex
labeling of GWA . The labeling induces a lexicographical total order of the edges, which in turn induces a
lexicographical total order on the set of all subgraphs of GWA . In order for the state I ∈ L2-Loop ∪ L2-Odd transit
to the J ∈ L2-Loop , it suffices that it updates, precisely once, those edges in I ⊕ J. In order to describe such
path, we first prove that there exist a injection from I ⊕ J to some unique disjoint partition I ⊕ J = ∪ki=0 Ci ,
where C0 is either a path or a cycle and C1 , · · · , Ck are cycles. Observe that since J ∈ L2-Loop , applying
symmetric difference with J does not change the parity of degrees of the vertices and I ⊕ J ∈ L2-Loop ∪ L2-Odd .
First consider the case when I ⊕ J ∈ L2-Odd . Then there exist a path between two odd-degree vertices in
I ⊕ J, since the sum of degrees over all vertices in a component is even. Among such paths, we pick C0 as the
path with the highest order according to the [n]-valued vertex labeling. Now observe that I ⊕ J\C0 ∈ L2-Loop
is Eulerian, which can be decomposed into disjoint set of cycles. We are able to choose a C1 , · · · , Ck uniquely
by recursively excluding a cycle with the highest order, i.e. we pick C1 as a cycle with highest order from
I ⊕ J\C0 , then pick C2 from I ⊕ J\C0 \C1 with the highes order, and so on. For the case when I ∈ L2-Loop ,
I ⊕ J ∈ L2-Loop is Eulerian and we can apply similar logic to obtain the unique decomposition into disjoint
cycles.
Now we are ready to describe γI,J , which updates the edges in I ⊕ J from C0 to Ck in order. If C0 is a
path, pick an endpoint with higher order of label and update the edges in the paths by it unwinding the edges
along the path until other endpoint is met. In the case of cycles, pick a vertex with highest order of label and
unwind the edges by a fixed orientation. Note that during the update of cycles, the number of odd-degree
vertices are at most 2, so the intermediate states are stil in L2-Loop ∪ L2-Odd . As a result, we have constructed a
path γI,F for each I ∈ L2-Loop ∪ L2-Odd and J ∈ L2-Loop where each edge correspond to an update on I ⊕ J
and |γI,F | = |I ⊕ J| ≤ m.
Next, we bound the corresponding Φ(Γ). First let L4-Odd denote the set of subgraphs with exactly 4
odd-degree vertices. We define a mapping ηT : HT → L2-Loop ∪ L2-Odd ∪ L4-Odd by the following:
ηT (I, J) := I ⊕ F ⊕ (A ∪ e),
where T = (A, A ⊕ e). Observe that ηT (I, J) agrees with I on the components that have already been
processed, and with J on the components that have not. We prove that ηT is an injection by reconstructing I
and J from ηT (I, J) given T = (A, A ⊕ e). To this end, observe that I ⊕ F = ηT (I, F ) ⊕ (A ∪ e) is uniquely
decided from ηT (I, F ) and (A ∪ e). Then given I ⊕ F , we are able to infer the decomposition C0 , C1 , · · · , Ck
of I ⊕ J by the rules defined previously. Moreover the updated edge e implies the current set Ci being updated.
Therefore we can infer the processed part of I ⊕ J. Then we can recover J by beginning in A and unwinding
the remaining edges in I ⊕ J that was not processed yet. Then we recover I via I = ηT (I, J) ⊕ (A ∪ e) ⊕ J
and therefore ηT is injective.
13

Next, we define a metric wWA such that given an edge set F ,
Y
wWA (F ) :=
|w(e)|.
e∈F

We complete the proof by showing that for any T = (A, A0 ) ∈ E, the following inequality holds:
(a)

X

Φ(Γ) ≤

I,J∈HT

π(I)π(J) (b)
≤
π(L2-Loop ) π(A)P (A, A0 )
1

X
I,J∈HT

(c)
2∆
Ψ(I)wWA (ηT (I, J)) ≤ ∆n4 .
wWA (L2-Loop )

First, (a) holds by definition of Φ. We prove (b) by the following chain of inequality:
1
π(I)π(J)
1
Ψ(I)wWA (I)nwWA (J)
=
π(L2-Loop ) π(A)P (A, A0 )
nwWA (L2-Loop ) Ψ(A)wWA (A)PWA (A, A0 )
(1)

2∆
1
Ψ(I)wWA (I)wWA (J)
wWA (L2-Loop )
wWA (A ∪ e)
2∆
(2)
=
Ψ(I)wWA (ηT (I, F )).
wWA (L2-Loop )
≤

In the above, (1) comes from the definition of the transition probability and (2) comes from the definition of
function wWA . Finally, we prove (c). First, we have
X

Ψ(Γ) ≤

(I,J)∈HT

2∆
Ψ(I)wWA (ηT (I, F ))
wWA (L2-Loop )

2∆
[wWA (L2-Loop ∪ L2-Odd ) + 2wWA (L2-Loop ∪ L2-Odd ∪ L4-Odd )]
wWA (L2-Loop )
(I,J)∈HT


wWA (L2-Odd )
wWA (L4-Odd )
= 2∆ (n + 2) + (n + 2) +
+2
,
wWA (L2-Loop )
wWA (L2-Loop )
≤

X

since ηT (I, J) is an injection on L2-Loop ∪ L2-Odd ∪ L4-Odd , and the set L2-Loop , L2-Odd , L4-Odd are disjoint. Now
we prove
 
 
wWA (L2-Odd )
n
wWA (L4-Odd )
n
≤
≤
,
wWA (L2-Loop )
2
wWA (L2-Loop )
4


4
which completes the proof of Lemma 1 since (n + 2) + (n + 2) + n2 + 2 n4 ≤ n2 . To this end, we let
LOdd (W ) denote the set of generalized loops having W as the set of odd degree vertices. Now observe the
following inequality:
X

(a)

wWA (F ) =

Y
Y
1 X Y
|w(e)|
(1 + (−1)dF (v) )
(1 + (−1)dF (v)+1 )
n
2
F ∈L e∈F

F ∈LOdd (W )

1
= n
2
=
(b)

≥

X

|w(e)|

σ∈{−1,1}V F ∈L e∈F

X

Y

σ∈{−1,+1}V

e=(u,v)∈E

X

=

X

Y

Y

σvdF (v)

s∈V

v∈W

(1 + |w(e)|σu σv )

Y

Y

v∈W

(1 + |w(e)|σu σv )

σ∈{−1,+1}V e=(u,v)∈E
(c)

s∈W

s∈V \W

X Y

wWA (F ).

F ∈LL2-Loop

14

σv
σv

In the above, (a) comes from the fact that 1 + (−1)dv (F ) = 2 if dv (F ) is even and 0 otherwise, so only the
terms corresponding to 2-regular loop becomes non-zero. For (b), the inequality comes from the fact that
1 + |w(e)|σ
u σv ≥ 0 and σ
v ≤ 1. For (c), the equality is from the fact that L2-Loop = LOdd (∅). Therefore we
P
P
have F ∈L(∅) |w(F )| ≥ F ∈L(W ) |w(F )|, leading to
wWA (L2-Odd )
=
wWA (L2-Loop )

P

W ⊆V,|W |=2

P

F ∈LOdd (W )

|wWA (F )|

 
n
≤
,
2

wWA (L2-Loop )

and the case for L4-Odd is done similarly. This completes the proof of Lemma 1.
B.2 Proof of Lemma 2
Given W ⊆ V , we let LOdd (W ) denote the set of generalized loops having W as the set of odd degree vertices.
where Odd(F ) is the set of odd-degree vertices in F . Now observe the following inequality:
X

Y
Y
1 X Y
|w(e)|
(1 + (−1)dF (v) )
(1 + (−1)dF (v)+1 )
n
2

(a)

wWA (F ) =

F ∈L e∈F

F ∈LOdd (W )

=

1
2n

X

Y

X

Y

|w(e)|

σ∈{−1,1}V F ∈L e∈F

=

s∈W

s∈V \W

X Y

σvdF (v)

s∈V

v∈W

(1 + |w(e)|σu σv )

Y

X

≥

Y

σv
σv

v∈W

σ∈{−1,+1}V e=(u,v)∈E
(b)

Y

(1 + |w(e)|σu σv )

σ∈{−1,+1}V e=(u,v)∈E
(c)

X

=

wWA (F ).

F ∈LL2-Loop

In the above, (a) comes from the fact that 1 + (−1)dv (F ) = 2 if dv (F ) is even and 0 otherwise, so only the
terms corresponding to 2-regular loop becomes non-zero. For (b), the inequality comes from the fact that
1 + |w(e)|σ
u σv ≥ 0 and σ
v ≤ 1. For (c), the equality is from the fact that L2-Loop = LOdd (∅). Therefore we
P
P
have F ∈L(∅) |w(F )| ≥ F ∈L(W ) |w(F )|, leading to
P

F ∈L2-Loop

P

πWA (F )

F ∈L2-Loop ∪L2-Odd

πWA (F )

n
=
≥

n

P

F ∈L2-Loop

P

|wWA (F )| +

F ∈L2-Loop

P

|wWA (F )|
P

WWA ⊆V,|W |=2

F ∈LOdd (W )

|wWA (F )|

n
1

,
n =
n
n+2 2

which completes the proof of Lemma 2.

Appendix C. Proof of Theorem 2
First, we quantify how much samples from Algorithm 1 are necessary for estimating some non-negative real
valued function f on L2-Loop . To this, we state the following lemma which is a straightforward application of
the known result in [8].

15

Lemma 3. Let f be a non-negative real-valued function defined on L2-Loop and bounded above by fmax ≥ 0.
Given 0 < ξ ≤ 1 and 0 < η ≤ 1/2, choose
504ξ −2 dlog η −1 efmax
Eπ2-Loop [f ]

s≥

N ≥ 1.2n log

T ≥ (m − n + 1) log 2 + 4∆mn4 log

24fmax
,
ξEπ2-Loop [f ]

8fmax
,
ξEπ2-Loop [f ]

and generate 2-regular loops F1 , F2 , · · · Fs using Algorithm 1 with inputs N and T . Then, it follows that

 1P
| s i |w(Fi )| − Eπ2-Loop (f )|
≤ ξ ≤ 1 − η.
P
Eπ2-Loop (f )
namely, samples of Algorithm 1 estimates Eπ2-Loop (f ) within approximation ratio 1 ± ξ with probability at least
1 − η.
First, recall that during each stage of simulated annealing, we approximate the expectation of the function
w(F )1/n with respect to the distribution π2-Loop (β), i.e.,
h
i
†
†
Eπ2-Loop (β) |w(F )|1/n = Z2-Loop
(βi+1 )/Z2-Loop
(βi ).


Hence, to apply Lemma 3, we bound maxF |w(F )|1/n and Eπ2-Loop (β) |w(F )|1/n as follows:
i
h
Eπ2-Loop (β) |w(F )|1/n ≥ wmin ,

|w(F )|1/n ≤ 1

where the first inequality is due to w(e) ≤ 1 for any e ∈ E and the second one is from |F | ≤ n for any
2-regular loop F . Thus, from Lemma 3 with parameters
−1
s ≥ 18144n2 ε−2 wmin
dlog(6nν −1 )e,

−1
N ≥ 1.2n log(144nε−1 wmin
),

−1
T ≥ (m − n + 1) log 2 + 4∆mn4 log(48nε−1 wmin
),

on each stage, we obtain

P

This implies that the product

†
†
|Hi − Z2-Loop
(βi+1 )/Z2-Loop
(βi )|
†
†
Z2-Loop
(βi+1 )/Z2-Loop
(βi )

Q

i Hi estimates

†
Z2-Loop
m−n+1
2


ε
ν
≤
≥1−
.
6n
6n

within approximation ratio in

[((1 − ε/6n)n , (1 + ε/6n)n ] ⊆ [1 − ε/3, 1 + ε/3]
with probability at least (1 − ν/6n)n ≥ 1 − ν/3, i.e.,

P

|2m−n+1

Q

i Hi −
†
Z2-Loop

†
Z2-Loop
|


ε
ν
≤
≥1− .
3
3

Next we define a non-negative real-valued random function g on L2-Loop as
(
1
if w(F ) < 0
g(F ) =
,
0
otherwise

16

namely, Eπ2-Loop [g(F )] = Pπ2-Loop [w(F ) < 0]. Since maxF g(F ) = 1, one can apply Lemma 3 with parameters
s ≥ 18144ζ(1 − 2ζ)−2 ε−2 dlog(3ν −1 )e,

N ≥ 1.2n log(144ε−1 (1 − 2ζ)−1 ),

T ≥ (m − n + 1) log 2 + 4∆mn4 log(48ε−1 (1 − 2ζ)−1 )
and have


P


|κ − Pπ2-Loop [w(F ) < 0]|
(1 − 2Pπ2-Loop [w(F ) < 0])ε
ν
≤
≥1− ,
Pπ2-Loop [w(F ) < 0]
6Pπ2-Loop [w(F ) < 0]
3

since ζ = Pπ2-Loop [w(F ) < 0]. Furthermore, after some algebraic calculations, one can obtain

|(1 − 2κ) − (1 − 2Pπ2-Loop [w(F ) < 0])|
ε
ν
≤
≥1− .
P
1 − 2Pπ2-Loop [w(F ) < 0]
3
3


†
The rest of the proof is straightforward since we estimate Z2-Loop = (1 − 2Pπ2-Loop [w(F ) < 0])Z2-Loop
by
Q
m−n+1
2
2
(1 − 2κ)2
i Hi , the approximation ratio is in [(1 − ε/3) , (1 + ε/3) ] ⊆ [1 − ε, 1 + ε] with probability
at least (1 − ν/3)2 ≥ 1 − ν.

Appendix D. Proof of Theorem 3
Given F ∈ L, we let the odd-degree vertices in F (i.e., dF (·) is odd) by v1 , v2 , · · · v2` for some integer ` ≥ 0.
Since we assume G is connected, there exist a set of paths P1 , P2 , · · · P` such that Pi is a path from v2i−1 to
v2i . Note that given any set of edges D ⊆ E, D ⊕ Pi changes the parities of dD (v2i−1 ), dD (v2i ), while others
remain same. Therefore, all degrees in F ⊕ P1 ⊕ · · · ⊕ P` become even. Then, due to the definition of cycle
basis, there exist some C1 , C2 , · · · Ck ∈ C such that
C1 ⊕ C2 · · · ⊕ Ck = F ⊕ P1 ⊕ · · · ⊕ P` ,
namely,
F = C1 ⊕ C2 · · · ⊕ Ck ⊕ P1 ⊕ · · · ⊕ P` .
This completes the proof of Theorem 3.

17

