Structured Stochastic Variational Inference

Introduction

Hierarchical Bayesian modeling is a powerful framework for learning from rich data sources. Unfortunately, the intractability of the posteriors of rich models drives practitioners to resort to approximate inference algorithms such as mean-field variational inference or Markov chain Monte Carlo (MCMC). These classes of methods have complementary strengths and weaknesses—MCMC methods have strong asymptotic guarantees of unbiasedness but are often slow, while mean-field variational inference is often faster but tends to misrepresent important qualities of the posterior of interest and is more vulnerable to local optima. Incremental versions of both methods based on stochastic optimization have been developed that are applicable to large datasets (Welling and Teh,, 2011; Hoffman et al.,, 2013).

In this paper we focus on variational inference. In particular, we are interested in using variational inference to estimate the parameters of high-dimensional Bayesian models with highly multimodal posteriors, such as mixture models, topic models, and factor models. We will be less concerned with uncertainty estimates, since they are difficult to trust and interpret in this setting.

Mean-field variational inference approximates the intractable posterior distribution implied by the model and data with a factorized approximating distribution in which all parameters are independent. This mean-field distribution is then tuned to minimize its Kullback-Leibler divergence to the posterior, which is equivalent to maximizing a lower bound on the marginal probability of the data. The restriction to factorized distributions makes the problem tractable, but reduces the fidelity of the approximation and introduces local optima (Wainwright and Jordan,, 2008).

A partial remedy is to weaken the mean-field factorization by restoring some dependencies, resulting in “structured” mean-field approximations (Saul and Jordan,, 1996). The applicability, speed, effectiveness, and ease-of-implementation of standard structured mean-field algorithms is limited because the lower bound implied by the structured distribution must be available in closed form.

More recent work manages these intractable variational lower bounds using stochastic optimization, which allows one to optimize functions that can only be computed approximately. For example, Ji et al., (2010) use mean-field approximations to the posteriors of “collapsed” models where some parameters have been analytically marginalized out, Salimans and Knowles, (2013) apply a structured approximation to the posterior of a stochastic volatility model, and Mimno et al., (2012) use a structured approximation to the posterior of a collapsed model.

In parallel, Hoffman et al., (2013) proposed the stochastic variational inference (SVI) framework, which uses stochastic optimization to apply mean-field variational inference to massive datasets. SVI splits the unobserved variables in a hierarchical model into global parameters $\beta$ (which are shared across all observations) and groups of local hidden variables $z_{1},\ldots,z_{N}$ (each of which is specific to a small group of observations $y_{n}$). The goal is to minimize the Kullback-Leibler (KL) divergence between a tractable approximating distribution $q(z,\beta)$ over the local and global parameters and the true posterior $p(z,\beta|y)$ over those parameters. SVI approximates posteriors much more quickly than traditional batch variational inference algorithms.

Like batch variational inference, SVI relies on the mean-field approximation, which requires that $q$ factorize as $(\prod_{k}q(\beta_{k}))\prod_{n,m}q(z_{n,m})$; that is, SVI approximates the joint posterior $p(z,\beta|y)$ with a distribution $q$ that cannot represent any dependencies between random variables. Mimno et al., (2012) proposed a variant for latent Dirichlet allocation that restores dependencies between sets of local hidden variables $z_{n,1:M}$ so that the approximating distribution has the form $q(z,\beta)=(\prod_{k}q(\beta_{k}))\prod_{n}q(z_{n,1:M})$, improving $q$’s ability to approximate the posterior $p(z,\beta|y)$. But their method still breaks the dependence between the global and local variables $\beta$ and $z$.

In this paper, we introduce structured stochastic variational inference (SSVI), a generalization of the SVI framework that can restore the dependence between global and local variables. SSVI approximates the posteriors $p(z,\beta|y)$ of a wide class of models with distributions of the form $q(z,\beta)=(\prod_{k}q(\beta_{k}))\prod_{n}q(z_{n}|\beta)$, allowing for arbitrary dependencies between $\beta$ and $z_{n}$.

In experiments on three models and datasets, we find that in all cases restoring these dependencies allows SSVI to find qualitatively and quantitatively better parameter estimates, avoiding the local optima and sensitivity to hyperparameters that plague mean-field variational inference.

Structured Stochastic Variational Inference

In this section, we will present two SSVI algorithms. We first review the class of models to which SSVI can be applied and the variational distributions that it employs.

As in SVI (Hoffman et al.,, 2013), we assume we have $N$ groups of observations $y_{1:N}$ and a probability model that factorizes as $\textstyle p(y,z,\beta)=p(\beta)\prod_{n}p(y_{n},z_{n}|\beta).$ The independence structure of such a model is visualized in figure 1. The global parameters $\beta$ are shared across all observations, and the local hidden variables $z_{1:N}$ are conditionally independent of one another given the global parameters $\beta$.

We will restrict our attention to conditionally conjugate models. We assume that the prior $p(\beta)$ is in a tractable exponential family $p(\beta)=h(\beta)\exp\{\eta\cdot t(\beta)-A(\eta)\}$. The base measure $h$ and log-normalizer $A$ are scalar-valued functions, $\eta$ is a vector of natural parameters, and $t(\beta)$ is a vector-valued sufficient statistic function. We futher assume that the joint likelihood of the local variables $y_{n}$ and $z_{n}$ given $\beta$ is of the form $p(y_{n},z_{n}|\beta)=\exp\{t(\beta)\cdot\eta_{n}(y_{n},z_{n})+g_{n}(y_{n},z_{n})\}$ (where $g_{n}$ is a real-valued function and $\eta_{n}$ is a vector-valued function).

This form for $p(y_{n},z_{n},\beta)$ includes all conjugate pairs of distributions $p(\beta)$, $p(y_{n},z_{n}|\beta)$ (Gelman et al.,, 2013); that is, it is the most general family of distributions for which the conditional $p(\beta|y,z)$ is in the same family as the prior $p(\beta)$. This conditional is

These restrictions are a weaker version of those imposed by Hoffman et al., (2013); the difference is that we make no assumptions about the tractability of the conditional distributions $p(z_{n}|y_{n},\beta)$ or $p(z_{n,m}|y_{n},z_{n,\backslash m},\beta)$. This work is therefore applicable to any model that fits in the SVI framework, including mixture models, LDA, hidden Markov models (HMMs), factorial HMMs, Kalman filters, factor analyzers, probabilistic matrix factorizations, hierarchical linear regression, hierarchical probit regression, and many other hierarchical models. Unlike SVI, it can also address models without tractable local conditionals, such as multilevel logistic regressions (Gelman and Hill,, 2007) or the correlated topic model (Blei and Lafferty,, 2006).

2 Approximating Distribution

Our goal is to approximate the intractable posterior $p(z,\beta|y)$ with a distribution $q(z,\beta)$ in some restricted, tractable family. We will choose a $q$ distribution from this family by solving an optimization problem, minimizing the Kullback-Leibler (KL) divergence between $q(z,\beta)$ and the posterior $p(z,\beta|y)$.

The simplest approach is to make the mean-field approximation, restricting $q$ to factorize so that $q(z,\beta)=q(\beta)\prod_{n}\prod_{m}q(z_{n,m}).$ This restriction dramatically simplifies the form of the KL-divergence between $q$ and the posterior, but this simplicity comes at a price. Every dependence that we break to make $q$ easier to work with makes $q$ less able to closely approximate the posterior $p(z,\beta|y)$. Breaking dependencies may also introduce additional local minima into the KL divergence between $q$ and the posterior; imposing independence assumptions places nonconvex constraints on the dual of the solution space, which may block the path from a bad solution to a good one (Wainwright and Jordan,, 2008).

Structured mean-field partially relaxes the mean-field independence restriction (Saul and Jordan,, 1996). Traditional structured mean-field algorithms require the practitioner to identify and exploit some model-specific structure; for example, Ghahramani and Jordan, (1997) exploited the availability of dynamic programming algorithms for HMMs to derive a structured mean-field algorithm for factorial HMMs.

Mimno et al., (2012) proposed a structured stochastic variational inference algorithm for latent Dirichlet allocation that depends less on model-specific structure, placing no restrictions on the joint distributions $q(z_{n,1:M})$ so that $q(z,\beta)$ factorizes as $q(z,\beta)=(\prod_{k}q(\beta_{k}))\prod_{n}q(z_{n}).$ The optimal $q(z_{n})$ may not be tractable to normalize, but it can still be sampled from using Markov chain Monte Carlo (MCMC), which is all that is necessary to generate a stochastic natural gradient for a stochastic variational inference algorithm. The result was a significant improvement in the quality of the inference algorithm’s ability to obtain high-quality estimates of model parameters. However, the approximate posterior of Mimno et al., (2012) still breaks the dependence between global parameters $\beta$ and local hidden variables $z$.

We introduce a framework for structured stochastic variational inference (SSVI) algorithms that restore the dependence between $\beta$ and $z$. Our variational distribution $q$ is of the form

The only remaining factorization we impose is between the elements of $\beta$; the conditional independence between the $z_{n}$s given $\beta$ is implied by the model in equation 1.

We will restrict $q(\beta)$ to be in the same exponential family as the prior $p(\beta)$, so that $q(\beta)=h(\beta)\exp\{\lambda\cdot t(\beta)-A(\lambda)\}.$ $\lambda$ is a vector of free parameters that controls $q(\beta)$. We also require that any dependence under $q$ between $z_{n}$ and $\beta$ be mediated by some vector-valued function $\gamma_{n}(\beta)$, so that we may write $q(z_{n}|\beta)=q(z_{n}|\gamma_{n}(\beta))$.

This form for $q$ allows for rich dependencies between nearly all variables in the model. This comes at a cost, however. In mean-field variational inference, we proceed by maximizing a lower bound on the marginal probability of the data; this is equivalent to minimizing the KL divergence from $q$ to the posterior (Bishop,, 2006). However, this lower bound contains expectations that become impossible to compute when we allow $z_{n}$ to depend on $\beta$ in $q$. This issue may seem insurmountable, but even though we cannot compute the variational lower bound, we can still optimize it using stochastic optimization.

3 The Structured Variational Objective

Our goal is to find a distribution $q(\beta,z)$ that has low KL divergence to the posterior $p(\beta,z|y)$. The KL divergence between $q$ and the full posterior is

Because the KL divergence must be non-negative, this yields the evidence lower bound (ELBO)

We used the conditional independence structure assumed in section 2.1 to break $\log p(y,z|\beta)$ into a sum over $n$. Our goal is to maximize the ELBO subject to some restrictions on $q$.

Before describing the form we choose for $\gamma_{n}(\beta)$, we first note that the second integral in equation 2.3 is itself a lower bound on the marginal probability of the $n$th group of observations:

Thus, for any particular value of $\beta$ we can maximize the global ELBO over $q(z_{n}|\beta)$ by minimizing the KL divergence between $q(z_{n}|\beta)$ and $p(z_{n}|y_{n},\beta)$. We will assume that the function $\gamma_{n}(\beta)$ that controls $q(z_{n}|\beta)=q(z_{n}|\gamma_{n}(\beta))$ is defined to do just that, so that $\gamma_{n}(\beta)$ is at a local maximum of this “local ELBO”, i.e.,

The function $\gamma_{n}(\beta)$ may be implicit; e.g., it might be evaluated by solving an optimization problem.

4 Algorithms

where $V(\beta,\lambda)$ is a matrix that is defined in terms of the cumulative distribution functions (CDFs) quantile functions (inverse-CDFs) of $q(\beta)$. Defining the CDF $Q_{k}(\beta_{k})\equiv\int_{-\infty}^{\beta_{k}}q(\beta_{k}^{\prime})d\beta_{k}^{\prime}$ and the quantile function $R_{k}(Q_{k}(\beta_{k}))\equiv\beta_{k}$, $V(\beta,\lambda)$ is defined as the product of two matrices: the inverse of the second derivative of the log-normalizer $A$ of $q$, and the Jacobian of $t(R(Q(\beta)))$ with respect to $\lambda$:

which avoids the difficulty and (modest) expense of computing derivatives of quantiles and log-normalizers.

5 A matrix of approaches

We are free to choose any form for $q(z_{n}|\beta)$ and any unbiased estimator $\hat{\eta}$ that we want; different choices have different properties:

Mean-field:

Non-bound-preserving approaches:

Although we derived SSVI assuming that $q(z_{n}|\beta)$ is chosen to maximize the local ELBO $\mathcal{L}_{n}$, one could obtain a distribution over $z_{n}$ in other ways. For example, for latent Dirichlet allocation one could adapt the CVB0 method of Asuncion et al., (2009) to use a fixed value of $\beta$, resulting in an algorithm akin to that of Foulds et al., (2013).

All of these choices of local variational distribution could also be used in a traditional mean-field setup where $q(\beta,z)=q(\beta)q(z)$, or as part of a variational maximum a posteriori (MAP) estimation algorithm. So for any model that enjoys conditional conjugacy we have a matrix of possible variational inference algorithms: we can match any “E-step” (e.g. mean-field or sampling from the exact conditional) used to approximate $p(z_{n}|y_{n},\beta)$ with any “M-step” (e.g. MAP, mean-field, SSVI, SSVI-A) used to update our approximation to $p(\beta|y)$.

6 Extensions

There are several ways in which the basic algorithms presented above can be extended:

As in SVI, we can compute an unbiased estimate of the sum over $n$ in equation 7 by only computing $\hat{\eta}_{n}$ for some randomly sampled subset of $S$ observations, resulting in the update

For large datasets, the reduced computational effort of only looking at a fraction of the data far outweighs the noise that this subsampling introduces. Taking a cue from the recent work of Broderick et al., (2013) and Wang and Blei, (2012), we suggest gradually ramping up the multiplier $N$ over the course of the first sweep over the dataset.

Hyperparameter updates and parameter hierarchies.

As in the mean-field stochastic variational inference framework of Hoffman et al., (2013), we can optimize any hyperparameters in our model by taking steps in the direction of the gradient of the ELBO with respect to those hyperparameters. We can also extend the framework developed in this paper to models with hierarchies of global parameters as in appendix A of (Hoffman et al.,, 2013).

Related Work

The idea of sampling from global variational distributions to optimize intractable variational inference problems has been proposed previously in several contexts. Ji et al., (2010); Nott et al., (2012); Gerrish, (2013); Paisley et al., (2012), and Ranganath et al., (2014) proposed sampling without a change of variables as a way of coping with non-conjugacy. Kingma and Welling, (2014) and Titsias and Lázaro-Gredilla, (2014) proposed methods that do use a change of variables, although their methods focus more on speed and/or dealing with nonconjugacy than on improving the accuracy of the variational approximation. Salimans and Knowles, (2013) also suggest using a change of variables as a way of dramatically reducing the variance of a stochastic gradient estimator.

Although some of the above methods use stochastic optimization to improve the quality of the mean-field approximation, there are major differences between these methods and SSVI. Ji et al., (2010) apply their method to models where some parameters have been analytically marginalized out, but do not consider explicitly structured variational distributions. Also, in our informal experiments we found that the variance of their gradient estimator was unacceptably high for high-dimensional problems.

Salimans and Knowles, (2013) apply their stochastic linear regression method to structured variational distributions in which the natural parameters of lower-level variational distributions are explicit functions of draws from higher-level variational distributions. In comparison, the implicit form we choose for the conditional variational distribution $q(z|\beta)$ allows for more complicated dependencies between global and local parameters. Also, their regression-based approach requires storing and multiplying matrices that become impractically large for high-dimensional problems.

A final difference is that the papers mentioned above do not exploit conjugacy relationships, which are central to the ease of implementation and efficiency of SSVI-A.

Two other related algorithms are due to Mimno et al., (2012) (which SSVI generalizes) and Wang and Blei, (2012). The algorithm of Wang and Blei, (2012) also uses sampling and conjugacy relationships to attempt to restore the dependencies broken in mean-field algorithms, although their method lacks guarantees of convergence or correctness.

Experiments

In this section we empirically evaluate SSVI and SSVI-A’s ability to estimate parameters for three hierarchical Bayesian models. In each case, we find that relaxing the mean-field approximation allows SSVI and SSVI-A to find significantly better parameter estimates than mean-field. We also find evidence suggesting that this superior performance is primarily due to SSVI/SSVI-A’s ability to avoid local optima.

We evaluated the quality of parameter estimates from SSVI and SSVI-A on the latent Dirichlet allocation (LDA) topic model fit to the 3,800,000-document Wikipedia dataset from (Hoffman et al.,, 2013). We compared with full mean-field stochastic variational inference (Blei et al.,, 2003; Hoffman et al., 2010a, ), a mean-field M-step with a Gibbs sampling E-step (Mimno et al.,, 2012), SSVI with Gibbs, and SSVI-A with Gibbs. Results for other E-step/M-step combinations are in appendix D.

To speed up learning, each update we subsample a minibatch of 1,000 documents rather than analyzing the whole dataset each iteration. We also experimented with various settings of the hyperparameters $\alpha$ and $\eta$, which mean-field variational inference for LDA is known to be quite sensitive to (Asuncion et al.,, 2009). For all algorithms we used a step size schedule $\rho^{(t)}=t^{-0.75}$.

We held out a test set of 10,000 documents, and periodically evaluated the average per-word marginal log probability assigned by the model to each test document, using the expected value under the variational distribution as a point estimate of the topics. We estimated marginal log probabilities with a Chib-style estimator (Wallach et al.,, 2009).

Figure 2 summarizes the results for $\alpha=0.1$, which yielded the best results for all algorithms. The method of Mimno et al., (2012) outperforms the online LDA algorithm of Hoffman et al., 2010a , but both methods are very sensitive to hyperparameter selection. SSVI achieves good results regardless of hyperparameter choice. SSVI-A’s performance is very slightly worse than that of SSVI.

We also evaluated the stochastic gradient Riemannian Langevin dynamics (SGRLD) algorithm for LDA, which Patterson and Teh, (2013) found outperformed the Gibbs-within-SVI method of Mimno et al., (2012). We experimented with various hyperparameter settings, including the optimal values reported by Patterson and Teh, (2013). The best per-word marginal log probability achieved by SGRLD was $-6.83$. SGRLD thus outperforms standard SVI regardless of hyperparameters, but only outperforms Gibbs-within-SVI for some hyperparameter settings, and achieves performance comparable to that of SSVI and SSVI-A. This suggests that the poor performance of Gibbs-within-SVI in the experiments of Patterson and Teh, (2013) may be due to their setting $\eta=0.1$.

With mini-batches of 1000 documents the computational costs of SSVI and SSVI-A were comparable to those of Gibbs-within-SVI: SSVI took about 30% longer than SSVI to analyze the same number of documents, and SSVI-A’s speed is very close to SSVI. Using Gibbs sampling in the “E-step” rather than variational inference costs about 50% more than a mean-field E-step. We found that SGLRD was significantly slower per document than either SSVI or SSVI-A, although this may be due to implementation differences; in principle the methods should have comparable costs.

Local optima:

The improved performance of SSVI and SSVI-A over SVI might be because SSVI/SSVI-A optimize an objective function that more closely approximates the KL divergence between $q(\beta)$ and $p(\beta|y)$ than the MF objective does. But it could also be because the MF objective includes undesirable local optima that the SSVI objective does not, and so SVI cannot help getting stuck in these local optima. Our experiments above show that SSVI and SSVI-A consistently find better parameter estimates than MF even with multiple restarts, but this does not rule out local optima as an explanation—even with many restarts it might be very difficult for MF to find a good local optimum.

To test this question, we intialized mean-field SVI with the variational parameters found by SVI, Gibbs-within-SVI, SSVI, and SSVI-A, and did another sweep through the Wikipedia dataset. Figure 3 plots the mean-field ELBOs obtained by running SVI from each initialization. SVI initialized with the result from the structured algorithms finds a much better local optimum of the ELBO on the training set than SVI initialized randomly.

This result suggests that the main weakness of mean-field methods may not be the inability of factorized distributions to adequately approximate the posterior, but the difficulty of finding a good local optimum of the ELBO (at least in high-dimensional, multimodal problems). Conversely, the improved performance of the structured methods may be due to the relative lack of nasty local optima in the structured ELBOs. That the structured methods find variational distributions that are also good in the mean-field setting suggests that the structured ELBO may resemble a smoother version of the mean-field ELBO.

2 Dirichlet process mixture of Bernoullis

SSVI-A’s performance closely mirrored that of CGS; both methods were more accurate than mean-field. Mean-field only discovered 17 of the 56 mixture components; the rest were not significantly associated with data. By contrast, SSVI-A and CGS used 54 and 55 components respectively. We also estimated (using Monte Carlo) the KL divergence between the true data-generating distribution $p(y|\pi,\phi)$ and $p(y|\hat{\pi},\hat{\phi})$, where $\hat{\pi}$ and $\hat{\phi}$ are the estimates of the posterior means of $\pi$ and $\phi$ obtained by mean-field, SSVI-A, and CGS. Mean-field achieved a KL divergence of 5.23, while CGS and SSVI-A achieved much lower KL divergences of 1.9 and 1.94, respectively.

3 Bayesian nonparametric nonnegative matrix factorization

We also evaluated the ability of SSVI-A to determine an appropriate number of active components in a Bayesian nonparametric model of audio magnitude spectrograms proposed by Nakano et al., (2011) as a variant on the GaP-NMF model of Hoffman et al., 2010b . The model, which we will call GaP-KL-NMF, assumes that a quantized magnitude spectrogram matrix $Y$ is sampled according to the following generative process:

Discussion

We have presented stochastic structured variational inference (SSVI), an algorithmic framework that uses stochastic variational inference to restore the dependencies between global and local unobserved variables that mean-field variational inference breaks. Experiments suggest that both SSVI and SSVI-A can fit models in a way that outperforms previously existing variational inference algorithms.

Appendix A Algorithm derivation

The convergence proofs that our stochastic variational inference algorithm relies on require that the ELBO be twice continuously differentiable in $\lambda$ (Hoffman et al.,, 2013). We also require that $q(\beta)$ and $p(y,z,\beta)$ be continuously differentiable in $\beta$ for any $y$ and $z$. Finally, we require that $q(z_{n}|\gamma_{n})$ be continuously differentiable in $\gamma_{n}$, that $\gamma_{n}(\beta)$ be continuously differentiable in $\beta$, and that it be possible to compute an unbiased estimate $\hat{\eta}_{n}$ of the expectation $\int_{z_{n}}q(z_{n}|\beta)\eta_{n}(y_{n},z_{n})dz_{n}$ for any $\beta$ (for example, using Markov chain Monte Carlo to sample from $q(z_{n}|\beta)$).

We begin by writing the bound in equation 2.3 as a function of $\lambda$. We define $\mathcal{L}(\lambda)$ as the ELBO achieved by setting $q(z_{n}|\beta)=q(z_{n}|\gamma_{n}(\beta))$ for all values of $\beta$:

where we define $\mathcal{L}_{n}(\beta,\gamma_{n})$ as a shorthand for the part of the ELBO that depends on $y_{n}$ and $z_{n}$.

We will optimize this bound using stochastic optimization. To do so, we need to consider the derivative of $\mathcal{L}(\lambda)$ with respect to $\lambda$. The derivative of the first term in the expectation simplifies:

Now, taking the derivative of $\mathcal{L}_{n}$ with respect to $\lambda$ using the chain rule yields

The second term vanishes because $\nabla_{\gamma_{n}}\mathcal{L}_{n}(\beta,\gamma_{n}(\beta))=0$ by the definition in equation 6, and so the derivative of equation 17 simplifies to

Now, if we sample the global parameters $\beta\sim q_{\beta}$ and compute an unbiased estimate $\hat{\eta}_{n}$ of the expectation $\int_{z_{n}}q(z_{n}|\gamma_{n}(\beta))\eta_{n}(y_{n},z_{n})dz_{n}$, then we can compute a random vector $g$ whose expectation is the true gradient:

A stochastic natural gradient can be obtained by preconditioning this stochastic gradient with the inverse of the Fisher information matrix of $q(\beta)$:

where $V(\beta,\lambda)$ is defined as in equation 8 as

Appendix B Derivatives of Quantile Functions

One definition of the quantile function is as the inverse of the cumulative distribution function (CDF) $Q_{k}(\beta_{k},\lambda)=q(\beta_{k}<\beta_{k})$. Writing down the definition of an inverse function and Differentiating both sides of this definition shows that

where we use the identities that the derivative of a function’s inverse is one over the derivative of that function and that the derivative of a CDF with respect to the random variable is the corresponding probability distribution function (PDF). The derivative of $Q_{k}$ with respect to $\lambda_{k}$ can be obtained numerically using finite differences or automatic differentiation. (The same is true of $R_{k}$, but CDFs are often much cheaper to compute than quantile functions.) For multivariate distributions defined as in equation 7 (main text) we can compute $\frac{\partial{R}}{\partial{\lambda}}$ as

where $\hat{q}_{k}$ is the PDF of the $k$th random variable obtained via the $k$th univariate quantile function $\hat{R}_{k}$ and $\hat{Q}_{k}$ is the CDF that is the inverse of $\hat{R}_{k}$.

Appendix C SSVI for Latent Dirichlet Allocation

In this section we demonstrate how to use SSVI to do approximate posterior inference on the popular topic model latent Dirichlet allocation (LDA) (Blei et al.,, 2003). LDA is a generative model of text that assumes that the words in a corpus of documents are generated according to the process

where $w_{n,m}\in\{1,\ldots,V\}$ is the index into the vocabulary of the $m$th word in the $n$th document, $z_{n,m}\in\{1,\ldots,K\}$ indicates which topic is responsible for $w_{n,m}$, $\theta_{n,k}$ is the prior probability of a word in document $n$ coming from topic $k$, and $\beta_{k,v}$ is the probability of drawing the word index $v$ from topic $k$. For simplicity we use symmetric Dirichlet priors.

LDA fits into the SSVI framework; the random variables $\beta$, $\theta$, $z$, and $w$ can be broken into global variables ($\beta$) and $N$ sets of local variables ($\theta_{n}$, $z_{n}$, and $w_{n}$) that are conditionally independent given the global variables, and the posterior over $\beta$ given $w$, $z$, and $\theta$ is in the same tractable exponential family as the prior (i.e., a Dirichlet):

So we could sample from $q_{\beta}$ by sampling $KV$ independent uniform random variables $u_{k,v}$, passing each through the gamma quantile function $\hat{R}$ to get $\beta_{k,v}^{\prime}\equiv\hat{R}(u_{k,v},\lambda_{k,v},1)$, and letting $\beta_{k,v}=\beta^{\prime}_{k,v}/\sum_{i}\beta^{\prime}_{k,i}$ so that we have $R(u_{k,v},\lambda_{k,v})\equiv\hat{R}(u_{k,v},\lambda_{k,v},1)/\sum_{i}\hat{R}(u_{k,i},\lambda_{k,i},1).$

To compute the update for $\lambda$ in algorithm 1 we need to know $(\frac{\partial^{2}A}{\partial\lambda\partial\lambda^{\top}})^{-1}(\frac{\partial{t}}{\partial{\beta}}|_{\beta^{(t)}}\frac{\partial{R}}{\partial{\lambda}}|_{u^{(t)},\lambda^{(t)}})^{\top}\eta_{n}(w_{n^{(t)}},z^{(t)})$. Since each $\beta_{k}$ is independent of all of the other topic vectors under $q$, we need only consider a single $\beta_{k}$ at a time. The sufficient statistic vector for the Dirichlet distribution is $t(\beta_{k})=\log\beta_{k}$, so we have

We now have everything we need to apply algorithm 1 to LDA.

Appendix D Full Matrix of LDA Results

We tested various combinations of E-steps and M-steps for latent Dirichlet allocation with 100 topics on the 3,800,000-document Wikipedia dataset from (Hoffman et al.,, 2013). To update the global variational distributions, we used traditional mean-field updates, SSVI updates, and SSVI-A updates. For the local variational distributions, we used the traditional mean-field approximation (Blei et al.,, 2003), the CVB0 algorithm of Asuncion et al., (2009), and Gibbs sampling as in (Mimno et al.,, 2012). We also experimented with various settings of the hyperparameters $\alpha$ and $\eta$, which mean-field variational inference for LDA is known to be quite sensitive to (Asuncion et al.,, 2009). For all algorithms we used mini-batches of 1000 documents and a step size schedule $\rho^{(t)}=t^{-0.75}$.

Figure 6 summarizes the results for $\alpha=0.1$, which yielded the best results for all variational algorithms. Using traditional mean-field inference (bottom row) to approximate $p(z_{n}|y_{n},\beta)$ degrades performance, but the CVB0 approximation (top row) works almost as well as Gibbs sampling (middle row) for the two SSVI algorithms. CVB0 is outperformed by Gibbs when using the mean-field M-step. The two SSVI algorithms perform comparably well, but the mean-field M-step (left column) is very sensitive to hyperparameter selection compared to SSVI and SSVI-A.

References