Asymptotically Exact, Embarrassingly Parallel MCMC
Communication costs, resulting from synchronization requirements during learning, can greatly slow down many parallel machine learning algorithms. In this paper, we present a parallel Markov chain Monte Carlo (MCMC) algorithm in which subsets of data are processed independently, with very little communication. First, we arbitrarily partition data onto multiple machines. Then, on each machine, any classical MCMC method (e.g., Gibbs sampling) may be used to draw samples from a posterior distribution given the data subset. Finally, the samples from each machine are combined to form samples from the full posterior. This embarrassingly parallel algorithm allows each machine to act independently on a subset of the data (without communication) until the final combination stage. We prove that our algorithm generates asymptotically exact samples and empirically demonstrate its ability to parallelize burn-in and sampling in several models.
Markov chain Monte Carlo (MCMC) methods are popular tools for performing approximate Bayesian inference via posterior sampling. One major benefit of these techniques is that they guarantee asymptotically exact recovery of the posterior distribution as the number of posterior samples grows. However, MCMC methods may take a prohibitively long time, since for data points, most methods must perform operations to draw a sample. Furthermore, MCMC methods might require a large number of “burn-in” steps before beginning to generate representative samples. Further complicating matters is the issue that, for many big data applications, it is necessary to store and process data on multiple machines, and so MCMC must be adapted to run in these data-distributed settings.
Researchers currently tackle these problems independently, in two primary ways. To speed up sampling, multiple independent chains of MCMC can be run in parallel [20, 11, 13]; however, each chain is still run on the entire dataset, and there is no speed-up of the burn-in process (as each chain must still complete the full burn-in before generating samples). To run MCMC when data is partitioned among multiple machines, each machine can perform computation that involves a subset of the data and exchange information at each iteration to draw a sample [10, 14, 18]; however, this requires a significant amount of communication between machines, which can greatly increase computation time when machines wait for external information [1, 7].
We aim to develop a procedure to tackle both problems simultaneously, to allow for quicker burn-in and sampling in settings where data are partitioned among machines. To accomplish this, we propose the following: on each machine, run MCMC on only a subset of the data (independently, without communication between machines), and then combine the samples from each machine to algorithmically construct samples from the full-data posterior distribution. We’d like our procedure to satisfy the following four criteria:
Each machine only has access to a portion of the data.
Each machine performs MCMC independently, without communicating (i.e. the procedure is “embarrassingly parallel”).
Each machine can use any type of MCMC to generate samples.
The combination procedure yields provably asymptotically exact samples from the full-data posterior.
The third criterion allows existing MCMC algorithms or software packages to be run directly on subsets of the data—the combination procedure then acts as a post-processing step to transform the samples to the correct distribution. Note that this procedure is particularly suitable for use in a MapReduce  framework. Also note that, unlike current strategies, this procedure does not involve multiple “duplicate” chains (as each chain uses a different portion of the data and samples from a different posterior distribution), nor does it involve parallelizing a single chain (as there are multiple chains operating independently). We will show how this allows our method to, in fact, parallelize and greatly reduce the time required for burn-in.
In this paper we will (1) introduce and define the subposterior density—a modified posterior given a subset of the data—which will be used heavily, (2) present methods for the embarrassingly parallel MCMC and combination procedure, (3) prove theoretical guarantees about the samples generated from our algorithm, (4) describe the current scope of the presented method (i.e. where and when it can be applied), and (5) show empirical results demonstrating that we can achieve speed-ups for burn-in and sampling while meeting the above four criteria.
2. Embarrassingly Parallel MCMC
The basic idea behind our method is to partition a set of i.i.d. data points into subsets, sample from the subposterior—the posterior given a data subset with an underweighted prior—in parallel, and then combine the resulting samples to form samples from the full-data posterior , where and .
More formally, given data partitioned into subsets , the procedure is:
For (in parallel):
Sample from the subposterior , where
Combine the subposterior samples to produce samples from an estimate of the subposterior density product , which is proportional to the full-data posterior, i.e. .
We want to emphasize that we do not need to iterate over these steps and the combination stage (step 3) is the only step that requires communication between machines. Also note that sampling from each subposterior (step 2) can typically be done in the same way as one would sample from the full-data posterior. For example, when using the Metropolis-Hastings algorithm, one would compute the likelihood ratio as instead of , where is the proposed move. In the next section, we show how the combination stage (step 3) is carried out to generate samples from the full-data posterior using the subposterior samples.
3. Combining Subposterior Samples
Our general idea is to combine the subposterior samples in such a way that we are implicitly sampling from an estimate of the subposterior density product function . If our density product estimator is consistent, then we can show that we are drawing asymptotically exact samples from the full posterior. Further, by studying the estimator error rate, we can explicitly analyze how quickly the distribution from which we are drawing samples is converging to the true posterior (and thus compare different combination algorithms).
In the following three sections we present procedures that yield samples from different estimates of the density product. Our first example is based on a simple parametric estimator motivated by the Bernstein-von Mises theorem ; this procedure generates approximate (asymptotically biased) samples from the full posterior. Our second example is based on a nonparametric estimator, and produces asymptotically exact samples from the full posterior. Our third example is based on a semiparametric estimator, which combines beneficial aspects from the previous two estimators while also generating asymptotically exact samples.
3.1. Approximate posterior sampling with a parametric estimator
The first method for forming samples from the full posterior given subposterior samples involves using an approximation based on the Bernstein-von Mises (Bayesian central limit) theorem, an important result in Bayesian asymptotic theory. Assuming that a unique, true data-generating model exists and is denoted , this theorem states that the posterior tends to a normal distribution concentrated around as the number of observations grows. In particular, under suitable regularity conditions, the posterior is well approximated by (where is the fisher information of the data) when is large . Since we aim to perform posterior sampling when the number of observations is large, a normal parametric form often serves as a good posterior approximation. A similar approximation was used in  in order to facilitate fast, approximately correct sampling. We therefore estimate each subposterior density with where and are the sample mean and covariance, respectively, of the subposterior samples. The product of the subposterior densities will be proportional to a Gaussian pdf, and our estimate of the density product function is
where the parameters of this distribution are
These parameters can be computed quickly and, if desired, online (as new subposterior samples arrive).
3.2. Asymptotically exact posterior sampling with nonparametric density product estimation
In the previous method we made a parametric assumption based on the Bernstein-von Mises theorem, which allows us to generate approximate samples from the full posterior. Although this parametric estimate has quick convergence, it generates asymptotically biased samples, especially in cases where the posterior is particularly non-Gaussian. In this section, we develop a procedure that implicitly samples from the product of nonparametric density estimates, which allows us to produce asymptotically exact samples from the full posterior. By constructing a consistent density product estimator from which we can generate samples, we ensure that the distribution from which we are sampling converges to the full posterior.
Given samples111For ease of description, we assume each machine generates the same number of samples, . In practice, they do not have to be the same. from a subposterior , we can write the kernel density estimator as,
where we have used a Gaussian kernel with bandwidth parameter . After we have obtained the kernel density estimator for subposteriors, we define our nonparametric density product estimator for the full posterior as
This estimate is the probability density function (pdf) of a mixture of Gaussians with unnormalized mixture weights . Here, we use to denote the set of indices for the samples (each from a separate machine) associated with a given mixture component, and we define
Although there are possible mixture components, we can efficiently generate samples from this mixture by first sampling a mixture component (based on its unnormalized component weight ) and then sampling from this (Gaussian) component. In order to sample mixture components, we use an independent Metropolis within Gibbs (IMG) sampler. This is a form of MCMC, where at each step in the Markov chain, a single dimension of the current state is proposed (i.e. sampled) independently of its current value (while keeping the other dimensions fixed) and then is accepted or rejected. In our case, at each step, a new mixture component is proposed by redrawing one of the current sample indices associated with the component uniformly and then accepting or rejecting the resulting proposed component based on its mixture weight. We give the IMG algorithm for combining subposterior samples in Algorithm 1.222Again for simplicity, we assume that we generate samples to represent the full posterior, where is the number of subposterior samples from each machine.
In certain situations, Algorithm 1 may have a low acceptance rate and therefore may mix slowly. One way to remedy this is to perform the IMG combination algorithm multiple times, by first applying it to groups of subposteriors and then applying the algorithm again to the output samples from each initial application. For example, one could begin by applying the algorithm to all pairs (leaving one subposterior alone if is odd), then repeating this process—forming pairs and applying the combination algorithm to pairs only—until there is only one set of samples remaining, which are samples from the density product estimate.
3.3. Asymptotically exact posterior sampling with semiparametric density product estimation
Our first example made use of a parametric estimator, which has quick convergence, but may be asymptotically biased, while our second example made use of a nonparametric estimator, which is asymptotically exact, but may converge slowly when the number of dimensions is large. In this example, we implicitly sample from a semiparametric density product estimate, which allows us to leverage the fact that the full posterior has a near-Gaussian form when the number of observations is large, while still providing an asymptotically unbiased estimate of the posterior density, as the number of subposterior samples .
We make use of a semiparametric density estimator for that consists of the product of a parametric estimator and a nonparametric estimator of the correction function . This estimator gives a near-Gaussian estimate when the number of samples is small, and converges to the true density as the number of samples grows. Given samples from a subposterior , we can write the estimator as
where we have used a Gaussian kernel with bandwidth parameter for the nonparametric component of this estimator. Therefore, we define our semiparametric density product estimator to be
This estimate is proportional to the pdf of a mixture of Gaussians with unnormalized mixture weights,
where and are given by Eq. 3.1 and 3.2. We can sample from this semiparametric estimate using the IMG procedure outlined in Algorithm 1, replacing the component weights with and the component parameters and with and .
We also have a second semiparametric procedure that may give higher acceptance rates in the IMG algorithm. We follow the above semiparametric procedure, where each component is a normal distribution with parameters and , but we use the nonparametric component weights instead of . This procedure is also asymptotically exact, since the semiparametric component parameters and approach the nonparametric component parameters and as , and thus this procedure tends to the nonparametric procedure given in Algorithm 1.
4. Method Complexity
Given data subsets, to produce samples in dimensions with the nonparametric or semiparametric asymptotically exact procedures (Algorithm 1) requires operations. The variation on this algorithm that performs this procedure times on pairs of subposteriors (to increase the acceptance rate; detailed in Section 3.2) instead requires only operations.
We have presented our method as a two step procedure, where first parallel MCMC is run to completion, and then the combination algorithm is applied to the sets of samples. We can instead perform an online version of our algorithm: as each machine generates a sample, it immediately sends it to a master machine, which combines the incoming samples333For the semiparametric method, this will involve an online update of mean and variance Gaussian parameters. and performs the accept or reject step (Algorithm 1, lines 3-12). This allows the parallel MCMC phase and the combination phase to be performed in parallel, and does not require transfering large volumes of data, as only a single sample is ever transferred at a time.
The total communication required by our method is transferring scalars ( samples from each of machines), and as stated above, this can be done online as MCMC is being carried out. Further, the communication is unidirectional, and each machine does not pause and wait for any information from other machines during the parallel sampling procedure.
5. Theoretical Results
Our second and third procedures aim to draw asymptotically exact samples by sampling from (fully or partially) nonparametric estimates of the density product. We prove the asymptotic correctness of our estimators, and bound their rate of convergence. This will ensure that we are generating asymptotically correct samples from the full posterior as the number of samples from each subposterior grows.
5.1. Density product estimate convergence and risk analysis
To prove (mean-square) consistency of our estimator, we give a bound on the mean-squared error (MSE), and show that it tends to zero as we increase the number of samples drawn from each subposterior. To prove this, we first bound the bias and variance of the estimator. The following proofs make use of similar bounds on the bias and variance of the nonparametric and semiparametric density estimators, and therefore the theory applies to both the nonparametric and semiparametric density product estimators.
Throughout this analysis, we assume that we have samples from each subposterior (), and that denotes the bandwidth of the nonparametric density product estimator (which is annealed to zero as in Algorithm 1). Let Hölder class on be defined as the set of all times differentiable functions whose derivative satisfies
We also define the class of densities to be
We also assume that all subposterior densities are bounded, i.e. that there exists some such that for all and .
First, we bound the bias of our estimator. This shows that the bias tends to zero as the bandwidth shrinks.
The bias of the estimator satisfies
for some .
For all ,
where we have used the fact that for some . ∎
Next, we bound the variance of our estimator. This shows that the variance tends to zero as the number of samples grows large and the bandwidth shrinks.
The variance of the estimator satisfies
for some and .
For all ,
where we have used the facts that for some and for some . ∎
Finally, we use the bias and variance bounds to bound the MSE, which shows that our estimator is consistent.
If , the mean-squared error of the estimator satisfies
for some and .
For all , using the fact that the mean-squared error is equal to the variance plus the bias squared, we have that
for some and . ∎
6. Method Scope
The theoretical results and algorithms in this paper hold for posterior distributions over finite-dimensional real spaces. These include generalized linear models (e.g. linear, logistic, or Poisson regression), mixture models with known weights, hierarchical models, and (more generally) finite-dimensional graphical models with unconstrained variables. This also includes both unimodal and multimodal posterior densities (such as in Section 8.2). However, the methods and theory presented here do not yet extend to cases such as infinite dimensional models (e.g. nonparametric Bayesian models ) nor to distributions over the simplex (e.g. topics in latent Dirichlet allocation ). In the future, we hope to extend this work to these domains.
7. Related Work
In [19, 2, 16], the authors develop a way to sample approximately from a posterior distribution when only a small randomized mini-batch of data is used at each step. In , the authors used a hypothesis test to decide whether to accept or reject proposals using a small set of data (adaptively) as opposed to the exact Metropolis-Hastings rule. This reduces the amount of time required to compute the acceptance ratio. Since all of these algorithms are still sequential, they can be directly used in our algorithm to generate subposterior samples to further speed up the entire sampling process.
Several parallel MCMC algorithms have been designed for specific models, such as for topic models [18, 14] and nonparametric mixture models . These approaches still require synchronization to be correct (or approximately correct), while ours aims for more general model settings and does not need synchronization until the final combination stage.
Consensus Monte Carlo  is perhaps the most relevant work to ours. In this algorithm, data is also portioned into different machines and MCMC is performed independently on each machine. Thus, it roughly has the same time complexity as our algorithm. However, the prior is not explicitly reweighted during sampling as we do in Eq 2.1, and final samples for the full posterior are generated by averaging subposterior samples. Furthermore, this algorithm has few theoretical guarantees. We find that this algorithm can be viewed as a relaxation of our nonparametric, asymptotically exact sampling procedure, where samples are generated from an evenly weighted mixture (instead of each component having weight ) and where each sample is set to instead of being drawn from . This algorithm is one of our experimental baselines.
8. Empirical Study
In the following sections, we demonstrate empirically that our method allows for quicker, MCMC-based estimation of a posterior distribution, and that our consistent-estimator-based procedures yield asymptotically exact results. We show our method on a few Bayesian models using both synthetic and real data. In each experiment, we compare the following strategies for parallel, communication-free sampling:444We did not directly compare with the algorithms that require synchronization since the setup of these experiments can be rather different. We plan to explore these comparisons in the extended version of this paper.
Single chain full-data posterior samples (regularChain)—Typical, single-chain MCMC for sampling from the full-data posterior.
Parametric subposterior density product estimate (parametric)—For sets of subposterior samples, the combination yielding samples from the parametric density product estimate.
Nonparametric subposterior density product estimate (nonparametric)—For sets of subposterior samples, the combination yielding samples from the nonparametric density product estimate.
Semiparametric subposterior density product estimate (semiparametric)—For sets of subposterior samples, the combination yielding samples from the semiparametric density product estimate.
Subposterior sample average (subpostAvg)—For sets of subposterior samples, the average of samples consisting of one sample taken from each subposterior.
Subposterior sample pooling (subpostPool)—For sets of subposterior samples, the union of all sets of samples.
Duplicate chains full-data posterior sample pooling (duplicateChainsPool)—For sets of samples from the full-data posterior, the union of all sets of samples.
To assess the performance of our sampling and combination strategies, we ran a single chain of MCMC on the full data for 500,000 iterations, removed the first half as burn-in, and considered the remaining samples the “groundtruth” samples for the true posterior density. We then needed a general method to compare the distance between two densities given samples from each, which holds for general densities (including multimodal densities, where it is ineffective to compare moments such as the mean and variance555In these cases, dissimilar densities might have similar low-order moments. See Section 8.2 for an example.). Following work in density-based regression , we use an estimate of the distance, , between the groundtruth posterior density and a proposed posterior density , where .
In the following experiments involving timing, to compute the posterior error at each time point, we collected all samples generated before a given number of seconds, and added the time taken to transfer the samples and combine them using one of the proposed methods. In all experiments and methods, we followed a fixed rule of removing the first samples for burn-in (which, in the case of combination procedures, was applied to each set of subposterior samples before the combination was performed).
Experiments were conducted with a standard cluster system. We obtained subposterior samples by submitting batch jobs to each worker since these jobs are all independent. We then saved the results to the disk of each worker and transferred them to the same machine which performed the final combination.
8.1. Generalized Linear Models
Generalized linear models are widely used for many regression and classification problems. Here we conduct experiments, using logistic regression as a test case, on both synthetic and real data to demonstrate the speed of our parallel MCMC algorithm compared with typical MCMC strategies.
8.1.1. Synthetic data
Our synthetic dataset contains 50,000 observations in dimensions. To generate the data, we drew each element of the model parameter and data matrix from a standard normal distribution, and then drew each outcome as (where denotes the row of )666Note that we did not explicitly include the intercept term in our logistic regression model.. We use Stan, an automated Hamiltonian Monte Carlo (HMC) software package,777http://mc-stan.org to perform sampling for both the true posterior (for groundtruth and comparison methods) and for the subposteriors on each machine. One advantage of Stan is that it is implemented with C++ and uses the No-U-Turn sampler for HMC, which does not require any user-provided parameters .
In Figure 1, we illustrate results for logistic regression, showing the subposterior densities, the subposterior density product, the subposterior sample average, and the true posterior density, for the number of subsets set to 10 (left) and 20 (right). Samples generated by our approach (where we draw samples from the subposterior density product via the parametric procedure) overlap with the true posterior much better than those generated via the subpostAvg (subposterior sample average) procedure— averaging of samples appears to create systematic biases. Futher, the error in averaging appears to increase as grows. In Figure 2 (left) we show the posterior error vs time. A regular full-data chain takes much longer to converge to low error compared with our combination methods, and simple averaging and pooling of subposterior samples gives biased solutions.
We next compare our combination methods with multiple independent “duplicate” chains each run on the full dataset. Even though our methods only require a fraction of the data storage on each machine, we are still able to achieve a significant speed-up over the full-data chains. This is primarily because the duplicate chains cannot parallelize burn-in (i.e. each chain must still take some steps before generating reasonable samples, and the time taken to reach these steps does not decrease as more machines are added). However, in our method, each subposterior sampler can take each step more quickly, effectively allowing us to decrease the time needed for burn-in as we increase . We show this empirically in Figure 2 (right), where we plot the posterior error vs time, and compare with full duplicate chains as is increased.
Using a Matlab implementation of our combination algorithms, all (batch) combination procedures take under twenty seconds to complete on a 2.5GHz Intel Core i5 with 16GB memory.
8.1.2. Real-world data
Here, we use the covtype (predicting forest cover types)888http://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets dataset, containing 581,012 observations in dimensions. A single chain of HMC running on this entire dataset takes an average of 15.76 minutes per sample; hence, it is infeasible to generate groundtruth samples for this dataset. Instead we show classification accuracy vs time. For a given set of samples, we perform classification using a sample estimate of the posterior predictive distribution for a new label with associated datapoint , i.e.
where and denote the observations, and . Figure 3 (left) shows the results for this task, where we use =50 splits. The parallel methods achieve a higher accuracy much faster than the single-chain MCMC algorithm.
8.1.3. Scalability with dimension
We investigate how the errors of our methods scale with dimensionality, to compare the different estimators implicit in the combination procedures. In Figure 3 (right) we show the relative posterior error (taken at 1000 seconds) vs dimension, for the synthetic data with =10 splits. The errors at each dimension are normalized so that the regularChain error is equal to 1. Here, the parametric (asymptotically biased) procedure scales best with dimension, and the semiparametric (asymptotically exact) procedure is a close second. These results also demonstrate that, although the nonparametric method can be viewed as implicitly sampling from a nonparametric density estimate (which is usually restricted to low-dimensional densities), the performance of our method does not suffer greatly when we perform parallel MCMC on posterior distributions in much higher-dimensional spaces.
8.2. Gaussian mixture models
In this experiment, we aim to show correct posterior sampling in cases where the full-data posterior, as well as the subposteriors, are multimodal. We will see that the combination procedures that are asymptotically biased suffer greatly in these scenarios. To demonstrate this, we perform sampling in a Gaussian mixture model. We generate 50,000 samples from a ten component mixture of 2-d Gaussians. The resulting posterior is multimodal; this can be seen by the fact that the component labels can be arbitrarily permuted and achieve the same posterior value. For example, we find after sampling that the posterior distribution over each component mean has ten modes. To sample from this multimodal posterior, we used the Metropolis-Hastings algorithm, where the component labels were permuted before each step (note that this permutation results in a move between two points in the posterior distribution with equal probability).
In Figure 4 we show results for splits, showing samples from the true posterior, overlaid samples from all five subposteriors, results from averaging the subposterior samples, and the results after applying our three subposterior combination procedures. This figure shows the 2-d marginal of the posterior corresponding to the posterior over a single mean component. The subpostAvg and parametric procedures both give biased results, and cannot capture the multimodality of the posterior. We show the posterior error vs time in Figure 5 (left), and see that our asymptotically exact methods yield quick convergence to low posterior error.
8.3. Hierarchical models
We show results on a hierarchical Poisson-gamma model of the following form
for =50,000 observations. We draw from the above generative process (after fixing values for , , , and ), and use =10 splits. We again perform MCMC using the Stan software package.
In Figure 5 (right) we show the posterior error vs time, and see that our combination methods complete burn-in and converge to a low posterior error very quickly relative to the subpostAvg and subpostPool procedures and full-data chains.
9. Discussion and Future Work
In this paper, we present an embarrassingly parallel MCMC algorithm and provide theoretical guarantees about the samples it yields. Experimental results demonstrate our method’s potential to speed up burn-in and perform faster asymptotically correct sampling. Further, it can be used in settings where data are partitioned onto multiple machines that have little intercommunication—this is ideal for use in a MapReduce setting. Currently, our algorithm works primarily when the posterior samples are real, unconstrained values and we plan to extend our algorithm to more general settings in future work.
-  Alekh Agarwal and John C Duchi, Distributed delayed stochastic optimization, Decision and Control (CDC), 2012 IEEE 51st Annual Conference on, IEEE, 2012, pp. 5451–5452.
-  Sungjin Ahn, Anoop Korattikara, and Max Welling, Bayesian posterior sampling via stochastic gradient fisher scoring, Proceedings of the 29th International Conference on Machine Learning, 2012, pp. 1591–1598.
-  David M Blei, Andrew Y Ng, and Michael I Jordan, Latent dirichlet allocation, The Journal of Machine Learning Research 3 (2003), 993–1022.
-  Jeffrey Dean and Sanjay Ghemawat, Mapreduce: simplified data processing on large clusters, Communications of the ACM 51 (2008), no. 1, 107–113.
-  Samuel J Gershman and David M Blei, A tutorial on bayesian nonparametric models, Journal of Mathematical Psychology 56 (2012), no. 1, 1–12.
-  Nils Lid Hjort and Ingrid K Glad, Nonparametric density estimation with a parametric start, The Annals of Statistics (1995), 882–904.
-  Qirong Ho, James Cipar, Henggang Cui, Seunghak Lee, Jin Kyu Kim, Phillip B. Gibbons, Gregory R. Ganger, Garth Gibson, and Eric P. Xing, More effective distributed ml via a stale synchronous parallel parameter server, Advances in Neural Information Processing Systems, 2013.
-  Matthew D Hoffman and Andrew Gelman, The no-u-turn sampler: Adaptively setting path lengths in hamiltonian monte carlo, arXiv preprint arXiv:1111.4246 (2011).
-  Anoop Korattikara, Yutian Chen, and Max Welling, Austerity in MCMC land: Cutting the Metropolis-Hastings budget, arXiv preprint arXiv:1304.5299 (2013).
-  John Langford, Alex J Smola, and Martin Zinkevich, Slow learners are fast, Advances in Neural Information Processing Systems, 2009.
-  Kathryn Blackmond Laskey and James W Myers, Population Markov chain Monte Carlo, Machine Learning 50 (2003), no. 1-2, 175–196.
-  Lucien Le Cam, Asymptotic methods in statistical decision theory, New York (1986).
-  Lawrence Murray, Distributed Markov chain Monte Carlo, Proceedings of Neural Information Processing Systems Workshop on Learning on Cores, Clusters and Clouds, vol. 11, 2010.
-  David Newman, Arthur Asuncion, Padhraic Smyth, and Max Welling, Distributed algorithms for topic models, The Journal of Machine Learning Research 10 (2009), 1801–1828.
-  Junier Oliva, Barnabás Póczos, and Jeff Schneider, Distribution to distribution regression, Proceedings of The 30th International Conference on Machine Learning, 2013, pp. 1049–1057.
-  Sam Patterson and Yee Whye Teh, Stochastic gradient riemannian langevin dynamics on the probability simplex, Advances in Neural Information Processing Systems, 2013.
-  Steven L. Scott, Alexander W. Blocker, and Fernando V. Bonassi, Bayes and big data: The consensus monte carlo algorithm, Bayes 250, 2013.
-  Alexander Smola and Shravan Narayanamurthy, An architecture for parallel topic models, Proceedings of the VLDB Endowment 3 (2010), no. 1-2, 703–710.
-  Max Welling and Yee W Teh, Bayesian learning via stochastic gradient Langevin dynamics, Proceedings of the 28th International Conference on Machine Learning, 2011, pp. 681–688.
-  Darren J Wilkinson, Parallel Bayesian computation, Statistics Textbooks and Monographs 184 (2006), 477.
-  Sinead Williamson, Avinava Dubey, and Eric P Xing, Parallel Markov chain Monte Carlo for nonparametric mixture models, Proceedings of the 30th International Conference on Machine Learning, 2013, pp. 98–106.