Variational consensus Monte Carlo

06/09/2015 ∙ by Maxim Rabinovich, et al. ∙ berkeley college 0

Practitioners of Bayesian statistics have long depended on Markov chain Monte Carlo (MCMC) to obtain samples from intractable posterior distributions. Unfortunately, MCMC algorithms are typically serial, and do not scale to the large datasets typical of modern machine learning. The recently proposed consensus Monte Carlo algorithm removes this limitation by partitioning the data and drawing samples conditional on each partition in parallel (Scott et al, 2013). A fixed aggregation function then combines these samples, yielding approximate posterior samples. We introduce variational consensus Monte Carlo (VCMC), a variational Bayes algorithm that optimizes over aggregation functions to obtain samples from a distribution that better approximates the target. The resulting objective contains an intractable entropy term; we therefore derive a relaxation of the objective and show that the relaxed problem is blockwise concave under mild conditions. We illustrate the advantages of our algorithm on three inference tasks from the literature, demonstrating both the superior quality of the posterior approximation and the moderate overhead of the optimization step. Our algorithm achieves a relative error reduction (measured against serial MCMC) of up to 39 of estimating 300-dimensional probit regression parameter expectations; similarly, it achieves an error reduction of 92 cluster comembership probabilities in a Gaussian mixture model with 8 components in 8 dimensions. Furthermore, these gains come at moderate cost compared to the runtime of serial MCMC, achieving near-ideal speedup in some instances.



There are no comments yet.


page 1

page 2

page 3

page 4

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Modern statistical inference demands scalability to massive datasets and high-dimensional models. Innovation in distributed and stochastic optimization has enabled parameter estimation in this setting, e.g. via stochastic [3] and asynchronous [20]

variants of gradient descent. Achieving similar success in Bayesian inference – where the target is a posterior distribution over parameter values, rather than a point estimate – remains computationally challenging.

Two dominant approaches to Bayesian computation are variational Bayes and Markov chain Monte Carlo (MCMC). Within the former, scalable algorithms like stochastic variational inference [11] and streaming variational Bayes [4] have successfully imported ideas from optimization. Within MCMC, adaptive subsampling procedures [2, 14], stochastic gradient Langevin dynamics [25], and Firefly Monte Carlo [16] have applied similar ideas, achieving computational gains by operating only on data subsets. These algorithms are serial, however, and thus cannot take advantage of multicore and multi-machine architectures. This motivates data-parallel MCMC algorithms such as asynchronous variants of Gibbs sampling [1, 8, 12].

Our work belongs to a class of communication-avoiding data-parallel MCMC algorithms. These algorithms partition the full dataset  into  disjoint subsets  where  denotes the data associated with core . Each core samples from a subposterior distribution,


and then a centralized procedure combines the samples into an approximation of the full posterior. Due to their efficiency, such procedures have recently received substantial attention [18, 22, 24].

One of these algorithms, consensus Monte Carlo (CMC), requires communication only at the start and end of sampling [22]. CMC proceeds from the intuition that subposterior samples, when aggregated correctly, can approximate full posterior samples. This is formally backed by the factorization


If one can approximate the subposterior densities

, using kernel density estimates for instance 

[18], it is therefore possible to recombine them into an estimate of the full posterior.

Unfortunately, the factorization does not make it immediately clear how to aggregate on the level of samples without first having to obtain an estimate of the densities themselves. CMC alters (2) to untie the parameters across partitions and plug in a deterministic link from the  to :


This approximation and an aggregation function motivated by a Gaussian approximation lie at the core of the CMC algorithm [22].

The introduction of CMC raises numerous interesting questions whose answers are essential to its wider application. Two among these stand out as particularly vital. First, how should the aggregation function be chosen to achieve the closest possible approximation to the target posterior? Second, when model parameters exhibit structure or must conform to constraints — if they are, for example, positive semidefinite covariance matrices or labeled centers of clusters — how can the weighted averaging strategy of Scott et al. [22] be modified to account for this structure?

In this paper, we propose variational consensus Monte Carlo (VCMC), a novel class of data-parallel MCMC algorithms that allow both questions to be addressed. By formulating the choice of aggregation function as a variational Bayes problem, VCMC makes it possible to adaptively choose the aggregation function to achieve a closer approximation to the true posterior. The flexibility of VCMC likewise supports nonlinear aggregation functions, and we exploit this versatility to introduce structured aggregation functions applicable to inference problems that are not purely vectorial.

An appealing benefit of the VCMC point of view is a clarification of the untying step leading to (3

). In VCMC, the approximate factorization corresponds to a variational approximation to the true posterior. This approximation can be viewed as the joint distribution of

and  in an augmented model that assumes conditional independence between the data partitions and posits a deterministic mapping from partition-level parameters to the single global parameter. The added flexibility of this point-of-view makes it possible to move beyond subposteriors and include alternative forms of (3) within the CMC framework. In particular, it is possible to define , using partial posteriors in place of subposteriors (cf. [23]). Although extensive investigation of this issue is beyond the scope of this paper, we provide some evidence in Section 6 that partial posteriors are a better choice in some circumstances and demonstrate that VCMC can provide substantial gains in both the partial posterior and subposterior settings.

Before proceeding, we outline the remainder of this paper. Below, in §2, we review CMC and related data-parallel MCMC algorithms. Next, we cast CMC as a variational Bayes problem in §3. We define the variational optimization objective in §4, addressing the challenging entropy term by relaxing it to a concave lower bound, and give conditions for which this leads to a blockwise concave maximization problem. In §5, we define several classes of aggregation functions, and design novel classes that enable structured aggregation of constrained and structured samples—e.g. positive semidefinite matrices and mixture model parameters. In §6, we evaluate the performance of VCMC and CMC relative to serial MCMC. We replicate experiments carried out by Scott et al. [22] and execute more challenging experiments in higher dimensions and with more data. Finally in §7, we summarize our approach and discuss several open problems generated by this work.

2 Related work

We focus on data-parallel MCMC algorithms for large-scale Bayesian posterior sampling. Several recent research threads propose schemes in the setting where the posterior factors as in (2). In general, these parallel strategies are approximate relative to serial procedures, and the specific algorithms differ in terms of the approximations employed and amount of communication required.

At one end of the communication spectrum are algorithms that fit into the MapReduce model [7]. First, parallel cores sample from subposteriors, defined in (1), via any Monte Carlo sampling procedure. The subposterior samples are then aggregated to obtain approximate samples from the full posterior. This leads to the challenge of designing proper and efficient aggregation procedures.

Scott et al. [22] propose consensus Monte Carlo (CMC), which constructs approximate posterior samples via weighted averages of subposterior samples; our algorithms are motivated by this work. Let denote the -th subposterior sample from core . In CMC, the aggregation function averages across each set of  samples  to produce one approximate posterior sample 

. Uniform averaging is a natural but naïve heuristic that can in fact be improved upon via a weighted average,


where in general, 

is a vector and 

can be a matrix. The authors derive weights motivated by the special case of a Gaussian posterior, where each subposterior is consequently also Gaussian. Let  be the covariance of the -th subposterior. This suggests weights  equal to the subposteriors’ inverse covariances. CMC treats arbitrary subpostertiors as Gaussians, aggregating with weights given by empirical estimates of  computed from the observed subposterior samples.

Neiswanger et al. [18] propose aggregation at the level of distributions rather than samples. Here, the idea is to form an approximate posterior via a product of density estimates fit to each subposterior, and then sample from this approximate posterior. The accuracy and computational requirements of this approach depend on the complexity of these density estimates. Wang and Dunson [24] develop alternate data-parallel MCMC methods based on applying a Weierstrass transform to each subposterior. These Weierstrass sampling procedures introduce auxiliary variables and additional communication between computational cores.

3 Consensus Monte Carlo as variational inference

Given the distributional form of the CMC framework (3), we would like to choose  so that the induced distribution on  is as close as possible to the true posterior. This is precisely the problem addressed by variational Bayes, which approximates an intractable posterior by the solution to the constrained optimization problem

where is the family of variational approximations to the distribution, usually chosen to make both optimization and evaluation of target expectations tractable. We thus view the aggregation problem in CMC as a variational inference problem, with the variational family given by all distributions , where each is in some function class and defines a density

In practice, we parameterize

by a finite dimensional set of parameters and optimize over it using projected stochastic gradient descent (SGD).

4 The variational optimization problem

Standard optimization of the variational Bayes objective uses the evidence lower bound (ELBO)


We can therefore recast the variational optimization problem in an equivalent form as

Unfortunately, the variational Bayes objective remains difficult to optimize. Indeed, by writing

we see that optimizing requires computing an entropy and its gradients. We can deal with this issue by deriving a lower bound on the entropy that relaxes the objective further.

Concretely, suppose that every can be decomposed as , with each a differentiable bijection. Since the come from subposteriors conditioning on different segments of the data, they are independent. The entropy power inequality [6] therefore implies


where denotes the Jacobian of the function evaluated at . The proof of this inequality can be found in the supplement.

This approach gives an explicit, easily computed approximation to the entropy—and this approximation is a lower bound, allowing us to interpret it simply as a further relaxation of the original inference problem. Furthermore, and crucially, it decouples and , thereby making it possible to optimize over without estimating the entropy of any . We note additionally that if we are willing to sacrifice concavity, we can use the tighter lower bound on the entropy given by (6).

Putting everything together, we can define our relaxed variational objective as


Maximizing this function is the variational Bayes problem we consider in the remainder of the paper.

Conditions for concavity

Under certain conditions, the problem posed above is blockwise concave. To see when this holds, we use the language of graphical models and exponential families. To derive the result in the greatest possible generality, we decompose the variational objective as

and prove concavity directly for , then treat our choice of relaxed entropy (7). We emphasize that while the entropy relaxation is only defined for decomposed aggregation functions, concavity of the partial objective holds for arbitrary aggregation functions. All proofs are in the supplement.

Suppose the model distribution is specified via a graphical model , so that , such that each conditional distribution is defined by an exponential family

If each of these log conditional density functions is log-concave in , we can guarantee that the log likelihood is concave in each individually.

Theorem 4.1 (Blockwise concavity of the variational cross-entropy).

Suppose that the model distribution is specified by a graphical model in which each conditional probability density is a log-concave exponential family. Suppose further that the variational aggregation function family satisfies such that we can decompose each aggregation function across nodes via

If each is a convex subset of some vector space , then the variational cross-entropy is concave in each individually.

Assuming that the aggregation function can be decomposed into a sum over functions of individual subposterior terms we can also prove concavity of our entropy relaxation (7).

Theorem 4.2 (Concavity of the relaxed entropy).

Suppose , with each function decomposing as for unique bijective . Then the relaxed entropy (7) is concave in .

As a result, we derive concavity of the variational objective in a broad range of settings.

Corollary 4.1 (Concavity of the variational objective).

Under the hypotheses of Theorems 4.1 and 4.2, the variational Bayes objective is concave in each individually.

5 Variational aggregation function families

The performance of our algorithm depends critically on the choice of aggregation function family . The family must be sufficiently simple to support efficient optimization, expressive to capture the complex transformation from the set of subposteriors to the full posterior, and structured to preserve structure in the parameters. We now illustrate some aggregation functions that meet these criteria.

5.1 Vector aggregation

In the simplest case, is an unconstrained vector. Then, a linear aggregation function makes sense, and it is natural to impose constraints to make this sum behave like a weighted average—i.e., each is a positive semidefinite (PSD) matrix and . For computational reasons, it is often desirable to restrict to diagonal .

5.2 Spectral aggregation

Cases involving structure exhibit more interesting behavior. Indeed, if our parameter is a PSD matrix , applying the vector aggregation function above to the flattened vector form of the parameter does not suffice. Denoting elementwise matrix product as , we note that this strategy would in general lead to .

We therefore introduce a more sophisticated aggregation function that preserves PSD structure. For this, given symmetric , define and to be orthogonal and diagonal matrices, respectively, such that . Impose further—and crucially—the canonical ordering . We can then define our spectral aggregation function by

Assuming , the output of this function is guaranteed to be PSD, as required. As above we restrict the set of  to the matrix simplex .

5.3 Combinatorial aggregation

Additional complexity arises with unidentifiable latent variables and, more generally, models with multimodal posteriors. Since this class encompasses many popular algorithms in machine learning, including factor analysis, mixtures of Gaussians and multinomials, and latent Dirichlet allocation (LDA), we now show how our framework can accommodate them.

For concreteness, suppose now that our model parameters are given by , where  denotes the number of global latent variables (e.g. cluster centers). We introduce discrete alignment parameters  that indicate how latent variables associated with partitions map to global latent variables. Each  is thus a one-to-one correspondence , with  denoting the index on worker core  of cluster center . For fixed , we then obtain the variational aggregation function

Optimization can then proceed in an alternating manner, switching between the alignments  and the weights , or in a greedy manner, fixing the alignments at the start and optimizing the weight matrices. In practice, we do the latter, aligning using a simple heuristic objective , where denotes the mean value of cluster center on partition . As  suggests, we set . This is permitted because the global order is only determined up to a permutation. We find that minimizing  via the Hungarian algorithm [15] leads to good alignments.

6 Empirical evaluation

We now evaluate VCMC on three inference problems, in a range of data and dimensionality conditions. In the vector parameter case, we compare directly to the simple weighting baselines corresponding to previous work on CMC [22]; in the other cases, we compare to structured analogues of these weighting schemes. Our experiments demonstrate the advantages of VCMC across the whole range of model dimensionality, data quantity, and availability of parallel resources.

Baseline weight settings

Scott et al. [22] studied linear aggregation functions with fixed weights,


corresponding to uniform averaging and Gaussian averaging, respectively, where denotes the standard empirical estimate of the covariance. These are our baselines for comparison.

Evaluation metrics

Since the goal of MCMC is usually to estimate event probabilities and function expectations, we evaluate algorithm accuracy for such estimates, relative to serial MCMC output. For each model, we consider a suite of test functions (e.g. low degree polynomials, cluster comembership indicators), and we assess the error of each algorithm  using the metric

In the body of the paper, we report median values of

, computed within each test function class. The supplement expands on this further, showing quartiles for the differences in

and .

6.1 Bayesian probit regression

Figure 1: High-dimensional probit regression

. Moment approximation error for the uniform and Gaussian averaging baselines and VCMC, relative to serial MCMC, for

(left) subposteriors and (right) partial posteriors. We assessed three groups of functions: first moments, with for ; pure second moments, with for ; and mixed second moments, with for . For brevity, results for pure second moments are relegated to Figure 6 in the supplement. Relative errors are truncated to in all cases.

We consider the nonconjugate probit regression model. In this case, we use linear aggregation functions as our function class. For computational efficiency, we also limit ourselves to diagonal . We use Gibbs sampling on the following augmented model:

This augmentation allows us to implement an efficient and rapidly mixing Gibbs sampler, where

We run two experiments: the first using a data generating distribution from Scott et al. [22], with data points and dimensions, and the second using data points and dimensions. As shown in Figure 1 and, in the supplement,111Due to space constraints, we relegate results for to the supplement. Figures 5 and 6, VCMC decreases the error of moment estimation compared to the baselines, with substantial gains starting at partitions (and increasing with ). We also run the high-dimensional experiment using partial posteriors [23] in place of subposteriors, and observe substantially lower errors in this case.

6.2 Normal-inverse Wishart model

Figure 2: High-dimensional normal-inverse Wishart model (). (Far left, left, right) Moment approximation error for the uniform and Gaussian averaging baselines and VCMC, relative to serial MCMC. Letting denote the

largest eigenvalue of

, we assessed three groups of functions: first moments, with for ; pure second moments, with for ; and mixed second moments, with for . (Far right) Graph of error in estimating as a function of (where ).

To compare directly to prior work [22], we consider the normal-inverse Wishart model

Here, we use spectral aggregation rules as our function class, restricting to diagonal  for computational efficiency. We run two sets of experiments: one using the covariance matrix from Scott et al. [22], with data points and dimensions, and one using a higher-dimensional covariance matrix designed to have a small spectral gap and a range of eigenvalues, with data points and dimensions. In both cases, we use a form of projected SGD, using  samples per iteration to estimate the variational gradients and running  iterations of optimization. We note that because the mean  is treated as a point-estimated parameter, one could sample  exactly using normal-inverse Wishart conjugacy [10]. As Figure 2 shows,222Due to space constraints, we compare to the experiment of Scott et al. [22] in the supplement. VCMC improves both first and second posterior moment estimation as compared to the baselines. Here, the greatest gains from VCMC appear at large numbers of partitions (

). We also note that uniform and Gaussian averaging perform similarly because the variances do not differ much across partitions.

6.3 Mixture of Gaussians

Figure 3: Mixture of Gaussians . Expectation approximation error for the uniform and Gaussian baselines and VCMC. We report the median error, relative to serial MCMC, for cluster comembership probabilities of pairs of test data points, for (left) and (right) , where we run the VCMC optimization procedure for 50 and 200 iterations, respectively. When , some comembership probabilities are estimated poorly by all methods; we therefore only use the 70% of comembership probabilities with the smallest errors across all the methods.

A substantial portion of Bayesian inference focuses on latent variable models and, in particular, mixture models. We therefore evaluate VCMC on the mixture of Gaussians model defined by

where the mixture weights and the prior and likelihood variances and are assumed known. We use the combinatorial aggregation functions defined in Section 5; we set , , , and uniform and generate data points in dimensions, using the model from Nishihara et al. [19]. The resulting inference problem is therefore -dimensional. All samples were drawn using the PyStan implementation of Hamiltonian Monte Carlo (HMC).

As Figure 3 shows, VCMC drastically improves moment estimation compared to the baseline Gaussian averaging (9). To assess how VCMC influences estimates in cluster membership probabilities, we generated  new test points from the model and analyzed cluster comembership probabilities for all pairs in the test set. Concretely, for each and in the test data, we estimated . Figure 3 shows the resulting boost in accuracy: when , VCMC delivers estimates close to those of serial MCMC, across all numbers of partitions; the errors are larger for . Unlike previous models, uniform averaging here outperforms Gaussian averaging, and indeed is competitive with VCMC.

6.4 Assessing computational efficiency

The efficiency of VCMC depends on that of the optimization step, which depends on factors including the step size schedule, number of samples used per iteration to estimate gradients, and size of data minibatches used per iteration. Extensively assessing the influence of all these factors is beyond the scope of this paper, and is an active area of research both in general and specifically in the context of variational inference [13, 17, 21]. Here, we provide an initial assessment of the computational efficiency of VCMC, taking the probit regression and Gaussian mixture models as our examples, using step sizes and sample numbers from above, and eschewing minibatching on data points.

Figure 4 shows timing results for both models. For the probit regression, while the optimization cost is not negligible, it is significantly smaller than that of serial sampling, which takes over seconds to produce effective samples.333We ran the sampler for iterations, including burnin steps, and kept every fifth sample. Across most numbers of partitions, approximately iterations—corresponding to less than seconds of wall clock time—suffices to give errors close to those at convergence. For the mixture, on the other hand, the computational cost of optimization is minimal compared to serial sampling. We can see this in the overall speedup of VCMC relative to serial MCMC: for sampling and optimization combined, low numbers of partitions achieve speedups close to the ideal value of , and large numbers still achieve good speedups of about . The cost of the VCMC optimization step is thus moderate—and, when the MCMC step is expensive, small enough to preserve the linear speedup of embarrassingly parallel sampling. Moreover, since the serial bottleneck is an optimization, we are optimistic that performance, both in terms of number of iterations and wall clock time, can be significantly increased by using techniques like data minibatching [9], adaptive step sizes [21], or asynchronous updates [20].

Figure 4: Error versus timing and speedup measurements. (Left) VCMC error as a function of number of seconds of optimization. The cost of optimization is nonnegligible, but still moderate compared to serial MCMC—particularly since our optimization scheme only needs small batches of samples and can therefore operate concurrently with the sampler. (Right) Error versus speedup relative to serial MCMC, for both CMC with Gaussian averaging (small markers) and VCMC (large markers).

7 Conclusion and future work

The flexibility of variational consensus Monte Carlo (VCMC) opens several avenues for further research. Following previous work on data-parallel MCMC, we used the subposterior factorization. Our variational framework can accomodate more general factorizations that might be more statistically or computationally efficient – e.g. the factorization used by Broderick et al. [4]. We also introduced structured sample aggregation, and analyzed some concrete instantiations. Complex latent variable models would require more sophisticated aggregation functions – e.g. ones that account for symmetries in the model [5] or lift the parameter to a higher dimensional space before aggregating. Finally, recall that our algorithm – again following previous work – aggregates in a sample-by-sample manner, cf. (4). Other aggregation paradigms may be useful in building approximations to multimodal posteriors or in boosting the statistical efficiency of the overall sampler.


We thank R.P. Adams, N. Altieri, T. Broderick, R. Giordano, M.J. Johnson, and S.L. Scott for for helpful discussions. Support for this project was provided by ONR under the Multidisciplinary University Research Initiative (MURI) program (N00014-11-1-0688). E.A. is supported by the Miller Institute for Basic Research in Science, University of California, Berkeley. M.R. is supported by a Fannie and John Hertz Foundation Fellowship and an NSF Graduate Research Fellowship. This research is supported in part by NSF CISE Expeditions Award CCF-1139158, LBNL Award 7076018, and DARPA XData Award FA8750-12-2-0331, and gifts from Amazon Web Services, Google, SAP, The Thomas and Stacey Siebel Foundation, Adatao, Adobe, Apple, Inc., Blue Goji, Bosch, C3Energy, Cisco, Cray, Cloudera, EMC2, Ericsson, Facebook, Guavus, HP, Huawei, Informatica, Intel, Microsoft, NetApp, Pivotal, Samsung, Schlumberger, Splunk, Virdata and VMware.


  • Asuncion et al. [2008] A. U. Asuncion, P. Smyth, and M. Welling. Asynchronous distributed learning of topic models. In Advances in Neural Information Processing Systems 21, NIPS ’08, pages 81–88, 2008.
  • Bardenet et al. [2014] R. Bardenet, A. Doucet, and C. Holmes. Towards scaling up Markov chain Monte Carlo: An adaptive subsampling approach. In Proceedings of the 31st International Conference on Machine Learning, ICML ’14, 2014.
  • Bertsekas [1990] D. P. Bertsekas. Nonlinear Programming. Athena Scientific, Belmont, MA, 2nd edition, 1990.
  • Broderick et al. [2013] T. Broderick, N. Boyd, A. Wibisono, A. C. Wilson, and M. I. Jordan. Streaming variational Bayes. In Advances in Neural Information Processing Systems 26, NIPS ’13, pages 1727–1735, 2013.
  • Campbell and How [2014] T. Campbell and J. P. How. Approximate decentralized Bayesian inference. In

    30th Conference on Uncertainty in Artificial Intelligence

    , UAI ’14, 2014.
  • Cover and Thomas [2006] T. M. Cover and J. A. Thomas. Elements of Information Theory (Wiley Series in Telecommunications and Signal Processing). Wiley-Interscience, 2006.
  • Dean and Ghemawat [2008] J. Dean and S. Ghemawat. MapReduce: Simplified data processing on large clusters. Communications of the ACM, 51(1):107–113, Jan. 2008.
  • Doshi-Velez et al. [2009] F. Doshi-Velez, D. A. Knowles, S. Mohamed, and Z. Ghahramani. Large scale nonparametric Bayesian inference: Data parallelisation in the Indian buffet process. In Advances in Neural Information Processing Systems 22, NIPS ’09, pages 1294–1302, 2009.
  • Duchi et al. [2011] J. C. Duchi, E. Hazan, and Y. Singer. Adaptive subgradient methods for online learning and stochastic optimization. Journal of Machine Learning Research, 12:2121–2159, 2011.
  • Gelman et al. [2013] A. Gelman, J. B. Carlin, H. S. Stern, D. B. Dunson, A. Vehtari, and D. B. Rubin. Bayesian Data Analysis, Third Edition. Chapman and Hall/CRC, 2013.
  • Hoffman et al. [2013] M. D. Hoffman, D. M. Blei, C. Wang, and J. Paisley. Stochastic variational inference. Journal of Machine Learning Research, 14(1):1303–1347, May 2013.
  • Johnson et al. [2013] M. Johnson, J. Saunderson, and A. Willsky. Analyzing Hogwild parallel Gaussian Gibbs sampling. In Advances in Neural Information Processing Systems 26, pages 2715–2723, 2013.
  • Johnson and Zhang [2013] R. Johnson and T. Zhang. Accelerating stochastic gradient descent using predictive variance reduction. In Advances in Neural Information Processing Systems 26, NIPS ’13, pages 315–323, 2013.
  • Korattikara et al. [2014] A. Korattikara, Y. Chen, and M. Welling. Austerity in MCMC land: Cutting the Metropolis-Hastings budget. In Proceedings of the 31st International Conference on Machine Learning, ICML ’14, 2014.
  • Kuhn [1955] H. W. Kuhn. The Hungarian method for the assignment problem. Naval Research Logistics Quarterly, 2(1-2):83–97, 1955.
  • Maclaurin and Adams [2014] D. Maclaurin and R. P. Adams. Firefly Monte Carlo: Exact MCMC with subsets of data. In Proceedings of 30th Conference on Uncertainty in Artificial Intelligence, UAI ’14, 2014.
  • Mandt and Blei [2014] S. Mandt and D. M. Blei. Smoothed gradients for stochastic variational inference. In Advances in Neural Information Processing Systems 27, NIPS ’14, pages 2438–2446, 2014.
  • Neiswanger et al. [2014] W. Neiswanger, C. Wang, and E. Xing. Asymptotically exact, embarrassingly parallel MCMC. In 30th Conference on Uncertainty in Artificial Intelligence, UAI ’14, 2014.
  • Nishihara et al. [2014] R. Nishihara, I. Murray, and R. P. Adams. Parallel MCMC with generalized elliptical slice sampling. Journal of Machine Learning Research, 15:2087–2112, 2014.
  • Niu et al. [2011] F. Niu, B. Recht, C. Ré, and S. Wright. Hogwild!: A lock-free approach to parallelizing stochastic gradient descent. In Advances in Neural Information Processing Systems 24, NIPS ’11, pages 693–701, 2011.
  • Ranganath et al. [2013] R. Ranganath, C. Wang, D. M. Blei, and E. P. Xing. An adaptive learning rate for stochastic variational inference. In Proceedings of the 30th International Conference on Machine Learning, ICML ’13, pages 298–306, 2013.
  • Scott et al. [2013] S. L. Scott, A. W. Blocker, and F. V. Bonassi. Bayes and big data: The consensus Monte Carlo algorithm. In Bayes 250, 2013.
  • Strathmann et al. [2015] H. Strathmann, D. Sejdinovic, and M. Girolami. Unbiased Bayes for big data: Paths of partial posteriors. arXiv preprint arXiv:1501.03326, 2015.
  • Wang and Dunson [2013] X. Wang and D. B. Dunson. Parallel MCMC via Weierstrass sampler. arXiv preprint arXiv:1312.4605, 2013.
  • Welling and Teh [2011] M. Welling and Y. W. Teh. Bayesian learning via stochastic gradient Langevin dynamics. In Proceedings of the 28th International Conference on Machine Learning, ICML ’11, 2011.

Appendix A Proofs

Proof of entropy relaxation.

We apply the entropy power inequality [6], which asserts that for independent -dimensional random vectors , the sum



where denotes differential entropy.

In our case, we have



equation (10) implies


we immediately see that

as required. ∎

Proof of Theorem 4.1.

We first define

Since , where the expectation is taken with respect to the subposteriors, which do not vary with , it suffices to show that is concave in each individually for each fixed . Furthermore, since is linear in by the definition of function addition, it actually suffices to show in each individually. To see why this holds, first observe that for each , we have


where is a function of that is constant in . By the log-concavity assumption, the sum of the first two terms of  in (12) is concave in . On the other hand, by basic properties of exponential families, each is convex in and hence in , making its negative concave. Since the remaining terms are linear or constant, is in fact concave in . The claim follows. ∎

Proof of Theorem 4.2.

Clearly it suffices to show that each is concave and for this it suffices to show that for fixed , is concave. This is immediate, however, since the Jacobian is a linear function and is a concave function. ∎

Appendix B Variational objective functions

We derive the variational objectives and gradients for the models we analyze. Throughout, we make the convention that for ,

denotes the trace inner product.

b.1 Bayesian probit regression

In this section, we compute the variational objective for the Bayesian probit regression model. For convenience, we define

In this notation, the variational objective takes the simple form

where .

This leads to the gradients

where we have additionally defined and

b.2 Normal-inverse Wishart model

The variational objective for the normal-inverse Wishart model takes the form


and we have compressed our notation by setting , , , and . As before, we have

where we have suppressed the constant depending on the since it does not vary with .

Recalling that is diagonal, we can obtain the gradients by first computing

where we have used to denote elementwise vector products. We then find

b.3 Mixture of Gaussians

Per the description of aggregation in Section 5, we define merged samples in the mixture of Gaussians model by the equations

where denotes the cluster index and denotes the alignment mapping indices on the master core to indices on worker core . Throughout this section, we treat the alignment variables as fixed.

Using this notation, we define



The variational objective then takes the form

with the usual equation

Some calculation then shows that the gradients with respect to the various are given by


This covers the case of general PSD matrices . When the matrices are restricted to be diagonal, we get the simplified gradient

where denotes elementwise multiplication of vectors.


this gives us all the information we need to implement an optimization procedure for the objective.

Appendix C Extended empirical evaluation

Figure 5: Five-dimensional probit regression . Moment approximation error for the uniform and Gaussian averaging baselines and VCMC, relative to serial MCMC. We assessed three groups of functions: (left) first moments, with for ; (center) pure second moments, with for ; and (right) mixed second moments, with for .
Figure 6: High-dimensional probit regression (). Moment approximation error for the uniform and Gaussian averaging baselines and VCMC, relative to serial MCMC, for subposteriors (left) and partial posteriors (right). Here we show the pure second moments.
Figure 7: Five-dimensional normal-inverse Wishart model (). Moment approximation error for the uniform and Gaussian averaging baselines and VCMC, relative to serial MCMC. Letting denote the largest eigenvalue of , we assessed three groups of functions: (left) first moments, with for ; (center) pure second moments, with for ; and (right) mixed second moments, with for .