Patterns of Scalable Bayesian Inference
Elaine Angelino, Matthew James Johnson, Ryan P. Adams
Chapter 1 Introduction
We have entered a new era of scientific discovery, in which computational insights are being integrated with large-scale statistical data analysis to enable researchers to ask both grander and more subtle questions about our natural world. This viewpoint asserts that we need not be limited to the narrow hypotheses that can be framed by traditional small-scale analysis techniques. Supporting new kinds of data-driven queries, however, requires that new methods be developed for statistical inference that can scale up along multiple axes — more samples, more dimensions, and greater model complexity — as well as scale out by taking advantage of modern parallel compute environments.
There are a variety of methodological frameworks for performing statistical inference, e.g., performing estimation and evaluating hypotheses; here we are concerned with the Bayesian formalism. In the Bayesian setting, queries about structure in data are framed as interrogations of the posterior distribution over parameters, missing data, and other unknowns; these unobserved quantities are treated as random variables. By conditioning on the data, the Bayesian hopes to not only perform point estimation, but also to understand the uncertainties associated with those estimates.
Accounting for uncertainty is central to Bayesian analysis, and so the computations associated with most common tasks – e.g., estimation, prediction, evaluation of hypotheses – are typically integrations. In some situations, it is possible to perform such integrations exactly, either by taking advantage of conjugate structure in the prior-likelihood pair, or by using dynamic programming when the dependencies between random variables are appropriately simple. Unfortunately, most real-world analysis problems are not amenable to these exact inference procedures and so most of the interest in Bayesian computation focuses on better methods of approximate inference.
There are two dominant paradigms for approximate inference in Bayesian models: Monte Carlo sampling methods and variational approximations. The Monte Carlo approach observes that integrations performed to query posterior distributions can be framed as expectations, and thus estimated with samples; such samples are most often generated via simulation from carefully designed Markov chains. Variational inference seeks to compute these integrals by approximating the posterior distribution with a more tractable alternative, where identification of the best approximation can then be performed using powerful optimization techniques.
In this paper, we examine how these techniques can be scaled up to larger problems and scaled out across parallel computational resources. This is not intended to be an exhaustive survey of a rapidly-evolving area of research; rather, we seek to identify the main ideas and themes that are emerging in this area, and articulate what we believe are some of the significant open questions and challenges.
The Bayesian paradigm is fundamentally about integration: integration computes posterior estimates and measures of uncertainty, eliminates nuisance variables or missing data, and averages models to compute predictions or perform model comparison. While some statistical methods, such as MAP estimation, can be described from a Bayesian perspective, in which case the prior might serves as a regularizer in an optimization problem, such methods are not inherently or exclusively Bayesian. Posterior integration is the distinguishing characteristic of Bayesian statistics, and so a defense of Bayesian ideas in the big data regime rests on the utility of integration.
But from a classical perspective, the big data setting might seem to be precisely where integration isn’t important: as the dataset grows, shouldn’t the posterior distribution concentrate towards a point mass? If big data means we end up making predictions with such concentrated posteriors, why not focus on point estimation and avoid the specification of priors and the burden of approximate integration?
These objections certainly apply to settings where the number of parameters is small and fixed (“tall data”). However, many models of interest have many parameters (“wide data”), or indeed have a number of parameters that grows along with the amount of data.
For example, an Internet company making inferences about its users’ viewing and buying habits may have terabytes of data in total but only a few observations for its newest customers, the ones most important to impress with personalized recommendations. Moreover, it may wish to adapt its model in an online way as data arrive, a task that benefits from calibrated posterior uncertainties (Stern et al., 2009). As another example, consider a healthcare company. As its dataset grows, it might hope to make more detailed and complex inferences about populations while also making careful predictions with calibrated uncertainty for each patient, even in the presence of massive missing data (Lawrence, 2015). These scaling issues also arise in astronomy, where hundreds of billions of light sources, such as stars, galaxies, and quasars, each have latent variables that must be estimated from very weak observations, and are coupled in a large hierarchical model (Regier et al., 2015). In Microsoft Bing’s sponsored search advertising, predictive probabilities inform the pricing in the keyword auction mechanism. This problem nevertheless must be solved at scale, with tens of millions of impressions per hour (Graepel et al., 2010).
These are the regimes where big data can be small (Lawrence, 2015) and the number and complexity of statistical hypotheses grows with the data. The Bayesian methods we survey in this paper may provide solutions to these challenges.
2 The fidelity of approximate integration
Bayesian inference may be important in some modern big data regimes, but exact integration in general is computationally out of reach. While decades of research in Bayesian inference in both statistics and machine learning have produced many powerful approximate inference algorithms, the big data setting poses some new challenges. Iterative algorithms that read the entire dataset before making each update become prohibitively expensive. Sequential computation that cannot leverage parallel and distributed computing resources is at a significant and growing disadvantage. Insisting on zero asymptotic bias from Monte Carlo estimates of expectations may leave us swamped in errors from high variance (Korattikara et al., 2014) or transient bias.
These challenges, and the tradeoffs that may be necessary to address them, can be viewed in terms of how accurate the integration in our approximate inference algorithms must be. Markov chain Monte Carlo (MCMC) algorithms that admit the exact posterior as a stationary distribution may be the gold standard for generically estimating posterior expectations, but if standard MCMC algorithms become intractable in the big data regime we must find alternatives and understand their tradeoffs. Indeed, someone using Bayesian methods for machine learning may be less constrained than a classical Bayesian statistician: if the ultimate goal is to form predictions that perform well according to a specific loss function, computational gains at the expense of the internal posterior representation may be worthwhile. The methods studied here cover a range of such approximate integration tradeoffs.
3 Outline
The remainder of this review is organized as five chapters. In Chapter 2, we provide relevant background material on exponential families, MCMC inference, mean field variational inference, and stochastic gradient optimization. The next three chapters survey recent algorithmic ideas for scaling Bayesian inference, highlighting theoretical results where possible. Each of these central technical chapters ends with a summary and discussion, identifying emergent themes and patterns as well as open questions. Chapters 3 and 4 focus on MCMC algorithms, which are inherently serial and often slow to converge; the algorithms in the first of these use various forms of data subsampling to scale up serial MCMC and in the second use a diverse array of strategies to scale out on parallel resources. In Chapter 5 we discuss two recent techniques for scaling variational mean field algorithms. Both process data in minibatches: the first applies stochastic gradient optimization methods and the second is based on incremental posterior updating. Finally, in Chapter 6 we provide an overarching discussion of the ideas we survey, focusing on challenges and open questions in large-scale Bayesian inference.
Chapter 2 Background
In this chapter we summarize background material on which the ideas in subsequent chapters are based. This chapter also serves to fix some common notation. Throughout the chapter, we avoid measure-theoretic definitions and instead assume that any density exists with respect either to Lebesgue measure or counting measure, depending on its context.
First, we cover some relevant aspects of exponential families. Second, we cover the foundations of Markov chain Monte Carlo (MCMC) algorithms, which are the workhorses of Bayesian statistics and are common in Bayesian machine learning. Indeed, the algorithms discussed in Chapters 3 and 4 are either MCMC algorithms or aim to approximate MCMC algorithms. Next, we describe the basics of mean field variational inference and stochastic gradient optimziation, both of which are used extensively in Chapter 5. Finally, we close the chapter with notes on computational architectures and useful notions for measuring performance.
Exponential families of densities play a key role in Bayesian analysis and many practical Bayesian methods. In particular, likelihoods that are exponential families yield natural conjugate prior families, which provide analytical and computational advantages in both MCMC and variational inference algorithms. Exponential families are also particularly relevant in the context of large datasets: in a precise sense, they are the only families of densities which admit a finite-dimensional sufficient statistic. Thus only exponential families allow arbitrarily large amounts of data to be summarized with a fixed-size description.
In this section we give basic definitions, notation, and results concerning exponential families. For perspectives from convex analysis see Wainwright and Jordan (2008), and for perspectives from differential geometry see Amari and Nagaoka (2007).
Exponential families are defined in terms of densities with respect to some underlying -finite measure, which we denote .
We say a parameterized family of densities is an exponential family if each density can be written as
where is an inner product on a finite-dimensional real vector space. We call the natural parameter vector, the statistic vector, the base density, and
We restrict our attention to families for which the support of the density does not depend on . When we say the family is written in natural parameters or natural coordinates, which we denote by writing . We say a family is regular if is open, and minimal if there is no nonzero such that is equal to a constant (-a.e.).
The statistic is sufficient in the sense of the Fisher-Neyman Factorization Theorem (Keener, 2010, Theorem 3.6) by construction
and hence contains all the information about that is relevant for the parameter . In the context of Bayesian analysis, in which is a random variable, this definition of sufficiency is equivalent to the conditional independence statement {\theta\mathchoice{\mathrel{\hbox to0.0pt{\displaystyle\perp\hss}\mkern 4.0mu{\displaystyle\perp}}}{\mathrel{\hbox to0.0pt{\textstyle\perp\hss}\mkern 4.0mu{\textstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptstyle\perp\hss}\mkern 4.0mu{\scriptstyle\perp}}}{\mathrel{\hbox to0.0pt{\scriptscriptstyle\perp\hss}\mkern 4.0mu{\scriptscriptstyle\perp}}}X\;|\;t(X)}. The Koopman-Pitman-Darmois Theorem shows that exponential families are the only families which provide this powerful summarization property, under some mild regularity conditions (Hipp, 1974).
Exponential families have many convenient analytical and computational properties. In particular, differentiating the log partition function generates cumulants:
For a regular exponential family of densities of the form (2.1) with , we have and
More generally, the moment generating function of can be written
and so derivatives of give cumulants of , where the first cumulant is the mean and the second and third cumulants are the second and third central moments, respectively.
To derive the form of the moment generating function, we write
For members of an exponential family, many quantities can be expressed generically in terms of the natural parameter, expected statistics under that parameter, and the log partition function.
When the family is regular, the score with respect to the natural parameter is
When the family is regular, the Fisher information with respect to the natural parameter is
Each follows from (2.1) and Proposition 2.1.2.
Below, we define a notion of conjugacy for pairs of families of distributions. Conjugate families are especially useful for Bayesian analysis and algorithms.
A parameterized (not necessarily exponential) family of densities is conjugate for a likelihood function if for every density in the posterior distribution
also belongs to , for some that may depend on and .
Conjugate pairs are particularly useful in Bayesian analysis because if we have a prior family and we observe data generated according to a likelihood then the posterior is in the same family as the prior. In the context of Bayesian updating, we call the hyperparameter and the posterior hyperparameter.
Given a regular exponential family likelihood, we can always define a conjugate prior, as shown in the next proposition.
then if we define the statistic and an exponential family of densities with respect to that statistic as
then the pair is a conjugate pair of families with
and hence we can write the posterior hyperparameter as
When the prior family is parameterized with its natural parameter, we have .
As a consequence of Proposition 2.1.7, if the prior family is written with natural parameters and we generate data according to the model
where the notation denotes that the random variables are independently and identically distributed, then has posterior hyperparameter . Therefore any tractable computations in the prior, such as simulation or computing expectations, are shared by the posterior. Furthermore, for inferences about , the entire dataset can be summarized by the statistic .
2 Markov Chain Monte Carlo inference
Markov chain Monte Carlo (MCMC) is a class of algorithms for estimating expectations with respect to distributions. These distributions may be intractable, such as most posterior distributions arising in Bayesian inference. Given a target distribution, a standard MCMC algorithm proceeds by simulating an ergodic random walk that admits the target distribution as its stationary distribution. As we develop in the following subsections, by collecting samples from the simulated trajectory and forming Monte Carlo estimates, expectations of many functions can be approximated to arbitrary accuracy. Thus MCMC is employed when samples or expectations from a distribution cannot be obtained directly, as is often the case with complex, high-dimensional systems arising across disciplines, such as estimating bulk material properties from molecular dynamics physics simulations or performing inference in Bayesian probabilistic models.
In this section, we first review the two underlying ideas behind MCMC algorithms: Monte Carlo methods and Markov chains. First we define the bias and variance of estimators. Next, we introduce Monte Carlo estimators based on independent and identically distributed samples. We then describe how Monte Carlo estimates can be formed using mutually dependent samples generated by a Markov chain simulation. Finally, we introduce two general MCMC algorithms commonly applied to Bayesian posterior inference, the Metropolis-Hastings and Gibbs sampling algorithms. Our exposition here mostly follows the standard treatment, such as in Brooks et al. (2011, Chapter 1), Geyer (1992), and Robert and Casella (2004).
Notions of bias and variance are fundamental to understanding and comparing estimator performance, and much of our discussion of MCMC methods is framed in these terms.
Consider using a scalar-valued random variable to estimate a fixed scalar quantity of interest . The bias and variance of the estimator are defined as
This decomposition provides a basic language for evaluating estimators and thinking about tradeoffs. Among unbiased estimators, those with lower variance are generally preferrable. However, when an unbiased estimator has high variance, a biased estimator that achieves low variance can have a lower overall mean squared error.
As we describe in the following sections, a substantial amount of the study of Bayesian statistical computation has focused on algorithms that produce asymptotically unbiased estimates of posterior expectations, in which the bias due to initialization is transient and is washed out relatively quickly. In this setting, the error is typically considered to be dominated by the variance term, which can be made as small as desired by increasing computation time without bound. When computation becomes expensive as in the big data setting, errors under a realistic computational budget may in fact be dominated by variance, as observed by Korattikara et al. (2014), or, as we argue in Chapter 6, transient bias. Several of the new algorithms we examine in Chapters 3 and 4 aim to adjust this tradeoff by allowing some asymptotic bias while effectively reducing the variance and transient bias contributions through more efficient computation.
2.2 Monte Carlo estimates from independent samples
where . In particular, as grows, the standard deviation of the sample average converges to zero at an asymptotic rate proportional to . More generally, for any real-valued measurable function , the Monte Carlo standard error (MCSE) in the estimate (2.29) asymptotically scales as regardless of the dimension of .
Monte Carlo estimators effectively reduce the problem of computing expectations to the problem of generating samples. However, the preceding statements require the samples used in the Monte Carlo estimate to be independent, and independent samples can be computationally difficult to generate. Instead of relying on independent samples, Markov chain Monte Carlo algorithms compute estimates using mutually dependent samples generated by simulating a Markov chain.
2.3 Markov chains
for all measurable sets . A Markov chain is memoryless in the sense that its future behavior depends only on the current state and is independent of its past history.
Given an initial density for , a Markov chain evolves this density from one time point to the next through iterative application of the transition operator. We write the application of the transition operator to a density to yield a new density as
Writing to denote repeated applications of the transition operator , the density of induced by and is then given by .
Markov chain simulation follows this iterative definition by iteratively sampling the next state using the current state and the transition operator. That is, after first sampling from , Markov chain simulation proceeds at time step by sampling according to the density induced by the fixed sample .
We are interested in Markov chains that converge in total variation to a unique stationary density in the sense that
for any initial distribution , where denotes the total variation norm on densities:
For a transition operator to admit as a stationary density, its application must leave invariant:
For a discussion of general conditions that guarantee a Markov chain converges to a unique stationary distribution, i.e., that the chain is ergodic, see Meyn and Tweedie (2009).
In some cases it is easy to show that a transition operator has a particular unique stationary distribution. In particular, it is clear that is the unique stationary distribution when a transition operator is reversible with respect to , i.e., it satisfies the detailed balance condition with respect to a density ,
which is a pointwise condition over . Integrating over on both sides gives:
which is precisely the required condition from (2.35). We can interpret (2.36) as stating that, for a reversible Markov chain starting from its stationary distribution, any transition is equilibrated by the corresponding reverse transition . Many MCMC methods are based on deriving reversible transition operators.
For a thorough introduction to Markov chains, see Robert and Casella (2004, Chapter 6) and Meyn and Tweedie (2009).
2.4 Markov chain Monte Carlo (MCMC)
Markov chain Monte Carlo (MCMC) methods simulate a Markov chain for which the stationary distribution is equal to a target distribution of interest, and use the simulated samples to form Monte Carlo estimates of expectations. That is, consider simulating a Markov chain with unique stationary density , as in Section 2.2.3, and collecting its trajectory into a set of samples . These collected samples can be used to form a Monte Carlo estimate for a function of a random variable with density via
Even though this Markov chain Monte Carlo estimate is not constructed from independent samples, it can asymptotically satisfy analogs of the Law of Large Numbers (LLN) and Central Limit Theorem (CLT) that were used to justify ordinary Monte Carlo methods in Section 2.2.2. We sketch these important results here.
where and where and denote the variance and covariance operators with the chain initialized in stationarity with . Thus standard error in the MCMC estimate also scales asymptotically as , with a constant that depends on the autocovariance function of the stationary version of the chain. See Meyn and Tweedie (2009, chapter 17) and Robert and Casella (2004, section 6.7) for precise statements of both the LLN and CLT for Markov chain Monte Carlo estimates and for conditions on the Markov chain which guarantee that these theorems hold.
These results show that the asymptotic behavior of MCMC estimates of the form (2.37) is generally comparable to that of ordinary Monte Carlo estimates as discussed in Section 2.2.2. However, in the non-asymptotic regime MCMC estimates differ from ordinary Monte Carlo estimates in an important respect: there is a transient bias due to initializing the Markov chain out of stationarity. That is, the initial distribution from which the first iterate is sampled is generally not the chain’s stationary distribution , since if it were then ordinary Monte Carlo could be performed directly. While the marginal distribution of each Markov chain iterate converges to the stationary distribution, the effects of initialization on the initial iterates of the chain contribute an error term to Eq. (2.37) in the form of a transient bias.
This transient bias does not factor into the asymptotic behavior described by the MCMC analogs of the LLN and the CLT; asymptotically, it decreases at a rate of at least and is hence dominated by the Monte Carlo standard error which decreases only at rate . However, its effects can be significant in practice, especially in machine learning. Whenever a sampled chain seems “unmixed” because its iterates are too dependent on the initialization, errors in MCMC estimates are dominated by this transient bias.
The simulation in Figure 2.1 illustrates these error terms in MCMC estimates and how they can behave as more Markov chain samples are collected. The LLN and CLT for MCMC describe the regime on the far right of the plot: the total error can be driven arbitrarily small because the MCMC estimates are asymptotically unbiased, and the total error is asymptotically dominated by the Monte Carlo standard error. However, before reaching the asymptotic regime, the error is often dominated by the transient initialization bias. Several of the new methods we survey can be understood as attempts to alter the traditional MCMC tradeoffs, as we discuss further in Chapter 6.
Transient bias can be traded off against Monte Carlo standard error by choosing different subsets of Markov chain samples in the MCMC estimator. As an extreme choice, instead of using the MCMC estimator (2.37) with the full set of Markov chain samples , transient bias can be minimized by forming estimates using only the last Markov chain sample:
However, this choice of MCMC estimator maximizes the Monte Carlo standard error, which asymptotically cannot be decreased below the posterior variance of the estimand. A practical choice is to form MCMC estimates using the last Monte Carlo samples, resulting in an estimator
With this choice, once the marginal distribution of the Markov chain iterates approaches the stationary distribution the error due to transient bias is reduced at up to exponential rates. See Figure 2.2 for an illustration. With any choice of MCMC estimator, transient bias can be asymptotically decreased at least as fast as , and potentially much faster, while MCSE can decrease only as fast as .
Using these ideas, MCMC algorithms provide a general means for estimating posterior expectations of interest: first construct an algorithm to simulate an ergodic Markov chain that admits the intended posterior density as its stationary distribution, and then simply run the simulation, collect samples, and form Monte Carlo estimates from the samples. The task then is to design an algorithm to simulate from such a Markov chain with the intended stationary distribution. In the following sections, we briefly review two canonical procedures for constructing such algorithms: Metropolis-Hastings and Gibbs sampling. For a thorough treatment, see Robert and Casella (2004) and Brooks et al. (2011, Chapter 1).
2.5 Metropolis-Hastings (MH) sampling
In the context of Bayesian posterior inference, the Metropolis-Hastings (MH) algorithm simulates a reversible Markov chain over a state space that admits the posterior density as its stationary distribution. The algorithm depends on a user-specified proposal density, , which can be evaluated numerically and sampled from efficiently, and also requires that the joint density can be evaluated (up to proportionality). The MH algorithm then generates a sequence of states according to Algorithm 1.
In each iteration, a proposal for the next state is drawn from the proposal distribution, conditioned on the current state . The proposal is stochastically accepted with probability given by the acceptance probability,
via comparison to a random variate drawn uniformly from the interval ${u<\alpha(\theta,\theta^{\prime})}{q(\theta^{\prime}\,|\,\theta)=q(\theta\,|\,\theta^{\prime})}$, in which case the acceptance probability is simply
Hastings (1970) later relaxed this by showing that the proposal distribution could be arbitrary.
One can show that the stationary distribution is indeed by showing that the MH transition operator satisfies detailed balance (2.36). The MH transition operator density is a two-component mixture corresponding to the ‘accept’ event and the ‘reject’ event:
To show detailed balance, it suffices to show the two balance conditions
and we can use the same manipulation as in (2.50) under the integral sign:
See Robert and Casella (2004, Section 7.3) for a more detailed treatment of the Metropolis-Hastings algorithm.
2.6 Gibbs sampling
Given a collection of random variables , the Gibbs sampling algorithm iteratively samples each variable conditioned on the sampled values of the others. When the random variables are Markov on a graph , the conditioning can be reduced to each variable’s respective Markov blanket, as in Algorithm 2. In the context of Bayesian inference, the posterior of interest may correspond to conditioning on some subset of the random variables, fixing them to observed values.
A variant of the systematic scan of Algorithm 2, in which nodes are traversed in a fixed order for each outer iteration, is the random scan, in which nodes are traversed according to a random permutation sampled for each outer iteration. An advantage of the random scan (and other variants) is that the chain becomes reversible and therefore simpler to analyze (Robert and Casella, 2004, Section 10.1.2). With the conditional independencies implied by a graph, some sampling steps may be performed in parallel.
The Gibbs sampling algorithm can be analyzed as a special case of the Metropolis-Hastings algorithm, where the proposal distribution is based on the conditional distributions and the acceptance probability is always one. If the Markov chain produced by a Gibbs sampling algorithm is ergodic, then the stationary distribution is the target distribution of (Robert and Casella, 2004, Theorem 10.6). The Markov chain for a Gibbs sampler can fail to be ergodic if, for example, the support of the target distribution is disconnected (Robert and Casella, 2004, Example 10.7). A sufficient condition for Gibbs sampling to be ergodic is that all conditional densities exist and are positive everywhere (Robert and Casella, 2004, Theorem 10.8).
For a more detailed treatment of Gibbs sampling theory, see Robert and Casella (2004, Chapters 6 and 10).
3 Mean field variational inference
In mean field, and variational inference more generally, the task is to approximate an intractable distribution, such as a complex posterior, with a distribution from a tractable family so that the posterior can be efficiently interrogated for estimations of interest. In this section we define the mean field optimization problem and derive the standard coordinate optimization algorithm. We also give some basic results on the relationship between mean field and both graphical model and exponential family structure. For concreteness and simpler notation, we work mostly with undirected graphical models; the results extend immediately to directed models.
For a probability density with respect to a base measure of the form
where is the unnormalized density, for all densities with respect to we have
To show the equality, with we write
The inequality follows from the property , known as Gibbs’s inequality, which follows from Jensen’s inequality and the fact that the logarithm is concave:
with equality if and only if (-a.e.).
We call the negative log of in (2.54) the energy and the variational lower bound, and say decomposes into the entropy minus the average energy as in (2.56). For two densities and with respect to the same base measure, is the Kullback-Leibler divergence from to , used as a measure of dissimilarity between pairs of densities (Amari and Nagaoka, 2007).
The variational inequality given in Proposition 2.3.1 is useful in inference because if we wish to approximate an intractable with a tractable by minimizing , we can equivalently choose to maximize , which is possible to evaluate since it does not include the partition function .
In the context of Bayesian inference, is usually an intractable posterior distribution of the form , is the unnormalized joint distribution , and is the marginal likelihood , which plays a central role in Bayesian model selection and the minimum description length (MDL) criterion (MacKay, 2002, Chapter 28) (Hastie et al., 2001, Chapter 7).
Given that graphical model structure can affect the complexity of probabilistic inference (Koller and Friedman, 2009) it is natural to consider families that factor according to tractable graphs.
Let be the density with respect to for a collection of random variables , and let
be a family of densities with respect to that factorize according to a graph with being the set of maximal cliques of . Then the mean field optimization problem is
where is defined as in (2.56).
Note that this optimization problem is not in general convex.In the sense of maximizing a concave objective over a convex set. However, when the model distribution is an exponential family the objective is concave in each individually, and hence an optimization procedure that updates each factor in turn, while holding the rest fixed, will converge to a local optimum (Wainwright and Jordan, 2008) (Bishop, 2006, Section 10.1.1) (Murphy, 2012, Section 22.3). We call such a coordinate ascent procedure on (2.63) a mean field algorithm.
For approximating families in a factored form, we can derive a generic update to be used in a mean field algorithm.
Given a mean field objective as in Definition 2.3.3, the optimal update to a factor fixing the other factors defined by is
where the expectation is over with
Dropping terms constant with respect to , we write
Finally, we note the simple form of updates for exponential family conjugate pairs.
If appears in only in an exponential family conjugate pair where
then the optimal factor is in the prior family with natural parameter
The result follows from substituting (2.69) and (2.70) into (2.64).
See Wainwright and Jordan (2008, Chapter 5) for a convex analysis perspective on mean field algorithms in graphical models composed of exponential families.
4 Stochastic gradient optimization
In this section we briefly review some basic ideas in stochastic gradient optimization. In particular, the basic algorithm we use in this paper is given in Algorithm 3 and sufficient conditions for its convergence to a local extreme point are given in Theorem 2.4.1.
Given a dataset , where each is a data minibatch, consider the optimization problem
where the objective function decomposes according to
In the context of variational Bayesian inference, the objective may be a variational lower bound on the log of the model evidence and may be the parameters of the variational family. In MAP inference, may be the log joint density and may be its parameters.
Using the decomposition of , we can compute unbiased Monte Carlo estimates of its gradient. In particular, if the random index is sampled from , denoting the probability of sampling index as , we have
Thus by considering a Monte Carlo approximation to the expectation over , we can generate stochastic approximate gradients of the objective using only a single at a time.
A stochastic gradient ascent algorithm uses these approximate gradients to perform updates and find a stationary point of the objective. At each iteration, such an algorithm samples a data minibatch, computes a gradient with respect to that minibatch, and takes a step in that direction. In particular, for a sequence of stepsizes and a sequence of positive definite matrices , a typical stochastic gradient ascent algorithm is given in Algorithm 3.
Stochastic gradient algorithms have very general convergence guarantees, requiring only weak conditions on the step size sequence and even the accuracy of the gradients themselves. We summarize a common set of sufficient conditions in Theorem 2.4.1. Proofs of this result, along with more general versions, can be found in Bertsekas and Tsitsiklis (1989) and Bottou (1998). Note also that while the construction here has assumed that the stochasticity in the gradients arises only from randomly subsampling a finite sum, more general versions allow for other sources of stochasticity, typically requiring only bounded variance and allowing some degree of bias (Bertsekas and Tsitsiklis, 1989, Section 7.8).
there exists a constant such that
there are positive constants and such that
and the stepsize sequence satisfies
then Algorithm 3 converges to a stationary point in the sense that
While stochastic optimization theory provides convergence guarantees, there is no general theory to analyze rates of convergence for nonconvex problems such as those that commonly arise in posterior inference. Indeed, the empirical rate of convergence often depends strongly on the variance of the stochastic gradient updates and on the choice of step size sequence. There are automatic methods to tune or adapt the sequence of stepsizes (Snoek et al., 2012; Ranganath et al., 2013), though we do not discuss them here. To make a single-pass algorithm, the minibatches can be sampled without replacement.
Chapter 3 MCMC with data subsets
In MCMC sampling for Bayesian inference, the task is to simulate a Markov chain that admits as its stationary distribution the posterior distribution of interest. While there are many standard procedures for constructing and simulating from such Markov chains, when the dataset is large many of these algorithms’ updates become computationally expensive. This growth in complexity naturally suggests the question of whether there are MCMC procedures that can generate approximate posterior samples without using the full dataset in each update. In this chapter, we focus on recent MCMC sampling schemes that scale Bayesian inference by operating on only subsets of data at a time.
In most Bayesian inference problems, the fundamental object of interest is the posterior density, which for fixed data is proportional to the product of the prior and the likelihood:
In this survey we are often concerned with posteriors where the data are conditionally independent given the model parameters , and hence the likelihood can be decomposed into a product of terms:
When is large, this factorization can be exploited to construct MCMC algorithms in which the updates depend only on subsets of the data.
In particular, we can use subsets of data to form an unbiased Monte Carlo estimate of the log likelihood and consequently the log joint density. The log likelihood is a sum of terms:
and we can approximate this sum using a random subset of terms
where is a uniformly random subset of . This approximation is an unbiased estimator and yields an unbiased estimate of the log joint density:
Several of the methods reviewed in this chapter exploit this estimator to perform MCMC updates.
2 Adaptive subsampling for Metropolis–Hastings
In traditional Metropolis–Hastings (MH), we evaluate the joint density to decide whether to accept or reject a proposal. As noted by Korattikara et al. (2014), because the value of the joint density depends on the full dataset, when is large this is an unappealing amount of computation to reach a binary decision. In this section, we survey ideas for using approximate MH tests that depend on only a subset of the full dataset. The resulting approximate MCMC algorithms proceed in each iteration by reading only as much data as required to satisfy some estimated error tolerance.
While there are several variations, the common idea is to model the probability that the outcome of such an approximate MH test differs from the exact MH test. This probability model allows us to construct an approximate MCMC sampler, outlined in Section 3.2.2, where the user specifies some tolerance for the error in an MH test and the amount of data evaluated is controlled by an adaptive stopping rule. Different models for the MH test error lead to different stopping rules. Korattikara et al. (2014) use a normal model to construct a -statistic hypothesis test, which we describe in Section 3.2.3. Bardenet et al. (2014) instead use concentration inequalities, which we describe in Section 3.2.4. Given an error model and resulting stopping rule, both schemes rely on an MH test based on a Monte Carlo estimate of the log joint density, which we summarize in Section 3.2.1. Our notation in this section follows Bardenet et al. (2014).
Bardenet et al. (2014) observe that similar ideas have been developed both in the context of simulated annealingSimulated annealing is a stochastic optimization heuristic that is operationally similar to MH. by the operations research community (Bulgak and Sanders, 1988; Alkhamis et al., 1999; Wang and Zhang, 2006), and in the context of MCMC inference for factor graphs (Singh et al., 2012).
In the Metropolis–Hastings algorithm (§2.2.5), the proposal is stochastically accepted when
where . Rearranging and using log probabilities gives
Scaling both sides by gives an equivalent threshold,
where on the left, is the average log likelihood ratio,
This subsampled average log likelihood ratio is an unbiased estimate of the average log likelihood ratio . However, an error is made in the event that the approximate test (3.12) disagrees with the exact test (3.8), and the probability of such an error event depends on the distribution of and not just its mean.
2.2 Approximate MH with an adaptive stopping rule
A nested sequence of data subsets, sampled without replacement, that converges to the complete dataset gives us a sequence of approximate MH tests that converges to the exact MH test. Modeling the error of such an approximate MH test gives us a mechanism for designing an approximate MH algorithm in which, at each iteration, we incrementally read more data until an adaptive stopping rule informs us that our error is less than some user-specified tolerance. Algorithm 4 outlines this approach. The function AvgLogLikeRatioEstimate computes according to an adaptive stopping rule that depends on an error model, i.e., a way to approximate or bound the probability that the approximate outcome disagrees with the full-data outcome:
We describe two possible error models in Sections 3.2.3 and 3.2.4.
A practical issue with adaptive subsampling is choosing the sizes of the data subsets. One approach, taken by Korattikara et al. (2014), is to use a fixed batch size and read more data points at a time. Bardenet et al. (2014) instead geometrically increase the total subsample size, and also discuss connections between adaptive stopping rules and related ideas such as bandit problems, racing algorithms and boosting.
2.3 Using a t𝑡t-statistic hypothesis test
Korattikara et al. (2014) propose an approximate MH acceptance probability that uses a parametric test of significance as its error model. By assuming a normal model for the log likelihood estimate , a -statistic hypothesis test then provides an estimate of whether the approximate outcome agrees with the full-data outcome, i.e., the expression in Equation (3.14). This leads to an adaptive framework as in Section 3.2.2 where, at each iteration, the data are processed incrementally until the -test satisfies some user-specified tolerance .
The mean estimate for based on the subset of size is equal to :
To obtain a confidence interval, we multiply this estimate by the finite population correction, giving:
If is large enough for the CLT to hold, the test statistic
follows a Student’s -distribution with degrees of freedom when . The tail probability for then gives the probability that the approximate and actual outcomes agree, and thus
is the probability that they disagree, where is the CDF of the Student’s -distribution with degrees of freedom. The -test thus gives an adaptive stopping rule, i.e., for any user-provided tolerance , we can incrementally increase until . We illustrate this approach in Algorithm 5.
2.4 Using concentration inequalities
Bardenet et al. (2014) propose an adaptive subsampling method that is mechanically similar to using a -test but instead uses concentration inequalities. In addition to a bound on the error (of the approximate acceptance probability) that is local to each iteration, concentration bounds yield a bound on the total variation distance between the approximate and true stationary distributions.
As in Section 3.2.3, we evaluate an approximate MH threshold based on a data subset of size , given in Equation (3.12). We bound the probability that the approximate binary outcome is incorrect via concentration inequalities that characterize the quality of as an estimate for . Such a concentration inequality is a probabilistic statement that, for and some constant ,
For example, in Hoeffding’s inequality without replacement (Serfling, 1974)
Bardenet et al. (2014) use a concentration bound to construct an adaptive stopping rule based on a strategy called empirical Bernstein stopping (Mnih et al., 2008). Let be a concentration bound as in Equation (3.23) or (3.25) and let be the associated error. This concentration bound states that with probability . If , then the approximate MH test agrees with the exact MH test with probability . We reproduce a helpful illustration of this scenario from Bardenet et al. (2014) in Figure 3.1. If instead , then we want to increase until this is no longer the case. Let be the stopping time, i.e., the number of data points evaluated using this criterion,
We can set according to a user-defined parameter so that gives an upper bound on the error of the approximate acceptance probability. Let and set
under sampling without replacement. Hence, with probability , the approximate MH test based on agrees with the exact MH test. In other words, the stopping rule for computing in Algorithm 4 is satisfied once we observe . We illustrate this approach in Algorithm 6, using Hoeffding’s inequality without replacement.
In their actual implementation, Bardenet et al. (2014) modify to reflect the number of batches processed instead of the subsample size . For example, suppose we use the concentration bound in Equation (3.23), i.e., Hoeffding’s inequality without replacement. Then after processing a subsample of size in batches, the adaptive stopping rule checks whether , where
Also, as mentioned in Section 3.2.2, Bardenet et al. (2014) geometrically increase the subsample size by a factor . In their experiments, they use the empirical Bernstein-Serfling bound (Bardenet and Maillard, 2015). For the hyperparameters, they set , , and , and remark that they empirically found their algorithm to be robust to the choice of .
2.5 Error bounds on the stationary distribution
In this and the next subsection, we reproduce some theoretical results from Korattikara et al. (2014) and Bardenet et al. (2014). After setting up some notation, we emphasize the most general aspects of these results, which apply to pairs of transition kernels whose differences are bounded, and thus are not specific to adaptive subsampling procedures. The central theorem is an upper bound on the difference between the stationary distributions of such pairs of kernels in the case of Metropolis–Hastings. Its proof depends on the ability to bound the difference in the acceptance probabilities, at each iteration, of the two MH transition kernels.
Let and be probability measures (distributions) with Radon–Nikodym derivatives (densities) and , respectively, and absolutely continuous with respect to measure . The total variation distance between and is
be the acceptance probability error of the approximate MH test, with respect to the exact test. Finally, let
be the worst case absolute acceptance probability error.
Let be uniformly geometrically ergodic, i.e., there exists an integer , probability measure on , and constant such that for all and ,
and thus there exists a constant such that for all and ,
It follows that there exists a constant such that for all and ,
The upper bound in Equation (3.38) depends on the worst case acceptance probability error. For adaptive subsampling schemes, this depends on the choice of adaptive procedure.
We briefly outline a proof from Korattikara et al. (2014) of a similar theorem that exploits a stronger assumption on . Specifically, assume satisfies the contraction condition,
For approximate MH with an adaptive stopping rule, , the maximum acceptance probability error, gives an upper bound on the one-step error. Korattikara et al. (2014) show how to calculate an upper bound on when using a -test. Using concentration inequalities leads to a simpler bound: by construction, the user-defined error tolerance, , directly gives an upper bound on (Bardenet et al., 2014).
Finally, we note that an adaptive subsampling schemes using a concentration inequality enables an upper bound on the stopping time (Bardenet et al., 2014).
3 Sub-selecting data via a lower bound on the likelihood
Maclaurin and Adams (2014) introduce Firefly Monte Carlo (FlyMC), an auxiliary variable MCMC sampling procedure that operates on only subsets of data in each iteration. At each iteration, the algorithm dynamically selects what data to evaluate based on the random indicators included in the Markov chain state. In addition, it generates samples from the exact target posterior rather than an approximation. However, FlyMC requires a lower bound on the likelihood with a particular “collapsible” structure (essentially an exponential family lower bound) and is therefore not as generally applicable. The algorithm’s performance depends on the tightness of the bound; it can achieve impressive gains in performance when model structure allows.
FlyMC samples from an augmented posterior that eliminates potentially many likelihood factors. Define
and let be a strictly positive lower bound on , i.e., . For each datum, we introduce a binary auxiliary variable conditionally distributed according to a Bernoulli distribution,
where the are independent for different . When the bound is tight, i.e., , then with probability 1. More generally, a tighter bound results in a higher probability that . Augmenting the density with gives:
Using Equations (3.40) and (3.41), we can now write:
Thus for any fixed configuration of we can evaluate the joint density using only the likelihood terms where and the bound values for each .
While Equation (3.43) still involves a product of terms, if the product of the bound terms can be evaluated without reading each corresponding data point then the joint density can be evaluated reading only the data for which . In particular, if the form of is an exponential family density, then the product can be evaluated using only a finite-dimensional sufficient statistic for the data . Thus by exploiting lower bounds in the exponential family, FlyMC can reduce the amount of data required at each iteration of the algorithm while maintaining the exact posterior as its stationary distribution. Maclaurin and Adams (2014) show an application of this methodology to Bayesian logistic regression.
FlyMC presents three main challenges. The first is constructing a collapsible lower bound, such as an exponential family, that is sufficiently tight. The second is designing an efficient implementation. Maclaurin and Adams (2014) discuss these issues and, in particular, design a cache-like data structure for managing the relationship between the indicator values and the data. Finally, it is likely that the inclusion of these auxiliary variables slows the mixing of the Markov chain, but Maclaurin and Adams (2014) only provide empirical evidence that this effect is small relative to the computational savings from using data subsets.
4 Stochastic gradients of the log joint density
In this section, we review recent efforts to develop MCMC algorithms inspired by stochastic optimization techniques. This is motivated by the existence of, first, MCMC algorithms that can be thought of as the sampling analogues of optimization algorithms, and second, scalable stochastic versions of these optimization algorithms.
Traditional gradient ascent or descent performs optimization by iteratively computing and following a local gradient (Dennis and Schnabel, 1983). In Bayesian MAP inference, the objective function is typically a log joint density and the update rule for gradient ascent is given by
for . As discussed in Section 2.4, stochastic gradient descent (SGD) is simple modification of gradient descent that exploits situations where the objective function decomposes into a sum of many terms. While the traditional gradient descent update depends on all the data, i.e.,
SGD forms an update based on only a data subset,
The iterates converge to a local extreme point of the log joint density in the sense that if the step size sequence satisfies
A common choice of step size sequence is for some and .
Welling and Teh (2011) propose stochastic gradient Langevin dynamics (SGLD), an approximate MCMC procedure that combines SGD with a simple kind of Langevin dynamics (Langevin Monte Carlo) (Neal, 1994). They extend the Metropolis-adjusted Langevin algorithm (MALA) that uses noisy gradient steps to generate proposals for a Metropolis–Hastings chain (Roberts and Tweedie, 1996). At iteration , the MH proposal is
where the injected noise is Gaussian. Notice that the scale of the noise is , i.e., is constant and set by the gradient step size parameter. The MALA proposal is thus a stochastic gradient step, constructed by adding noise to a step in the direction of the gradient.
SGLD modifies the Langevin dynamics in Equation (3.48) by using stochastic gradients based on data subsets, as in Equation (3.46), and requiring that the step size parameter satisfy Equation (3.47). Thus, at iteration , the proposal is
where . Notice that the injected noise decays with the gradient step size parameter, but at a slower rate. Specifically, if decays as , then decays as . As in MALA, the SGLD proposal is a stochastic gradient step, where the noise comes from subsampling as well as the injected noise.
An actual Metropolis–Hastings algorithm would accept or reject the proposal in Equation (3.49) by evaluating the full (log) joint density at and , but this is precisely the computation we wish to avoid. Welling and Teh (2011) observe that as , in both Equations (3.48) and (3.49). In this limit, the probability of accepting the proposal converges to , but the chain stops completely. The authors suggest that can be decayed to a value that is large enough for efficient sampling, yet small enough for the acceptance probability to essentially be . These assumptions lead to a scheme where , for all , and all proposals are accepted, therefore the acceptance probability is never evaluated. We show this scheme in Algorithm 7. Without the stochastic MH acceptance step, however, asymptotic samples are no longer guaranteed to represent the target distribution.
In more recent work, Patterson and Teh (2013) apply SGLD to Riemann manifold Langevin dynamics (Girolami and Calderhead, 2011) and Chen et al. (2014) combine the idea of SGD with Hamiltonian Monte Carlo (HMC), an improved generalization of Langevin dynamics (Neal, 1994, 2010). Finally, we note that all the methods in this section require gradient information that might not be readily computable.
5 Summary
In this chapter, we have surveyed three recent approaches to scaling MCMC that operate on subsets of data. Below and in Table 3.1, we summarize and compare adaptive subsampling approaches (§3.2), FlyMC (§3.3), and SGLD (§3.4) along several axes.
Adaptive subsampling approaches replace the Metropolis–Hastings (MH) test, a function of all the data, with an approximate test that depends on only a subset. FlyMC is an auxiliary variable method that stochastically replaces likelihood computations with a collapsible lower bound. Stochastic gradient Langevin dynamics (SGLD) replaces gradients in a Metropolis-adjusted Langevin algorithm (MALA) with stochastic gradients based on data subsets and eliminates the Metropolis–Hastings test.
Each of the methods exploits assumptions or additional problem structure. Adaptive subsampling methods require an error model that accurately represents the probability that an approximate MH test will disagree with the exact MH test. A normal model (Korattikara et al., 2014) or concentration bounds (Bardenet et al., 2014) represent natural choices; under certain conditions, tighter concentration bounds may apply. FlyMC requires a strictly positive collapsible lower bound on the likelihood, essentially an exponential family lower bound, which may not in general be available. SGLD requires the log gradients of the prior and likelihood.
While all the methods use subsets of data, their access patterns differ. Adaptive subsampling and SGLD require randomization to avoid issues of bias due to data order, but this randomization can be achieved by permuting the data before each pass and hence these algorithms allow data access that is mostly sequential. In contrast, FlyMC operates on random subsets of data determined by the Markov chain itself, leading to a random access pattern. However, subsets from one iteration to the next tend to be correlated, and motivate implementation details such as the proposed cache data structure.
FlyMC does not introduce additional hyperparameters that require tuning. Both adaptive subsampling methods and SGLD introduce hyperparameters that can significantly affect performance. Both are mini-batch methods, and thus have the batch size as a tuning parameter. In adaptive subsampling methods, the stopping criterion is evaluated potentially more than once before it is satisfied. This motivates schemes that geometrically increase the amount of data processed whenever the stopping criterion is not satisfied, which introduces additional hyperparameters. Adaptive subsampling methods additionally provide a single tuning parameter that allows the user to control the error at each iteration. Finally, since these adaptive methods define an approximate MH test, they implicitly also require that the user specify a proposal distribution. For SGLD, the user must specify an annealing schedule for the step size parameter; in particular, it should converge to a small positive value so that the injected noise term dominates, while not being too large compared to the scale of the posterior distribution.
FlyMC is exact in the sense that the target posterior distribution is a marginal of its augmented state space. The adaptive subsampling approaches and SGLD are approximate methods in that neither has a stationary distribution equal to the target posterior. The adaptive subsampling approaches bound the error of the MH test at each iteration, and for MH transition kernels with uniform ergodicity this one-step error bound leads to an upper bound on the total variation distance between the approximate stationary distribution and the target posterior distribution. The theoretical analysis of SGLD is less clear (Sato and Nakagawa, 2014).
6 Discussion
The methods surveyed in this chapter achieve computational gains by using data subsets in place of an entire dataset of interest. The adaptive subsampling algorithms (§3.2) are more successful when a small subsample leads to an accurate estimator for the exact MH test’s accept/reject decision. Intuitively, such an estimator is easier to construct when the log posterior values at the proposed and current states are significantly different. This tends to be true far away from the mode(s) of the posterior, e.g., in the tails of a distribution that decay exponentially fast, compared to the area around a mode, which is locally more flat. Thus, these algorithms tend to evaluate more data when the chain is in the vicinity of a mode, and less data when the chain is far away (which tends to be the case for an arbitrary initial condition). SGLD (§3.4) exhibits somewhat related behavior. Recall that SGLD behaves more like SGD when the update rule is dominated by the gradient term, which tends to be true during the initial execution phase. Similar to SGD, the chain progresses toward a mode at a rate that depends on the accuracy of the stochastic gradients. For a log posterior target, stochastic gradients tend to be more accurate estimators of true gradients far away from the mode(s). In contrast, the MAP-tuned version of FlyMC (§3.3) requires the fewest data evaluations when the chain is close to the MAP, since by design, the lower likelihood bounds are tightest there. Meanwhile, the untuned version of FlyMC tends to exhibit the opposite behavior.
The Metropolis-Hastings algorithm requires the user to specify a proposal distribution. Fixing proposal distribution can be problematic, because the behavior of MH is sensitive to the proposal distribution and can furthermore change as the chain converges. A common solution, employed e.g., by Bardenet et al. (2014), is to use an adaptive MH scheme (Haario et al., 2001; Andrieu and Moulines, 2006). These algorithms tune the proposal distribution during execution, using information from the samples as they are generated, in a way that provably converges asymptotically. Often, it is desirable for the proposal distribution to be close to the target. This motivates adaptive schemes that fit a distribution to the observed samples and use this fitted model as the proposal distribution. For example, a simple online procedure can update the mean and covariance of a multidimensional Gaussian model as follows:
where indexes the MH iterations and controls the speed with which the adaptation vanishes. An appropriate choice is for . The tutorial by Andrieu and Thoms (2008) provides a review of this and other, more sophisticated, adaptive MH algorithms.
The subsampling-based methods in this chapter are conceptually modular, and some may be combined. For example, it might be of interest to consider a ‘tunable’ version of FlyMC that achieves even greater computational efficiency at the cost of its original exactness. For example, we might use an adaptive subsampling scheme (§3.2) to evaluate only a subset of terms in Equation (3.43); this subset would need to represent terms corresponding to both possible values of . As another example, Korattikara et al. (2014) suggest using adaptive subsampling as a way to ‘fix up’ SGLD. Recall that the original SGLD algorithm completely eliminates the MH test and blindly accepts all proposals, in order to avoid evaluating the full posterior. A reasonable compromise is to instead evaluate a fraction of the data within the adaptive subsampling framework, since this bounds the per-iteration error.
The estimator in Equation (3.4) based on a data subset is an unbiased estimator for the log likelihood; to be explicit,
is not an unbiased estimate of the likelihood. While it is possible to transform Equation (3.50) into an unbiased likelihood estimate, e.g., using a Poisson estimator (Wagner, 1987; Papaspiliopoulos, 2009; Fearnhead et al., 2010), it is not necessarily non-negative, which is a requirement to incorporate the estimator into a Metropolis-Hastings algorithm. In general, we cannot derive estimators that are both unbiased and nonnegative (Jacob and Thiery, 2015; Lyne et al., 2015). Pseudo-marginal MCMC algorithms,Pseudo-marginal MCMC is also known as exact-approximate sampling. first introduced by Lin et al. (2000), rely on non-negative unbiased likelihood estimators to construct unbiased MCMC procedures (Andrieu and Roberts, 2009). In this context, methods for constructing unbiased non-negative likelihood estimators include importance sampling (Beaumont, 2003) and particle filters (Andrieu et al., 2010; Doucet et al., 2015).
Chapter 4 Parallel and distributed MCMC
MCMC procedures that take advantage of parallel computing resources form another broad approach to scaling Bayesian inference. Because the computational requirements of inference often scale with the amount of data involved, and because large datasets may not even fit on a single machine, these approaches often focus on data parallelism. In this chapter we consider several approaches to scaling MCMC by exploiting parallel computation, either by adapting classical MCMC algorithms or by defining new simulation dynamics that are inherently parallel.
One way to use parallel computing resources is to run multiple sequential MCMC algorithms at once. However, running identical chains in parallel does not reduce the transient bias in MCMC estimates of posterior expectations, though it would reduce their variance. Instead of using parallel computation only to collect more MCMC samples and thus reduce only estimator variance without improving transient bias, it is often preferable to use computational resources to speed up the simulation of the chain itself. Section 4.1 surveys several methods that use parallel computation to speed up the execution of MCMC procedures, including both basic methods and more recent ideas.
Alternatively, instead of adapting serial MCMC procedures to exploit parallel resources, another approach is to design new approximate algorithms that are inherently parallel. Section 4.2 summarizes some recent ideas for simulations that can be executed in a data-parallel manner and have their results aggregated or corrected to represent posterior samples.
An advantage to parallelizing standard MCMC algorithms is that they retain their theoretical guarantees and analyses. Indeed, a common goal is to produce identical samples under serial and parallel execution, so that parallel resources enable speedups without introducing new approximations. This section first summarizes some basic opportunities for parallelism in MCMC and then surveys the speculative execution framework for MH.
The MH algorithm has a straightforward opportunity for parallelism. In particular, if the target posterior can be written as
then when the number of likelihood terms is large it may be beneficial to parallize the evaluation of the product of likelihoods. The communication between processors is limited to transmitting the value of the parameter and the scalar values of likelihood products. This basic parallelization, which naturally fits in a bulk synchronous parallel (BSP) computational model, exploits conditional independence in the probabilistic model, namely that the data are indepdendent given the parameter.
Gibbs sampling algorithms can exploit more fine-grained conditional independence structure, and are thus a natural fit for graphical models which express such structure. Given a graphical model and a corresponding graph coloring with colors that partitions the set of random variables into groups, the random variables in each color group can be resampled in parallel while conditioning on the values in the other groups (Gonzalez et al., 2011). Thus graphical models provide a natural perspective on opportunities for parallelism. See Figure 4.1 for some examples.
These opportunities for parallelism, while powerful in some cases, are limited by the fact that they require frequent global synchronization and communication. Indeed, at each iteration it is often the case that every element of the dataset is read by some processor and many processors must mutually communicate. The methods we survey in the remainder of this chapter aim to mitigate these limitations by adjusting the allocation of parallel resources or by reducing communication.
1.2 Speculative execution and prefetching
Another class of parallel MCMC algorithms uses speculative parallel execution to accelerate individual chains. This idea is called prefetching in some of the literature and appears to have received only limited attention.
As shown in Algorithm 1, the body of a MH implementation is a loop containing a single conditional statement and two associated branches. We can thus view the possible execution paths as a binary tree, illustrated in Figure 4.2. The vanilla version of parallel prefetching speculatively evaluates all paths in this binary tree on parallel processors (Brockwell, 2006). The sampled path will be exactly one of these, so with processors this approach achieves a speedup of with respect to single core execution, ignoring communication and bookkeeping overheads.
Naïve prefetching can be improved by observing that the two branches in Algorithm 1 are not taken with equal probability. For typical algorithm tunings, the reject branch tends to be more probable; a classic result for the optimal MH acceptance rate in the Gaussian case is 0.234 (Roberts et al., 1997), so prefetching scheduling policies can be built around the expectation of rejection. Angelino et al. (2014) provides a thorough review of these strategies.
Parallel predictive prefetching makes more efficient use of parallel resources by dynamically predicting the outcome of each MH test (Angelino et al., 2014). In the case of Bayesian inference, these predictions can be constructed in the same manner as the approximate MH algorithms based on subsets of data, as discussed in Section 3.2.2. Furthermore, these predictions can be made in the context of an error model, e.g., with the concentration inequalities used by Bardenet et al. (2014). This yields a straightforward and rational mechanism for allocating parallel cores to computations most likely to fall along the true execution path.
Algorithms 8 and 9 sketch pseudocode for an implementation of parallel predictive prefetching that follows a master-worker pattern. See Angelino (2014) for a formal description of the algorithm and implementation details.
2 Defining new data-parallel dynamics
In this section we survey two ideas for performing inference using new data-parallel dynamics. These algorithms define new dynamics in the sense that their iterates do not form ergodic Markov chains which admit the posterior distribution as an invariant distribution, and thus they do not qualify as classical MCMC schemes. Instead, while some of the updates in these algorithms resemble standard MCMC updates, the overall dynamics are designed to exploit parallel and distributed computation. A unifying theme of these new methods is to perform local computation on data while controlling the amount of global synchronization or communication.
One such family of ideas involves the definition of subposteriors, defined using only subsets of the full dataset. Inference in the subposteriors can be performed in parallel, and the results are then globally aggregated into an approximate representation of the full posterior. Because the synchronization and communication costs—as well as the approximation quality—are determined by the aggregation step, several such aggregation procedures have been proposed. In Section 4.2.1 we summarize some of these proposals.
Another class of data-parallel dynamics does not define independent subposteriors but instead, motivated by Gibbs sampling, focuses on simulating from local conditional distributions with out-of-date information. In standard Gibbs sampling, updates can be parallelized in models with conditional independence structure (Section 4.1), but without such structure the Gibbs updates may depend on the full dataset and all latent variables, and thus must be performed sequentially. These sequential updates can be especially expensive with large or distributed datasets. A natural approximation to consider is to run the same local Gibbs updates in parallel with out-of-date global information and only infrequent communication. While such a procedure loses the theoretical guarantees provided by standard Gibbs sampling analysis, some empirical and theoretical results are promising. We refer to this broad class of methods as Hogwild Gibbs algorithms, and we survey some particular algorithms and analyses in Section 4.2.2.
Suppose we want to divide the evaluation of the posterior across parallel cores. We can divide the data into partition elements, , also called shards, and factor the posterior into corresponding subposteriors, as
The contribution from the original prior is down-weighted so that the posterior is equal to the product of the subposteriors, i.e., . Note that a subposterior is not the same as the posterior formed from the corresponding partition, i.e.,
Once a large dataset has been partitioned across multiple machines, a natural alternative is to try running MCMC inference on each partition element separately and in parallel. This yields samples from each subposterior in Equation 4.3, but there is no obvious choice for how to combine them in a coherent fashion to form approximate samples of the full posterior. In this section, we survey various proposals for forming such a consensus solution from the subposterior samples. Algorithm 10 outlines the structure of consensus strategies for embarrassingly parallel posterior sampling. This terminology, used by Huang and Gelman (2005) and Scott et al. (2013), invokes related notions of consensus, notably those that have existed for decades in the optimization literature on data-parallel algorithms in decentralized or distributed settings. We discuss this topic briefly in Section 4.2.1.
Below, we present two recent consensus strategies for combining subposterior samples, through weighted averaging and density estimation, respectively. The earlier report by Huang and Gelman (2005) proposes four consensus strategies, based either on normal approximations or importance resampling; the authors focus on Gibbs sampling for hierarchical models and do not evaluate any actual parallel implementations. Another consensus strategy is the recently proposed variational consensus Monte Carlo (VCMC) algorithm, which casts the consensus problem within a variational Bayes framework (Rabinovich et al., 2015).
Throughout this section, Gaussian densities provide a useful reference point and motivate some of the consensus strategies. Consider the jointly Gaussian model
To arrive at an expression for the subposteriors, we begin by factoring the joint distribution into an appropriate product:
Thus the subposteriors are also Gaussian:
The challenge now is to design an appropriate set of weights.
where is the posterior covariance in Equation 4.8, then
where is the posterior mean in Equation 4.9. A similar calculation shows that .
Thus for the Gaussian model, is distributed according to the posterior distribution, indicating that Equation 4.16 gives the appropriate weights. Each weight matrix is a function of , the prior covariance, and the subposterior covariances . We can form a Monte Carlo estimate of each using the empirical sample covariance . Algorithm 11 summarizes this consensus approach with weighted averaging. While this weighting is optimal in the Gaussian setting, Scott et al. (2013) shows it to be effective in some non-Gaussian models. Scott et al. (2013) also suggests weighting each dimension of a sample by the reciprocal of its marginal posterior variance, effectively restricting the weight matrices to be diagonal.
Finally, one can sample from this posterior density estimator using MCMC; ideally, this density is straightforward to obtain and sample. In general, however, density estimation can yield complex models that are not amenable to efficient sampling.
Neiswanger et al. (2014) explore three density estimation approaches of various complexities. Their first approach assumes a parametric model and is therefore approximate. Specifically, they fit a Gaussian to each set of subposterior samples, yielding
where and are the empirical mean and covariance, respectively, of the samples from the th subposterior. This product of Gaussians simplifies to a single Gaussian , where
These parameters are straightforward to compute and the overall density estimate can be sampled with reasonable efficiency and even in parallel, if desired. Algorithm 12 summarizes this consensus strategy based on fits to Gaussians.
In the case when the model is jointly Gaussian, the parametric density estimator we form is , with and given in Equations 4.21 and 4.20, respectively. In this special case, the estimator exactly represents the Gaussian posterior. However, recall that we could have instead written the exact posterior directly as , where and are in Equations 4.9 and 4.8, respectively. Thus, computing the exact posterior is more or less as expensive as computing the density estimator, i.e., local matrix inversions (or corresponding linear system solves).
The second approach proposed by Neiswanger et al. (2014) is to use a nonparametric kernel density estimate (KDE) for each subposterior. Suppose we obtain samples from the th subposterior, then its KDE with bandwidth parameter has the following functional form:
i.e., the KDE is a mixture of kernels, each centered at one of the samples. If we use samples from each subposterior, then the density estimator for the full posterior is a complicated function with terms, since it the a product of such mixtures, and is therefore very challenging to sample from. Neiswanger et al. (2014) use a Gaussian KDE for each subposterior, and from this derive a density estimator for the full posterior that is a mixture of Gaussians with unnormalized mixture weights. They also consider a third, semi-parametric approach to density estimation given by the product of a parametric (Gaussian) model and a nonparametric (Gaussian KDE) correction. As the number of samples , the nonparametric and semi-parametric density estimates exactly represent the subposterior densities and are therefore asymptotically exact. Unfortunately, their complex mixture representations grow exponentially in size, rendering them somewhat unwieldy in practice.
The consensus strategies surveyed so far are embarrassingly parallel. These methods obtain samples from each subposterior independently and in parallel, and from these attempt to construct samples that (approximately) represent the posterior post-hoc. The methods in this section proceed similarly, but introduce some amount of information sharing between the parallel samplers. This communication pattern is reminiscent of the alternating direction method of multipliers (ADMM) algorithm for data-parallel convex optimization; for a detailed treatment of ADMM, see the review by Boyd et al. (2011).
Weierstrass samplers (Wang and Dunson, 2013) are named for the Weierstrass transform:
which was introduced by Weierstrass (1885). The transformed function is the convolution of a one-dimensional function with a Gaussian density of standard deviation , and so converges pointwise to as ,
where is the Dirac delta function. For > 0, can be thought of as a smoothed approximation to . Equivalently, if is the density of a random variable , then is the density of a noisy measurement of , where the noise is an additive Gaussian with zero mean and standard deviation .
Wang and Dunson (2013) analyzes a more general class of Weierstrass transforms by defining a multivariate version and also allowing non-Gaussian kernels:
For simplicity, we restrict our attention to the one-dimensional Weierstrass transform.
Weierstrass samplers use Weierstrass transforms on subposterior densities to define an augmented model. Let denote the -th subposterior,
so that the full posterior can be approximated as
The integrand of (4.25) defines the joint density of an augmented model that includes the as auxiliary variables:
The posterior of interest can then be approximated by the marginal distribution of in the augmented model,
with pointwise equality in the limit as . Thus by running MCMC in the augmented model, producing Markov chain samples of both and , we can generate approximate samples of the posterior. Furthermore, the augmented model is more amenable to parallelization due to its conditional independence structure: conditioned on , the subposterior parameters are rendered independent.
The same augmented model construction can be motivated without explicit reference to the Weierstrass transform of densities. Consider the factor graph model of the posterior in Figure 4.3(a), which represents the definition of the posterior in terms of subposterior factors,
This model can be equivalently expressed as a model where each subposterior depends on an exact local copy of . That is, writing as the local copy of for subposterior , the posterior is the marginal of a new augmented model given by
This new model can be represented by the factor graph in Figure 4.3(b), with potentials . Finally, rather than taking the to be exact local copies of , we can instead relax them to be noisy Gaussian measurements of :
Thus the potentials enforce some consistency across the noisy local copies of the parameter but allow them to be decoupled, where the amount of decoupling depends on . With smaller values of the approximate model is more accurate, but the local copies are more coupled and hence sampling in the augmented model is less efficient.
We can construct a Gibbs sampler for the joint distribution in Equation 4.25 by alternately sampling from and , for . It follows from Equation 4.25 that
where . The remaining Gibbs updates follow from Equation 4.25, which directly yields
This Gibbs sampler allows for parallelism but requires communication at every round. A straightforward parallel implementation, shown in Algorithm 13, generates the updates for in parallel, but the update for depends on the most recent values of all the . Wang and Dunson (2013) describes an approximate variant of the full Gibbs procedure that avoids frequent communication by only occasionally updating . In other efforts to exploit parallelism while avoiding communication, the authors propose alternate Weierstrass samplers based on importance sampling and rejection sampling.
2.2 Hogwild Gibbs
Instead of designing new data-parallel algorithms from scratch, another approach is to take an existing MCMC algorithm and execute its updates in parallel at the expense of accuracy or theoretical guarantees. In particular, Hogwild Gibbs algorithms take a Gibbs sampling algorithm (§2.2.6) with interdependent sequential updates (e.g., due to collapsed parameters or lack of graphical model structure) and simply run the updates in parallel anyway, using only occasional communication and out-of-date (stale) information from other processors. Because these strategies take existing algorithms and let the updates run ‘hogwild’ in the spirit of Hogwild! stochastic gradient descent in convex optimization (Recht et al., 2011), we refer to these methods as Hogwild Gibbs.
Similar approaches have a long history. Indeed, Gonzalez et al. (2011) attributes a version of this strategy, Synchronous Gibbs, to the original Gibbs sampling paper (Geman and Geman, 1984). However, these strategies have seen renewed interest, particularly due to extensive empirical work on Approximate Distributed Latent Dirichlet Allocation (AD-LDA) (Newman et al., 2007, 2009; Asuncion et al., 2008; Liu et al., 2011; Ihler and Newman, 2012), which showed that running collapsed Gibbs sampling updates in parallel allowed for near-perfect parallelism without a loss in predictive likelihood performance. With the growing challenge of scaling MCMC both to not only big datasets but also big models, it is increasingly important to understand when and how these approaches may be useful.
In this section, we first define some variations of Hogwild Gibbs based on examples in the literature. Next, we survey the empirical results and summarize the current state of theoretical understanding.
Here we define some Hogwild Gibbs methods and related schemes, such as the stale synchronous parameter server. In particular, we consider bulk-synchronous parallel and asynchronous variations. We also fix some notation used for the remainder of the section.
For all of the Hogwild Gibbs algorithms, as with standard Gibbs sampling, we are given a collection of random variables, where , and we assume that we can sample from the conditional distributions , where denotes . For the Hogwild Gibbs algorithms, we also assume we have processors, each of which is assigned a set of variables on which to perform MCMC updates. We represent an assignment of variables to processors by fixing a partition of , so that the th processor performs updates on the state values indexed by .
A bulk-synchronous parallel (BSP) Hogwild Gibbs algorithm assigns variables to processors and alternates between performing parallel processor-local updates and global synchronization steps. During epoch , the th processor performs MCMC updates, such as Gibbs updates, on the variables without communicating with the other processors; in particular, these updates are computed using out-of-date values for all . After all processors have completed their local updates, all processors communicate the updated state values in a global synchronization step and the system advances to the next epoch. We summarize this Hogwild Gibbs variant in Algorithm 14, in which the local MCMC updates are taken to be Gibbs updates.
Several special cases of the BSP Hogwild Gibbs scheme have been of interest. The Synchronous Gibbs scheme of Gonzalez et al. (2011) associates one variable with each processor, so that for each (in which case we may take since no local iterations are needed with a single variable). One may also consider the case where the partition is arbitrary and is very large, in which case the local MCMC iterations may converge and exact block samples are drawn on each processor using old statistics from other processors for each outer iteration. Finally, note that setting and reduces to standard Gibbs sampling on a single processor.
Another Hogwild Gibbs pattern involves performing updates asynchronously. That is, processors might communicate only by sending messages to one another instead of by a global synchronization. Versions of this Hogwild Gibbs pattern has proven effective both for collapsed latent Dirichlet allocation topic model inference (Asuncion et al., 2008), and for Indian Buffet Process inference (Doshi-Velez et al., 2009). A version was also explored in the Gibbs sampler of the Stale Synchronous Parameter (SSP) server of Ho et al. (2013), which placed an upper bound on the staleness of the entries of the state vector on each processor.
There are many possible communication strategies in the asynchronous setting, and so we follow a version of the random communication strategy employed by Asuncion et al. (2008). In this approach, after performing some number of local updates, a processor sends its updated state information to a set of randomly-chosen processors and receives updates from other processors. The processor then updates its state representation and performs another round of local updates. A version of this asynchronous Hogwild Gibbs strategy is summarized in Algorithm 15.
Despite its empirical successes, theoretical understanding of Hogwild Gibbs algorithms is limited. There are two settings in which some analysis has been offered: first, in a variant of AD-LDA, i.e., Hogwild Gibbs applied to Latent Dirichlet Allocation models, and second in the jointly Gaussian case.
The work of Ihler and Newman (2012) provides some understanding of the effectiveness of a variant of AD-LDA by bounding in terms of run-time quantities the one-step error probability induced by proceeding with sampling steps in parallel, thereby allowing an AD-LDA user to inspect the computed error bound after inference (Ihler and Newman, 2012, Section 4.2). In experiments, the authors empirically demonstrate very small upper bounds on these one-step error probabilities, e.g., a value of their parameter meaning that at least of samples are expected to be drawn just as if they were sampled sequentially. However, this per-sample error does not necessarily provide a direct understanding of the effectiveness of the overall algorithm because errors might accumulate over sampling steps; indeed, understanding this potential error accumulation is of critical importance in iterative systems. Furthermore, the bound is in terms of empirical run-time quantities, and thus it does not provide guidance on which other models the Hogwild strategy may be effective. Ihler and Newman (2012, Section 4.3) also provides approximate scaling analysis by estimating the order of the one-step bound in terms of a Gaussian approximation and some distributional assumptions.
The jointly Gaussian case is more tractable for analysis (Johnson et al., 2013; Johnson, 2014). In particular, Johnson (2014, Theorem 7.6.6) shows that for the BSP Hogwild Gibbs process to be stable, i.e., to form an ergodic Markov chain and have a well-defined stationary distribution, for any variable partition and any iteration schedule it suffices for the model’s joint Gaussian precision matrix to satisfy a generalized diagonal dominance condition. Because the precision matrix contains the coefficients of the log potentials in a Gaussian graphical model, the diagonal dominance condition captures the intuition that Hogwild Gibbs should be stable when variables do not interact too strongly. Johnson (2014, Proposition 7.6.8) gives a more refined condition for the case where the number of processor-local Gibbs iterations is large.
When a bulk-synchronous parallel Gaussian Hogwild Gibbs process defines an ergodic Markov chain and has a stationary distribution, Johnson (2014, Chapter 7) also provides an understanding of how that stationary distribution relates to the model distribution. Because both the model distribution and the Hogwild Gibbs process stationary distribution are Gaussian, accuracy can be measured in terms of the mean vector and covariance matrix. Johnson (2014, Proposition 7.6.1) shows that the mean of a stable Gaussian Hogwild Gibbs process is always correct. Johnson (2014, Propositions 7.7.2 and 7.7.3) identify a tradeoff in the accuracy of the process covariance matrix as a function of the number of processor-local Gibbs iterations: at least when the processor interactions are sufficiently weak, more processor-local iterations between synchronization steps increase the accuracy of the covariances among variables within each processor but decrease the accuracy of the covariances between variables on different processors. Johnson (2014, Proposition 7.7.4) also gives a more refined error bound as well as an inexpensive way to correct covariance estimates for the case where the number of processor-local Gibbs iterations is large.
3 Summary
Many ideas for parallelizing MCMC have been proposed, exhibiting many tradeoffs. These ideas vary in generality, in faithfulness to the posterior, and in the parallel computation architectures for which they are best suited. Here we summarize the surveyed methods, emphasizing their relative strengths on these criteria. See Table 4.1 for an overview.
Independent instances of serial MCMC algorithms can be run in an embarrassingly parallel manner, requiring only minimal communication between processors to ensure distinct initializations and to collect samples. This approach can reduce Monte Carlo variance by increasing the number of samples collected in any time budget, achieving an ideal parallel speedup, but does nothing to accelerate the warm-up period of the chains during which the transient bias is eliminated (see Section 2.2.4 and Chapter 6). That is, using parallel resources to run independent chains does nothing to improve mixing unless there is some mechanism for information sharing as in Nishihara et al. (2014). In addition, running an independent MCMC chain on each processor requires each processor to access the full dataset, which may be problematic for especially large datasets. These considerations motivate both subposterior methods and Hogwild Gibbs.
Some MCMC algorithms applied to models with particular structure allow for straightforward parallel implementation. In particular, when the likelihood is factorized across data points, the computation of the Metropolis–Hastings acceptance probability can be parallelized. This strategy lends itself to a bulk-synchronous parallel (BSP) computational model. Parallelizing MH in this way yields exact MCMC updates and can be effective at reducing the mixing time required by serial MH, but it requires a simple likelihood function and its implementation requires frequent synchronization and communication, mitigating parallel speedups unless the likelihood function is very expensive.
Gibbs sampling also presents an opportunity for direct parallelization for particular graphical model structures. In particular, given a graph coloring of the graphical model, variables corresponding to nodes assigned a particular color are conditionally mutually independent and can be updated in parallel without communication. However, frequent synchronization and significant communication can be required to transmit sampled values to neighbors after each update. Relaxing both the strict conditional independence requirements and synchronization requirements motivates Hogwild Gibbs.
The prefetching algorithms studied in Section 4.1.2 use speculative execution to transform traditional (serial) Metropolis–Hastings into a parallel algorithm without incurring approximate updates or requiring any model structure. The implementation naturally follows a master-worker pattern, where the master allocates (possibly speculative) computational work, such as proposal generation or (partial) density evaluation, to worker processors. Ignoring overheads, basic prefetching algorithms achieve at least logarithmic speedup in the number of processors available. More sophisticated scheduling by the master, such as predictive prefetching (Angelino et al., 2014), can increase speedup significantly. While this method is very general and yields the same iterates as serial MH, the speedup can be limited.
Subposterior methods, such as the consensus Monte Carlo algorithms and the Weierstrass samplers of Section 4.2.1, allow for data parallelism and minimal communication because each subposterior Markov chain can be allocated to a processor and simulation can proceed independently. Communication is required only for final sample aggregation in consensus Monte Carlo or the periodic resampling of the global parameter in the Weierstrass sampler. In consensus Monte Carlo, the quality of the approximate inference depends on both the effectiveness of the aggregation strategy and the extent to which dependencies in the posterior can be factorized into subposteriors. The Weierstrass samplers directly trade off approximation quality and the amount of decoupling between subposteriors.
The consensus Monte Carlo approach originated at Google (Scott et al., 2013) and naturally fits the MapReduce programming model, allowing it to be executed on large computational clusters. Recent work has extended consensus Monte Carlo and provides tools for designing simple consensus strategies (Rabinovich et al., 2015), but the generality and approximation quality of subposterior methods remain unclear. The Weierstrass sampler fits well into a BSP model.
Hogwild Gibbs of Section 4.2.2 also allows for data parallelism but avoids factorizing the posterior as in consensus Monte Carlo or instantiating coupled copies of a global parameter as in the Weierstrass sampler. Instead, processor-local sampling steps (such as local Gibbs updates) are performed with each processor treating other processors’ states as fixed at stale values; processors can communicate updated states less frequently, either via synchronous or asynchronous communication. Hogwild Gibbs variants span a range of parallel computation paradigms from fully synchronous BSP to fully asynchronous message-passing. While Hogwild Gibbs has proven effective in practice for several models, its applicability and approximation tradeoffs remain unclear.
4 Discussion
The ideas surveyed in this chapter suggest several challenges and questions.
The new data-parallel methods surveyed here, namely consensus Monte Carlo, the Weierstrass samplers, and Hogwild Gibbs, do not generate samples that are asymptotically distributed according to the target posterior. Instead, each generates samples that are asymptotically distributed according to a distribution that is meant to approximate the target posterior. While the nature of these approximations differ, each has a tradeoff between parallelism and accuracy: increasing parallelism by using more processors decreases the faithfulness of the asymptotic posterior approximation. This tradeoff may be inherent to most data-parallel MCMC schemes, though parallel predictive prefetching strategies do not suffer the same drawback.
The performance of data-parallel methods may be significantly affected by the data partitioning that assigns data subsets to processors. In the case of Hogwild Gibbs, it is probably best to choose a data partition that minimizes the strength of cross-processor dependence. Similarly, in the case of subposterior methods, some factorizations may be more effective than others. Since data paritioning is likely to have significant practical effects for all of these methods, it may be fruitful to develop and analyze general heuristics for assigning data to processors.
The works surveyed in this chapter introduce several alternative approximations. While each is well motivated, it is unclear how to choose the most appropriate method for a given model, or how to think about and compare the various approximations and tradeoffs. A more unified perspective is necessary, through empirical comparison or through analyzing their application to simple models that are tractable for analysis.
Chapter 5 Scaling variational mean field algorithms
Variational inference is a standard paradigm for posterior inference in Bayesian models. Because variational methods pose inference as an optimization problem, ideas in scalable optimization can in principle yield scalable posterior inference algorithms. In this chapter, we consider such scalable algorithms mainly in the context of mean field variational inference, which is often called variational Bayes.
These scalable variational inference algorithms can be compared to the algorithms of Chapters 3 and 4 in the same way that variational methods are usually compared to MCMC. That is, because inference is typically performed in a family of distributions that does not include the exact posterior, it can be said that variational methods do not fully instantiate the Bayesian computation that MCMC methods do (at least, when given unbounded computation time). Indeed, MAP inference, in which the posterior is represented only as a single atom, is an extreme case of variational inference. More generally, mean field variational families typically provide only unimodal approximations, and additionally cannot represent some posterior correlations near particular modes. As a result, MCMC methods can provide better performance even when the Markov chain only explores a single mode in a reasonable number of iterations.
Despite these potential shortcomings, variational inference is widely used in machine learning because the computational advantage over MCMC can be significant. This computational advantage is particularly salient in the context of scaling inference to large datasets. The big data context may also inform the relative cost of performing inference in a constrained variational family rather than attempting to represent the posterior exactly: when the posterior is concentrated, a variational approximation may suffice. While such questions may ultimately need to be explored empirically on a case-by-case basis, the scalable variational inference methods surveyed in this chapter provide the tools for such an exploration.
In this chapter we summarize two patterns of scalable variational inference. First, in Section 5.1, we discuss the application of stochastic gradient optimization methods to mean field variational inference problems. Second, in Section 5.2, we describe an alternative approach that instead leverages the idea of incremental posterior updating to develop an inference algorithm with minibatch-based updates.
Stochastic gradient optimization is a powerful tool for scaling optimization algorithms to large datasets, and it has been applied to mean field variational inference problems to great effect. While many traditional algorithms for optimizing mean field objective functions, including both gradient-based and coordinate optimization methods, require re-reading the entire dataset in each iteration, the stochastic gradient framework allows each update to be computed with respect to minibatches of the dataset while providing very general asymptotic convergence guarantees.
In this section we first summarize the stochastic variational inference (SVI) framework of Hoffman et al. (2013), which applies to models with complete-data conjugacy. Next, we discuss alternatives and extensions which can handle more general models at the cost of updates with greater variance and, hence, slower convergence.
This section follows the development in Hoffman et al. (2013). It depends on results from stochastic gradient optimization theory; see Section 2.4 for a review. For notational simplicity we consider each minibatch to consist of only a single observation; the generalization to minibatches of arbitrary sizes is immediate.
Many common probabilistic models are hierarchical: they can be written in terms of global latent variables (or parameters), local latent variables, and observations. That is, many models can be written as
where denotes global latent variables, denotes local latent variables, and denotes observations. See Figure 5.1 for a graphical model. Given such a class of models, the mean field variational inference problem is to approximate the posterior for fixed data with a distribution of the form by finding a local minimum of the KL divergence from the approximating distribution to the posterior or, equivalently, finding a local maximum of the marginal likelihood lower bound
Hoffman et al. (2013) develops a stochastic gradient ascent algorithm for such models that leverages complete-data conjugacy. Gradients of with respect to the parameters of have a convenient form if we assume the prior and each complete-data likelihood are a conjugate pair of exponential family densities. That is, if we have
then conjugacy identifies the statistic of the prior with the natural parameter and log partition function of the likelihood via
Conjugacy implies the optimal variational factor has the same form as the prior; that is, without loss of generality we can write in the same form as (5.3),
for some variational parameter .
Given this conjugacy structure, we can find a simple expression for the gradient of with respect to the global variational parameter , optimizing out the local variational factor . That is, we write the variational objective over global parameters as
Writing the optimal parameters of as , note that when is partially optimized to a stationary point of , so that at , the chain rule implies that the gradient with respect to the global variational parameters simplifies:
Because the optimal local factor can be computed with local mean field updates for a fixed value of the global variational parameter , we need only find an expression for the gradient in terms of the optimized local factors.
To find an expression for the gradient that exploits conjugacy structure, using (5.6) we can substitute
into the definition of in (5.2). Using the optimal form of , we have
where the constant does not depend on . Using the identity for natural exponential families that
Thus we can compute the gradient of with respect to the global variational parameters as
where the first two terms come from applying the product rule.
The matrix is the Fisher information of the variational family, since
In the context of stochastic gradient ascent, we can cancel the multiplication by the matrix simply by choosing the sequence of positive definite matrices in Algorithm 3 to be . This choice yields a stochastic natural gradient ascent algorithm (Amari and Nagaoka, 2007), where the updates are stochastic approximations to the natural gradient
Natural gradients effectively include a second-order quasi-Newton correction for local curvature in the variational family, making the updates invariant to reparameterization of the variational family and thus often improving performance of the algorithm. More importantly, at least for the case of complete-data conjugate families considered here, natural gradient steps are in fact easier to compute than ‘flat’ gradient steps in either the natural parameterization or moment parameterization of the variational family .
Therefore a stochastic natural gradient ascent algorithm on the global variational parameter proceeds at iteration by sampling a minibatch and taking a step of some size in an approximate natural gradient direction via
where we have assumed the minibatches are of equal size to simplify notation. The local variational factor is computed using a local mean field update on the data minibatch and the global variational factor. That is, if is not further factorized in the mean field approximation, it is computed according to
We summarize the general SVI algorithm in Algorithm 16.
1.2 Stochastic gradients with general nonconjugate models
The development of SVI in the preceding section assumes that and are a conjugate pair of exponential families. This assumption led to a particularly convenient form for the natural gradient of the mean field variational objective and hence an efficient stochastic gradient ascent algorithm. However, when models do not have this conjugacy structure, more general algorithms are required.
In this section we review Black Box Variational Inference (BBVI), which is a stochastic gradient algorithm for variational inference that can be applied at scale (Ranganath et al., 2014). The “black box” name suggests its generality: while the stochastic variational inference of Section 5.1.1 requires particular model structure, BBVI only requires that the model’s log joint distribution can be evaluated. It also makes few demands of the variational family, since it only requires that the family can be sampled and that the gradient of its log joint with respect to the variational parameters can be computed efficiently. With these minimal requirements, BBVI is not only useful in the big-data setting but also a tool for handling nonconjugate variational inference more generally. Because BBVI uses Monte Carlo approximation to compute stochastic gradient updates, it fits naturally into a stochastic gradient optimization framework, and hence it has the additional benefit of yielding a scalable algorithm simply by adding minibatch sampling to its updates at the cost of increasing their variance. In this subsection we review the general BBVI algorithm and then compare it to the SVI algorithm of Section 5.1.1. For a review of Monte Carlo estimation, see Section 2.2.2.
We consider a general model including parameters and observations divided into minibatches. The distribution of interest is the posterior and we write the variational family as , where we suppress the particular mean field factorization structure of from the notation. The mean field variational lower bound is then
Taking the gradient with respect to the variational parameter and expanding the expectation into an integral, we have
where we have moved the gradient into the integrand and applied the product rule to yield two terms. The first term is identically zero:
where we have used . To write the second term of (5.23) in a form that allows convenient Monte Carlo approximation, we first note the identity
and hence we can write the second term of (5.23) as
where in the final line we have written the expectation as a Monte Carlo estimate using a set of samples , where for . Notice that the gradient is written as a weighted sum of gradients of the variational log density with respect to the variational parameters, where the weights depend on the model log joint density.
The BBVI algorithm uses the Monte Carlo estimate (5.30) to compute stochastic gradient updates. This gradient estimator is also known as the score function estimator (Kleijnen and Rubinstein, 1996; Gelman and Meng, 1998). The variance of these updates, and hence the convergence of the overall stochastic gradient algorithm, depends both on the sizes of the gradients of the variational log density and on the variance of . Large variance in the gradient estimates can lead to very slow optimization, and so Ranganath et al. (2014) proposes and evaluates two variance reduction schemes, including a control variate method as well as a Rao-Blackwellization method that can exploit factorization structure in the variational family.
To provide a scalable version of BBVI, gradients can be further approximated by subsampling minibatches of data. That is, using we write (5.29) and (5.30) as
with the minibatch index distributed uniformly over and the minibatches are assumed to be the same size for simpler notation. This subsampling over minibatches further increases the variance of the updates and thus may further limit the rate of convergence of the algorithm. We summarize this version of the BBVI algorithm in Algorithm 17.
It is instructive to compare the fully general BBVI algorithm applied to hierarchical models to the SVI algorithm of Section 5.1.1; this comparison not only shows the benefits of exploiting conjugacy structure but also suggests a potential Rao-Blackwellization scheme. Taking and and starting from (5.22) and (5.29), we can write the gradient as
where is a set of samples with for . Thus if the entropy of the local variational distribution and the expectations with respect to of the log density can be computed without resorting to Monte Carlo estimation, then the resulting update would likely have a lower variance than the BBVI update that requires sampling over both and .
This comparison also makes clear the advantages of exploiting conjugacy in SVI: when the updates of Section 5.1.1 can be used, neither nor needs to be sampled. Furthermore, while BBVI uses stochastic gradients in its updates, the SVI algorithm of Section 5.1.1 uses stochastic natural gradients, adapting to the local curvature of the variational family. Computing stochastic natural gradients in BBVI would require both computing the Fisher information matrix of the variational family and solving a linear system with it.
1.3 Exploiting reparameterization for some nonconjugate models
While the score function estimator developed for BBVI in Section 5.1.2 is sufficiently general to handle essentially any model, some nonconjugate models admit convenient stochastic gradient estimators that can have lower variance. In particular, in settings where the latent variables are continuous (or any discrete latent variables that can be marginalized efficiently) samples from some variational distributions can be reparameterized in a way that enables an alternative stochastic gradient estimator. This technique is related to non-centered reparameterizations (Papaspiliopoulos et al., 2007) and has recently been called the reparameterization trick (Kingma and Welling, 2014; Rezende et al., 2014).
The reparameterization trick applies when samples , where has parameter , can be written as
where is a random variable with a distribution that does not depend on and where can be computed efficiently for almost every value of . In this case, we can compute stochastic estimates of the gradient of the variational objective by first writing a Monte Carlo approximation of the objective function itself:
This Monte Carlo approximation is a differentiable unbiased estimate of as a function of the variational parameter , and so we can form a Monte Carlo estimate of the gradient of the variational objective simply by differentiating it:
This estimator often has lower variance than the fully general score function estimator (Kingma and Welling, 2014) and can be easier to compute.
2 Streaming variational Bayes (SVB)
Streaming variational Bayes (SVB) provides an alternative framework in which to derive minibatch-based scalable variational inference (Broderick et al., 2013). While the methods of Section 5.1 generally apply stochastic gradient optimization algorithms to a fixed variational mean field objective, SVB instead considers the streaming data setting, in which case there may be no fixed dataset size and hence no fixed variational objective. To handle streaming data, the SVB approach is based on the classical idea of Bayesian updating, in which a posterior is updated to reflect new data as they become available. This sequence of posteriors is approximated by a sequence of variational models, and each variational model is computed from the previous variational model via an incremental update on new data.
More concretely, given a prior over a parameter and a (possibly infinite) sequence of data minibatches , each distributed independently according to a likelihood distribution , we consider the sequence of posteriors
Given an approximation updating algorithm one can compute a corresponding sequence of approximations
with . This sequential updating view naturally suggests an online or one-pass algorithm in which the update (5.40) is applied successively to each of a sequence of minibatches.
A sequence of such updates may also exploit parallel or distributed computing resources. For example, the sequence of approximations may be computed as
where indexes a set of data minibatches for which each update is computed in parallel before being combined in the final update from to .
This combination of partial results is especially appealing when the prior and the family of approximating distributions are in the same exponential family,
for a prior natural parameter and a sequence of variational parameters . In the exponential family case, the updates (5.42) can be written
where we may take the algorithm to return an updated natural parameter, .
Finally, similar updates can be performed in an asynchronous distributed master-worker setting. Each worker can process a minibatch and send the corresponding natural parameter increment to a master process, which updates the global variational parameter and transmits back the updated variational parameter along with a new data minibatch. In symbols, we can write that a worker operating on minibatch for some minibatch index computes the update increment according to
where is the index of the global variational parameter used in the worker’s computation. Upon receiving an update, the master updates its global variational parameter synchronously according to
We summarize a version of this process in Algorithms 18 and 19.
A related algorithm, which we do not detail here, is Memoized Variational Inference (MVI) (Hughes and Sudderth, 2013; Hughes et al., 2015). While this algorithm is designed for the fixed dataset setting rather than the streaming setting, the updates can be similar to those of SVB. In particular, MVI optimizes the mean field objective in the conjugate exponential family setting using the mean field coordinate descent algorithm but with an atypical update order, in which only some local variational factors are updated at a time. This update order enables minibatch-based updating but, unlike the stochastic gradient algorithms, does not optimize out the other local variational factors not included in the minibatch and instead leaves them fixed.
Streaming variational inference algorithms similar to SVB have also recently been studied in some Bayesian nonparametric mixture models (Tank et al., 2015).
3 Summary
The methods of Section 5.1 apply stochastic optimization to variational mean field inference objectives. In optimization literature and practice, stochastic gradient methods have a large body of both theoretical and empirical support, and so such methods offer a compelling framework for scalable inference. The streaming ideas surveyed in Section 5.2 are less well understood, but by treating the streaming setting, rather than the setting of a large fixed-size dataset, they may extend the reach of Bayesian modeling and inference.
All of this chapter’s scalable approaches to mean field variational inference are based on processing minibatches of data. These algorithms arrive at this data access pattern via two routes: the first applies stochastic gradient optimization to mean field variational inference (§5.1) and the second considers the streaming data setting (§5.2). SVI (§5.1.1) and BBVI (§5.1.2) optimize the variational objective by replacing full gradient updates with stochastic gradient updates. In both SVI and BBVI, these approximate gradients arise from randomly sampling data minibatches, while in BBVI there is additional stochasticity due to the Monte Carlo approximation required to handle nonconjugate structure. In contrast to SVI and BBVI, SVB (§5.2) processes data minibatches to drive incremental posterior updates, constructing a sequence of approximate posterior distributions that correspond to classical sequential Bayesian updating without having a single fixed objective to optimize.
As with most approaches to scaling MCMC samplers for Bayesian inference, these minibatch-based variational inference methods depend on model structure. In SVI and scalable BBVI, minibatches map to terms in a factorization of the joint probability. In SVB, minibatches map to a sequence of likelihoods to be incorporated into the variational posterior. Some of these methods further depend on and exploit exponential family and conjugacy structure. SVI is based on complete-data conjugacy, while BBVI was specifically developed for nonconjugate models. SVB is a general framework, but in the conjugate exponential family case the updates can be written in terms of simple updates to natural parameters. A direction for future research might be to develop new methods based on identifying and exploiting some ‘middle ground’ between the structural requirements of SVI and BBVI, or similarly of SVB with and without exponential family structure.
4 Discussion
The minibatch-based variational inference methods developed in this chapter suggest parallel and asynchronous variants. In the case of SVB, distributed and asynchronous versions, such as the master-worker pattern depicted by Algorithms 18 and 19, have been empirically studied Broderick et al. (2013). However, we lack theoretical understanding about these procedures, and it is unclear how to define and track notions of convergence or stability. Methods based on stochastic gradients, such as SVI, can naturally be extended to exploit parallel and asynchronous (or “Hogwild”) variants of stochastic gradient ascent. In such parallel settings, these optimization-based techniques benefit from powerful gradient convergence results (Bertsekas and Tsitsiklis, 1989, Section 7.8), though tuning such algorithms is still a challenge. Other parallel versions of these ideas and algorithms have also been developed in Campbell and How (2014) and Campbell et al. (2015).
Chapter 6 Challenges and questions
In this review, we have examined a variety of different views on scaling Bayesian inference up to large datasets and greater model complexity and out to parallel compute resources. Several different themes have emerged, from techniques that exploit subsets of data for computational savings to proposals for distributing inference computations across multiple machines. Progress is being made, but there remain significant open questions and outstanding challenges to be tackled as this research programme moves forward.
One of the key insights underpinning much of the recent work on scaling Bayesian inference can be framed in terms of a kind of bias-variance tradeoff. Traditional MCMC theory provides asymptotically unbiased estimators for which the error can eventually be driven arbitrarily small. However, in practice, under limited computational budgets the error can be significant. This error has two components: transient bias, in which the samples produced are too dependent on the Markov chain’s initialization, and Monte Carlo standard error, in which the samples collected may be too few or too highly correlated to produce good estimates.
Figure 6.1 illustrates the error regimes and tradeoffs in traditional MCMC.See also Section 2.2.4 Asymptotic analysis describes the regime on the right of the plot, after the sampler has mixed sufficiently well. In this regime, the marginal distribution of each sample is essentially equal to the target distribution, and the transient bias from initialization, which affects only the early samples in the Monte Carlo sum, is washed out rapidly at least at a rate. The dominant source of error is due to Monte Carlo standard error, which diminishes only at a rate.
However, machine learning practitioners using MCMC often find themselves in another regime: in the middle of the plot, the error is decreasing but dominated instead by the transient bias. The challenge in practice is often to get through this regime, or even to get into it at all. When the underlying Markov chain does not mix sufficiently well or when the transitions cannot be computed sufficiently quickly, getting to this regime may be practically infeasible for a realistic computational budget.
Several of the new MCMC techniques we have studied aim to address this challenge. In particular, the parallel predictive prefetching method of Section 4.1.2 accelerates this phase of MCMC without affecting the stationary distribution. Other methods instead introduce approximate transition operators that can be executed more efficiently. For example, the adaptive subsampling methods of Section 3.2 and the stochastic gradient sampler of Section 3.4 can execute updates more efficiently by operating only on data subsets, while the Weierstrass and Hogwild Gibbs samplers of Sections 4.2.1 and 4.2.2, respectively, execute more quickly by leveraging data parallelism. These transition operators are approximate in that they do not admit the exact target distribution as a stationary distribution: instead, the stationary distribution is only intended to be close to the target. Framed in terms of Monte Carlo estimates, these approximations effectively accelerate the execution of the chain at the cost of introducing an asymptotic bias. Figure 6.2 illustrates this new tradeoff.
Allowing some asymptotic bias to reduce transient bias or even Monte Carlo variance is likely to enable MCMC inference at a new scale. However, both the amount of asymptotic bias introduced by these methods and the ways in which it depends on model and algorithm parameters remain unclear. More theoretical understanding and empirical study is necessary to guide machine learning practice.
Scalability in the context of Bayesian inference is ultimately about spending computational resources to better interrogate posterior distributions. It is therefore important to consider whether there are fundamental limits to what can be achieved by, e.g., spending more money on Amazon EC2, for either faster computers or more of them.
In parallel systems, linear scaling is ideal: twice as much computational power yields twice as much useful work. Unfortunately, even if this lofty parallel speedup goal is achieved, the asymptotic picture for MCMC is dim: in the asymptotic regime, doubling the number of samples collected can only reduce the Monte Carlo standard error by a factor of . This scaling means that there are diminishing returns to purchasing additional computational resources, even if those resources provide linear speedup in terms of accelerating the execution of the MCMC algorithm.
Interestingly, variational methods may not suffer from such intrinsic limits. In particular, the stochastic gradient variational inference methods surveyed in Section 5.1 can utilize optimization methods that, at least in smooth convex settings, can converge at least at rates (Bubeck, 2015). When these convergence properties are maintained for achieving local minima in non-convex problems, applying additional computational resources would not inherently suffer from the problem of diminishing marginal returns.
With all the ideas surveyed here, one thing is clear: there are many alternatives for how to scale Bayesian inference. How should we compare these alternative algorithms? Can we tell when any of these algorithms work well in an absolute sense?
One standard approach for evaluating MCMC procedures is to define a set of scalar-valued test functions (or estimands of interest) and compute effective sample size (Gelman et al., 2014, Section 11.5) as a function of wall-clock time. However, in complex models designing an appropriately comprehensive set of test functions may be difficult. Furthermore, many such measures require the Markov chain to mix and do not account for any asyptotic bias (Gorham and Mackey, 2015), hence limiting their applicability to measuring the performance of many of the new inference methods studied here.
To confront these challenges, one recently-proposed approach (Gorham and Mackey, 2015) draws on Stein’s method, classically used as an analytical tool, to design an efficiently-computable measure of discrepancy between a target distribution and a set of samples. A natural measure of discrepancy between a target density and a (weighted) sample distribution , where for some set of samples and weights , is to consider their largest absolute difference across a large class of test functions:
Such operators can be designed using infinitessimal generators from continuous-time ergodic Markov processes, and Gorham and Mackey (2015) suggest using the operator
which requires computing only the gradient of the target log density. Furthermore, while the optimization in (6.2) is infinite-dimensional in general and might have infinitely many smoothness constraints from , Gorham and Mackey (2015) shows that for the sample distribution the test function need only be evaluated at the finitely-many sample points and that only a small number of constraints must be enforced. This new performance metric does not require assumptions on whether the samples are generated from an unbiased, stationary Markov chain, and so it may provide clear ways to compare across a broad spectrum sampling-based approximate inference algorithms.
Another recently-proposed approach attempts to estimate or bound the KL divergence from an algorithm’s approximate posterior representation to the true posterior, at least when applied to synthetic data. This approach, called bidirectional Monte Carlo (BDMC) (Grosse et al., 2015), can be applied to measure the performance of both variational mean field algorithms as well as annealed importance sampling (AIS) and sequential Monte Carlo (SMC) algorithms. By rearranging the variational identity (2.55), we can write the KL divergence from an approximating distribution to a target posterior in terms of the log marginal likelihood and an expectation with respect to :
Because the expectation can be readily computed in a mean field setting or stochastically lower-bounded when using AIS (Grosse et al., 2015, Section 4.1), with a stochastic upper bound on we can use (6.4) to compute a stochastic upper bound on the KL divergence . BDMC provides a method to compute such stochastic upper bounds on for synthetic datasets , and so may enable new performance metrics that apply to both sampling-based algorithms as well as variational mean field algorithms. However, while MCMC transition operators are used to construct AIS algorithms, BDMC does not directly apply to evaluating the performance of such transition operators in standard MCMC inference.
Developing performance metrics and evaluation procedures is critical to making progress. As observed in Grosse et al. (2015),
In many application areas of machine learning, especially supervised learning, benchmark datasets have spurred rapid progress in developing new algorithms and clever refinements to existing algorithms. […] So far, the lack of quantitative performance evaluations in marginal likelihood estimation, and in sampling-based inference more generally, has left us fumbling around in the dark.
By developing better ways to measure the performance of these Bayesian inference algorithms, we will be much better equipped to compare, improve, and extend them.
Acknowledgements
This work was funded in part by NSF IIS-1421780 and the Alfred P. Sloan Foundation. E.A. is supported by the Miller Institute for Basic Research in Science, University of California, Berkeley. M.J. is supported by a fellowship from the Harvard/MIT Joint Grants program.