Ensuring Rapid Mixing and Low Bias for Asynchronous Gibbs Sampling

1 Introduction

Gibbs sampling is one of the most common Markov chain Monte Carlo methods used with graphical models [Koller and Friedman, 2009]. In this setting, Gibbs sampling (Algorithm 1) operates iteratively by choosing at random a variable from the model at each timestep, and updating it by sampling from its conditional distribution given the other variables in the model. Often, it is applied to inference problems, in which we are trying to estimate the marginal probabilities of some query events in a given distribution.

Algorithm 1 Gibbs sampling

Require: Variables xi for 1 ≤ i ≤ n, and distribution π.

for t = 1 to T do

Sample s uniformly from {1, . . . , n}.

Re-sample xs uniformly from Pπ(Xs|X{1,...,n}\{s}).

end for

Gibbs sampling is often applied to sparse graphical models, where each of these updates needs to read the values of only a small subset of the variables; therefore each update can be computed very quickly on modern hardware. Because of this and other useful properties of Gibbs sampling, many systems use Gibbs sampling to perform inference on big data [Newman et al., 2007; Lunn et al., 2009;

McCallum et al., 2009; Smola and Narayanamurthy, 2010; Theis et al., 2012; Zhang and R´e, 2014].

Since Gibbs sampling is such a ubiquitous algorithm, it is important to try to optimize its execution speed on modern hardware. Unfortunately, while modern computer hardware has been trending towards more parallel architectures [Sutter, 2005], traditional Gibbs sampling is an inherently sequential algorithm; that is, the loop in Algorithm 1 is not directly parallelizable. Furthermore, for sparse models, very little work happens within each iteration, meaning it is difficult to extract much parallelism from the body of this loop. Since traditional Gibbs sampling parallelizes so poorly, it is interesting to study variants of Gibbs sampling that can be parallelized. Several such variants have been proposed, including applications to latent Dirichlet allocation [Newman et al., 2007; Smola and Narayanamurthy, 2010] and distributed constraint optimization problems [Nguyen et al., 2013].

In one popular variant, multiple threads run the Gibbs sampling update rule in parallel without locks, a strategy called asynchronous or HOGWILD! execution-in this paper, we use these two terms interchangeably. This idea was previously proposed [Smola and Narayanamurthy, 2010], but not analyzed theoretically, and has been shown to give empirically better results on many models [Zhang and R´e, 2014]. But when can we be sure that HOGWILD! Gibbs sampling will produce accurate results? Except for the case of Gaussian random variables [Johnson et al., 2013], there is no existing analysis by which we can ensure that asynchronous Gibbs sampling will be useful for a particular application. Even the problems posed by HOGWILD!-Gibbs are poorly understood, and their solutions more so.

In the following sections, we will identify two main issues when analyzing asynchronous Gibbs sampling. Firstly, we will show by example that, surprisingly, HOGWILD!-Gibbs can be biased-unlike sequential Gibbs, it does not always produce samples that are arbitrarily close to the target distribution. Secondly, we will show that the mixing time (the time for the chain to become close to its stationary distribution) of asynchronous Gibbs sampling can be exponentially greater than that of the corresponding sequential chain.

To address the issue of bias, we need some way to describe the distance between the target distribution π and the distribution of the samples produced by HOGWILD!-Gibbs. The standard notion to use here is the total variation distance, but for the task of computing marginal probabilities, it gives an overestimate on the error caused by bias. To better describe the bias, we introduce a new notion of statistical distance, the sparse variation distance. While this relaxed notion of statistical distance is interesting in its own right, its main benefit here is that it uses a more local view of the chain to more tightly measure the effect of bias.

Our main goal is to identify conditions under which the bias and mixing time of asynchronous Gibbs can be bounded. One parameter that has been used to great effect in the analysis of Gibbs sampling is the total influence α of a model. The total influence measures the degree to which the marginal distribution of a variable can depend on the values of the other variables in the model-this parameter has appeared as part of a celebrated line of work on Dobrushin's condition (α < 1), which ensures the rapid mixing of spin statistics systems [Dobrushin, 1956; Dyer et al., 2006; Hayes, 2006]. It turns out that we can use this parameter to bound both the bias and mixing time of HOGWILD!-Gibbs, and so we make the following contributions:

• We describe a way to statistically model the asynchronicity in HOGWILD!-Gibbs sampling.
• To bound the bias, we prove that for classes of models with bounded total influence α = O(1), if sequential Gibbs sampling achieves small sparse variation distance to π in O(n) steps, where n is the number of variables, then HOGWILD!-Gibbs samples achieve the same distance in at most O(1) more steps.
• For models that satisfy Dobrushin's condition (that is, α < 1), we show that the mixing time bounds of sequential and HOGWILD!-Gibbs sampling differ only by a factor of 1 +O(n-1 ).

1.1 RelatedWork

Much work has been done on the analysis of parallel Gibbs samplers. One simple way to parallelize Gibbs sampling is to run multiple chains independently in parallel: this heuristic uses parallelism to produce more samples overall, but does not produce accurate samples more quickly. Additionally, this strategy is sometimes worse than other strategies on a systems level [Smola and Narayanamurthy, 2010; Zhang and R´e, 2014] because it requires additional memory.

The HOGWILD!-Gibbs sampling algorithm was inspired by a line of work on parallelizing stochastic gradient descent (SGD) by running it asynchronously [Niu et al., 2011; Liu et al., 2015; De Sa et al., 2015; Mania et al., 2015; Liu and Wright, 2015]. Several of these results suggest modeling the race conditions inherent in HOGWILD! SGD as noise in a stochastic process; in this paper, we will apply a similar stochastic process model to Gibbs sampling.

2 Modeling Asynchronicity

In this section, we describe a statistical model for asynchronous Gibbs sampling by adapting the hardware model outlined in [De Sa et al., 2015]. We will focus on the case where our target distribution π ranges over a discrete space.

Any HOGWILD!-Gibbs implementation involves some number of threads each repeatedly executing the Gibbs update rule on a single copy of the model (typically stored in RAM). We assume that this model serializes all writes, such that we can speak of the state of the system after t writes have occurred. We call this time t, and we will model the HOGWILD! system as a stochastic process indexed by t. Let xi,t denote the value of variable i at time t. For HOGWILD!-Gibbs sampling, the sampler does not get to use exactly the values of xi,t; rather it has access to a cache containing potentially stale values. To model this, we define ˜ vi,t = xi,t-˜ τi,t to be the potentially stale model state used to compute the update at time t, where ˜ τi,t ≥ 0 is a delay parameter that represents how old the currently-cached value for variable i could be. The distribution of these delays ˜ τi,t depends on the number of threads and the specifics of the hardware [Niu et al., 2011]. To make our results general, we will not assume any particular distribution on the delays. Instead, we require only a bound on the expected delay, E [˜ τi,t|Ft] ≤ τ for a hardware-dependent constant τ , and a very weak bound on its magnitude ˜ τi,t ≤ n, where n is the number of variables in the model.

3 The First Challenge: Bias

Perhaps the most basic result about sequential Gibbs sampling is the fact that, in the limit of large numbers of samples, it is unbiased. Unfortunately, this is not the case for HOGWILD! Gibbs sampling, because the race conditions from the asynchrony add bias to the samples. To understand why, we can look at a simple example. Consider running HOGWILD! Gibbs on a model with two variables X1 and X2 each taking on values in {0, 1}, and having distribution

Suppose that the state is currently (1, 1) and two threads, T1 and T2, simultaneously update X1 and X2 respectively. Since T1 reads state (1, 1) it will update X1 to 0 or 1 each with probability 0.5; the same will be true for T2 and X2. Therefore, after this happens, every state will have probability 0.25; this includes the state (0, 0) which should never occur! Over time, this race condition will produce samples with value (0, 0) with some non-zero frequency; this is an example of bias introduced by the HOGWILD! sampling. Worse, this bias is not just theoretical: Figure 1 illustrates how the measured distribution for this model is affected by two-thread asynchronous execution. To measure the amount of this bias, it is standard to use the total variation distance.

Definition 1 (Total Variation Distance). The total variation distance [Levin et al., 2009] between two probability measures µ and ν on probability space Ω is defined as

that is, the maximum difference between the probabilities that µ and ν assign to a single event A.

Here, we observe a total variation distance of 9.8% between the measured and the true distribution. Unlike in the sequential case, this bias doesn't disappear as the number of samples goes to infinity. This example shows that asynchronous Gibbs sampling will not necessarily samples that approach the target distribution. Instead, the samples may Distribution of Sequential vs. HOGWILD! Gibbs approach some other distribution, which we hope is sufficiently similar for some practical purpose. Often, the purpose of Gibbs sampling is to estimate the marginal distributions of individual variables or of events that each depend on only a small number of variables. To characterize the accuracy of these estimates, the total variation distance is too conservative: it depends on the difference over all the events in the space, when most of these are events that we do not care about. To address this, we introduce a new metric.

Definition 2 (Sparse Variation Distance). For any event A in a probability space Ω over a set of variables V , let |A| denote the number of variables upon which A depends. Then, for any two distributions µ and ν over Ω, we define the ω-sparse variation distance to be

For the wide variety of applications that use sampling for marginal estimation, the sparse variation distance measures the quantity we actually care about: the maximum possible bias in the marginal distribution of the samples. As we will show, asynchronous execution seems to have less effect on the sparse variation distance than the total variation distance, because sparse variation distance uses a more localized view of the chain. For example, in Figure 1, the total variation distance between the sequential and HOGWILD! distributions is 9.8%, while the 1-sparse variation distance is only 0.4%.

This definition suggests the question: how long do we have to run before our samples have low sparse variation distance from the target distribution? To answer this question, we introduce the following definition.

Definition 3 (Sparse Estimation Time). The ω-sparse estimation time of a stochastic sampler with distribution P (t)µ0 at time t and target distribution π is the first time t at which, for any initial distribution µ0, the estimated distribution is within sparse variation distance ϵ of π,

When analyzing Gibbs sampling, we can bound the sparse estimation time (and bias) using a parameter called the total influence. While we arrived at this condition independently, it has been studied before in the context of Dobrushin's condition, which ensures rapid mixing of Gibbs sampling.

Definition 4 (Total Influence). Let π be a distribution over some set of variables I. Let Bj be the set of state pairs (X, Y ) which differ only at variable j. Let πi(·|XI\{i}) denote the conditional distribution in π of variable i given the other variables in state X. Then define α, the total influence of π, as

We say the model satisfies Dobrushin's condition if α < 1.

One way to think of total influence for factor graphs is as a generalization of maximum degree; indeed, if a factor graph has maximum degree ∆, it can easily be shown that α ≤ ∆. It turns out that if we can bound both this parameter and the sparse estimation time of sequential Gibbs sampling, we can prove a simple bound on the asymptotic sparse estimation time for asynchronous Gibbs sampling.

Theorem 1. Assume that we have a class of distributions with bounded total influence α = O(1). For each distribution π in the class, let ¯ tSE-seq(ω)(π, ϵ) be an upper bound on the ω-sparse estimation time of its sequential Gibbs sampler, and assume that it is a convex, decreasing function of ϵ. Further assume that, for any fixed ϵ, across all distributions,

where n is the number of variables in the model. 1 Then, for any fixed ϵ, the sparse estimation time of HOGWILD!-Gibbs across all models is bounded by

Roughly, this means that HOGWILD!-Gibbs sampling "works" (in the sense of producing output with arbitrarily small bias ϵ) asymptotically on all classes of problems for which we know marginal estimation is "fast" and the total influence is bounded. Since the sparse estimation times here are measured in iterations, and the asynchronous sampler is able, due to parallelism, to run many more iterations in the same amount of wall clock time, this result implies that HOGWILD!-Gibbs can be much faster than sequential Gibbs for producing estimates of similar quality.

4 The Second Challenge: Mixing Times

Even though the HOGWILD!-Gibbs sampler produces biased estimates, it is still interesting to analyze how long we need to run it before the samples it produces are independent of its initial conditions. To measure the efficiency of a Markov chain, it is standard to use the mixing time.

Definition 5 (Mixing Time). The mixing time [Levin et al., 2009] of a stochastic process with transition matrix P (t) at time t and target distribution π is the first time t at which, for any initial distribution µ0, the estimated distribution is within TV-distance ϵ of P (t)π. That is,

1 In many practical systems [Neubig, 2014; Shin et al., 2015], Gibbs sampling is naively run a fixed number of passes through the dataset, which works for exactly those models for which accurate marginal estimates can be achieved after O(n) samples. This suggests that models with sparse estimation time O(n) are practical and commonly occurring.

Unfortunately, as was the case with bias, the mixing time can also be greatly affected by running asynchronously. We can prove that there are classes of models for which asynchronous execution disastrously increases the mixing time.

Lemma 1. There exist classes of distributions for which the mixing time of sequential Gibbs sampling is O(n logn) but the mixing time of HOGWILD!-Gibbs sampling, even with delay τ = O(1), can be exp(Ω(n)).

This Lemma shows that fast mixing of the sequential sampler alone is not sufficient to guarantee fast mixing of the HOGWILD! chain. Consequently, we look for classes of models for which we can say something about the mixing time of both sequential and HOGWILD!-Gibbs. Dobrushin's condition is well known to imply rapid mixing of sequential Gibbs, and it turns out that we can leverage it again here to bound the mixing time of HOGWILD!-Gibbs.

Theorem 2. Assume that we run Gibbs sampling on a distribution that satisfies Dobrushin's condition, α < 1. Then the mixing time of sequential Gibbs will be bounded by

Under the same conditions, the mixing time of HOGWILD!- Gibbs will be bounded by

Dobrushin's condition holds for a wide variety of problems, including the Ising model at high temperatures. We can compare these two mixing time results as

the bounds on the mixing times differ by a negligible factor of 1+O(n-1 ). This result shows that, for problems that satisfy Dobrushin's condition, HOGWILD!-Gibbs sampling mixes in about the same time as sequential Gibbs sampling, and is therefore a practical choice for generating samples.

5 Experiments

Now that we have derived a theoretical characterization of the behavior of HOGWILD!-Gibbs sampling, we examine whether this characterization holds up under experimental evaluation. Specifically, we want to check whether increasing the expected delay parameter τ actually increases the mixing time as predicted by Equation 1.

To do this, we simulated HOGWILD!-Gibbs sampling running on a random synthetic Ising model graph of order n = 1000, degree ∆ = 3, and inverse temperature β = 0.2. This model has total influence α ≤ 0.6, and Theorem 2 guarantees that it will mix rapidly. To estimate the mixing time, which is difficult to calculate experimentally, we use a technique called coupling to the future, which produces asymptotically tight upper-bound estimates of the mixing time. The blue series in Figure 2 displays the estimate of tmix(1/4) across a range of τ 2 . The red line in Figure 2 shows the mixing time that would be predicted by naively applying Equation 1 using the estimate of the sequential mixing time as a starting point-we can see that it is a very good match for the experimental results. This experiment shows that, at least for one archetypal model, our theory accurately characterizes the behavior of HOGWILD! Gibbs sampling as the delay parameter τ is changed, and that using HOGWILD!-Gibbs doesn't cause the model to catastrophically fail to mix.

2 We sampled the delays ˜ τi,t to be i.i.d. according to the maximum-entropy distribution supported on {0, 1, . . . , 200} consistent with a particular assignment of τ .

6 Conclusion

We analyzed HOGWILD!-Gibbs sampling, a heuristic for parallelized MCMC sampling, on discrete-valued graphical models. First, we constructed a statistical model for HOGWILD!-Gibbs by adapting a model already used for the analysis of asynchronous SGD. Next, we illustrated a major issue with HOGWILD!-Gibbs sampling: that it produces biased samples. To address this, we proved that if for some class of models with bounded total influence, only O(n) sequential Gibbs samples are necessary to produce good marginal estimates, then HOGWILD!-Gibbs sampling produces equally good estimates after only O(1) additional steps. Additionally, for models that satisfy Dobrushin's condition (α < 1), we proved mixing time bounds for sequential and asynchronous Gibbs sampling that differ by only a factor of 1 +O(n-1 ).

Acknowledgments

The authors acknowledge the support of: DARPA FA8750-12-2-0335; NSF IIS-1247701; NSF CCF-1111943; DOE 108845; NSF CCF-1337375; DARPA FA8750-132-0039; NSF IIS-1353606; ONR N000141210041 and N000141310129; NIH U54EB020405; Oracle; NVIDIA; Huawei; SAP Labs; Sloan Research Fellowship; Moore Foundation; American Family Insurance; Google; and Toshiba. "The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of DARPA, AFRL, NSF, ONR, NIH, or the U.S. Government."