Probabilistic Duality for Parallel Gibbs Sampling
without Graph Coloring
Lars Mescheder1 , Sebastian Nowozin2 , and Andreas Geiger1

arXiv:1611.06684v1 [cs.LG] 21 Nov 2016

1

Autonomous Vision Group, MPI Tübingen2
Microsoft Research, Cambridge
November 23, 2016

Abstract
We present a new notion of probabilistic duality for random variables involving mixture distributions. Using this notion, we show how
to implement a highly-parallelizable Gibbs sampler for weakly coupleddiscrete pairwise graphical models with strictly positive factors that
requires almost no preprocessing and is easy to implement. Moreover,
we show how our method can be combined with blocking to improvemixing. Even though our method leads to inferior mixing times compared to a sequential Gibbs sampler, we argue that our method is still
very useful for large dynamic networks, where factors are added andremoved on a continuous basis, as it is hard to maintain a graph coloring in this setup. Similarly, our method is useful for parallelizing Gibbs
sampling in graphical models that do not allow for graph colorings witha small number of colors such as densely connected graphs.

1

Introduction

Inference in general discrete graphical models is hard. Besides variational
methods, the main approach for inference in such models is through runninga Markov chain that in the limit draws samples from the true posterior
distribution.
One such Markov-Chain is given by the so-called Gibbs-sampler thatwas first introduced introduced by Geman and Geman [4]. In each step,
one random variable is resampled given all the others, using the conditional
probability distributions in the graphical model. Under mild hypotheses, theGibbs sampler produces an ergodic Markov chain that converges to the target
distribution. The main appeal of the Gibbs sampler lies in its simplicity andease of implementation. Unfortunately, for highly coupled random variables,
mixing of this Markov chain can be prohibitively slow. Moreover, in order1

to achieve ergodicity, we have to sample one variable after another, yielding
an inherently sequential algorithm. However, with the advent of affordable
parallel computing hardware in the form of GPUs, it is desirable to havea parallel sampling algorithm. Early attempts tried running all the update
steps in parallel [4], however this update schedule usually does not converge
to the target distribution [8].
Another common approach to parallel Gibbs-sampling is to compute agraph coloring of the underlying graph and then perform Gibbs-sampling
blockwise [5]. However, it is not always straightforward to find an appropriate graph coloring and it is hard to maintain such a graph coloring in a
dynamic setup, i.e. when factors are added and removed on a continuous
basis.
In this paper, we show a simple method of parallelizing a Gibbs-samplerthat does not require a graph coloring. This is particularly useful for situations in which a graph coloring is hard to obtain or the graph topology
changes frequently, which requires to maintain or to recompute the graphcoloring.

2

Related WorkAn early attempt at parallelization for Ising models was described by Swendsen and Wang [12]. Their method was generalized to arbitrary probabilistic graphical models in [1]. However, while the Swendsen Wang-algorithm
mixes fast for the Ising model with no unary potentials, this needs not be thecase for general probabilistic graphical models [6, 7]. Higdon [6] presents amethod for performing partial Swendsen Wang-updates. However, sampling
using Higdon’s method requires sampling from a coarser graphical model
as a subproblem, which Higdon tackles using conventional sampling methods. Our dualization strategy allows to circumvent this step, so that only
standard clusterwise sampling as in [12] is required.
[5] describes two ways of parallelizing Gibbs sampling in discrete Markov
random fields. The first method relies on computing graph colorings, thesecond on decomposing the graph into blocks (called splashes) consistingof subgraphs with limited tree width. Both methods are complimentary inthat graph colorings work well for loosely-coupled graphical models whereas
splash sampling works well in the strongly coupled case. However, computing a minimal graph coloring is a NP-hard problem [3] and the number ofcolorings necessary depends on the graph. Moreover, it is hard to maintaina graph coloring in a dynamic setting in which the graph topology is notconstant anymore. Our approach does not suffer from these issues and requires almost no preprocessing. Moreover, our approach can be combinedwith splash sampling. Whereas the approach in [5] requires the splashes tobe induced subgraphs of the graphical model, our approach allows to select2arbitrary subgraphs of the graphical model as splashes, making it possible
to use splashes containing many variables.
Schmidt et al. show in [10] how a Gaussian scale mixtures (GSMs) canbe used for efficiently sampling from fields of experts [9]. Our approach isvery similar to theirs, but deals with discrete graphical models. Moreover,whereas Schmidt et al. started with a model in a primal-dual formulationand trained it on data, we focus on decomposing existing graphical models
and mainly use duality for inference. In fact, both techniques can be subsumed in the framework of exponential family harmoniums [13], which makes
it possible to deal with models consisting of both discrete and continuous
random variables.
Martens et al. [7] show how to sample from a discrete graphical modelusing auxiliary variables. Their approach is similar to ours, but relies on
computing the (sparse) Cholesky decomposition of a large matrix beforehand. Our approach does not have this issue.
A variant of our algorithm that computes expectations instead of performing sampling corresponds to the mean-field-algorithm for junction treeapproximations in [14]. We can show that our algorithm minimizes an upperbound to the true mean-field objective. However, whereas the mean-fieldalgorithm in [14] has to recalibrate the tree for each update of a potential,
our algorithm updates all the potentials at once with only one run of thejunction-tree-algorithm.
Schwing et al. [11] designed a system that performs belief propagation
in a distributed way. Our approach has a similar goal but follows a different
strategy: while Schwing et al. achieve parallelism through a convex formulation, we augment the probabilistic model with additional random variables.

3

Probabilistic duality

We first define the notion of duality of random variables. Note that thisdefinition is very similar to the Lagrange functional of a convex optimization
problem in convex analysis.
Definition 1. Let x ∈ Ωx and θ ∈ Ωθ denote random variables. We callfunctions s : Ωx → V and r : Ωθ → V to a common vector space with somebilinear form link functions. We say x and θ are dual to each other via (s, r)
if the joint distributions can be written asp(x, θ) = h(x)g(θ)ehs(x),r(θ)iwith positive functions h : Ωx → R and g : Ωθ → R.
In the language of [13], a dual pair of random-variables is simply an
exponential family harmonium. For a pair of link functions (s, r) and h :Ωx → R, g : Ωθ → R some real valued functions we now define the conceptof an (s, r)-transform.
3

Definition 2. The (s, r)-transforms of h : Ωx → R, g : Ωθ → R is definedasX
h(x)ehs(x),r(θ)i
H(θ) :=
x

G(x) :=

Xg(θ)ehs(x),r(θ)i .

θNote the resemblance to the notion convex conjugacy. The following
simple lemma is central to the rest of the theory:
Lemma 1. Let be x and θ two dually related random variables with joint
probability density as above. Thenp(x) = h(x)G(x)p(θ) = H(θ)g(θ)h(x) hs(x),r(θ)ip(x | θ) =eH(θ)g(θ) hs(x),r(θ)ip(θ | x) =e.
G(x)Formally, this is similar to the notion of duality in convex optimization,
we call the problem of sampling from p(x) the primal and the corresponding
problem of sampling from p(θ) the dual sampling problem. The primal anddual problems are linked to each other via the conditional densities p(x | θ)and p(θ | x). Note that both p(x | θ) and p(θ | x) are in the exponential
family.
Another view is that we represent p(x) as a mixture of probability distributions p(x | θ) from a specified exponential family determined by s(x) andh(x). The density of the mixture parameters is then given by p(θ).
To obtain a dual formulation of a sampling problem, we have to decompose p(x) asXp(x) = h(x)G(x) = h(x)g(θ)ehs(x),r(θ)iθwith some functions g(θ) and h(x). p(x) is then the marginal ofp(x, θ) = h(x)g(θ)ehs(x),r(θ)i .
We show how this can be done for the discrete case in Section 4.1.
Further evidence can be incorporated by replacing s(x) and h(x) withs̃(x) = s(x, xe ) and h̃(x, xe ) = h(x, xe ).
The lemma already shows a possible strategy how to sample from p(x)using a simple Gibbs sampler: first sample from p(x | θ), then from p(θ | x)4

and so on. As both x and θ are generally high-dimensional, sampling fromp(x | θ) and from p(θ | x) could potentially cause problems. However, weshow that p(x | θ) and p(θ | x) both factorize in Markov random fields forappropriate choices s and r which yields algorithms that are very easy toparallelize.

4

Duality in MRFs

We now take a closer look at sampling in Markov random fields. The
Hammersley-Clifford theorem states that the joint probability density ofnodes in the MRF can be decomposed asN

1 Y
pi (x)
Zp(x) =

i=1

with some probability measures pi and normalization constant Z. In
general the pi only depend on a small subset of the components of x (e.g.
the cliques in the MRF). Assume that for every pi we have a random variableθi , so that x and θi are dual to each other via (s, ri ) with respect to pi . We
can then writepi (x, θi ) = hi (x)gi (θi )ehs(x),ri (θi )i .
P
Q
Theorem 1.QLet θ := (θ1 , . . . , θN ) and r(θ) := i ri (θi ), h(x) := i hi (x)and g(θ) := i gi (θi ).
Then p(x) is the marginal ofp(x, θ) ∝ h(x)g(θ)ehs(x),r(θ)i .
Thus if p(x, θ) ∝ h(x)g(θ)ehs(x),r(θ)i , x and θ are dual to each other. The
marginal distribution of θ is given byYp(θ) ∝ H(θ)
gi (θ).
i

Proof. We havep(x) ∝Yhi (x)Gi (x) =iYhi (x)

i

=
=

X

θ1

···

X Y
θN

hi (x)#"
Yih(x)g(θ)egi (θi )ehs(x),ri (θi )i 

θi

"
X


X

i

hs(x),r(θ)i

θ

5

.

#
gi (θi ) ehs(x),Pi ri (θi )iCorollary 1. p(x | θ) and p(θ | x) are given byPYp(x | θ) ∝hi (x)ehs(x), i ri (θi )iip(θ | x) ∝Ygi (θi )ehs(x),Pi ri (θi )i.

i

In particular, p(θ | x) factorizes and if hi ≡ 1 and every component ofs(x) only depends on one component of x, p(x | θ) factorizes as well. In
particular, this is true for the standard choice s(x) = x.

4.1

Binary pairwise MRFs

We now want to turn to the special case of a pairwise binary MRF. It turns
out that finding a dual representation is equivalent to finding an appropriate
factorization of the probability table.
Theorem 2. Let P be proportional to the probability table of two binary
random variables x1 and x2 . Assume we are given a factorization P = BC |with B, C ∈ R2×2 , where both B and C have strictly positive entries. Let
B2,1
B1,1
C2,1
α2 = log
C1,1
B1,2 C1,2
q = log
B1,1 C1,1
B2,2 B1,1
β1 = log
B1,2 B2,1
C2,2 C1,1
β2 = log.
C1,2 C2,1

α1 = logThenp(x1 , x2 ) ∝X

h(x)g(θ)ehx,r(θ)iθ∈{0,1}with
r(θ) = θ

β1
β2


.

h(x) = eα1 x1 eα2 x2g(θ)= eqθ .
Proof. This follows fromP =X  B1,i 
i=1,2

B2,i

after a simple calculation.
6

C1,i C2,iWe will now show how to find such a factorization:
Lemma 2. If P is symmetric and det P ≥ 0, P can be factored in the form
P = B B | , where
 √



√
π 1
p12
p11 cos(ϕ)
p11 sin(ϕ)√√
B=.
with ϕ = − arccos √p22 sin(ϕ)
p22 cos(ϕ)
4 2p11 p22
Proof. B is well defined because 0 ≤ detP = p11 p22 − p212 and positive
because ϕ ∈ 0, π4 . Now letB=
We have
B B| =



b̃|1
b̃|2


.

kb̃1 k2 hb̃1 , b̃2 ihb̃1 , b̃2 i kb̃2 k2

.

Due to trigonometric considerations


√
π

p11
c p11 p22|√
with c = cos− 2ϕ .
BB =c p11 p22
p222This shows B B | = P as required.
Remark 1. For ϕ = π4 − 21 arccos(a), we have√1 √cos(ϕ) = ( 1 + a + 1 − a)2√
1 √sin(ϕ) = ( 1 + a − 1 − a).
2Lemma 3. For any P , then


0
p−1
12
0 p−1
21


P

is symmetric.
Lemma 4. If det Pi < 0 , then


0 1
1 0


P

has positive determinant.
In summary, we have found a strictly positive factorization of any strictly
positive 2 × 2 matrix. Together with Theorem 1 this yields a dual representation for every binary pairwise MRF.
7

4.2

General discrete MRFs

When the variables are allowed to have multiple states, we can convert any
discrete pairwise MRF into a binary MRF using 0−1-encoding and additional
hard-constraints that ensure that exactly one binary variable belonging to arandom variable in the original MRF has value 1. All inference algorithms
in this paper therefore generalize to this situation.
Dualizing a n × m factor in this way introduces nm auxiliary binary
random variables to the model. Note however, that no new random variables need to be introduced for 1-entries in the factor. For example, for a
Potts-factor of order n, only n auxiliary binary random variables have to be
introduced per factor.
Arbitrary discrete MRFs with higher-order factors work as well, as longas we can find an appropriate positive tensor factorization of the probability
table. Moreover, it is also possible to perform inference approximately byfitting a mixture of Bernoulli / mixture of Dirichlet distributions to thefactors using expectation maximization.

4.3

Swendsen-Wang and local constraints

As it turns out, the Swendsen-Wang algorithm can be seen as a degenerate
special case of our formalism for a particular choice of s(x). Moreover,more general local constraint models can potentially be derived from thisformalism.
More explicitly, for the Ising model lets(x) := (−I(xe1 = xe2 ))e∈E
where E denotes the set of edges and(
0 if xe1 = xe2
I(xe1 = xe2 ) =.
∞ elseThe Ising factor of the form  −w
e i
1
e−wi
Pi ∝
=
e−wi
e−wi
1

e−wi
e−wi




+

1 − e−wi0

can then be decomposed as
Pi (xe1 , xe2 ) =X

g(θi )e−θi I(xe1 =xe2 ) .

θi ∈{0,1}

whereg(0) = e−wig(1) = 1 − e−wi .
80
1 − e−wiThe primal dual sampling algorithm then proceeds as follows:p(θi | x) ∝ g(θi )e−θi I(xe1 =xe2 )p(x | θi ) ∝ e−Pi θi I(xe1 =xe2 ),which are just the update rules for the Swendsen-Wang algorithm.
The partial Swendsen-Wang method by Higdon [6] can be regarded as adecomposition of the form 
1 − α e−wi
α 0
Pi ∝
+
.
e−wi 1 − α
0 αThis leads to the method described by Higdon, where we are left with sampling from coarser Ising-model. By applying a factorization as in Section
4.1 to the first term enables us to circumvent this step, so that all clusters
can be sampled independently (the latent variables θ then have 3 different
states).
Similarly, the generalization of the Swendson-Wang algorithm in [1] canbe regarded as a multiplicative decomposition of the form

 w
e i 1
? P̃i ,
Pi ∝1 ewi
where we used the ?-operator to indicate componentwise multiplication. Applying the decomposition above to the first factor, yields (a variant of) the
method in [1]. We can use our method to further decompose P̃i , allowing to
update all clusters in parallel.

5

Inference

5.1

Sampling

Having a primal-dual representation of p(x) of the form
p(x, θ) ∝ h(x)g(θ)ehs(x),r(θ)i
we can sample from p(x, θ) (and thereby from p(x)) by blockwise Gibbssampling, i.e.
x(t+1) ∼ p(x | θ(t) ) ∝ h(x)ehs(x),r(θθ(t+1) ∼ p(θ | x(t+1) ) ∝ g(θ)ehs(x(t) )i(t+1) ),r(θ)i.

For discrete pairwise MRFs with the dual representations as described above,
both distributions factor, so that sampling can be done in parallel, e.g. on
the GPU. Effectively, we have converted our model to a restricted Boltzmann
machine.
95.2

Estimation of the log-partition-function

The logarithm normalization constant of an unnormalized probability distribution, the so-called log-partition-function, is important for model-selection
and related tasks. In this section, we provide a simple estimator for thisquantity that can be used for any dual pair of random variables.
The following defines an unbiased estimator for the partition function Z:
V (x, θ) =p̃(x)p̃(θ)p(x)p(θ)=Z,
p̃(x, θ)p(x, θ)where p̃(x) and p̃(θ) are the unnormalized probability distributions. Indeed,
we haveZp(x)p(θ)
E [V (x, θ)] = Z
p(x, θ)dxdθ = Z.
p(x, θ)Written in terms of G and H, V (x, θ) can be written as
V (x, θ) = G(x)H(θ)e−hs(x),r(θ)i .
Note thatE [− log V (x, θ)] − [− log E [V (x, θ)]] = I(x, θ),
where I(x, θ) is the mutual information between x and θ. This is a measurefor the uncertainty of V (x, θ), as it can be interpeted as a generalized variancefor the convex function x 7→ − log(x).
This also implies that the expectation oflog V (x, θ)yields a lower bound to the log-partition function. In practice, V (x, θ) has
too much variance to be useful. Therefore, we estimate the expectation oflog V (x, θ), which yields a lower bound to the log-partition function.
Example 1. For Q
the Swendsen-Wang-duality for the Ising-model, we haveh(x) = 1, G(x) = e Pe (xe1 , xe2 ) and therefore
XY
H(θ) =e−θe I(xe1 =xe2 ) = 2C(θ) ,xewhere C(θ) is the number of clusters defined by θ. Therefore
log V (x, θ) = log 2 · C(θ) + log p̃(x),where p̃(x) is the unnormalized distribution of the Ising-model.
A natural question to ask is how this estimator relates to the estimate
obtained by running naive mean-field on the primal distribution p(x) only.
The following negative result shows that in most cases the estimate obtainedby mean-field approximations is preferable:
10Lemma 5. We haveI(x, θ) = Eθ KL(p(x | θ), p(x)) ≥ min KL(p(x | ξ), p(x)),ξ(1)

where we define p(x | ξ) ∝ h(x) exp(hs(x), ξi).
Proof. The equality in (1) can be obtained by a straightforward calculation.
The inequality is a simple application of the fact that the expectation over
θ is always bigger than the minimum over θ.
Note, however, that it is not always straightforward to find ξ that minimizes KL(p(x | ξ), p(x)). This is for example the case for the SwendsenWang-representation from Example 1.

5.3

MAP- and mean-field inference

The concept of probabilistic duality that we introduced in Section 3 is also
useful to derive parallel MAP- and mean-field inference algorithms.
By applying EM to p(x, θ) we can also compute local MAP-assignmentsto p(x) in parallel. The updates readx(t+1) = argmaxx h(x)ehs(x),ξ
ξ (t+1) = E r(θ) | x(t+1) .

(t) iSimilar, we can compute mean-field assignments to p(x, θ) using the updatesη (t+1) = E s(x) | ξ (t)ξ (t+1) = E r(θ) | η (t+1) ,where the expectations are taken over the distributions
p(x | ξ) ∝ h(x)ehs(x),ξip(θ | η) ∝ g(θ)ehη,r(θ)i .
Note that these updates have the advantage over ICM and standard naive
mean field that they can again run in parallel and still have convergenceguarantees.
Using this algorithm, we can show that we minimize an upper bound to
the true mean-field objective KL(p(x | ξ), p(x)):
Lemma 6. We havemin KL (p(x | ξ)p(θ, η), p(x, θ)) ≥ KL(p(x | ξ), p(x)).
η

11

(2)

x4
x3
θ4,7
x7

x5
x2

θ1,7
x1

x6

θ1,6

Figure 1: Example of blocking using a tree as a subgraph. Only the dual
variables on the dashed edges are used. Sampling then alternatingly samples all θe variables and all the x variables on the tree. Similarly, the EM
algorithm alternatingly maximizes all x over the tree using max-productbelief-propagation and takes conditional expectations over the θe variables
to update the unary potentials. Our tree-mean-field algorithm does the same,
except that it uses sum-product belief propagation.
Proof. We havep(x, θ)KL(p(x | ξ)p(θ, η), p(x, θ)) = Ex|ξ Eθ|η − logp(x | ξ)p(θ | η)p(x, θ)p(x)≥ Ex|ξ − log Eθ|η= Ex|ξ − logp(x | ξ)p(θ | η)p(x | ξ)= KL(p(x | ξ), p(x)).

Lemma 6 implies that traditional mean-field updates are still preferable
to the ones from our method. Indeed, in practice we found that our method
can lead to poor approximations in presence of many factors. However, it isstill possible to first run our fast parallel algorithm and then fine-tune theresult using traditional mean-field updates.

5.4

Blocking

Gibbs sampling as described in Section 5.1 can still be prohibitively slow inpresence of strongly correlated random variables. Similarly, EM and meanfield updates tend to get stuck in local optima in this situation. This problem
12also occurs for standard Gibbs sampling, ICM and naive mean-field updates.
A common way out is to introduce blocking, e.g. as in [5] for Gibbs sampling.
Unfortunately, blocking is only possible with respect to induced subgraphs for traditional algorithms. As it turns out, our primal-dual decomposition allows to perform blocking with respect to arbitrary subgraphs. Thisis illustrated in Figure 1. The idea is to decompose the dual variables θ into
two subsets θ0 and θ1 , so that p(x, θ0 | θ1 ) is tractable. This is the case,when p(x | θ1 ) is tractable, because thenp(x, θ0 | θ1 ) = p(θ0 | x)p(x | θ1 ).
Note that p(x | θ1 ) is tractable, if the graph obtained by removing all
the factors belonging to θ1 has low tree-width.
(t)
For blocked Gibbs sampling, we sample in each step p(x | θ1 ) and then(t)p(θ1 | x(t+1) ). As a variation of this process, we can vary the decompositionof θ into θ0 and θ1 in each step.
When we perform max-product-belief propagation for x and take expec(t)tations for θ1 we obtain a new inference algorithm for MAP-inference. Notethat in each step, we maximize over all x variables at once. Similarly, we obtain a probabilistic inference algorithm by changing the max-product-belief
propagation step for x by sum-product-belief propagation.
Note also that both the EM, as well as the mean-field algorithm are
guaranteed to increase the objective function in each step.
In this framework, the standard sequential Gibbs sampler can also beinterpreted as a blocked Gibbs sampler, where blocking is performed withrespect to one primal and all neighboring dual variables. Unfortunately,
as blocking generally improves mixing of a Gibbs chain [5], this implies thatthe standard sequential Gibbs sampler has better mixing properties than theparallel primal-dual sampling algorithm. Still, the primal dual formulation
allows for more flexible blocking schemes, potentially making it possible toimprove on the mixing properties of the standard sequential Gibbs sampler
in some situations.

6

Experimental Results

We tested our method on 3 synthetic graphical models. The first model
consists of a 50 × 50-Ising grid and coupling strenghts ranging from β = 0.1
to β = 0.5. Even though the Ising grid is two-colorable and it is therefore
trivial to implement a parallel Gibbs sampler in this setting, this is no longer
possible when the Graph topology is dynamic, i.e. we remove and add factors
from time to time. Maintaining a coloring in this setting is itself a hard
problem. The second model is given by a random graph with N = 1000
variables and F = k N factors, where k ∈ {2, 4, 8, 16, 32, 64}. Both theunitary and pairwise log-potentials were sampled from a normal distribution
13

(a) For the Ising model.

(b) For the fully connected Ising model.

Figure 2: Number of iterations needed to achieve a PSRF below 1.01.
with mean 0 and a standard deviation of 1. The last model consists of a fully
connected Ising model of N = 100 variables and coupling strengths ranging
from β = 0.01 to β = 0.015. Note that for such models, there is an algorithm
that computes the partition function and the marginals in polynomial time[2]. However, this is not longer the case, when the potentials have varying
coupling strengths.
For all these models we compute the potential scale reduction factor(PSRF) for both a sequential Gibbs sampler and our primal-dual-sampler
by running 10 Markov chains in parallel. From the PSRF we compute an
estimate of the mixing time of the Markov chain by taking the first index,
so that the PSRF remains below some specified threshhold afterwards.
The result for the Ising grid is shown in Figure 2a. For both the primaldual sampler and the sequential Gibbs sampler, we plotted the number ofseeps over the whole grid to achieve a PSRF below 1.01. As expected, boththe sequential Gibbs sampler and the primal-dual sampler mix slower as we
increase the coupling strength. Moreover, even though the primal-dual sampler mixes slower than the sequential Gibbs sampler, in our experiments theratio of the mixing times was between 2 and 7 for all the coupling strength.
Therefore, even though a Gibbs sampler based on a two coloring is preferable
in the static setting, our primal-dual sampler becomes a viable alternativein the dynamic setting.
Similar results where obtained for the random graphs. As expected,
mixing of the primal dual sampler became worse as the number of factors
per vertix increases. While our primal-dual-sampler can be an interestingalternative when the factor-to-vertex ratio is low (e.g. k ≈ 2), we do notrecommend our method for models with many more factors than variables ifthese factors are not very weak.
The result for the fully connected Ising-model is shown in Figure 2b.
As there is no coloring available for a fully connected graphical model, we

14compare the number of full sweeps of our primal-dual sampler against the
number of single-site updates of the sequential Gibbs sampler. We see that
our method leads to improved mixing in this setting.

7

Conclusion

We have introduced a new concept of duality to random variales and showed
its usefulness in performing inference in probabilistic graphical models. In
particular, we demonstrated how to obtain a highly-parallizing Gibbs sampler. Even though this parallel Gibbs sampler has inferior mixing properties
compared to the sequential Gibbs sampler, we believe that it can still bevery useful in settings, where a good graph coloring is hard to obtain or
the graph topology changes frequently. Possible extensions of our approach
include good algorithms for selecting appropriate subgraphs for blocking.
Moreover, as primal-dual representations are not unique, we believe thatfurther progress can be made by deriving new decompositions. Another lineof research is to generalize our ideas to higher order factors, both exactly andin approximate ways. We believe that this is possible and allows to apply
the methods in this paper to arbitrary discrete graphical models.

Acknowledgements
This work was supported by Microsoft Research through its PhD Scholarship
Programme.

References[1] Adrian Barbu and Song-Chun Zhu. Generalizing swendsen-wang to sampling arbitrary posterior probabilities. IEEE Transactions on PatternAnalysis and Machine Intelligence, 27(8):1239–1253, 2005.
[2] Boris Flach. A class of random fields on complete graphs with tractablepartition function. IEEE transactions on pattern analysis and machine
intelligence, 35(9):2304–2306, 2013.
[3] Michael R Garey, David S Johnson, and Larry Stockmeyer. Some simplified np-complete problems. In Proceedings of the sixth annual ACM
symposium on Theory of computing, pages 47–63. ACM, 1974.
[4] Stuart Geman and Donald Geman. Stochastic relaxation, gibbs distributions, and the bayesian restoration of images. IEEE Transactions on
pattern analysis and machine intelligence, (6):721–741, 1984.

15

[5] Joseph Gonzalez, Yucheng Low, Arthur Gretton, and Carlos Guestrin.
Parallel gibbs sampling: From colored fields to thin junction trees. In
AISTATS, volume 15, pages 324–332, 2011.
[6] David M Higdon. Auxiliary variable methods for markov chain montecarlo with applications. Journal of the American Statistical Association,
93(442):585–595, 1998.
[7] James Martens and Ilya Sutskever. Parallelizable sampling of markov
random fields. In AISTATS, pages 517–524, 2010.
[8] David Newman, Padhraic Smyth, Max Welling, and Arthur U Asuncion. Distributed inference for latent dirichlet allocation. In Advancesin neural information processing systems, pages 1081–1088, 2007.
[9] Stefan Roth and Michael J Black. Fields of experts. International
Journal of Computer Vision, 82(2):205–229, 2009.
[10] Uwe Schmidt, Qi Gao, and Stefan Roth. A generative perspective on
mrfs in low-level vision. In Computer Vision and Pattern Recognition(CVPR), 2010 IEEE Conference on, pages 1751–1758. IEEE, 2010.
[11] Alexander Schwing, Tamir Hazan, Marc Pollefeys, and Raquel Urtasun.
Distributed message passing for large scale graphical models. In Computer vision and pattern recognition (CVPR), 2011 IEEE conference on,
pages 1833–1840. IEEE, 2011.
[12] Robert H Swendsen and Jian-Sheng Wang. Nonuniversal critical dynamics in monte carlo simulations. Physical review letters, 58(2):86,
1987.
[13] Max Welling, Michal Rosen-Zvi, and Geoffrey E Hinton. Exponential
family harmoniums with an application to information retrieval. In
Advances in neural information processing systems, pages 1481–1488,
2004.
[14] Wim Wiegerinck. Variational approximations between mean field theory
and the junction tree algorithm. In Proceedings of the Sixteenth conference on Uncertainty in artificial intelligence, pages 626–633. Morgan
Kaufmann Publishers Inc., 2000.

16