Streaming Variational Bayes

Introduction

Large, streaming data sets are increasingly the norm in science and technology. Simple descriptive statistics can often be readily computed with a constant number of operations for each data point in the streaming setting, without the need to revisit past data or have advance knowledge of future data. But these time and memory restrictions are not generally available for the complex, hierarchical models that practitioners often have in mind when they collect large data sets. Significant progress on scalable learning procedures has been made in recent years [e.g., 1, 2]. But the underlying models remain simple, and the inferential framework is generally non-Bayesian. The advantages of the Bayesian paradigm (e.g., hierarchical modeling, coherent treatment of uncertainty) currently seem out of reach in the Big Data setting.

An exception to this statement is provided by , who have shown that a class of approximation methods known as variational Bayes (VB) can be usefully deployed for large-scale data sets. They have applied their approach, referred to as stochastic variational inference (SVI), to the domain of topic modeling of document collections, an area with a major need for scalable inference algorithms. VB traditionally uses the variational lower bound on the marginal likelihood as an objective function, and the idea of SVI is to apply a variant of stochastic gradient descent to this objective. Notably, this objective is based on the conceptual existence of a full data set involving $D$ data points (i.e., documents in the topic model setting), for a fixed value of $D$. Although the stochastic gradient is computed for a single, small subset of data points (documents) at a time, the posterior being targeted is a posterior for $D$ data points. This value of $D$ must be specified in advance and is used by the algorithm at each step. Posteriors for $D^{\prime}$ data points, for $D^{\prime}\neq D$, are not obtained as part of the analysis.

We view this lack of a link between the number of documents that have been processed thus far and the posterior that is being targeted as undesirable in many settings involving streaming data. In this paper we aim at an approximate Bayesian inference algorithm that is scalable like SVI but is also truly a streaming procedure, in that it yields an approximate posterior for each processed collection of $D^{\prime}$ data points—and not just a pre-specified “final” number of data points $D$. To that end, we return to the classical perspective of Bayesian updating, where the recursive application of Bayes theorem provides a sequence of posteriors, not a sequence of approximations to a fixed posterior. To this classical recursive perspective we bring the VB framework; our updates need not be exact Bayesian updates but rather may be approximations such as VB. This approach is similar in spirit to assumed density filtering or expectation propagation , but each step of those methods involves a moment-matching step that can be computationally costly for models such as topic models. We are able to avoid the moment-matching step via the use of VB. We also note other related work in this general vein: MCMC approximations have been explored by , and VB or VB-like approximations have also been explored by .

Although the empirical success of SVI is the main motivation for our work, we are also motivated by recent developments in computer architectures, which permit distributed and asynchronous computations in addition to streaming computations. As we will show, a streaming VB algorithm naturally lends itself to distributed and asynchronous implementations.

Streaming, distributed, asynchronous Bayesian updating

Streaming Bayesian updating. Consider data $x_{1},x_{2},\ldots$ generated iid according to a distribution $p(x\>|\>\Theta)$ given parameter(s) $\Theta$. Assume that a prior $p(\Theta)$ has also been specified. Then Bayes theorem gives us the posterior distribution of $\Theta$ given a collection of $S$ data points, $C_{1}:=(x_{1},\ldots,x_{S})$:

where $p(C_{1}\>|\>\Theta)=p(x_{1},\ldots,x_{S}\>|\>\Theta)=\prod_{s=1}^{S}p(x_{s}\>|\>\Theta).$

Suppose we have seen and processed $b-1$ collections, sometimes called minibatches, of data. Given the posterior $p(\Theta\>|\>C_{1},\ldots,C_{b-1})$, we can calculate the posterior after the $b$th minibatch:

That is, we treat the posterior after $b-1$ minibatches as the new prior for the incoming data points. If we can save the posterior from $b-1$ minibatches and calculate the normalizing constant for the $b$th posterior, repeated application of Eq. (1) is streaming; it automatically gives us the new posterior without needing to revisit old data points.

In complex models, it is often infeasible to calculate the posterior exactly, and an approximation must be used. Suppose that, given a prior $p(\Theta)$ and data minibatch $C$, we have an approximation algorithm $\mathcal{A}$ that calculates an approximate posterior $q$: $q(\Theta)=\mathcal{A}(C,p(\Theta))$. Then, setting $q_{0}(\Theta)=p(\Theta)$, one way to recursively calculate an approximation to the posterior is

When $\mathcal{A}$ yields the posterior from Bayes theorem, this calculation is exact. This approach already differs from that of , which we will see (Sec. 3.2) directly approximates $p(\Theta\mid C_{1},\ldots,C_{B})$ for fixed $B$ without making intermediate approximations for $b$ strictly between $1$ and $B$.

Distributed Bayesian updating. The sequential updates in Eq. (2) handle streaming data in theory, but in practice, the $\mathcal{A}$ calculation might take longer than the time interval between minibatch arrivals or simply take longer than desired. Parallelizing computations increases algorithm throughput. And posterior calculations need not be sequential. Indeed, Bayes theorem yields

That is, we can calculate the individual minibatch posteriors $p(\Theta\mid C_{b})$, perhaps in parallel, and then combine them to find the full posterior $p(\Theta\mid C_{1},\ldots,C_{B})$.

Given an approximating algorithm $\mathcal{A}$ as above, the corresponding approximate update would be

for some approximating distribution $q$, provided the normalizing constant for the right-hand side of Eq. (4) can be computed.

Variational inference methods are generally based on exponential family representations , and we will make that assumption here. In particular, we suppose $p(\Theta)\propto\exp\{\xi_{0}\cdot T(\Theta)\};$ that is, $p(\Theta)$ is an exponential family distribution for $\Theta$ with sufficient statistic $T(\Theta)$ and natural parameter $\xi_{0}$. We suppose further that $\mathcal{A}$ always returns a distribution in the same exponential family; in particular, we suppose that there exists some parameter $\xi_{b}$ such that

When we make these two assumptions, the update in Eq. (4) becomes

where the normalizing constant is readily obtained from the exponential family form. In what follows we use the shorthand $\xi\leftarrow\mathcal{A}(C,\xi_{0})$ to denote that $\mathcal{A}$ takes as input a minibatch $C$ and a prior with exponential family parameter $\xi_{0}$ and that it returns a distribution in the same exponential family with parameter $\xi$.

So, to approximate $p(\Theta\>|\>C_{1},\ldots,C_{B})$, we first calculate $\xi_{b}$ via the approximation primitive $\mathcal{A}$ for each minibatch $C_{b}$; note that these calculations may be performed in parallel. Then we sum together the quantities $\xi_{b}-\xi_{0}$ across $b$, along with the initial $\xi_{0}$ from the prior, to find the final exponential family parameter to the full posterior approximation $q$. We previously saw that the general Bayes sequential update can be made streaming by iterating with the old posterior as the new prior (Eq. (2)). Similarly, here we see that the full posterior approximation $q$ is in the same exponential family as the prior, so one may iterate these parallel computations to arrive at a parallelized algorithm for streaming posterior computation.

We emphasize that while these updates are reminiscent of prior-posterior conjugacy, it is actually the approximate posteriors and single, original prior that we assume belong to the same exponential family. It is not necessary to assume any conjugacy in the generative model itself nor that any true intermediate or final posterior take any particular limited form.

Asynchronous Bayesian updating. Performing $B$ computations in parallel can in theory speed up algorithm running time by a factor of $B$, but in practice it is often the case that a single computation thread takes longer than the rest. Waiting for this thread to finish diminishes potential gains from distributing the computations. This problem can be ameliorated by making computations asynchronous. In this case, processors known as workers each solve a subproblem. When a worker finishes, it reports its solution to a single master processor. If the master gives the worker a new subproblem without waiting for the other workers to finish, it can decrease downtime in the system.

Our asynchronous algorithm is in the spirit of Hogwild! . To present the algorithm we first describe an asynchronous computation that we will not use in practice, but which will serve as a conceptual stepping stone. Note in particular that the following scheme makes the computations in Eq. (6) asynchronous. Have each worker continuously iterate between three steps: (1) collect a new minibatch $C$, (2) compute the local approximate posterior $\xi\leftarrow\mathcal{A}(C,\xi_{0})$, and (3) return $\Delta\xi:=\xi-\xi_{0}$ to the master. The master, in turn, starts by assigning the posterior to equal the prior: $\xi^{(\textrm{post})}\leftarrow\xi_{0}$. Each time the master receives a quantity $\Delta\xi$ from any worker, it updates the posterior synchronously: $\xi^{(\textrm{post})}\leftarrow\xi^{(\textrm{post})}+\Delta\xi$. If $\mathcal{A}$ returns the exponential family parameter of the true posterior (rather than an approximation), then the posterior at the master is exact by Eq. (4).

A preferred asynchronous computation works as follows. The master initializes its posterior estimate to the prior: $\xi^{(\textrm{post})}\leftarrow\xi_{0}$. Each worker continuously iterates between four steps: (1) collect a new minibatch $C$, (2) copy the master posterior value locally $\xi^{(\textrm{local})}\leftarrow\xi^{(\textrm{post})}$, (3) compute the local approximate posterior $\xi\leftarrow\mathcal{A}(C,\xi^{(\textrm{local})})$, and (4) return $\Delta\xi:=\xi-\xi^{(\textrm{local})}$ to the master. Each time the master receives a quantity $\Delta\xi$ from any worker, it updates the posterior synchronously: $\xi^{(\textrm{post})}\leftarrow\xi^{(\textrm{post})}+\Delta\xi$.

The key difference between the first and second frameworks proposed above is that, in the second, the latest posterior is used as a prior. This latter framework is more in line with the streaming update of Eq. (2) but introduces a new layer of approximation. Since $\xi^{(\textrm{post})}$ might change at the master while the worker is computing $\Delta\xi$, it is no longer the case that the posterior at the master is exact when $\mathcal{A}$ returns the exponential family parameter of the true posterior. Nonetheless we find that the latter framework performs better in practice, so we focus on it exclusively in what follows.

We refer to our overall framework as SDA-Bayes, which stands for (S)treaming, (D)istributed, (A)synchronous Bayes. The framework is intended to be general enough to allow a variety of local approximations $\mathcal{A}$. Indeed, SDA-Bayes works out of the box once an implementation of $\mathcal{A}$—and a prior on the global parameter(s) $\Theta$—is provided. In the current paper our preferred local approximation will be VB.

Case study: latent Dirichlet allocation

In what follows, we consider examples of the choices for the $\Theta$ prior and primitive $\mathcal{A}$ in the context of latent Dirichlet allocation (LDA) . LDA models the content of $D$ documents in a corpus. Themes potentially shared by multiple documents are described by topics. The unsupervised learning problem is to learn the topics as well as discover which topics occur in which documents.

More formally, each topic (of $K$ total topics) is a distribution over the $V$ words in the vocabulary: $\beta_{k}=(\beta_{kv})_{v=1}^{V}$. Each document is an admixture of topics. The words in document $d$ are assumed to be exchangeable. Each word $w_{dn}$ belongs to a latent topic $z_{dn}$ chosen according to a document-specific distribution of topics $\theta_{d}=(\theta_{dk})_{k=1}^{K}$. The full generative model, with Dirichlet priors for $\beta_{k}$ and $\theta_{d}$ conditioned on respective parameters $\eta_{k}$ and $\alpha$, appears in .

To see that this model fits our specification in Sec. 2, consider the set of global parameters $\Theta=\beta$. Each document $w_{d}=(w_{dn})_{n=1}^{N_{d}}$ is distributed iid conditioned on the global topics. The full collection of data is a corpus $C=w=(w_{d})_{d=1}^{D}$ of documents. The posterior for LDA, $p(\beta,\theta,z\mid C,\eta,\alpha)$, is equal to the following expression up to proportionality:

The posterior for just the global parameters $p(\beta\>|\>C,\eta,\alpha)$ can be obtained from $p(\beta,\theta,z\>|\>C,\eta,\alpha)$ by integrating out the local, document-specific parameters $\theta,z$. As is common in complex models, the normalizing constant for Eq. (7) is intractable to compute, so the posterior must be approximated.

To apply SDA-Bayes to LDA, we use the prior specified by the generative model. It remains to choose a posterior-approximation algorithm $\mathcal{A}$. We consider two possibilities here: variational Bayes (VB) and expectation propagation (EP). Both primitives take Dirichlet distributions as priors for $\beta$ and both return Dirichlet distributions for the approximate posterior of the topic parameters $\beta$; thus the prior and approximate posterior are in the same exponential family. Hence both VB and EP can be utilized as a choice for $\mathcal{A}$ in the SDA-Bayes framework.

Mean-field variational Bayes. We use the shorthand $p_{D}$ for Eq. (7), the posterior given $D$ documents. We assume the approximating distribution, written $q_{D}$ for shorthand, takes the form

The problem of VB is to find the best approximating $q_{D}$, defined as the collection of variational parameters $\lambda,\gamma,\phi$ that minimize the KL divergence from the true posterior: $\textrm{KL}\left(q_{D}\>\|\>p_{D}\right)$. Even finding the minimizing parameters is a difficult optimization problem. Typically the solution is approximated by coordinate descent in each parameter as in Alg. 1. The derivation of VB for LDA can be found in and Sup. Mat. A.1.

Expectation propagation. An EP algorithm for approximating the LDA posterior appears in Alg. 6 of Sup. Mat. B. Alg. 6 differs from , which does not provide an approximate posterior for the topic parameters, and is instead our own derivation. Our version of EP, like VB, learns factorized Dirichlet distributions over topics.

2 Other single-pass algorithms for approximate LDA posteriors

The algorithms in Sec. 3.1 pass through the data multiple times and require storing the data set in memory—but are useful as primitives for SDA-Bayes in the context of the processing of minibatches of data. Next, we consider two algorithms that can pass through a data set just one time (single pass) and to which we compare in the evaluations (Sec. 4).

Stochastic variational inference. VB uses coordinate descent to find a value of $q_{D}$, Eq. (8), that locally minimizes the KL divergence, $\textrm{KL}\left(q_{D}\>\|\>p_{D}\right)$. Stochastic variational inference (SVI) is exactly the application of a particular version of stochastic gradient descent to the same optimization problem. While stochastic gradient descent can often be viewed as a streaming algorithm, the optimization problem itself here depends on $D$ via $p_{D}$, the posterior on $D$ data points. We see that, as a result, $D$ must be specified in advance, appears in each step of SVI (see Alg. 2), and is independent of the number of data points actually processed by the algorithm. Nonetheless, while one may choose to visit $D^{\prime}\neq D$ data points or revisit data points when using SVI to estimate $p_{D}$ , SVI can be made single-pass by visiting each of $D$ data points exactly once and then has constant memory requirements. We also note that two new parameters, $\tau_{0}>0$ and $\kappa\in(0.5,1]$, appear in SVI, beyond those in VB, to determine a learning rate $\rho_{t}$ as a function of iteration $t$: $\rho_{t}:=(\tau_{0}+t)^{-\kappa}$.

Sufficient statistics. On each round of VB (Alg. 1), we update the local parameters for all documents and then compute $\lambda_{kv}\leftarrow\eta_{kv}+\sum_{d=1}^{D}\phi_{dvk}n_{dv}$. An alternative single-pass (and indeed streaming) option would be to update the local parameters for each minibatch of documents as they arrive and then add the corresponding terms $\phi_{dvk}n_{dv}$ to the current estimate of $\lambda$ for each document $d$ in the minibatch. This essential idea has been proposed previously for models other than LDA by and forms the basis of what we call the sufficient statistics update algorithm (SSU): Alg. 3. This algorithm is equivalent to SDA-Bayes with $\mathcal{A}$ chosen to be a single iteration over the global variable $\lambda$ of VB (i.e., updating $\lambda$ exactly once instead of iterating until convergence).

Evaluation

We follow (and further ) in evaluating our algorithms by computing (approximate) predictive probability. Under this metric, a higher score is better, as a better model will assign a higher probability to the held-out words.

We calculate predictive probability by first setting aside held-out testing documents $C^{(\textrm{test})}$ from the full corpus and then further setting aside a subset of held-out testing words $W_{d,\textrm{test}}$ in each testing document $d$. The remaining (training) documents $C^{(\textrm{train})}$ are used to estimate the global parameter posterior $q(\beta)$, and the remaining (training) words $W_{d,\textrm{train}}$ within the $d$th testing document are used to estimate the document-specific parameter posterior $q(\theta_{d})$. In all cases, we estimate $q(\theta_{d})$ for evaluative purposes using VB since direct EP estimation takes prohibitively long. To calculate predictive probability, an approximation is necessary since we do not know the predictive distribution—just as we seek to learn the posterior distribution. Specifically, we calculate the normalized predictive distribution and report “log predictive probability” as

To facilitate comparison with SVI, we use the Wikipedia and Nature corpora of in our experiments. These two corpora represent a range of sizes (3,611,558 training documents for Wikipedia and 351,525 for Nature) as well as different types of topics. We expect words in Wikipedia to represent an extremely broad range of topics whereas we expect words in Nature to focus more on the sciences. We further use the vocabularies of and SVI code available online at . We hold out $10{,}000$ Wikipedia documents and $1{,}024$ Nature documents (not included in the counts above) for testing. In the results presented in the main text, we follow in fitting an LDA model with $K=100$ topics and hyperparameters chosen as: $\forall k,\alpha_{k}=1/K$, $\forall(k,v),\eta_{kv}=0.01$. For both Wikipedia and Nature, we set the parameters in SVI according to the optimal values of the parameters described in Table 1 of (number of documents $D$ correctly set in advance, step size parameters $\kappa=0.5$ and $\tau_{0}=64$).

Figs. 3(a) and 3(b) demonstrate that both SVI and SDA are sensitive to minibatch size when $\eta_{kv}=0.01$, with generally superior performance at larger batch sizes. Interestingly, both SVI and SDA performance improve and are steady across batch size when $\eta_{kv}=1$ (Figs. 3(a) and 3(b)). Nonetheless, we use $\eta_{kv}=0.01$ in what follows in the interest of consistency with . Moreover, in the remaining experiments, we use a large minibatch size of $2^{15}=32{,}768$. This size is the largest before SVI performance degrades in the Nature data set (Fig. 3(b)).

Performance and timing results are shown in Table 1. One would expect that with additional streaming capabilities, SDA-Bayes should show a performance loss relative to SVI. We see from Table 1 that such loss is small in the single-thread case, while SSU performs much worse. SVI is faster than single-thread SDA-Bayes in this single-pass setting.

Full SDA-Bayes improves run time with no performance cost. We handicap SDA-Bayes in the above comparisons by utilizing just a single thread. In Table 1, we also report performance of SDA-Bayes with 32 threads and the same minibatch size. In the synchronous case, we consider minibatch size to equal the total number of data points processed per round; therefore, the minibatch size equals the number of data points sent to each thread per round times the total number of threads. In the asynchronous case, we analogously report minibatch size as this product.

Fig. 2 shows the performance of SDA-Bayes when we run with $\{1,2,4,8,16,32\}$ threads while keeping the minibatch size constant. The goal in such a distributed context is to improve run time while not hurting performance. Indeed, we see dramatic run time improvement as the number of threads grows and in fact some slight performance improvement as well. We tried both a parallel version and a full distributed, asynchronous version of the algorithm; Fig. 2 indicates that the speedup and performance improvements we see here come from parallelizing—which is theoretically justified by Eq. (3) when $\mathcal{A}$ is Bayes rule. Our experiments indicate that our Hogwild!-style asynchrony does not hurt performance. In our experiments, the processing time at each thread seems to be approximately equal across threads and dominate any communication time at the master, so synchronous and asynchronous performance and running time are essentially identical. In general, a practitioner might prefer asynchrony since it is more robust to node failures.

SVI is sensitive to the choice of total data size $D$. The evaluations above are for a single posterior over $D$ data points. Of greater concern to us in this work is the evaluation of algorithms in the streaming setting. We have seen that SVI is designed to find the posterior for a particular, pre-chosen number of data points $D$. In practice, when we run SVI on the full data set but change the input value of $D$ in the algorithm, we can see degradations in performance. In particular, we try values of $D$ equal to $\{0.01,0.1,1,10,100\}$ times the true $D$ in Fig. 3(c) for the Wikipedia data set and in Fig. 3(d) for the Nature data set.

A practitioner in the streaming setting will typically not know $D$ in advance, or multiple values of $D$ may be of interest. Figs. 3(c) and 3(d) illustrate that an estimate may not be sufficient. Even in the case where $D$ is known in advance, it is reasonable to imagine a new influx of further data. One might need to run SVI again from the start (and, in so doing, revisit the first data set) to obtain the desired performance.

SVI is sensitive to learning step size. use cross-validation to tune step-size parameters $(\tau_{0},\kappa)$ in the stochastic gradient descent component of the SVI algorithm. This cross-validation requires multiple runs over the data and thus is not suited to the streaming setting. Figs. 3(e) and 3(f) demonstrate that the parameter choice does indeed affect algorithm performance. In these figures, we keep $D$ at the true training data size.

have observed that the optimal $(\tau_{0},\kappa)$ may interact with minibatch size, and we further observe that the optimal values may vary with $D$ as well. We also note that recent work has suggested a way to update $(\tau_{0},\kappa)$ adaptively during an SVI run .

EP is not suited to LDA. Earlier attempts to apply EP to the LDA model in the non-streaming setting have had mixed success, with in particular finding that EP performance can be poor for LDA and, moreover, that EP requires “unrealistic intermediate storage requirements.” We found this to also be true in the streaming setting. We were not able to obtain competitive results with EP; based on an 8-thread implementation of SDA-Bayes with an EP primitiveWe chose 8 threads since any fewer was too slow to get results and anything larger created too high of a memory demand on our system., after over 91 hours on Wikipedia (and $6.7\times 10^{4}$ data points), log predictive probability had stabilized at around $-7.95$ and, after over 97 hours on Nature (and $9.7\times 10^{4}$ data points), log predictive probability had stabilized at around $-8.02$. Although SDA-Bayes with the EP primitive is not effective for LDA, it remains to be seen whether this combination may be useful in other domains where EP is known to be effective.

Discussion

We have introduced SDA-Bayes, a framework for streaming, distributed, asynchronous computation of an approximate Bayesian posterior. Our framework makes streaming updates to the estimated posterior according to a user-specified approximation primitive. We have demonstrated the usefulness of our framework, with variational Bayes as the primitive, by fitting the latent Dirichlet allocation topic model to the Wikipedia and Nature corpora. We have demonstrated the advantages of our algorithm over stochastic variational inference and the sufficient statistics update algorithm, particularly with respect to the key issue of obtaining approximations to posterior probabilities based on the number of documents seen thus far, not posterior probabilities for a fixed number of documents.

We thank M. Hoffman, C. Wang, and J. Paisley for discussions, code, and data and our reviewers for helpful comments. TB is supported by the Berkeley Fellowship, NB by a Hertz Foundation Fellowship, and ACW by the Chancellor’s Fellowship at UC Berkeley. This research is supported in part by NSF CISE Expeditions award CCF-1139158, DARPA XData Award FA8750-12-2-0331, and AMPLab sponsor donations from Amazon Web Services, Google, SAP, Blue Goji, Cisco, Clearstory Data, Cloudera, Ericsson, Facebook, General Electric, Hortonworks, Intel, Microsoft, NetApp, Oracle, Samsung, Splunk, VMware and Yahoo!. This material is based upon work supported in part by the Office of Naval Research under contract/grant number N00014-11-1-0688.

References

Appendix A Variational Bayes

As described in the main text, the idea of VB is to find the distribution $q_{D}$ that best approximates the true posterior, $p_{D}$. More specifically, the optimization problem of VB is defined as finding a $q_{D}$ to minimize the KL divergence between the approximating distribution and the posterior:

Typically $q_{D}$ takes a particular, constrained form, and finding the optimal $q_{D}$ amounts to finding the optimal parameters for $q_{D}$. Moreover, the optimal parameters usually cannot be expressed in closed form, so often a coordinate descent algorithm is used.

For the LDA model, we have $q_{D}$ in the form of Eq. (8) and $p_{D}$ defined by Eq. (7). We wish to find the following variational parameters (i.e., parameters to $q_{D}$): $\lambda$ (describing each topic), $\gamma$ (describing the topic proportions in each document), and $\phi$ (describing the assignment of each word in each document to a topic).

Finding $q_{D}$ to minimize the KL divergence between $q_{D}$ and $p_{D}$ is equivalent to finding $q_{D}$ to maximize the evidence lower bound (ELBO),

since $p(x_{1:D})$ is constant in $q_{D}$. The VB optimization problem is often phrased in terms of the ELBO instead of the KL divergence.

The ELBO for LDA can be written as follows, where the model parameters are $\beta,\theta,z$ and the data is $w$; $\eta$ and $\alpha$ are fixed hyperparameters.

The expectations in $q$ in the previous equation can be evaluated as follows. The equations below make use of the digamma function $\psi$ and trigamma function $\psi_{1}$. Here,

A.1.2 Coordinate ascent

We maximize the ELBO via coordinate ascent in each dimension of the variational parameters: $\lambda$, $\gamma$, and $\phi$.

Choose a topic index $k$. Fix $\gamma$, $\phi$, and each $\lambda_{j}$ for $j\neq k$. Then we can write the ELBO’s functional dependence on $\lambda_{k}$ as follows, where “const” is a constant in $\lambda_{k}$.

The partial derivative of $\textrm{ELBO}(\lambda_{k})$ with respect to one of the dimensions of $\lambda_{k}$, say $\lambda_{kv}$, is

From the last line of the previous equation, we see that one can set the gradient of $\textrm{ELBO}(\lambda_{k})$ to zero by setting

Equivalently, if $n_{dv}$ is the number of occurrences (tokens) of word type $v$ in document $d$, then the update may be written

Now choose a document $d$. Fix $\lambda$, $\phi$, and $\gamma_{c}$ for $c\neq d$. Then we can express the functional dependence of the ELBO on $\gamma_{d}$ as follows.

The partial derivative of $\textrm{ELBO}(\gamma_{d})$ with respect to one of the dimensions of $\gamma_{d}$, say $\gamma_{dk}$, is

As for the $\lambda$ case above, one obvious way to achieve a gradient of $\textrm{ELBO}(\gamma_{d})$ equal to zero is to set

Finally, consider fixing $\lambda$, $\gamma$, and $\phi_{cu}$ for $(c,u)\neq(d,v)$. In this case, the dependence of the ELBO on $\phi_{dv}$ can be written as follows.

The partial derivative of $\textrm{ELBO}(\phi_{dv})$ with respect to one of the dimensions of $\phi_{dv}$, say $\phi_{dvk}$, is

Using the method of Lagrange multipliers to incorporate the constraint that $\sum_{k=1}^{K}\phi_{dvk}=1$, we wish to find $\rho$ and $\phi_{dvk}$ such that

achieves the desired outcome in Eq. (9). Here, $\propto_{k}$ indicates that the proportionality is across $k$. The optimal choice of $\rho$ is expressed via this proportionality. The above assignment may also be written as

The coordinate-ascent algorithm iteratively updates the parameters $\lambda$, $\gamma$, and $\phi$. In practice, we usually iterate the updates for the “local” parameters $\phi$ and $\gamma$ until they converge, then update the “global” parameter $\lambda$, and repeat. The resulting batch variational Bayes algorithm is presented in Alg. 1.

A.2 SDA-Bayes VB

For a fixed hyperparameter $\alpha$, we can think of BatchVB as an algorithm that takes input in the form of a prior on topic parameters $\beta$ and a minibatch of documents. In particular, let $C_{b}$ be the $b$th minibatch of documents; for documents with indices in $\mathcal{D}_{b}$, these documents can be summarized by the word counts $(n_{d})_{d\in\mathcal{D}_{b}}$. Then, in the notation of Eq. (2), we have $\Theta=\beta$, $\mathcal{A}=\textrm{BatchVB}$, and

In general, the $b$th posterior takes the same form and therefore can be summarized by its parameters $\lambda^{(b)}$:

In this case, if we set the prior parameters to $\lambda^{(0)}_{k}:=\eta_{k}$, Eq. (2) becomes the following algorithm.

Next, we apply the asynchronous, distributed updates described in the “Asynchronous Bayesian updating” portion of Sec. 2 to the batch VB primitive and LDA model. In this case, $\lambda^{(\textrm{post})}$ is the posterior parameter estimate maintained at the master, and each worker updates this value after a local computation. The posterior after seeing a collection of minibatches is $q(\beta)=\prod_{k=1}^{K}\textrm{Dirichlet}(\beta_{k}|\lambda^{(\textrm{post})}_{k})$.

Appendix B Expectation Propagation

Our batch expectation propagation (EP) algorithm for LDA learns a posterior for both the document-specific topic mixing proportions $(\theta_{d})_{d=1}^{D}$ and the topic distributions over words $(\beta_{k})_{k=1}^{K}$. By contrast, the algorithm in learns only the former and so is not appropriate to the model in Sec. 3.

For consistency, we also follow in making a distinction between token and type word updates, where a token refers to a particular word instance and a type refers to all words with the same vocabulary value. Let $C=(w_{d})_{d=1}^{D}$ denote the set of documents that we observe, and for each word $v$ in the vocabulary, let $n_{dv}$ denote the number of times $v$ appears in document $d$.

We begin by collapsing (i.e., integrating out) the word assignments $z$ in the posterior (7) of LDA. We can express the collapsed posterior as

For each document-word pair $(d,v)$, consider approximating the term $\sum_{k=1}^{K}\theta_{dk}\beta_{kv}$ above by

where $q(\beta_{k}\mid\lambda_{k})=\textrm{Dirichlet}(\beta_{k}\mid\lambda_{k})$ and $q(\theta_{d}\mid\gamma_{d})=\textrm{Dirichlet}(\theta_{d}\mid\gamma_{d})$, with the parameters

Optimization problem. We seek to find the optimal parameters $(\chi,\zeta)$ by minimizing the (reverse) KL divergence:

This joint minimization problem is not tractable, and the idea of EP is to proceed iteratively by fixing most of the factors in Eq. (10) and minimizing the KL divergence over the parameters related to a single word.

More formally, suppose we already have a set of parameters $(\chi,\zeta)$. Consider a document $d$ and word $v$ that occurs in document $d$ (i.e., $n_{dv}\geq 1$). We start by removing the component of $q$ related to $(d,v)$ in Eq. (10). Following , we subtract out the effect of one occurrence of word $v$ in document $d$, but at the end of this process we update the distribution on the type level. In doing so, we use the following shorthand for the remaining global parameters:

We obtain an improved estimate of the posterior $q$ by updating the parameters from $(\lambda,\gamma)$ to $(\hat{\lambda},\hat{\gamma})$, where

Solution to the optimization problem. First, note that for $d^{\prime}:d^{\prime}\neq d$, we have $\hat{\gamma}_{d^{\prime}}=\gamma_{d^{\prime}}$.

Now consider the index $d$ chosen on this iteration. Since $\beta$ and $\theta$ are Dirichlet-distributed under $q$, the minimization problem in Eq. (12) reduces to solving the moment-matching equations

These can be solved via Newton’s method though recommends solving exactly for the first and “average second” moments of $\beta_{ku}$ and $\theta_{dk}$, respectively, instead. We choose the latter approach for consistency with ; our own experiments also suggested taking the approach of was faster than Newton’s method with no noticeable performance loss. The resulting moment updates are

We then set $(\chi_{kdv})_{k=1}^{K}$ and $\zeta_{dv}$ such that the new global parameters $(\lambda_{k})_{k=1}^{K}$ and $\gamma_{d}$ are equal to the optimal parameters $(\hat{\lambda}_{k})_{k=1}^{K}$ and $\hat{\gamma}_{d}$. The resulting algorithm is presented below (Alg. 6).

The results in the main text (Sec. 4) are reported for Alg. 6. We also tried a slightly modified EP algorithm that makes token-level updates to parameter values, rather than type-level updates. This modified version iterates through each word placeholder in document $d$; that is, through pairs $(d,n)$ rather than pairs $(d,v)$ corresponding to word values. Since there are always at least as many $(d,n)$ pairs as $(d,v)$ pairs with $n_{dv}\geq 1$ (and usually many more of the former), the modified algorithm requires many more iterations. In practice, we find better experimental performance for the modified EP algorithm in terms of log predictive probability as a function of number of data points in the training set seen so far: e.g., leveling off at about $-7.96$ for Nature vs. $-8.02$. However, the modified algorithm is also much slower, and still returns much worse results than SDA-Bayes or SVI, so we do not report these results in the main text.Here and in the main text we run EP with $\eta=1$. We also tried EP with $\eta=0.01$, but the positivity check for $\lambda^{\setminus(d,v)}_{ku}$ and $\gamma^{\setminus(d,v)}_{dk}$ on line $(\dagger)$ in Algorithm 6 always failed and as a result none of the parameters were updated.

B.2 SDA-Bayes EP

Putting a batch EP algorithm for LDA into the SDA-Bayes framework is almost identical to putting a batch VB algorithm for LDA into the SDA-Bayes framework. This similarity is to be expected since SDA-Bayes works out of the box with a batch approximation algorithm in the correct form.

For a fixed hyperparameter $\alpha$, we can think of BatchEP as an algorithm (just like BatchVB) that takes input in the form of a prior on topic parameters $\beta$ and a minibatch of documents. The same setup and notation from Sup. Mat. A.2 applies. In this case, Eq. (2) becomes the following algorithm.

This algorithm is exactly the same as Alg. 4 but with a batch EP primitive instead of a batch VB primitive.

Next, we apply the asynchronous, distributed updates described in the “Asynchronous Bayesian updating” portion of Sec. 2 to the batch EP primitive and LDA model. Again, the setup and notation from Sup. Mat. A.2 applies, and we find the following algorithm.

Indeed, the recipe outlined here applies more generally to other primitives besides EP and VB.