Variational Inference: A Review for Statisticians

Introduction

One of the core problems of modern statistics is to approximate difficult-to-compute probability densities. This problem is especially important in Bayesian statistics, which frames all inference about unknown quantities as a calculation about the posterior. Modern Bayesian statistics relies on models for which the posterior is not easy to compute and corresponding algorithms for approximating them.

In this paper, we review variational inference (vi), a method from machine learning for approximating probability densities (Jordan et al., 1999; Wainwright and Jordan, 2008). Variational inference is widely used to approximate posterior densities for Bayesian models, an alternative strategy to Markov chain Monte Carlo (mcmc) sampling. Compared to mcmc, variational inference tends to be faster and easier to scale to large data—it has been applied to problems such as large-scale document analysis, computational neuroscience, and computer vision. But variational inference has been studied less rigorously than mcmc, and its statistical properties are less well understood. In writing this paper, our hope is to catalyze statistical research on variational inference.

First, we set up the general problem. Consider a joint density of latent variables $\mathbf{z}=z_{1:m}$ and observations $\mathbf{x}=x_{1:n}$,

In Bayesian models, the latent variables help govern the distribution of the data. A Bayesian model draws the latent variables from a prior density $p(\mathbf{z})$ and then relates them to the observations through the likelihood $p(\mathbf{x}\,|\,\mathbf{z})$. Inference in a Bayesian model amounts to conditioning on data and computing the posterior $p(\mathbf{z}\,|\,\mathbf{x})$. In complex Bayesian models, this computation often requires approximate inference.

For decades, the dominant paradigm for approximate inference has been mcmc (Hastings, 1970; Gelfand and Smith, 1990). In mcmc, we first construct an ergodic Markov chain on $\mathbf{z}$ whose stationary distribution is the posterior $p(\mathbf{z}\,|\,\mathbf{x})$. Then, we sample from the chain to collect samples from the stationary distribution. Finally, we approximate the posterior with an empirical estimate constructed from (a subset of) the collected samples.

mcmc sampling has evolved into an indispensable tool to the modern Bayesian statistician. Landmark developments include the Metropolis-Hastings algorithm (Metropolis et al., 1953; Hastings, 1970), the Gibbs sampler (Geman and Geman, 1984) and its application to Bayesian statistics (Gelfand and Smith, 1990). mcmc algorithms are under active investigation. They have been widely studied, extended, and applied; see Robert and Casella (2004) for a perspective.

However, there are problems for which we cannot easily use this approach. These arise particularly when we need an approximate conditional faster than a simple mcmc algorithm can produce, such as when data sets are large or models are very complex. In these settings, variational inference provides a good alternative approach to approximate Bayesian inference.

Rather than use sampling, the main idea behind variational inference is to use optimization. First, we posit a family of approximate densities $\mathcal{Q}$. This is a set of densities over the latent variables. Then, we try to find the member of that family that minimizes the Kullback-Leibler (kl) divergence to the exact posterior,

Finally, we approximate the posterior with the optimized member of the family $q^{*}(\cdot)$.

We emphasize that mcmc and variational inference are different approaches to solving the same problem. mcmc algorithms sample a Markov chain; variational algorithms solve an optimization problem. mcmc algorithms approximate the posterior with samples from the chain; variational algorithms approximate the posterior with the result of the optimization.

Comparing variational inference and mcmc. When should a statistician use mcmc and when should she use variational inference? We will offer some guidance. mcmc methods tend to be more computationally intensive than variational inference but they also provide guarantees of producing (asymptotically) exact samples from the target density (Robert and Casella, 2004). Variational inference does not enjoy such guarantees—it can only find a density close to the target—but tends to be faster than mcmc. Because it rests on optimization, variational inference easily takes advantage of methods like stochastic optimization (Robbins and Monro, 1951; Kushner and Yin, 1997) and distributed optimization (though some mcmc methods can also exploit these innovations (Welling and Teh, 2011; Ahmed et al., 2012)).

Thus, variational inference is suited to large data sets and scenarios where we want to quickly explore many models; mcmc is suited to smaller data sets and scenarios where we happily pay a heavier computational cost for more precise samples. For example, we might use mcmc in a setting where we spent 20 years collecting a small but expensive data set, where we are confident that our model is appropriate, and where we require precise inferences. We might use variational inference when fitting a probabilistic model of text to one billion text documents and where the inferences will be used to serve search results to a large population of users. In this scenario, we can use distributed computation and stochastic optimization to scale and speed up inference, and we can easily explore many different models of the data.

Data set size is not the only consideration. Another factor is the geometry of the posterior distribution. For example, the posterior of a mixture model admits multiple modes, each corresponding label permutations of the components. Gibbs sampling, if the model permits, is a powerful approach to sampling from such target distributions; it quickly focuses on one of the modes. For mixture models where Gibbs sampling is not an option, variational inference may perform better than a more general mcmc technique (e.g., Hamiltonian Monte Carlo), even for small datasets (Kucukelbir et al., 2015). Exploring the interplay between model complexity and inference (and between variational inference and mcmc) is an exciting avenue for future research (see Section 5.4).

The relative accuracy of variational inference and mcmc is still unknown. We do know that variational inference generally underestimates the variance of the posterior density; this is a consequence of its objective function. But, depending on the task at hand, underestimating the variance may be acceptable. Several lines of empirical research have shown that variational inference does not necessarily suffer in accuracy, e.g., in terms of posterior predictive densities (Blei and Jordan, 2006; Braun and McAuliffe, 2010; Kucukelbir et al., 2016); other research focuses on where variational inference falls short, especially around the posterior variance, and tries to more closely match the inferences made by mcmc (Giordano et al., 2015). In general, a statistical theory and understanding around variational inference is an important open area of research (see Section 5.2). We can envision future results that outline which classes of models are particularly suited to each algorithm and perhaps even theory that bounds their accuracy. More broadly, variational inference is a valuable tool, alongside mcmc, in the statistician’s toolbox.

It might appear to the reader that variational inference is only relevant to Bayesian analysis. Indeed, both variational inference and mcmc have had a significant impact on applied Bayesian computation and we will be focusing on Bayesian models here. We emphasize, however, that these techniques also apply more generally to computation about intractable densities. mcmc is a tool for simulating from densities and variational inference is a tool for approximating densities. One need not be a Bayesian to have use for variational inference.

Research on variational inference. The development of variational techniques for Bayesian inference followed two parallel, yet separate, tracks. Peterson and Anderson (1987) is arguably the first variational procedure for a particular model: a neural network. This paper, along with insights from statistical mechanics (Parisi, 1988), led to a flurry of variational inference procedures for a wide class of models (Saul et al., 1996; Jaakkola and Jordan, 1996; Jaakkola and Jordan, 1997; Ghahramani and Jordan, 1997; Jordan et al., 1999). In parallel, Hinton and Van Camp (1993) proposed a variational algorithm for a similar neural network model. Neal and Hinton (1999) (first published in 1993) made important connections to the expectation maximization (em) algorithm (Dempster et al., 1977), which then led to a variety of variational inference algorithms for other types of models (Waterhouse et al., 1996; MacKay, 1997).

Modern research on variational inference focuses on several aspects: tackling Bayesian inference problems that involve massive data; using improved optimization methods for solving Equation (1) (which is usually subject to local minima); developing generic variational inference, algorithms that are easy to apply to a wide class of models; and increasing the accuracy of variational inference, e.g., by stretching the boundaries of $\mathcal{Q}$ while managing complexity in optimization.

Organization of this paper. Section 2 describes the basic ideas behind the simplest approach to variational inference: mean-field inference and coordinate-ascent optimization. Section 3 works out the details for a Bayesian mixture of Gaussians, an example model familiar to many readers. Sections 4.1 and 4.2 describe variational inference for the class of models where the joint density of the latent and observed variables are in the exponential family—this includes many intractable models from modern Bayesian statistics and reveals deep connections between variational inference and the Gibbs sampler of Gelfand and Smith (1990). Section 4.3 expands on this algorithm to describe stochastic variational inference (Hoffman et al., 2013), which scales variational inference to massive data using stochastic optimization (Robbins and Monro, 1951). Finally, with these foundations in place, Section 5 gives a perspective on the field—applications in the research literature, a survey of theoretical results, and an overview of some open problems.

Variational inference

The goal of variational inference is to approximate a conditional density of latent variables given observed variables. The key idea is to solve this problem with optimization. We use a family of densities over the latent variables, parameterized by free “variational parameters.” The optimization finds the member of this family, i.e., the setting of the parameters, that is closest in kl divergence to the conditional of interest. The fitted variational density then serves as a proxy for the exact conditional density. (All vectors defined below are column vectors, unless stated otherwise.)

Let $\mathbf{x}=x_{1:n}$ be a set of observed variables and $\mathbf{z}=z_{1:m}$ be a set of latent variables, with joint density $p(\mathbf{z},\mathbf{x})$. We omit constants, such as hyperparameters, from the notation.

The inference problem is to compute the conditional density of the latent variables given the observations, $p(\mathbf{z}\,|\,\mathbf{x})$. This conditional can be used to produce point or interval estimates of the latent variables, form predictive densities of new data, and more.

The denominator contains the marginal density of the observations, also called the evidence. We calculate it by marginalizing out the latent variables from the joint density,

For many models, this evidence integral is unavailable in closed form or requires exponential time to compute. The evidence is what we need to compute the conditional from the joint; this is why inference in such models is hard.

Note we assume that all unknown quantities of interest are represented as latent random variables. This includes parameters that might govern all the data, as found in Bayesian models, and latent variables that are “local” to individual data points.

Bayesian mixture of Gaussians. Consider a Bayesian mixture of unit-variance univariate Gaussians. There are $K$ mixture components, corresponding to $K$ Gaussian distributions with means $\boldsymbol{\mathbf{\mu}}=\{\mu_{1},\ldots,\mu_{K}\}$. The mean parameters are drawn independently from a common prior $p(\mu_{k})$, which we assume to be a Gaussian $\mathcal{N}(0,\sigma^{2})$; the prior variance $\sigma^{2}$ is a hyperparameter. To generate an observation $x_{i}$ from the model, we first choose a cluster assignment $c_{i}$. It indicates which latent cluster $x_{i}$ comes from and is drawn from a categorical distribution over $\{1,\dots,K\}$. (We encode $c_{i}$ as an indicator $K$-vector, all zeros except for a one in the position corresponding to $x_{i}$’s cluster.) We then draw $x_{i}$ from the corresponding Gaussian $\mathcal{N}(c_{i}^{\top}\boldsymbol{\mathbf{\mu}},1)$.

For a sample of size $n$, the joint density of latent and observed variables is

The latent variables are $\mathbf{z}=\{\boldsymbol{\mathbf{\mu}},\mathbf{c}\}$, the $K$ class means and $n$ class assignments.

The integrand in Equation (8) does not contain a separate factor for each $\mu_{k}$. (Indeed, each $\mu_{k}$ appears in all $n$ factors of the integrand.) Thus, the integral in Equation (8) does not reduce to a product of one-dimensional integrals over the $\mu_{k}$’s. The time complexity of numerically evaluating the $K$-dimensional integral is $\mathcal{O}(K^{n})$.

If we distribute the product over the sum in (8) and rearrange, we can write the evidence as a sum over all possible configurations $\mathbf{c}$ of cluster assignments,

Here each individual integral is computable, thanks to the conjugacy between the Gaussian prior on the components and the Gaussian likelihood. But there are $K^{n}$ of them, one for each configuration of the cluster assignments. Computing the evidence remains exponential in $K$, hence intractable.

2 The evidence lower bound

In variational inference, we specify a family $\mathcal{Q}$ of densities over the latent variables. Each $q(\mathbf{z})\in\mathcal{Q}$ is a candidate approximation to the exact conditional. Our goal is to find the best candidate, the one closest in kl divergence to the exact conditional. The kl divergence is an information-theoretic measure of proximity between two densities. It is asymmetric—that is, $\textsc{kl}\left(q\|p\right)\neq\textsc{kl}\left(p\|q\right)$—and nonnegative. It is minimized when $q(\cdot)=p(\cdot)$. Inference now amounts to solving the following optimization problem,

Once found, $q^{*}(\cdot)$ is the best approximation of the conditional, within the family $\mathcal{Q}$. The complexity of the family determines the complexity of this optimization.

However, this objective is not computable because it requires computing the evidence $\log p(\mathbf{x})$ in Equation (3). (That the evidence is hard to compute is why we appeal to approximate inference in the first place.) To see why, recall that kl divergence is where all expectations are taken with respect to $q(\mathbf{z})$. Expand the conditional,

This reveals its dependence on $\log p(\mathbf{x})$.

Because we cannot compute the kl, we optimize an alternative objective that is equivalent to the kl up to an added constant,

This function is called the evidence lower bound (elbo). The elbo is the negative kl divergence of Equation (11) plus $\log p(\mathbf{x})$, which is a constant with respect to $q(\mathbf{z})$. Maximizing the elbo is equivalent to minimizing the kl divergence.

Examining the elbo gives intuitions about the optimal variational density. We rewrite the elbo as a sum of the expected log likelihood of the data and the kl divergence between the prior $p(\mathbf{z})$ and $q(\mathbf{z})$,

Which values of $\mathbf{z}$ will this objective encourage $q(\mathbf{z})$ to place its mass on? The first term is an expected likelihood; it encourages densities that place their mass on configurations of the latent variables that explain the observed data. The second term is the negative divergence between the variational density and the prior; it encourages densities close to the prior. Thus the variational objective mirrors the usual balance between likelihood and prior.

Another property of the elbo is that it lower-bounds the (log) evidence, $\log p(\mathbf{x})\geq\textsc{elbo}(q)$ for any $q(\mathbf{z})$. This explains the name. To see this notice that Equations (11) and (12) give the following expression of the evidence,

The bound then follows from the fact that $\textsc{kl}\left(\cdot\right)\geq 0$ (Kullback and Leibler, 1951). In the original literature on variational inference, this was derived through Jensen’s inequality (Jordan et al., 1999).

The relationship between the elbo and $\log p(\mathbf{x})$ has led to using the variational bound as a model selection criterion. This has been explored for mixture models (Ueda and Ghahramani, 2002; McGrory and Titterington, 2007) and more generally (Beal and Ghahramani, 2003). The premise is that the bound is a good approximation of the marginal likelihood, which provides a basis for selecting a model. Though this sometimes works in practice, selecting based on a bound is not justified in theory. Other research has used variational approximations in the log predictive density to use vi in cross-validation based model selection (Nott et al., 2012).

Finally, many readers will notice that the first term of the elbo in Equation (12) is the expected complete log-likelihood, which is optimized by the em algorithm (Dempster et al., 1977). The em algorithm was designed for finding maximum likelihood estimates in models with latent variables. It uses the fact that the elbo is equal to the log likelihood $\log p(\mathbf{x})$ (i.e., the log evidence) when $q(\mathbf{z})=p(\mathbf{z}\,|\,\mathbf{x})$. em alternates between computing the expected complete log likelihood according to $p(\mathbf{z}\,|\,\mathbf{x})$ (the E step) and optimizing it with respect to the model parameters (the M step). Unlike variational inference, em assumes the expectation under $p(\mathbf{z}\,|\,\mathbf{x})$ is computable and uses it in otherwise difficult parameter estimation problems. Unlike em, variational inference does not estimate fixed model parameters—it is often used in a Bayesian setting where classical parameters are treated as latent variables. Variational inference applies to models where we cannot compute the exact conditional of the latent variables.Two notes: (a) Variational em is the em algorithm with a variational E-step, i.e., a computation of an approximate conditional. (b) The coordinate ascent algorithm of Section 2.4 can look like the em algorithm. The “E step” computes approximate conditionals of local latent variables; the “M step” computes a conditional of the global latent variables.

3 The mean-field variational family

We described the elbo, the variational objective function in the optimization of Equation (10). We now describe a variational family $\mathcal{Q}$, to complete the specification of the optimization problem. The complexity of the family determines the complexity of the optimization; it is more difficult to optimize over a complex family than a simple family.

In this review we focus on the mean-field variational family, where the latent variables are mutually independent and each governed by a distinct factor in the variational density. A generic member of the mean-field variational family is

Each latent variable $z_{j}$ is governed by its own variational factor, the density $q_{j}(z_{j})$. In optimization, these variational factors are chosen to maximize the elbo of Equation (12).

We emphasize that the variational family is not a model of the observed data—indeed, the data $\mathbf{x}$ does not appear in Equation (14). Instead, it is the elbo, and the corresponding kl minimization problem, that connects the fitted variational density to the data and model.

Notice we have not specified the parametric form of the individual variational factors. In principle, each can take on any parametric form appropriate to the corresponding random variable. For example, a continuous variable might have a Gaussian factor; a categorical variable will typically have a categorical factor. We will see in Sections 4, 4.1 and 4.2 that there are many models for which properties of the model determine optimal forms of the mean-field variational factors $q_{j}(z_{j})$.

Finally, though we focus on mean-field inference in this review, researchers have also studied more complex families. One way to expand the family is to add dependencies between the variables (Saul and Jordan, 1996; Barber and Wiegerinck, 1999); this is called structured variational inference. Another way to expand the family is to consider mixtures of variational densities, i.e., additional latent variables within the variational family (Bishop et al., 1998). Both of these methods potentially improve the fidelity of the approximation, but there is a trade off. Structured and mixture-based variational families come with a more difficult-to-solve variational optimization problem.

Bayesian mixture of Gaussians (continued). Consider again the Bayesian mixture of Gaussians. The mean-field variational family contains approximate posterior densities of the form

Following the mean-field recipe, each latent variable is governed by its own variational factor. The factor $q(\mu_{k};m_{k},s^{2}_{k})$ is a Gaussian distribution on the $k$th mixture component’s mean parameter; its mean is $m_{k}$ and its variance is $s^{2}_{k}$. The factor $q(c_{i};\varphi_{i})$ is a distribution on the $i$th observation’s mixture assignment; its assignment probabilities are a $K$-vector $\varphi_{i}$.

Here we have asserted parametric forms for these factors: the mixture components are Gaussian with variational parameters (mean and variance) specific to the $k$th cluster; the cluster assignments are categorical with variational parameters (cluster probabilities) specific to the $i$th data point. In fact, these are the optimal forms of the mean-field variational density for the mixture of Gaussians.

With the variational family in place, we have completely specified the variational inference problem for the mixture of Gaussians. The elbo is defined by the model definition in Equation (7) and the mean-field family in Equation (15). The corresponding variational optimization problem maximizes the elbo with respect to the variational parameters, i.e., the Gaussian parameters for each mixture component and the categorical parameters for each cluster assignment. We will see this example through in Section 3.

Visualizing the mean-field approximation. The mean-field family is expressive because it can capture any marginal density of the latent variables. However, it cannot capture correlation between them. Seeing this in action reveals some of the intuitions and limitations of mean-field variational inference.

Consider a two dimensional Gaussian distribution, shown in violet in Figure 1. This density is highly correlated, which defines its elongated shape.

The optimal mean-field variational approximation to this posterior is a product of two Gaussian distributions. Figure 1 shows the mean-field variational density after maximizing the elbo. While the variational approximation has the same mean as the original density, its covariance structure is, by construction, decoupled.

Further, the marginal variances of the approximation under-represent those of the target density. This is a common effect in mean-field variational inference and, with this example, we can see why. The kl divergence from the approximation to the posterior is in Equation (LABEL:eq:kl-simple). It penalizes placing mass in $q(\cdot)$ on areas where $p(\cdot)$ has little mass, but penalizes less the reverse. In this example, in order to successfully match the marginal variances, the circular $q(\cdot)$ would have to expand into territory where $p(\cdot)$ has little mass.

4 Coordinate ascent mean-field variational inference

Using the elbo and the mean-field family, we have cast approximate conditional inference as an optimization problem. In this section, we describe one of the most commonly used algorithms for solving this optimization problem, coordinate ascent variational inference (cavi) (Bishop, 2006). cavi iteratively optimizes each factor of the mean-field variational density, while holding the others fixed. It climbs the elbo to a local optimum.

Because of the mean-field family assumption—that all the latent variables are independent—the expectations on the right hand side do not involve the $j$th variational factor. Thus this is a valid coordinate update.

cavi can also be seen as a “message passing” algorithm (Winn and Bishop, 2005), iteratively updating each random variable’s variational parameters based on the variational parameters of the variables in its Markov blanket. This perspective enabled the design of automated software for a large class of models (Wand et al., 2011; Minka et al., 2014). Variational message passing connects variational inference to the classical theories of graphical models and probabilistic inference (Pearl, 1988; Lauritzen and Spiegelhalter, 1988). It has been extended to nonconjugate models (Knowles and Minka, 2011) and generalized via factor graphs (Minka, 2005).

Finally, cavi is closely related to Gibbs sampling (Geman and Geman, 1984; Gelfand and Smith, 1990), the classical workhorse of approximate inference. The Gibbs sampler maintains a realization of the latent variables and iteratively samples from each variable’s complete conditional. Equation (17) uses the same complete conditional. It takes the expected log, and uses this quantity to iteratively set each variable’s variational factor.Many readers will know that we can significantly speed up the Gibbs sampler by marginalizing out some of the latent variables; this is called collapsed Gibbs sampling. We can speed up variational inference with similar reasoning; this is called collapsed variational inference. It has been developed for the same class of models described here (Sung et al., 2008; Hensman et al., 2012). These ideas are outside the scope of our review.

Derivation. We now derive the coordinate update in Equation (17). The idea appears in Bishop (2006), but the argument there uses gradients, which we do not. Rewrite the elbo of Equation (12) as a function of the $j$th variational factor $q_{j}(z_{j})$, absorbing into a constant the terms that do not depend on it,

We have rewritten the first term of the elbo using iterated expectation. The second term we have decomposed, using the independence of the variables (i.e., the mean-field assumption) and retaining only the term that depends on $q_{j}(z_{j})$.

Up to an added constant, the objective function in Equation (18) is equal to the negative kl divergence between $q_{j}(z_{j})$ and $q^{*}_{j}(z_{j})$ from Equation (17). Thus we maximize the elbo with respect to $q_{j}$ when we set $q_{j}(z_{j})=q^{*}_{j}(z_{j})$.

5 Practicalities

Here, we highlight a few things to keep in mind when implementing and using variational inference in practice.

This is not always a disadvantage. Some models, such as the mixture of Gaussians (Section 3 and appendix B) and mixed-membership model (Appendix C), exhibit many posterior modes due to label switching: swapping cluster assignment labels induces many symmetric posterior modes. Representing one of these modes is sufficient for exploring latent clusters or predicting new observations.

Assessing convergence. Monitoring the elbo in cavi is simple; we typically assess convergence once the change in elbo has fallen below some small threshold. However, computing the elbo of the full dataset may be undesirable. Instead, we suggest computing the average log predictive of a small held-out dataset. Monitoring changes here is a proxy to monitoring the elbo of the full data. (Unlike the full elbo, held-out predictive probability is not guaranteed to monotonically increase across iterations of cavi.)

Numerical stability. Probabilities are constrained to live within $$. Precisely manipulating and performing arithmetic of small numbers requires additional care. When possible, we recommend working with logarithms of probabilities. One useful identity is the “log-sum-exp” trick,

The constant $\alpha$ is typically set to $\max_{i}x_{i}$. This provides numerical stability to common computations in variational inference procedures.

A complete example: Bayesian mixture of Gaussians

As an example, we return to the simple mixture of Gaussians model of Section 2.1. To review, consider $K$ mixture components and $n$ real-valued data points $x_{1:n}$. The latent variables are $K$ real-valued mean parameters $\boldsymbol{\mathbf{\mu}}=\mu_{1:K}$ and $n$ latent-class assignments $\mathbf{c}=c_{1:n}$. The assignment $c_{i}$ indicates which latent cluster $x_{i}$ comes from. In detail, $c_{i}$ is an indicator $K$-vector, all zeros except for a one in the position corresponding to $x_{i}$’s cluster. There is a fixed hyperparameter $\sigma^{2}$, the variance of the normal prior on the $\mu_{k}$’s. We assume the observation variance is one and take a uniform prior over the mixture components.

The joint density of the latent and observed variables is in Equation (7). The variational family is in Equation (15). Recall that there are two types of variational parameters—categorical parameters $\varphi_{i}$ for approximating the posterior cluster assignment of the $i$th data point and Gaussian parameters $m_{k}$ and $s^{2}_{k}$ for approximating the posterior of the $k$th mixture component.

We combine the joint and the mean-field family to form the elbo for the mixture of Gaussians. It is a function of the variational parameters $\mathbf{m}$, $\mathbf{s}^{2}$, and $\boldsymbol{\mathbf{\varphi}}$,

In each term, we have made explicit the dependence on the variational parameters. Each expectation can be computed in closed form.

The cavi algorithm updates each variational parameter in turn. We first derive the update for the variational cluster assignment factor; we then derive the update for the variational mixture component factor.

We first derive the variational update for the cluster assignment $c_{i}$. Using Equation (17),

The terms in the exponent are the components of the joint density that depend on $c_{i}$. The expectation in the second term is over the mixture components $\boldsymbol{\mathbf{\mu}}$.

The first term of Equation (21) is the log prior of $c_{i}$. It is the same for all possible values of $c_{i}$, $\log p(c_{i})=-\log K$. The second term is the expected log of the $c_{i}$th Gaussian density. Recalling that $c_{i}$ is an indicator vector, we can write

We use this to compute the expected log probability,

Thus the variational update for the $i$th cluster assignment is

Notice it is only a function of the variational parameters for the mixture components.

2 The variational density of the mixture-component means

We turn to the variational density $q(\mu_{k};m_{k},s_{k}^{2})$ of the $k$th mixture component. Again we use Equation (17) and write down the joint density up to a normalizing constant,

This calculation reveals that the coordinate-optimal variational density of $\mu_{k}$ is an exponential family with sufficient statistics $\{\mu_{k},\mu_{k}^{2}\}$ and natural parameters $\{\sum_{i=1}^{n}\varphi_{ik}x_{i},-1/2\sigma^{2}-\textstyle\sum_{i=1}^{n}\varphi_{ik}/2\}$, i.e., a Gaussian. Expressed in terms of the variational mean and variance, the updates for $q(\mu_{k})$ are

These updates relate closely to the complete conditional density of the $k$th component in the mixture model. The complete conditional is a posterior Gaussian given the data assigned to the $k$th component. The variational update is a weighted complete conditional, where each data point is weighted by its variational probability of being assigned to component $k$.

3 CAVI for the mixture of Gaussians

Algorithm 2 presents coordinate-ascent variational inference for the Bayesian mixture of Gaussians. It combines the variational updates in Equation (21) and Equation (33). The algorithm requires computing the elbo of Equation (20). We use the elbo to track the progress of the algorithm and assess when it has converged.

Once we have a fitted variational density, we can use it as we would use the posterior. For example, we can obtain a posterior decomposition of the data. We assign points to their most likely mixture assignment $\hat{c}_{i}=\operatorname*{arg\,max}_{k}\varphi_{ik}$ and estimate cluster means with their variational means $m_{k}$.

We can also use the fitted variational density to approximate the predictive density of new data. This approximate predictive is a mixture of Gaussians,

where $p(x_{\textrm{new}}\,|\,m_{k})$ is a Gaussian with mean $m_{k}$ and unit variance.

4 Empirical study

We present two analyses to demonstrate the mixture of Gaussians algorithm in action. The first is a simulation study; the second is an analysis of a data set of natural images.

Simulation study. Consider two-dimensional real-valued data $\mathbf{x}$. We simulate $K=5$ Gaussians with random means, covariances, and mixture assignments. Figure 3 shows the data; each point is colored according to its true cluster. Figure 3 also illustrates the initial variational density of the mixture components—each is a Gaussian, nearly centered, and with a wide variance; the subpanels plot the variational density of the components as the cavi algorithm progresses.

The progression of the elbo tells a story. We highlight key points where the elbo develops “elbows”, phases of the maximization where the variational approximation changes its shape. These “elbows” arise because the elbo is not a convex function in terms of the variational parameters; cavi iteratively reaches better plateaus.

Finally, we plot the logarithm of the Bayesian predictive density as approximated by the variational density. Here we report the average across held-out data. Note this plot is smoother than the elbo.

Image analysis. We now turn to an experimental study. Consider the task of grouping images according to their color profiles. One approach is to compute the color histogram of the images. Figure 4 shows the red, green, and blue channel histograms of two images from the imageclef data (Villegas et al., 2013). Each histogram is a vector of length 192; concatenating the three color histograms gives a 576-dimensional representation of each image, regardless of its original size in pixel-space.

We use cavi to fit a Gaussian mixture model with thirty clusters to image histograms. We randomly select two sets of ten thousand images from the imageclef collection to serve as training and testing datasets. Figure 5 shows similarly colored images assigned to four randomly chosen clusters. Figure 6 shows the average log predictive accuracy of the testing set as a function of time. We compare cavi to an implementation in Stan (Stan Development Team, 2015), which uses a Hamiltonian Monte Carlo-based sampler (Hoffman and Gelman, 2014). (Details are in Appendix B.) cavi is orders of magnitude faster than this sampling algorithm.This is not a definitive comparison between variational inference and mcmc. Other samplers, such as a collapsed Gibbs sampler, may perform better than Hamiltonian Monte Carlo sampling.

Variational inference with exponential families

We described mean-field variational inference and derived cavi, a general coordinate-ascent algorithm for optimizing the elbo. We demonstrated this approach on a simple mixture of Gaussians, where each coordinate update was available in closed form.

The mixture of Gaussians is one member of the important class of models where each complete conditional is in the exponential family. This includes a number of widely used models, such as Bayesian mixtures of exponential families, factorial mixture models, matrix factorization models, certain hierarchical regression models (e.g., linear regression, probit regression), stochastic blockmodels of networks, hierarchical mixtures of experts, and a variety of mixed-membership models (which we will discuss below).

Working in this family simplifies variational inference: it is easier to derive the corresponding cavi algorithm, and it enables variational inference to scale up to massive data. In Section 4.1, we develop the general case. In Section 4.2, we discuss conditionally conjugate models, i.e., the common Bayesian application where some latent variables are “local” to a data point and others, usually identified with parameters, are “global” to the entire data set. Finally, in Section 4.3, we describe stochastic variational inference (Hoffman et al., 2013), a stochastic optimization algorithm that scales up variational inference in this setting.

Consider the generic model $p(\mathbf{z},\mathbf{x})$ of Section 2.1 and suppose each complete conditional is in the exponential family:

where $z_{j}$ is its own sufficient statistic, $h(\cdot)$ is a base measure, and $a(\cdot)$ is the log normalizer (Brown, 1986). Because this is a conditional density, the parameter $\eta_{j}(\mathbf{z}_{-j},\mathbf{x})$ is a function of the conditioning set.

Consider mean-field variational inference for this class of models, where we fit $q(\mathbf{z})=\prod_{j}q_{j}(z_{j})$. The exponential family assumption simplifies the coordinate update of Equation (16),

This update reveals the parametric form of the optimal variational factors. Each one is in the same exponential family as its corresponding complete conditional. Its parameter has the same dimension and it has the same base measure $h(\cdot)$ and log normalizer $a(\cdot)$.

Having established their parametric forms, let $\nu_{j}$ denote the variational parameter for the $j$th variational factor. When we update each factor, we set its parameter equal to the expected parameter of the complete conditional,

This expression facilitates deriving cavi algorithms for many complex models.

2 Conditional conjugacy and Bayesian models

One important special case of exponential family models are conditionally conjugate models with local and global variables. Models like this come up frequently in Bayesian statistics and statistical machine learning, where the global variables are the “parameters” and the local variables are per-data-point latent variables.

Conditionally conjugate models. Let $\beta$ be a vector of global latent variables, which potentially govern any of the data. Let $\mathbf{z}$ be a vector of local latent variables, whose $i$th component only governs data in the $i$th “context.” The joint density is

The mixture of Gaussians of Section 3 is an example. The global variables are the mixture components; the $i$th local variable is the cluster assignment for data point $x_{i}$.

We will assume that the modeling terms of Equation (40) are chosen to ensure each complete conditional is in the exponential family. In detail, we first assume the joint density of each $(x_{i},z_{i})$ pair, conditional on $\beta$, has an exponential family form,

where $t(\cdot,\cdot)$ is the sufficient statistic.

Next, we take the prior on the global variables to be the corresponding conjugate prior (Diaconis et al., 1979; Bernardo and Smith, 1994),

This prior has natural (hyper)parameter $\alpha=[\alpha_{1},\alpha_{2}]^{\top}$, a column vector, and sufficient statistics that concatenate the global variable and its log normalizer in the density of the local variables.

With the conjugate prior, the complete conditional of the global variables is in the same family. Its natural parameter is

Turn now to the complete conditional of the local variable $z_{i}$. Given $\beta$ and $x_{i}$, the local variable $z_{i}$ is conditionally independent of the other local variables $\mathbf{z}_{-i}$ and other data $\mathbf{x}_{-i}$. This follows from the form of the joint density in Equation (40). Thus

We further assume that this density is in an exponential family,

This is a property of the local likelihood term $p(z_{i},x_{i}\,|\,\beta)$ from Equation (41). For example, in the mixture of Gaussians, the complete conditional of the local variable is a categorical.

Variational inference in conditionally conjugate models. We now describe cavi for this general class of models. Write $q(\beta\,|\,\lambda)$ for the variational posterior approximation on $\beta$; we call $\lambda$ the “global variational parameter”. It indexes the same exponential family density as the prior. Similarly, let the variational posterior $q(z_{i}\,|\,\varphi_{i})$ on each local variable $z_{i}$ be governed by a “local variational parameter” $\varphi_{i}$. It indexes the same exponential family density as the local complete conditional. cavi iterates between updating each local variational parameter and updating the global variational parameter.

This is an application of Equation (39), where we take the expectation of the natural parameter of the complete conditional in Equation (44).

The global variational update applies the same technique. It is

Here we take the expectation of the natural parameter in Equation (43).

cavi optimizes the elbo by iterating between local updates of each local parameter and global updates of the global parameters. To assess convergence we can compute the elbo at each iteration (or at some lag), up to a constant that does not depend on the variational parameters,

This is the elbo in Equation (12) applied to the joint in Equation (40) and the corresponding mean-field variational density; we have omitted terms that do not depend on the variational parameters. The last term is

cavi for the mixture of Gaussians model (Algorithm 2) is an instance of this method. Appendix C presents another example of cavi for latent Dirichlet allocation (lda), a probabilistic topic model.

3 Stochastic variational inference

Modern applications of probability models often require analyzing massive data. However, most posterior inference algorithms do not easily scale. cavi is no exception, particularly in the conditionally conjugate setting of Section 4.2. The reason is that the coordinate ascent structure of the algorithm requires iterating through the entire data set at each iteration. As the data set size grows, each iteration becomes more computationally expensive.

An alternative to coordinate ascent is gradient-based optimization, which climbs the elbo by computing and following its gradient at each iteration. This perspective is the key to scaling up variational inference using stochastic variational inference (svi) (Hoffman et al., 2013), a method that combines natural gradients (Amari, 1998) and stochastic optimization (Robbins and Monro, 1951).

svi focuses on optimizing the global variational parameters $\lambda$ of a conditionally conjugate model. The flow of computation is simple. The algorithm maintains a current estimate of the global variational parameters. It repeatedly (a) subsamples a data point from the full data set; (b) uses the current global parameters to compute the optimal local parameters for the subsampled data point; and (c) adjusts the current global parameters in an appropriate way. svi is detailed in Algorithm 3. We now show why it is a valid algorithm for optimizing the elbo.

The natural gradient of the elbo. In gradient-based optimization, the natural gradient accounts for the geometric structure of probability parameters (Amari, 1982, 1998). Specifically, natural gradients warp the parameter space in a sensible way, so that moving the same distance in different directions amounts to equal change in symmetrized kl divergence. The usual Euclidean gradient does not enjoy this property.

In exponential families, we find the natural gradient with respect to the parameter by premultiplying the usual gradient by the inverse covariance of the sufficient statistic, $a^{\prime\prime}(\lambda)^{-1}$. This is the inverse Riemannian metric and the inverse Fisher information matrix (Amari, 1982).

Conditionally conjugate models enjoy simple natural gradients of the elbo. We focus on gradients with respect to the global parameter $\lambda$. Hoffman et al. (2013) derive the Euclidean gradient of the elbo,

We can use this natural gradient in a gradient-based optimization algorithm. At each iteration, we update the global parameters,

Substituting Equation (51) into the second term reveals a special structure,

Notice this does not require additional types of calculations other than those for coordinate ascent updates. At each iteration, we first compute the coordinate update. We then adjust the current estimate to be a weighted combination of the update and the current variational parameter.

Though easy to compute, using the natural gradient has the same cost as the coordinate update in Equation (47); it requires summing over the entire data set and computing the optimal local variational parameters for each data point. With massive data, this is prohibitively expensive.

Stochastic optimization of the elbo. Stochastic variational inference solves this problem by using the natural gradient in a stochastic optimization algorithm. Stochastic optimization algorithms follow noisy but cheap-to-compute gradients to reach the optimum of an objective function. (In the case of the elbo, stochastic optimization will reach a local optimum.) In their seminal paper, Robbins and Monro (1951) proved results implying that optimization algorithms can successfully use noisy, unbiased gradients, as long as the step size sequence satisfies certain conditions. This idea has blossomed (Spall, 2003; Kushner and Yin, 1997). Stochastic optimization has enabled modern machine learning to scale to massive data (Le Cun and Bottou, 2004).

Our aim is to construct a cheaply computed, noisy, unbiased natural gradient. We expand the natural gradient in Equation (51) using Equation (43)

where $\varphi^{*}_{i}$ indicates that we consider the optimized local variational parameters (at fixed global parameters $\lambda$) in Equation (46). We construct a noisy natural gradient by sampling an index from the data and then rescaling the second term,

Finally, we set the step size sequence. It must follow the conditions of Robbins and Monro (1951),

Many sequences will satisfy these conditions, for example $\epsilon_{t}=t^{-\kappa}$ for $\kappa\in(0.5,1]$. The full svi algorithm is in Algorithm 3.

We emphasize that svi requires no new derivation beyond what is needed for cavi. Any implementation of cavi can be immediately scaled up to a stochastic algorithm.

Probabilistic topic models. We demonstrate svi with a probabilistic topic model. Probabilistic topic models are mixed-membership models of text, used to uncover the latent “topics” that run through a collection of documents. Topic models have become a popular technique for exploratory data analysis of large collections (Blei, 2012).

In detail, each latent topic is a distribution over terms in a vocabulary and each document is a collection of words that comes from a mixture of the topics. The topics are shared across the collection, but each document mixes them with different proportions. (This is the hallmark of a mixed-membership model.) Thus topic modeling casts topic discovery as a posterior inference problem. Posterior estimates of the topics and topic proportions can be used to summarize, visualize, explore, and form predictions about the documents.

One motivation for topic modeling is to get a handle on massive collections of documents. Early inference algorithms were based on coordinate ascent variational inference (Blei et al., 2003) and analyzed collections in the thousands or tens of thousands of documents. (Appendix C presents this algorithm). With svi, topic models scale up to millions of documents; the details of the algorithm are in Hoffman et al. (2013). Figure 7 illustrates topics inferred using the latent Dirichlet allocation model (Blei et al., 2003) from 1.8M articles from the New York Times. This analysis would not have been possible without svi.

Discussion

We described variational inference, a method that uses optimization to make probabilistic computations. The goal is to approximate the conditional density of latent variables $\mathbf{z}$ given observed variables $\mathbf{x}$, $p(\mathbf{z}\,|\,\mathbf{x})$. The idea is to posit a family of densities $\mathcal{Q}$ and then to find the member $q^{*}(\cdot)$ that is closest in kl divergence to the conditional of interest. Minimizing the kl divergence is the optimization problem, and its complexity is governed by the complexity of the approximating family.

We then described the mean-field family, i.e., the family of fully factorized densities of the latent variables. Using this family, variational inference is particularly amenable to coordinate-ascent optimization, which iteratively optimizes each factor. (This approach closely connects to the classical Gibbs sampler (Geman and Geman, 1984; Gelfand and Smith, 1990).) We showed how to use mean-field vi to approximate the posterior density of a Bayesian mixture of Gaussians, discussed the special case of exponential families and conditional conjugacy, and described the extension to stochastic variational inference (Hoffman et al., 2013), which scales mean-field variational inference to massive data.

Researchers in many fields have used variational inference to solve real problems. Here we focus on example applications of mean-field variational inference and structured variational inference based on the kl divergence. This discussion is not exhaustive; our intention is to outline the diversity of applications of variational inference.

Computational biology. vi is widely used in computational biology, where probabilistic models provide important building blocks for analyzing genetic data. For example, vi has been used in genome-wide association studies (Carbonetto and Stephens, 2012; Logsdon et al., 2010), regulatory network analysis (Sanguinetti et al., 2006), motif detection (Xing et al., 2004), phylogenetic hidden Markov models (Jojic et al., 2004), population genetics (Raj et al., 2014), and gene expression analysis (Stegle et al., 2010).

Computer vision and robotics. Since its inception, variational inference has been important to computer vision. Vision researchers frequently analyze large and high-dimensional data sets of images, and fast inference is important to successfully deploy a vision system. Some of the earliest examples included inferring non-linear image manifolds (Bishop and Winn, 2000) and finding layers of images in videos (Jojic and Frey, 2001). As other examples, variational inference is important to probabilistic models of videos (Chan and Vasconcelos, 2009; Wang and Mori, 2009), image denoising (Likas and Galatsanos, 2004), tracking (Vermaak et al., 2003; Yu and Wu, 2005), place recognition and mapping for robotics (Cummins and Newman, 2008; Ramos et al., 2012), and image segmentation with Bayesian nonparametrics (Sudderth and Jordan, 2009). Du et al. (2009) uses variational inference in a probabilistic model to combine the tasks of segmentation, clustering, and annotation.

Computational neuroscience. Modern neuroscience research also requires analyzing very large and high-dimensional data sets, such as high-frequency time series data or high-resolution functional magnetic imaging data. There have been many applications of variational inference to neuroscience, especially for autoregressive processes (Roberts and Penny, 2002; Penny et al., 2003, 2005; Flandin and Penny, 2007; Harrison and Green, 2010). Other applications of variational inference to neuroscience include hierarchical models of multiple subjects (Woolrich et al., 2004), spatial models (Sato et al., 2004; Zumer et al., 2007; Kiebel et al., 2008; Wipf and Nagarajan, 2009; Lashkari et al., 2012; Nathoo et al., 2014), brain-computer interfaces (Sykacek et al., 2004), and factor models (Manning et al., 2014; Gershman et al., 2014). There is a software toolbox that uses variational methods for solving neuroscience and psychology research problems (Daunizeau et al., 2014).

Natural language processing and speech recognition. In natural language processing, variational inference has been used for solving problems such as parsing (Liang et al., 2007, 2009), grammar induction (Kurihara and Sato, 2006; Naseem et al., 2010; Cohen and Smith, 2010), models of streaming text (Yogatama et al., 2014), topic modeling (Blei et al., 2003), and hidden Markov models and part-of-speech tagging (Wang and Blunsom, 2013). In speech recognition, variational inference has been used to fit complex coupled hidden Markov models (Reyes-Gomez et al., 2004) and switching dynamic systems (Deng, 2004).

Other applications. There have been many other applications of variational inference. Fields in which it has been used include marketing (Braun and McAuliffe, 2010), optimal control and reinforcement learning (Van Den Broek et al., 2008; Furmston and Barber, 2010), statistical network analysis (Wiggins and Hofman, 2008; Airoldi et al., 2008), astrophysics (Regier et al., 2015), and the social sciences (Erosheva et al., 2007; Grimmer, 2011). General variational inference algorithms have been developed for a variety of classes of models, including shrinkage models (Armagan et al., 2011; Armagan and Dunson, 2011; Neville et al., 2014), general time-series models (Roberts et al., 2004; Barber and Chiappa, 2006; Archambeau et al., 2007b, a; Johnson and Willsky, 2014; Foti et al., 2014), robust models (Tipping and Lawrence, 2005; Wang and Blei, 2015), and Gaussian process models (Titsias and Lawrence, 2010; Damianou et al., 2011; Hensman et al., 2014).

2 Theory

Though researchers have not developed much theory around variational inference, there are several threads of research about theoretical guarantees of variational approximations. As we mentioned in the introduction, one of our purposes for writing this paper is to catalyze research on the statistical theory around variational inference.

Below, we summarize a variety of results. In general, they are all of the following type: treat vi posterior means as point estimates (or use m-step estimates from variational em) and confirm that they have the usual frequentist asymptotics. (Sometimes the research finds that they do not enjoy the same asymptotics.) Each result revolves around a single model and a single family of variational approximations.

You et al. (2014) study the variational posterior for a classical Bayesian linear model. They put a normal prior on the coefficients and an inverse gamma prior on the response variance. They find that, under standard regularity conditions, the mean-field variational posterior mean of the parameters are consistent in the frequentist sense. Ormerod et al. (2014) build on their earlier work with a spike-and-slab prior on the coefficients and find similar consistency results.

Hall et al. (2011a, b) examine a simple Poisson mixed-effects model, one with a single predictor and a random intercept. They use a Gaussian variational approximation and estimate parameters with variational em. They prove consistency of these estimates at the parametric rate and show asymptotic normality with asymptotically valid standard errors.

Celisse et al. (2012) and Bickel et al. (2013) analyze network data using stochastic blockmodels. They show asymptotic normality of parameter estimates obtained using a mean-field variational approximation. They highlight the computational advantages and theoretical guarantees of the variational approach over maximum likelihood for dense, sparse, and restricted variants of the stochastic blockmodel.

Westling and McCormick (2015) study the consistency of vi through a connection to M-estimation. They focus on a broader class of models (with posterior support in real coordinate space) and analyze an automated vi technique that uses a Gaussian variational approximation (Kucukelbir et al., 2015). They derive an asymptotic covariance matrix estimator of the variational approximation and show its robustness to model misspecification.

Finally, Wang and Titterington (2006) analyze variational approximations to mixtures of Gaussians. Specifically, they consider Bayesian mixtures with conjugate priors, the mean-field variational approximation, and an estimator that is the variational posterior mean. They confirm that cavi converges to a local optimum, that the vi estimator is consistent, and that the vi estimate and maximum likelihood estimate (mle) approach each other at a rate of $\mathcal{O}(\nicefrac{{1}}{{n}})$. In Wang and Titterington (2005), they show that the asymptotic variational posterior covariance matrix is “too small”—it differs from the mle covariance (i.e., the inverse Fisher information) by a positive-definite matrix.

3 Beyond conditional conjugacy

We focused on models where the complete conditional is in the exponential family. Many models, however, do not enjoy this property. A simple example is Bayesian logistic regression,

where $\sigma(\cdot)$ is the logistic function. The posterior density of the coefficients is not in an exponential family and we cannot apply the variational inference methods we discussed above. Specifically, we cannot compute the expectations in the first term of the elbo in Equation (12) or the coordinate update in Equation (17).

Exploring variational methods for such models has been a fruitful area of research. An early example is Jaakkola and Jordan (1997); Jaakkola and Jordan (2000), who developed a variational bound tailored to logistic regression. Blei and Lafferty (2007) later adapted their idea to nonconjugate topic models, and researchers have continued to improve the original bound (Khan et al., 2010; Marlin et al., 2011; Ermis and Bouchard, 2014). In other work, Braun and McAuliffe (2010) derived a variational inference algorithm for the discrete choice model, which also lies outside of the class of conditionally conjugate models. They developed a delta method to approximate the difficult-to-compute expectations. Finally, Wand et al. (2011) use auxiliary variable methods, quadrature, and mixture approximations to handle a variety of likelihood terms that fall outside of the exponential family.

More recently, researchers have generalized nonconjugate inference, seeking recipes that can be used across many models. Wang and Blei (2013) adapted Laplace approximations and the delta method to this end, improving inference in nonconjugate generalized linear models and topic models; this approach is also used by Bugbee et al. (2016) for semi-parametric regression. Knowles and Minka (2011) generalized the Jaakkola and Jordan (1997); Jaakkola and Jordan (2000) bound in a message-passing algorithm and Wand (2014) further simplified and extended their approach. Tan and Nott (2013, 2014) applied these message-passing methods to generalized linear mixed models (and also combined them with svi). Rohde and Wand (2015) unified many of these algorithmic developments and provided practical insights into their numerical implementations.

Finally, there has been a flurry of research on optimizing difficult variational objectives with Monte Carlo (mc) estimates of the gradient. The idea is to write the gradient of the elbo as an expectation, compute mc estimates of it, and then use stochastic optimization with repeated mc gradients. This first appeared independently in several papers (Ji et al., 2010; Nott et al., 2012; Paisley et al., 2012; Wingate and Weber, 2013). The newest approaches avoid any model-specific derivations, and are termed “black box” inference methods. As examples, see Kingma and Welling (2014); Rezende et al. (2014); Ranganath et al. (2014, 2016); Salimans and Knowles (2014); Titsias and Lázaro-Gredilla (2014), and Tran et al. (2016). Kucukelbir et al. (2016) leverage these ideas toward an automatic vi technique that works on any model written in the probabilistic programming system Stan (Stan Development Team, 2015). This is a step towards a derivation-free, easy-to-use vi algorithm.

4 Open problems

There are many open avenues for statistical research in variational inference.

We focused on optimizing $\textsc{kl}\left(q(\mathbf{z})||p(\mathbf{z}\,|\,\mathbf{x})\right)$ as the variational objective function. A promising avenue of research is to develop variational inference methods that optimize other measures, such as $\alpha$-divergence measures. As one example, expectation propagation (Minka, 2001) is inspired by the kl divergence “in the other direction,” between $p(\mathbf{z}\,|\,\mathbf{x})$ and $q(\mathbf{z})$. Other work has developed divergences based on lower bounds that are tighter than the elbo (Barber and de van Laar, 1999; Leisink and Kappen, 2001). While alternative divergences may be difficult to optimize, they may give better approximations (Minka, 2005; Opper and Winther, 2005).

Though it is flexible, the mean-field family makes strong independence assumptions. These assumptions help with scalable optimization, but they limit the expressibility of the variational family. Further, they can exacerbate issues with local optima of the objective and underestimating posterior variances; see Figure 1. A second avenue of research is to develop better approximations while maintaining efficient optimization.

As we mentioned above, structured variational inference has its roots in the early days of the method (Saul and Jordan, 1996; Barber and Wiegerinck, 1999). More recently, Hoffman and Blei (2015) use generic structured variational inference in a stochastic optimization algorithm; Kucukelbir et al. (2016), Challis and Barber (2013), and Tan and Nott (2016) take advantage of Gaussian variational families with non-diagonal covariance; Giordano et al. (2015) post-process the mean-field parameters to correct for underestimating the variance; and Ranganath et al. (2016) embed the mean-field parameters themselves in a hierarchical model to induce variational dependencies between latent variables.

The interface between variational inference and mcmc remains relatively unexplored. Freitas et al. (2001) used fitted variational distributions as a component of a proposal distribution for Metropolis-Hastings. Mimno et al. (2012) and Hoffman and Blei (2015) studied mcmc as a method of approximating coordinate updates, e.g., to include structure in the variational family. Salimans et al. (2015) propose a variational approximation to the mcmc chain; their method enables an explicit trade off between computational accuracy and speed. Understanding how to combine these two strategies for approximate inference is a ripe area for future research. A principled analysis of when to use (and combine) variational inference and mcmc would have both theoretical and practical impact in the field.

Finally, the statistical properties of variational inference are not yet well understood, especially in contrast to the wealth of analysis of mcmc techniques. There has been some progress; see Section 5.2. A final open research problem is to understand variational inference as an estimator and to understand its statistical profile relative to the exact posterior.

Appendix A Bayesian Linear Regression with Automatic RelevanceDetermination

Automatic relevance determination (ard) assigns a separate prior for each regression coefficient; the idea is to automatically shrink irrelevant coefficients during inference (MacKay, 1992; Neal, 2012; Tipping, 2001; Wipf and Nagarajan, 2008). ard works by pairing the prior precision of each regression coefficient with its own latent variable $\alpha_{d}$. The hyper-prior on these relevance variables $\alpha$ encourages small values; this, in turn, selects relevant regression coefficients.

Here we present a conditionally conjugate Bayesian linear regression model with an ard prior, based on Drugowitsch (2013). All Gaussian distributions below follow the precision parameterization.

Define a Gaussian likelihood with precision parameter $\tau$ as

ard then posits the following hierarchical prior

where $\alpha$ is a $D$-dimensional relevance variable

Here $a_{0},b_{0},c_{0},$ and $d_{0}$ are fixed hyper-parameters. The latent variable $\alpha$ determines the relevance of each regression coefficient.

The posterior $p(\beta,\tau,\alpha\mid\mathbf{y}\;;\;\mathbf{x})$ is not available in closed form. A simpler model that does not include $\alpha$ admits a closed form posterior because the normal-gamma distribution is conjugate to a normal likelihood with unknown mean and precision.

We derive a cavi algorithm for this model. First, factorize the variational approximation as

Here we consider $\beta$ and $\tau$ in a single factor.

Begin by applying Equation (17) to identify the optimal form of $q(\beta,\tau)$ as

The optimal variational approximation to the regression coefficients and the precision is thus a normal-gamma with the following parameters:

Next consider the optimal form of the relevance variables $\alpha$. Again, apply Equation (17) to identify the optimal form of $q(\alpha)=\prod_{d=1}^{D}q(\alpha_{d})$ as

The optimal variational approximation to the relevance variable is thus a Gamma with the following parameters:

where $[\cdot]_{d}$ indicates the $d$th diagonal entry of a matrix.

Iteratively updating $a_{*},b_{*},c_{*},d_{*},V_{*}^{-1},$ and $\beta_{*}$ defines cavi for this model. These quantities also define the elbo; Drugowitsch (2013) presents the additional algebra that computes the elbo.

Appendix B Gaussian Mixture Model of Image Histograms

The joint density of the model factorizes as

The likelihood is Gaussian with precision parameterization

The marginal over cluster assignments is a categorical distribution,

The prior over the mixing vector is a Dirichlet distribution with fixed hyperparameters $a_{0}$,

The prior over mean and precision parameters is a normal-gamma distribution with hyperparameters $m_{0}$, $b_{0}$, $\alpha_{0}$, $\beta_{0}$,

We use the following values for the hyperparameters

Bishop (2006, Chapter 10.2) derives a cavi algorithm for this model.

Figure 8 presents Stan code that implements this model. Since Stan does not support discrete latent variables, we marginalize over the assignment variables.

Appendix C Latent Dirichlet Allocation

Probabilistic topic models are mixed-membership models of text, used to uncover the latent “topics” that run through a collection of documents. Topic models have become a popular technique for exploratory data analysis of large collections (Blei, 2012).

Latent Dirichlet allocation (lda) (Blei et al., 2003) is a conditionally conjugate topic model (Section 4.2). It treats documents as containing multiple topics, where a topic is a distribution over words in a vocabulary.

Following the notation of Hoffman et al. (2013), let $K$ be a specific number of topics and $V$ the size of the vocabulary. lda defines the following generative process:

draw a distribution over words $\beta_{k}\sim\textrm{Dir}_{V}(\eta)$.

draw a vector of topic proportions $\theta_{d}\sim\textrm{Dir}_{K}(\alpha)$.

draw a topic assignment $z_{dn}\sim\textrm{Mult}(\theta_{d})$, then

draw a word $w_{dn}\sim\textrm{Mult}(\beta_{z_{dn}})$.

The posterior $p(\beta,\theta,z\mid w)$ is not available in closed form. While the topic assignments $z$ and their proportions $\theta$ enjoy a conjugate relationship, the introduction of the topics $\beta$ renders this posterior analytically intractable.

We derive a cavi algorithm for this model, based on Hoffman et al. (2013). Posit a mean-field variational family

Since lda is a conditionally conjugate model, we can directly identify the family of each factor (Section 4.2).

Begin with the complete conditional of the topic assignment. This is a multinomial,

The variational approximation to the topic assignments is also a multinomial distribution, with parameters $\phi_{dn}$.

Follow with the complete conditional of the topic proportions. This is a $K$-dimensional Dirichlet

The variational approximation to the topic proportions is also a $K$-dimensional Dirichlet with parameters $\gamma_{d}$.

End with the complete conditional of the topics. This is a $V$-dimensional Dirichlet

In words, the $v$th element of the $k$th topic is a Dirichlet with parameter equal to the sum of the fixed scalar $\eta$ and the number of times term $v$ (denoted by $w_{dn}$) was assigned to topic $k$ (denoted by $z_{dn}^{k}$). The variational approximation to the topic proportions is a $V$-dimensional Dirichlet with parameters $\lambda_{k}$.

The cavi updates for the topic assignment and topic proportions require iterating over the following for each word within each document until convergence:

This is a direct application of Equation (46) to the complete conditionals above.

The updates for $\phi$ and $\gamma$ depend on the variational parameters for the topics $\lambda$. The update for the topics, in turn, depends on the variational parameters for the topic proportions. That update is

This update only depends on the variational parameter for the topic assignments $\phi_{dn}$.

Algorithm 4 presents the full cavi algorithm for lda. A similar computation defines the elbo for lda; Hoffman et al. (2013) present the additional algebra for the elbo.

References