SAME but Different: Fast and HighQuality Gibbs Parameter Estimation
Abstract
Gibbs sampling is a workhorse for Bayesian inference but has several limitations when used for parameter estimation, and is often much slower than nonsampling inference methods. SAME (State Augmentation for Marginal Estimation) [14, 8] is an approach to MAP parameter estimation which gives improved parameter estimates over direct Gibbs sampling. SAME can be viewed as cooling the posterior parameter distribution and allows annealed search for the MAP parameters, often yielding very high quality (lower loss) estimates. But it does so at the expense of additional samples per iteration and generally slower performance. On the other hand, SAME dramatically increases the parallelism in the sampling schedule, and is an excellent match for modern (SIMD) hardware. In this paper we explore the application of SAME to graphical model inference on modern hardware. We show that combining SAME with factored sample representation (or approximation) gives throughput competitive with the fastest symbolic methods, but with potentially better quality. We describe experiments on Latent Dirichlet Allocation, achieving speeds similar to the fastest reported methods (online Variational Bayes) and lower crossvalidated loss than other LDA implementations. The method is simple to implement and should be applicable to many other models.
SAME but Different: Fast and HighQuality Gibbs Parameter Estimation
Huasha Zhao Computer Science Division UC Berkeley Berkeley, CA 94720 hzhao@cs.berkeley.edu Biye Jiang Computer Science Division UC Berkeley Berkeley, CA 94720 bjiang@cs.berkeley.edu John Canny Computer Science Division UC Berkeley Berkeley, CA 94720 jfc@cs.berkeley.edu
1 Introduction
Many machine learning problems can be formulated as inference on a joint distribution where represents observed data, a set of parameters, and represents latent variables. Both and are latent in general, but it is useful to distinguish between them  encodes which of a class of models represents the current situation, while represent local labels or missing data. One generally wants to optimize over while marginalizing over . And the output of the algorithm is a value or distribution over while the are often ignored.
Gibbs sampling is a very general approach to posterior estimation for , but it provides samples only rather than MAP estimates. But therein lies a problem: sampling is a sensible approach to marginal estimation, but can be a very inefficient approach to optimization. This is particularly true when the dimension of is large compared to (which is true e.g. in Latent Dirichlet Allocation and probabilistic recommendation algorithms). Such models have been observed to require many samples (thousands to hundreds of thousands) to provide good parameter estimates. Hybrid approaches such as MonteCarlo EM have been developed to address this issue  a MonteCarlo method such as Gibbs sampling is used to estimate the expected values in the Estep while an optimization method is applied to the parameters in the Mstep. But this requires a separate optimization strategy (usually gradientbased), a way to compute the dependence on the parameters symbolically, and analysis of the accuracy of the Estep estimates.
SAME (State Augmentation for Marginal Estimation) [14, 8] is a simple approach to MAP parameter estimation that remains within the Gibbs framework^{1}^{1}1SAME is a general approach to MCMC MAP estimation, but in this paper we will focus on its realization in Gibbs samplers. SAME replicates the latent state with additional states. This has the effect of “cooling” the marginal distribution on , which sharpens its peaks and causes samples to approach local optima. The conditional distribution remains the same, so we are still marginalizing over a full distribution on . By making the temperature a controllable parameter, the parameter estimates can be annealed to reach better local optima. In both [14, 8] and the present paper we find that this approach gives better estimates than competing approaches. The novelty of the present paper is showing that SAME estimation can be very fast, and competitive with the fastest symbolic meethods. Thus it holds the potential to be the method of choice for many inference problems.
Specifically, we define a new joint distribution
(1) 
which models copies of the original system with tied parameters and independent latent variable blocks . The marginalized conditional . And
(2) 
where . So which is up to a constant factor equal to , a power of the original marginal parameter distribution. Thus it has the same optima, including the global optimum, but its peaks are considerably sharpened. In what follows we will often demote to a subscript since it fixed, writing as etc.
This new distribution can be written as a Gibbs distribution on , as , from which we see that is an inverse temperature parameter ( is Boltzmann’s constant). Increasing amounts to cooling the distribution.
Gibbs sampling from the new distribution is usually straightforward given a sampler for the original. It is perhaps not obvious why sampling from a more complex system could improve performance, but we have added considerable parallelism since we can sample various “copies” of the system concurrently. It will turn out this approach is complementary to using a factored form for the posterior. Together these methods gives us ordersofmagnitude speedup over other samplers for LDA.
The rest of the paper is organized as follows. Section 2 summarizes related work on parameter inference for probabilistic models and their limitations. Section 3 introduces the SAME sampler. We discuss in Section 4 a factored approximation that considerably accelerates sampling. A hardwareoptimized implementation of the algorithm for LDA is described in Section 5. Section 6 presents the experimental results and finally Section 7 concludes the paper.
2 Related Work
2.1 EM and related Algorithms
The ExpectationMaximization (EM) algorithm [7] is a popular method for parameter estimation of graphical models of the form we are interested in. The EM algorithm alternates between updates in the expectation (E) step and maximization (M) step. The Estep marginalizes out the latent variables and computes the expectation of the likelihood as a function of the parameters. The E step computes a function to be optimized in the Mstep.
For EM to work, one has to compute the expectation of the sufficient statistics of the likelihood function, where the expectation is over . It also requires a method to optimize the function. In practice, the iterative update equations can be hard to derive. Moreover, the EM algorithm is a gradientbased method, and therefore is only able to find locallyoptimal solutions.
Variational Bayes (VB) [12] is an EMlike algorithm that uses a parametric approximation to the posterior distribution of both parameters and other latent variables, and attempts to optimize the fit (e.g. using KLdivergence) to the observed data. Parameter estimates are usually taken from the means of the parameter approximations (the hyperparameters), unlike EM where point parameter estimates are used. It is common to assume a coordinatefactored form for this approximate posterior. The factored form simplifies inference, but makes strong assumptions about the distribution (effectively eliminating interactions). It makes most sense for parameter estimates, but is a strong constraint to apply to other latent variables. VB has been criticized for lack of quality (higher loss) on certain problems because of this. Nevertheless it works well for many problems (e.g. on LDA). This motivated us to introduce a similar factored approximation in our method described laster.
2.2 Gibbs Sampling
Gibbs samplers [9] work by stepping through individual (or blocks of) latent variables, and in effect perform a simulation of the stochastic system described by the joint distribution. The method is simple to apply on graphical models since each variable is conditioned locally by its Markov blanket. They are easy to define for conjugate distributions and can generalize to nonconjugate graphical models using e.g. slice sampling. Gibbs samplers are often the first method to be derived, and sometimes the only practical learning strategy for complex graphical models.
Gibbs samplers only require that the joint density is strictly positive over the sample space which is known as the Gibbs distribution. However, Gibbs sampling only gives samples from the distribution of , it can be difficult to find ML or MAP values from those samples. Furthermore, it can be very slow, especially for models with high dimensions. Our results suggest this slow convergence is often due to large variance in the parameters in standard samplers.
2.3 Monte Carlo EM
Monte Carlo EM [15] is a hybrid approach that uses MCMC (e.g. Gibbs sampling) to approximate the expected value with the mean of the loglikelihood of the samples. The method has to optimize using a numerical method (conjugate gradient etc.). Like standard EM, it can suffer from convergence problems, and may only find a local optimum of likelihood.
2.4 MessagePassing Methods
Belief propagation [16] and Expectation propagation [13] use local (nodewise) updates to infer posterior parameters in graphical models in a manner reminiscent of Gibbs sampling. But they are exact only for a limited class of models. Recently variational messagepassing [17] has extended the class of models for which parametric inference is possible to conjugateexponential family graphical models. However similar to standard VB, the method uses a coordinatefactored approximation to the posterior which effectively eliminates interactions (although they can be added at high computational cost by using a factor graph). It also finds only local optima of the posterior parameters.
3 SAME Parameter Estimation
SAME estimation involves sampling multiple ’s independently and inferring using the aggregate of ’s.
3.1 Method
We use the notation and similarly for .
Sampling in the SAME sampler is exactly the same as for the standard sampler. Since the groups are independent of each other, we can use the sampling function for the original distribution, conditioned only on the other components of the same group: .
Sampling is only slightly more complicated. We want to sample from
(3) 
where and if we ignore the normalizing constants:
(4) 
which is now expressed as a product of conditionals from the original sampler . Inference in the new model will be tractable if we are able to sample from a product of the distributions . This will be true for many distributions. e.g. for exponential family distributions in canonical form, the product is still an exponential family member. A product of Dirichlet distributions is Dirichlet etc., and in general this distribution represents the parameter estimate obtained by combining evidence from independent observations. The normalizing constant will usually be implied from the closedform parameters of this distribution.
Adjusting sample number at different iterations allows annealing of the estimate.
4 CoordinateFactored Approximation
Certain distributions (including LDA) have the property that the latent variables are independent given and . That is . Therefore the ’s can be sampled (without approximation) in parallel. Furthermore, rather than a single sample from (e.g. a categorical sample for a discrete ) we can construct a SAME Gibbs sampler by taking samples. These samples will now have a multinomial distribution with count and probability vector . Let denote the count for among the samples, and denote the conditional probability that .
We can introduce still more parallelism by randomizing the order in which we choose which from which to sample. The count for variable is then replaced by random variable and the coordinatewise distributions of become independent Poisson variables:
(5) 
when the are independent given , the counts fully capture the results of taking the (independent) samples. These samples can be generated very fast, and completely in parallel. is no longer constrained to be an integer, or even to be . Indeed, each sample no longer corresponds to execution of a block of code, but is simply an increment in the value of the Poisson random number generator. In LDA, it is almost as fast to generate all samples for a single word for a large (up to about 100) as it is to generate one sample. This is a source of considerable speedup in our LDA implementation.
For distributions where the are not independent, i.e. when , we can still perform independent sampling as an approximation. This approach is quite similar to the coordinatefactored approximation often used in Variational Bayes. We leave the details to a forthcoming paper.
5 Implementation of SAME Gibbs LDA
Our SAME LDA sampler implementation is described in Algorithm 3. Samples are taken directly from Poisson distributions (line of the algorithm) as described earlier.
5.1 Latent Dirichlet Allocation
LDA is a generative process for modeling a collection of documents in a corpus. In the LDA model, the words are observed, the topics are latent variables and parameters are where is document topic distributions and is word topic distributions. Subscript denotes document and its word. Given and , are independent, which also implies that we can sample from without any information about from previous samples. A similar result holds for the parameters, which can be sampled given only the current counts from the . We can therefore alternate parameter and latent variable inference using a pair of samplers: and . Such blocked sampling maximizes parallelism and works very well with modern SIMD hardware.
The sampling of ’s uses the Poisson formula derived earlier. The conditional distributions are multiple independent Dirichlet’s. In practice, we collapse out , so we in effect sample a new given the from a previous iteration. The update sampler follows a Dirichlet compound multinomial distribution (DCM) and can be derived as,
(6) 
where and are numbers of words and topics respectively, are counts across all samples and all documents which are defined as,
(7) 
where is the number of documents and superscript of denotes the sample of the hidden topic variable. In Equation 6, dot(s) in the subscript of denotes integrating (summing) over that dimension(s), for example, . As shown in Equation 6, sample counts are sufficient statistics for the update formula.
5.2 Minibatch Processing
For scalability (to process datasets that will not fit in memory), we implement LDA using a minibatch update strategy. In minibatch processing, data is read from a dataset and processed in blocks. Minibatch algorithms have been shown to be very efficient for approximate inference [4, 3].
In Algorithm 3, is a sparse (nwords x ndocs) matrix that encodes the subset of documents (minibatch) to process at period . The global model (across minibatches) is the wordtopic matrix . It is updated as the weighted average of the current model and the new update, see line . The weight is determined by according to [11]. We do not explicitly denote passes over the dataset. The data are treated instead as an infinite stream and we examine the crossvalidated likelihood as a function of the number of minibatch steps, up to .
5.3 GPU optimizations
GPUs are extremely wellsuited to Gibbs sampling by virtue of their large degree of parallelism, and also because of their extremely high throughput in nonuniform random number generation (thanks to hardwareaccelerated transcendental function evaluation). For best performance we must identify and accelerate the bottleneck steps in the algorithm. First, computing the normalizing factor in equation 6 is a bottleneck step. It involves evaluating the product of two dense matrix at only nonzeros of a sparse matrix . Earlier we developed a kernel (SDDMM) for this step which is a bottleneck for several other factor models including our online Variational Bayes LDA and Sparse Factor Analysis [6].
Second, line 7 of the algorithm is another dominant step, and has the same operation count as the SDDMM step. We wrote a custom kernel that implements lines 411 of the algorithm with almost the same speed as SDDMM. Finally, we use a matrixcaching strategy [5] to eliminate the need for memory allocation on the GPU in each iteration.
6 Experiments
6.1 Lda
In this section, we evaluate the performance of the SAME/Factored Gibbs sampler against several other algorithms and systems. We implemented LDA using our SAME approach, VB online and VB batch, all using GPU acceleration. The code is open source and distributed as part of the BIDMach project [5, 6] on github ^{2}^{2}2https://github.com/BIDData/BIDMach/wiki. We compare our systems with four other systems: 1) Blei et al’s VB batch implementation for LDA [2], 2) Griffiths et al’s collapsed Gibbs sampling (CGS) for LDA [10], 3) Yan et al’s GPU parallel CGS for LDA [18], 4) Ahmed et al’s cluster CGS for LDA [1].
1) and 2) are both C/C++ implementations of the algorithm ^{3}^{3}3 2) provides a Matlab interface for CGS LDA on CPUs. 3) is the stateoftheart GPU implementation of parallel CGS. To our best knowledge, 3) is the fastest singlemachine LDA implementation to date, and 4) is the fastest cluster implementation of LDA.
All the systems/algorithms are evaluated on a single PC equipped with a single 8core CPU (Intel E52660) and a dualcore GPU (Nvidia GTX690), except GPU CGS and cluster CGS. Each GPU core comes with 2GB memory. Only one core of GPU is used in the benchmark. The benchmark for GPU CGS is reported in [18]. They run the algorithm on a machine with a GeForce GTX280 GPU. The benchmarks we use for cluster CGS were reported on 100 and 1000 Yahoo cluster nodes [1].
Two datasets are used for evaluation. 1) NYTimes. The dataset has approximately 300k New York Times news articles. There are 100k unique words and 100 million tokens in the corpus. 2) PubMed. The dataset contains about 8.2 million abstracts from collections of US National Library of Medicine. There are 140k unique tokens, and 730 million words in the corpus.
6.1.1 Convergence
We first compare the convergence of SAME Gibbs LDA with other methods. We choose for the SAME sampler, because convergence speed only improves marginally beyond , however, runtime per pass over the dataset starts to increase noticeably beyond while being relatively flat for . The minibatch size is set to be of the total number of examples. We use the per word log likelihood as a measure of convergence.
Figure 1 shows the crossvalidated likelihood as a function of the number of passes through the data for both datasets, up to 20 passes. As we can see, the SAME Gibbs sampler converges to a higher quality result than VB online and VB batch on both datasets. The SAME sampler and VB online converged after 45 passes for the NYTimes dataset and 2 passes for the PubMed dataset. VB batch converges in about 20 passes.
All three methods above are able to produce higher quality than CGS within 20 passes over the dataset. It usually takes 10002000 passes for CGS to converge for typical datasets such as NYTimes, as can be seen in Figure 0(c).
Figure 0(c) plots loglikelihood against the number of samples taken per word. We compare the SAME sampler with , and CGS. To reach the same number of samples CGS and the SAME sampler with need to run 100 times as many passes over the data as the SAME sampler with . That is Figure 0(c) shows 20 passes of SAME sampler with and 2000 passes of the other two methods. At the beginning, CGS converges faster. But in the long run, SAME Gibbs leads to better convergence. Notice that SAME with is not identical to CGS in our implementation, because of minibatching and the moving average estimate for .
6.1.2 Runtime Performance
We further measure the runtime performance of different methods with different implementations for LDA. Again we fix the sample size for SAME GS at . All the runtime number are generated by running each method/system on the NYTimes dataset with , except that Yan’s parallel CGS [18] reports their benchmark for . We also report time to convergence, which is defined as the time to reach the loglikelihood for standard CGS at the iteration.
SAME GS BIDMach  VB online BIDMach  VB batch BIDMach  VB batch Blei et al  CGS Griffiths et al  GPU CGS Yan et al  

Runtime per pass (s)  30  20  20  12600  225  5.4 
Time to converge (s)  90  40  400  252000  225000  5400 
Results are illustrated in the Table 1. The SAME Gibbs sampler takes 90 seconds to converge on the NYTimes dataset. As comparison, the other CPU implementations take around 6070 hours to converge, and Yan’s GPU implementation takes 5400 second to converge for . Our system demonstrates two orders of magnitude improvement over the state of the art. Our implementation of online VB takes about 40 seconds to converge on the NYTimes dataset. The Gibbs sampling method is close in performance (less than a factor of 3) to that.
Finally, we compare our system with the stateofart cluster implementation of CGS LDA [1], using time to convergence. Ahmed et al. [1] processed 200 million news articles on 100 machines and 1000 iterations in 2mins/iteration = 2000 minutes overall = 120k seconds. We constructed a repeating stream of news articles (as did [1]) and ran for two iterations  having found that this was sufficient for news datasets of comparable size. This took 30k seconds, which is 4x faster, on a single GPU node. Ahmed et al. also processed 2 billion articles on 1000 machines to convergence in 240k seconds, which is slightly less than linear scaling. Our system (which is sequential on minibatches) simply scales linearly to 300k seconds on this problem. Thus singlemachine performance of GPUaccelerated SAME GS is almost as fast as a custom cluster CGS on 1000 nodes, and 4x faster than 100 nodes.
6.1.3 Multisampling and Annealing
The effect of sample number is studied in Figure 2 and Table 2. As expected, more samples yields better convergences given fixed number of passes through the data. And the benefit comes almost free thanks to the SIMD parallelism offered by the GPU. As shown in Table 2, is only modestly slower than . However, the runtime for is much longer than .
We next studied the effects of dynamic adjustment of , i.e. annealing. We compare constant scheduling (fixed ) with

linear scheduling, ,

logarithmic scheduling, ,

invlinear scheduling, .
is the total number of iterations. The average sample size per iteration is fixed to for all 4 configurations. As shown in Figure 1(b), we cannot identify any particular annealing schedule that is significantly better then fixed sample size. Invlinear is slighter faster at the beginning but has the highest initial sample number.
m=500  m=100  m=20  m=5  m=1  

Runtime per iteration (s)  50  30  25  23  20 
7 Conclusions and Future Work
This paper described hardwareaccelerated SAME Gibbs Parameter
estimation  a method to improve and accelerate parameter estimation
via Gibbs sampling. This approach reduces the number of passes over
the dataset while introducing more parallelism into the sampling
process. We showed that the approach meshes very well with SIMD
hardware, and that a GPUaccelerated implementation of cooled GS for
LDA is faster than other sequential systems, and comparable with the
fastest cluster implementation of CGS on 1000 nodes. The code is at
https://github.com/BIDData/BIDMach
SAME GS is applicable to many other problems, and we are currently exploring the method for inference on general graphical models. The coordinatefactored sampling approximation is also being applied to this problem in conjunction with full sampling (the approximation reduces the number of full samples required) and we anticipate similar large improvements in speed. The method provides the quality advantages of sampling and annealed estimation, while delivering performance with custom (symbolic) inference methods. We believe it will be the method of choice for many inference problems in future.
References
 [1] A. Ahmed, M. Aly, J. Gonzalez, S. Narayanamurthy, and A. J. Smola. Scalable inference in latent variable models. In Proceedings of the fifth ACM international conference on Web search and data mining, pages 123–132. ACM, 2012.
 [2] D. M. Blei, A. Y. Ng, and M. I. Jordan. Latent dirichlet allocation. the Journal of machine Learning research, 3:993–1022, 2003.
 [3] L. Bottou. Largescale machine learning with stochastic gradient descent. In Proceedings of COMPSTAT’2010, pages 177–186. Springer, 2010.
 [4] L. Bottou and O. Bousquet. The tradeoffs of large scale learning. In NIPS, volume 4, page 2, 2007.
 [5] J. Canny and H. Zhao. Bidmach: Largescale learning with zero memory allocation. In BigLearn Workshop, NIPS, 2013.
 [6] J. Canny and H. Zhao. Big data analytics with small footprint: Squaring the cloud. In Proceedings of the 19th ACM SIGKDD international conference on Knowledge discovery and data mining, pages 95–103. ACM, 2013.
 [7] A. P. Dempster, N. M. Laird, D. B. Rubin, et al. Maximum likelihood from incomplete data via the em algorithm. Journal of the Royal statistical Society, 39(1):1–38, 1977.
 [8] A. Doucet, S. Godsill, and C. Robert. Marginal maximum a posteriori estimation using markov chain monte carlo. Statistics and Computing, 12:77–84, 2002.
 [9] S. Geman and D. Geman. Stochastic relaxation, gibbs distributions, and the bayesian restoration of images. Pattern Analysis and Machine Intelligence, IEEE Transactions on, (6):721–741, 1984.
 [10] T. L. Griffiths and M. Steyvers. Finding scientific topics. Proceedings of the National academy of Sciences of the United States of America, 101(Suppl 1):5228–5235, 2004.
 [11] M. D. Hoffman, D. M. Blei, and F. R. Bach. Online learning for latent dirichlet allocation. In NIPS, volume 2, page 5, 2010.
 [12] M. I. Jordan, Z. Ghahramani, T. S. Jaakkola, and L. K. Saul. An introduction to variational methods for graphical models. Machine learning, 37(2):183–233, 1999.
 [13] T. Minka. Expectation propagation for approximate bayesian inference. 2001.
 [14] C. Robert, A. Doucet, and S. Godsill. Marginal MAP estimation using markov chain montecarlo. In IEEE Int. Conf. on Acoustics, Speech and Signal Processing, volume 3, pages 1753–1756. IEEE, 1999.
 [15] G. C. Wei and M. A. Tanner. A monte carlo implementation of the em algorithm and the poor man’s data augmentation algorithms. Journal of the American Statistical Association, 85(411):699–704, 1990.
 [16] Y. Weiss. Correctness of local probability propagation in graphical models with loops. Neural Computation, 12(1):1–41, 2000.
 [17] J. Winn and C. Bishop. Variational message passing. Journal of Machine Learning Research, 6:661¡V694, 2005.
 [18] F. Yan, N. Xu, and Y. Qi. Parallel inference for latent dirichlet allocation on graphics processing units. In NIPS, volume 9, pages 2134–2142, 2009.