Bayesian inference has become a crucial component of machine learning. It allows us to systematically reason about parameter uncertainty. The central object of interest in Bayesian inference is the posterior distribution of model parameters given observations. This review focuses on variational inference (VI): a methodology that makes Bayesian inference computationally efficient and scalable to large data sets.
Bayesian machine learning frequently relies on probabilistic latent variable models, such as Gaussian mixture models, Hidden Markov models, Latent Dirichlet Allocation, stochastic block models, and Bayesian deep learning architectures. Computing the exact Bayesian posterior requires to sum or integrate over all latent variables, which can be in the millions or billions for complex models and large-scale applications. Exact inference is therefore typically intractable in these models, and approximations are needed.
The central idea of VI is to approximate the model posterior by a simpler distribution. To this end, one minimizes the Kullback-Leibler divergence between the posterior and the approximating distribution. This approach circumvents computing intractable normalization constants. It only requires knowledge of the joint distribution of the observations and the latent variables. This methodology along with its recent refinements will be reviewed in this paper.
Within the field of approximate Bayesian inference, VI falls into the class of optimization-based approaches [14, 62]. This class also contains methods such as loopy belief propagation  and expectation propagation (EP) 
. On the contrary, Markov Chain Monte Carlo (MCMC) approaches rely on sampling[61, 151, 22]. By construction, MCMC is often unbiased, and thus converges to the true posterior in the limit, but it can be slow to converge. Optimization-based methods, on the other hand, are often faster but may suffer from oversimplified posterior approximations [14, 205]. In recent years, there has been considerable progress in both fields [15, 7], and in particular on bridging the gap between these methods [1, 90, 169, 154, 113]. In fact, recent progress in scalable VI partly relies on fusing optimization-based and sampling-based methods. While this review focuses on VI, readers interested in EP and MCMC are referred to, e.g.,  and .
The origins of VI date back to the 1980s. Mean field methods, for instance, have their origins in statistical physics, where they played a prominent role in the statistical mechanics of spin glasses [121, 147]
. Early applications of variational methods also include the study of neural networks[144, 149]. The latter work inspired the computer science community of the 1990s to adopt variational methods in the context of probabilistic graphical models [172, 73, 79, 143] .
In recent years, several factors have driven a renewed interest in variational methods. The modern versions of VI differ significantly from earlier formulations. Firstly, the availability of large datasets triggered the interest in scalable
approaches, e.g., based on stochastic gradient descent[18, 67]. Secondly, classical VI is limited to conditionally conjugate exponential family models, a restricted class of models described in [205, 67]. In contrast, black box VI algorithms [74, 79, 154] and probabilistic programs facilitate generic VI, making it applicable to a range of complicated models. Thirdly, this generalization has spurred research on more accurate variational approximations, such as alternative divergence measures [128, 103, 221] and structured variational families . Finally, amortized
inference employs complex functions such as neural networks to predict variational distributions conditioned on data points, rendering VI an important component of modern Bayesian deep learning architectures such as variational autoencoders. In this work, we discuss important papers concerned with each of these four aspects.
While several excellent reviews of VI exist, we believe that our focus on recent developments in scalable, generic, accurate and amortized VI goes beyond those efforts. Both  and  date back to the early s and do not cover the developments of recent years. Similarly,  is an excellent resource, especially regarding structured approximations and the information geometrical aspects of VI. However, it was published prior to the widespread use of stochastic methods in VI. Among recent introductions,  contains many examples, empirical comparisons, and explicit model calculations but focuses less on recent developments while  mainly focuses on scalable MCMC. Our review concentrates on the advances of the last 10 years prior to the publication of this paper. Complementing previous reviews, we skip example calculations to focus on a more exhaustive survey of the recent literature. For readers who are new to the field, we recommend to read Chapter 10 on approximate inference in  as a preparation.
We survey the trends and developments in VI in a self-contained manner. Section 2 covers basic concepts, such as variational distributions and the evidence lower bound. In the succeeding sections, we concentrate on recent advances and identify four main research directions: scalable VI (Section 3), generic VI (Section 4), accurate VI (Section 5), and amortized VI (Section 6). We finalize the review with a discussion (Section 7) and concluding remarks (Section 8).
2 Variational Inference
We begin this review with a brief tutorial on variational inference, presenting the mathematical foundations of this procedure and explaining the basic mean-field approximation.
The generative process is specified by observations , as well as latent variables and a joint distribution . We use bold font to explicitly indicate sets of variables, i.e. , where is the total number of latent variables and , where is the total number of observations in the dataset. The variational distribution is defined over the latent variables and has variational parameters . 111Note that the number variational parameter is not necessary the same as the number of latent variables. Latent variables can be shared among multiple data points, and individual data points can have multiple latent variables.
2.1 Inference as Optimization
The central object of interest in Bayesian statistics is the posterior distribution of latent variables given observations:
For most models, this integral is high dimensional, thus computing the normalization term is intractable.
Instead of computing the posterior normalization, the basic idea of VI is to approximate the posterior with a simpler distribution. This involves a variational distribution , characterized by a set of variational parameters . These parameters are tuned to obtain the best matching. Finally, the optimized variational distribution is taken as a proxy for the posterior. In this way, VI turns Bayesian inference into an optimization problem over variational parameters.
For two distributions and , a divergence measures the difference between the distributions, such that and only for . VI amounts to minimizing a divergence between the variational distribution and the posterior. We show below that this does not require knowing the posterior normalization.
While various divergence measures exist [4, 128, 6, 181], the most commonly used divergence is the Kullback-Leibler (KL) divergence [93, 14], which is also referred to as relative entropy or information gain:
As seen in Eq. 2, the KL divergence is asymmetric; . Depending on the ordering, we obtain two different approximate inference methods. As we show below, VI employs . In contrast and as an aside, expectation propagation (EP)  optimizes
for local moment matching, which is not reviewed in this paper222We refer the readers to the EP roadmap for more information about advancements of EP. https://tminka.github.io/papers/ep/roadmap.html.
2.2 Variational Objective
Classical VI aims at determining a variational distribution that is as close as possible to the posterior , measured in terms of the KL divergence. Minimizing the KL divergence to zero would guarantee that the variational distribution matches the exact posterior. However, in practice this is rarely possible: the variational distribution is usually under-parameterized and thus not sufficiently flexible to capture the full complexity of the true posterior.
As follows, we will show that minimizing the KL divergence is equivalent to maximizing a related quantity, the Evidence Lower BOund (ELBO)
. The ELBO is a lower bound on the log marginal probability of the data and can be derived fromusing Jensen’s inequality as follows:
It can be shown (see Appendix A.1) that the difference between the true log marginal probability of the data and the ELBO is the KL divergence between the variational distribution and the posterior:
Thus, maximizing the ELBO is equivalent to minimizing the KL divergence between and , where and replace and
for the sake of brevity. Since the ELBO is a conservative estimate of this marginal, it is sometimes taken as an estimate of how well the model fits the data. The ELBO can also be used for model selection.
In traditional VI, computing the ELBO amounts to analytically solving the expectations over . This restricts the class of tractable models to the so-called conditionally conjugate exponential family (see Appendix A.2 and ). For an example calculation to derive the ELBO analytically for a mixture of Gaussians, we refer to . Section 4 introduces modern alternatives to computing these expectations.
2.3 Mean Field Variational Inference
There is a tradeoff in choosing expressive enough to approximate the posterior well, and simple enough to lead to a tractable approximation . A common choice is a fully factorized distribution, also called mean field distribution. A mean field approximation assumes that all latent variables are independent, which simplifies derivations. However, this independence assumption also leads to less accurate results especially when the true posterior variables are highly dependent. Section 5 discusses a more expressive class of variational distributions.
Mean Field Variational Inference (MFVI) has its origins in the mean field theory of physics . In this approximation, the variational distribution factorizes, and each factor is governed by its own variational parameter:
For notational simplicity, we omit the variational parameters for the remainder of this section. We now review how to maximize the ELBO , defined in Eq. 3, under a mean field assumption.
A fully factorized variational distribution allows one to optimize via simple iterative updates. To see this, we focus on updating the variational parameter associated with latent variable . Inserting the mean field distribution into Eq. 3 allows us to express the ELBO as follows:
Above, indicates the set excluding . The constant contains all terms that are constant with respect to , such as the entropy term associated with . We have thus separated the full expectation into an inner expectation over , and an outer expectation over .
Eq. 6 assumes the form of a negative KL divergence, which is maximized for variable by
Exponentiating and normalizing this result yields:
Using Eq. 8, the variational distribution can be updated iteratively for each latent variable until convergence. Similar updates also form the basis for the variational message passing algorithm  (Appendix A.3).
2.4 Beyond Vanilla Variational Inference
Classical MFVI has historically played an important role, however, it is limited in multiple ways when it comes to modern applications. One of the challenges is to scale VI to big datasets. This will be addressed in Section 3, where we show that VI can be combined with stochastic optimization and distributed computing to achieve this goal. Big datasets and fast algorithms allow for more sophisticated models. In order to make VI tractable for this modern class of models (in particular for so-called non-conjugate ones), Section 4 describes methods that make VI both easier to use and more generic. Furthermore, certain models and applications require more accurate inference techniques, such as improved variational approximations and tighter bounds. A popular stream of research is concerned with alternative divergence measures beyond the KL divergence, and will be reviewed in Section 5, where we also review non-mean field variational approximations. Finally, we describe in Section 6 how neural networks can be used to amortize the estimiation of certain local latent variables. This leads to a significant speedup for many models and bridges the gap between Bayesian inference and modern representation learning.
3 Scalable Variational Inference
In this section, we survey scalable VI. Big datasets raise new challenges for the computational feasibility of Bayesian algorithms, making scalable inference techniques essential. We begin by reviewing stochastic variational inference (SVI) in Section 3.1, which uses stochastic gradient descent (SGD) to scale VI to large datasets. Section 3.2
discusses practical aspects of SVI, such as adaptive learning rates and variance reduction. Further approaches to improve on the scalability of VI are discussed in Section3.3; these include sparse inference, collapsed inference, and distributed inference.
This section follows the general model structure of global and local hidden variables, assumed in . Fig. 1 depicts the corresponding graphical model where the latent variable includes local (per data point) variables and global variable . Similarly, the variational parameters are given by , where the variational parameter corresponds to the global latent variable, and
denotes the set of local variational parameters. Furthermore, the model depends on hyperparameters. The mini-batch size is denoted by .
3.1 Stochastic Variational Inference
We showed that VI frames Bayesian inference as an optimization problem. For many models of interest, the variational objective has a special structure, namely, it is the sum over contributions from all individual data points. Problems of this type can be solved efficiently using stochastic optimization [161, 18]. Stochastic Variational Inference amounts to applying stochastic optimization to the objective function encountered in VI [71, 67, 65, 209], thereby scaling VI to very large datasets. Using stochastic optimization in the context of VI was proposed in [171, 71, 67]. We follow the conventions of  which presents SVI for models of the conditionally conjugate exponential family class.
The ELBO of the general graphical model shown in Fig. 1 has the following form:
We assume that the variational distribution is given by . Here, we also assume that the expectations in Eq. 9 are analytically tractable, yielding a closed-form objective.
Eq. 9 could be optimized by coordinate descent (Section 2), or gradient descent on the ELBO. In both cases, every iteration or gradient step scales with , and is therefore expensive for large data. In contrast, SVI solves this problem in the spirit of stochastic gradient descent . In every iteration, one randomly selects mini-batches of size to obtain a stochastic estimate of the ELBO ,
where is the variable index from the mini-batch. Then, the gradient of Eq. 10 is computed, which gives a noisy estimator of the direction of steepest ascent of the true ELBO.
An important result of  is that using natural gradients instead of standard gradients in SVI simplifies the variational updates for models in the conditionally conjugate exponential family. Natural gradients, first studied in  and introduced to VI in [171, 70, 69], take the information geometry of the model into account. They are obtained by pre-multiplying the gradient with the inverse Fisher information matrix. While we skip further discussions for brevity, interested readers are referred to Appendix A.4 and . Recent advances in this direction include [64, 106, 118].
SVI requires the same conditions for convergence as regular SGD. The minibatch indices must be drawn uniformly at random. The size of the minibatch satisfies . Larger values of reduce the variance of the stochastic natural gradient. When , SVI reduces to traditional batch VI when the learning rate is set to . However, computational savings are only obtained for . The optimal choice of emerges from a trade-off between the computational overhead associated with processing a mini-batch, such as performing inference over global parameters (favoring larger mini-batches which have lower gradient noise, allowing larger learning rates), and the cost of iterating over local parameters in the mini-batch (favoring small mini-batches). Additionally, this tradeoff is also affected by memory structures in modern hardware such as GPUs. The learning rate needs to decrease with iterations , satisfying the Robbins-Monro conditions and . This guarantees that every point in the parameter space can be reached, while the gradient noise decreases quickly enough to ensure convergence .
Sometimes SVI is referred to as online VI [65, 209]. These methods are equivalent under the assumptions that the volume of the data is known. In streaming applications, the mini-batches arrive sequentially from a data source, but the SVI updates are the same. However, when is unknown, it is unclear how to set the scale parameter in Eq. 10. To this end,  introduces the concept of the population posterior which depends on the unknown size of the dataset. This concept allows one to apply online VI with respect to the expected ELBO over the population.
3.2 Tricks of the Trade for SVI
The convergence speed of SGD, forming the basis of SVI, depends on the variance of the gradient estimates. Smaller gradient noise allows for larger learning rates and leads to faster convergence. This section covers tricks of the trade in the context of SVI, such as adaptive learning rates and variance reduction. Some of these approaches are generally applicable in SGD.
Adaptive Learning Rate and Mini-batch Size
. Due to the law of large numbers, increasing the mini-batch size reduces the stochastic gradient noise, allowing larger learning rates. To accelerate the learning procedure, one can either optimally adapt the mini-batch size for a given learning rate, or optimally adjust the learning rate to a fixed mini-batch size.
We begin by discussing learning rate adaptation. In each iteration, the empirical gradient variance can guide the adaptation of the learning rate which is inversely proportional to the gradient noise. Popular optimization methods that make use of this idea include RMSProp, AdaGrad , AdaDelta  and Adam . These methods are not specific to SVI but are frequently used in this context; for more details we refer interested readers to .
 first introduced adaptive learning rates for the global variational parameter in SVI, where the optimal learning rate was shown to satisfy
Above, denotes the optimal global variational parameter, and the current estimate. is the covariance matrix of the variational parameter in this mini-batch. Since is unknown,  showed how to estimate the optimal learning rate in an online fashion.
Instead of adapting the learning rate, the mini-batch size can be adapted while keeping the learning rate fixed. This achieves similar effects [10, 26, 37, 184]. In order to decrease the SGD variance,  proposed to choose the mini-batch size proportionally to the value of the objective function relative to its optimum. In practice, the estimated gradient noise covariance and the magnitude of the gradient are used to estimate the optimal mini-batch size.
In addition to controlling the optimization path through the learning rate and mini-batch size, we can reduce the variance, thereby enabling larger gradient steps. Variance reduction is often employed in SVI to achieve faster convergence. As follows, we summarize the literature on how to accomplish this goal via control variates, non-uniform sampling, and other approaches.
A control variate is a stochastic term that can be added to the stochastic gradient such that its expectation remains the same, but its variance is reduced . A control variate needs to be correlated with the stochastic gradient, and easy to compute. Using control variates for variance reduction is common in Monte Carlo simulation and stochastic optimization [165, 208]. Several authors have suggested the use of control variates in the context of SVI [146, 208, 78, 154].
As a prominent example, we discuss the stochastic variance reduced gradient (SVRG) method . In SVRG, one constructs a control variate which takes advantage of previous gradients from all data points, and one exploits that gradients along the optimization path are correlated. The standard stochastic gradient update is replaced by
indicates the estimated objective (here the negative ELBO) based on the current set of mini-batch indices, is a snapshot of after every iterations, and is the batch gradient computed over all the data points, . Since has expectation zero, it is a control variate.
SVRG requires a full pass through the dataset every iteration to compute the full gradients, even though a full pass can be relaxed to a very large mini-batch for large data sets. For smooth but not strongly convex objectives, SVRG was shown to achieve the asymptotic convergence rate , compared to of SGD. Many other control variates are used in practice [146, 203, 140]. We present another popular type of a control variate, the score function control variate, in Section 4.2.
Instead of subsampling data points with equal probability, non-uniform sampling can be used to select mini-batches with a lower gradient variance. Several authors suggested variants of importance sampling in the context of mini-batch selection [32, 148, 226, 55]. Although effective, these methods are not always practical, as the computational complexity of the sampling mechanism relates to the dimensionality of model parameters . Alternative methods aim at de-correlating similar points and sampling diversified mini-batches. These methods include stratified sampling , where one samples data from pre-defined subgroups based on meta-data or labels, clustering-based sampling 
, which amounts to clustering the data using k-means and then sampling data from every cluster with adjusted probabilities, and diversified mini-batch sampling[223, 224] using repulsive point processes to suppress the probability of data points with similar features in the same mini-batch. All of these methods have been shown to reduce variance and can also be used for learning on imbalanced data.
A number of alternative methods have been developed that contribute to variance reduction for SVI. A popular approach relies on Rao-Blackwellization, which is used in . The Rao-Blackwellization theorem (see Appendix A.5) generally states that a conditional estimation has lower variance if there exists a valid statistic that it can be conditioned on. Inspired by Rao-Blackwellization, the local expectation gradients method  has been proposed. This method splits the computation of the gradient of the ELBO into a Monte Carlo estimation and an exact expectation so that the contribution of each latent dimension to the gradient variance is optimally taken into account. Another related approach has been developed for SVI, which averages expected sufficient statistics over a sliding window of mini-batches to obtain a natural gradient with smaller mean squared error .
3.3 Collapsed, Sparse, and Distributed VI
In contrast to using stochastic optimization for faster convergence, this section presents methods that leverage the structure of certain models to achieve the same goal. In particular, we focus on collapsed, sparse, parallel, and distributed inference.
Collapsed variational inference (CVI) relies on the idea of analytically integrating out certain model parameters [188, 94, 182, 197, 97, 64, 83]. Due to the reduced number of parameters to be estimated, inference is typically faster. Collapsed inference is commonly constrained in the traditional conjugate exponential families, where the ELBO assumes an analytical form during marginalization. For these models, one can either marginalize out these latent variables before the ELBO is derived, or eliminate them afterwards [83, 64].
Several authors have proposed CVI for topic models [188, 94] where one can either collapse the topic proportions  or the topic assignments . In addition to these model specific derivations,  unifies existing model-specific CVI approaches and presents a general collapsed inference method for models in the conjugate exponential family class.
The computational benefit of CVI depends strongly on the statistics of the collapsed variables. Additionally, collapsing latent random variables can make other inference techniques tractable. For models such as topic models, we can collapse the discrete variables and only infer the continuous ones. This allows the usage of inference networks (Section6) [122, 180].
More generally, CVI does not solve all problems. On the one side, integrating out certain model variables makes the ELBO tighter, since the marginal likelihood does not have to get lower-bounded in these variables. On the other hand, besides mathematical challenges, marginalizing variables can introduce additional dependencies between variables. For example, collapsing the global variables in Latent Dirichlet Allocation introduces non-local dependencies between the assignment variables, making distributed inference harder.
Sparse inference introduces additional low-rank approximations into the variational approach, enabling more scalable inference [177, 195, 63]. Sparse inference can be either interpreted as a modeling choice or as an inference scheme .
Sparse inference methods are often encountered in the Gaussian Process (GPs) literature. The computational cost of learning GPs is , where is the number of data points. This cost is caused by the inversion of the kernel matrix of size , which hinders the application of GPs to big data sets. The idea of sparse inference in GPs  is to introduce inducing points. Inducing points can be interpreted as pseudo-inputs that reflect the original data, but yield a more sparse representation since . With inducing points, only a sized matrix needs to be inverted, and consequently the computational complexity of this method is .  collapses the distribution of inducing points, and  further extends this work to a stochastic version with a computational complexity of . Additionally, sparse inducing points make inference in Deep GPs tractable .
Parallel and Distributed Inference
Variational inference can be adjusted to distributed computing scenarios, where subsets of the data or parameters are distributed among several machines. [135, 219, 49, 138, 21]. Distributed inference schemes are often required in large scale scenarios, where data and computations are distributed across several machines. Independent latent variable models are trivially parallelizable. However, model specific designs such as reparametrizations might be required to enable efficient distributed inference . Current computing resources make VI applicable to web-scale data analysis .
4 Generic VI: Beyond the conjugate exponential family
In this section, we review techniques which aim at making VI more generic. This includes making VI applicable to a broader class of models, and also to make VI more automatic, eliminating the need for model-specific calculations. This makes VI more accessible and easier to use.
Variational inference was originally limited to conditionally conjugate models, for which the ELBO could be computed analytically before it was optimized [67, 220]. In this section, we introduce methods that relax this requirement and simplify inference. Central to this section are stochastic gradient estimators of the ELBO that can be computed for a broader class of models.
We start with the Laplace approximation in Section 4.1 and illustrate its limitations. We will then introduce black box variational inference (BBVI) which removes the need for analytic solutions. We discuss BBVI methods that rely on the REINFORCE or score function gradient in Section 4.2 and a different form of BBVI, which uses reparameterization gradients, in Section 4.3. Other approaches for non-conjugate VI are discussed in Section 4.4.
4.1 Laplace’s Method and Its Limitations
While not being the main focus of this survey, we briefly review the Laplace approximation as an alternative to non-conjugate inference. The Laplace (or Gaussian) approximation approximates the posterior by a Gaussian distribution
. To this end, one seeks the maximum of the posterior and computes the inverse of its Hessian. These two entities are taken as the mean and covariance of the Gaussian posterior approximation. To make this approach feasible, the log posterior needs to be twice-differentiable. According to the Bernstein von Mises theorem (a.k.a. Bayesian central limit theorem), the posterior approaches a Gaussian asymptotically in the limit of large data, and the Laplace approximation becomes exact (provided that the model is under-parameterized). The approach can be applied to approximate the maximum a posteriori (MAP) mean and covariance, predictive densities, and marginal posterior densities . The Laplace method has also been extended to more complex models such as belief networks with continuous variables .
This approximation suffers mainly from being purely local and depending only on the curvature of the posterior around the optimum; KL minimization typically approximates the posterior covariance more accurately. Additionally, the Laplace approximation is limited to the Gaussian variational family and does not apply to discrete variables . Computationally, the method requires the inversion of a potentially large Hessian, which can be costly in high dimensions. This makes this approach intractable in setups with a large number of parameters.
4.2 REINFORCE Gradients
In classical VI, the ELBO is first derived analytically, and then optimized. This procedure is usually restricted to models in the conditionally conjugate exponential family . For many models, including Bayesian deep learning architectures or complex hierarchical models, the ELBO contains intractable expectations with no known or simple analytical solution. Even if an analytic solution is available, the analytical derivation of the ELBO often requires time and mathematical expertise. In contrast, BBVI proposes a generic inference algorithm for which only the generative process of the data has to be specified. The main idea is to represent the gradient as an expectation, and to use Monte Carlo techniques to estimate this expectation.
As discussed in Section 2, in general VI aims at maximizing the ELBO, which is equivalent to minimizing the KL divergence between the variational posterior and target distribution. To maximize the ELBO, one needs to follow the gradient or stochastic gradient of the variational parameters. The key insight of BBVI is that one can obtain an unbiased gradient estimator by sampling from the variational distribution without having to compute the ELBO analytically[146, 154].
For a broad class of models, the gradient of the ELBO can be expressed as an expectation with respect to the variational distribution :
The full gradient , involving the expectation over , can now be approximated by a stochastic gradient estimator by sampling from :
where . Thus, BBVI provides black box gradient estimators for VI. Moreover, it only requires the practitioner to provide the joint distribution of observations and latent variables without the need to derive the gradient of the ELBO explicitly. The quantity is also known as the score function and is part of the REINFORCE algorithm .
The derivation of Eq. 13 requires the log derivative trick which can be applied to any bound. While the ELBO in combination with the KL results in Eq. 13, other divergence measures lead to additional terms in the REINFORCE gradients (B.3 in ).
A direct implementation of stochastic gradient ascent based on Eq. 14 suffers from high variances of the estimated gradients. Much of the success of BBVI can be attributed to variance reduction through Rao-Blackwellization and control variates . As one of the most important advancements of modern approximate inference, BBVI as been extended and made amortized inference feasible, see Section 6.1.
Variance Reduction for BBVI
BBVI requires a different set of techniques for variance reduction than those reviewed for SVI in Section 3.2. In contrast to SVI where the noise resulted from subsampling from a finite set of data points, the BBVI noise originates from random variables with possibly infinite support. In this case, techniques such as SVRG are not applicable, since the full gradient is not a sum over finitely many terms and cannot be kept in memory. Hence, BBVI involves a different set of control variates and other methods which shall briefly be reviewed here.
The arguably most important control variate in BBVI is the score function control variate , where one subtracts a Monte Carlo expectation of the score function from the gradient estimator:
As required, the score function control variate has expectation zero under the variational distribution. The weight is selected such that it minimizes the variance of the gradient.
While the original BBVI paper introduces both Rao-Blackwellization and control variates,  points out that good choices for control variates might be model-dependent. They further elaborate on local expectation gradients, which take only the Markov blanket of each variable into account. A different approach is presented by , which introduces overdispersed importance sampling. By sampling from a proposal distribution that belongs to an overdispersed exponential family and that places high mass on the tails of the variational distribution, the variance of the gradient is reduced.
4.3 Reparameterization Gradient VI
An alternative to the REINFORCE gradients introduced in Section 4.2 are the so-called reparameterization gradients. These gradients are obtained by representing the variational distribution as a deterministic parametric transformation of a noise distribution. Empirically, reparameterization gradients are often found to have lower variance than REINFORCE gradients.
The reparameterization trick allows to estimate the gradient of the ELBO by Monte Carlo samples by representing random variables as deterministic functions of noise distributions. This gives low-variance stochastic gradients for a large class of models without having to compute analytic expectations.
In more detail, the trick states that a random variable with a distribution can be expressed as a transformation of a random variable that comes from a noise distribution, such as uniform or Gaussian. For example, if , then where [85, 160].
More generally, the random variable is given by a parameterized, deterministic function of random noise, . Importantly, the noise distribution is considered independent of the parameters of , and therefore and share the same parameters . This allows to compute any expectation over as an expectation over by the theory behind the change of variables in integrals. We can now build a stochastic gradient estimator of the ELBO by pulling the gradient into the expectation, and approximating it by samples from the noise distribution:
Often, the entropy term can be computed analytically, which can lead to a lower gradient variance .
Note that the gradient of the log joint distribution enters the expectation. This is in contrast to the REINFORCE gradient, where the gradient of the variational distribution is taken (Eq. 14). The advantage of taking the gradient of the log joint is that this term is more informed about the direction of the maximum posterior mode. The lower variance of the reparameterization gradient may be attributed to this property.
While the variance of this estimator (Eq. 16) is often lower than the variance of the score function gradient (Eq. 14), a theoretical analysis shows that this is not guaranteed, see Chapter 3 in .  shows that the reparameterization gradient can be divided into a path derivative and the score function. Omitting the score function in the vicinity of the optimum can result in an unbiased gradient estimator with lower variance. Reparameterization gradients are also the key to variational autoencoders [85, 160] which we discuss in detail in Section 6.2.
The reparameterization trick does not trivially extend to many distributions, in particular to discrete ones. Even if a reparameterization function exists, it may not be differentiable. In order for the reparameterization trick to apply to discrete distributions, the variational distributions require further approximations. Several groups have addressed this problem. In [75, 111], the categorical distribution is approximated with the help of the Gumbel-Max trick and by replacing the argmax operation with a softmax operator. A temperature parameter controls the degree to which the softmax can approximate the categorical distribution. The closer it resembles a categorical distribution, the higher the variance of the gradient. The authors propose annealing strategies to improve convergence. Similarly, a stick-breaking process is used in 
to approximate the Beta distribution with the Kumaraswamy distribution.
As many of these approaches rely on approximations of individual distributions, there is growing interest in more general methods that are applicable without specialized approximations. The generalized reparameterization gradient  achieves this by finding an invertible transformation between the noise and the latent variable of interest. The authors derive the gradient of the ELBO which decomposes the expected likelihood into the standard reparameterization gradient and a correction term. The correction term is only needed when the transformation weakly depends on the variational parameters. A similar division is derived by  which proposes an accept-reject sampling algorithm for reparameterization gradients that allows one to sample from expressive posteriors. While reparameterization gradients often demonstrate lower variance than the score function, the use of Monte Carlo estimates still suffers from the injected noise. The variance can be further reduced with control variates [123, 162] or Quasi-Monte Carlo methods .
4.4 Other Generalizations
Finally, we survey a number of approaches that consider VI in non-conjugate models but do not follow the BBVI principle. Since the ELBO for non-conjugate models contains intractable integrals, these integrals have to be approximated somehow, either using some form of Taylor approximations (including Laplace approximations), lower-bounding the ELBO further such that the resulting integrals are tractable, or using some form of Monte Carlo estimators. Approximation methods which involve inner optimization routines [16, 206, 222] often become prohibitively slow for practical inference tasks. In contrast, approaches based on additional lower bounds with closed-form updates [88, 207, 82] can be computationally more efficient. Examples include extensions of the variational message passing algorithm  to non-conjugate models [88, 207]. Furthermore, 
proposed a technique based on stochastic linear regression to estimate the parameters of a fixed variational distribution based on Monte Carlo approximations of certain sufficient statistics. Recently, proposed a hybrid approach, where inference is split into a conjugate and a non-conjugate part.
5 Accurate VI:
Beyond KL and Mean Field
In this section, we present various methods that aim at improving the accuracy of standard VI. Previous sections dealt with making VI scalable and applicable to non-conjugate exponential family models. Most of the work in those areas, however, still addresses the standard setup of MFVI and employs the KL divergence as a measure of distance between distributions. Here we review recent developments that go beyond this setup, with the goal of avoiding poor local optima and increasing the accuracy of VI. Inference networks, normalizing flows, and related methods may also be considered as non-standard VI, but are discussed in Section 6.
We start by reviewing the origins of MFVI in statistical physics and describe its limitations (Section 5.1). We then discuss alternative divergence measures in Section 5.2. Structured variational approximations beyond mean field are discussed in Section 5.3, followed by alternative methods that do not fall into the previous two classes (Section 5.4).
5.1 Origins and Limitations of Mean Field VI
Variational methods have a long tradition in statistical physics. The mean field method was originally applied to model spin glasses, which are certain types of disordered magnets where the magnetic spins of the atoms are not aligned in a regular pattern 
. A simple example for such a spin glass model is the Ising model, a model of binary variables on a lattice with pairwise couplings. To estimate the resulting statistical distribution of spin states, a simpler, factorized distribution is used as a proxy. This is done with the goal of approximating the marginal probabilities of the spins pointing up or down (also called ’magnetization’) as well as possible, while ignoring all correlations between the spins. The many interactions of a given spin with its neighbors are replaced by a single interaction between a spin and the effective magnetic field (a.k.a.mean field) of all other spins. This explains the name origin. Physicists typically denote the negative log posterior as an energy or Hamiltonian function. This language has been adopted by the machine learning community for approximate inference in both directed and undirected models, summarized in Appendix A.6 for the reader’s reference.
Mean field methods were first adopted in neural networks by Anderson and Peterson in 1987 , and later gained popularity in the machine learning community [172, 79, 143]. The main limitation of mean field approximations is that they explicitly ignore correlations between different variables e.g., between the spins in the Ising model. Furthermore,  showed that the more possible dependencies are broken by the variational distribution, the more non-convex the optimization problem becomes. Conversely, if the variational distribution contains more structure, certain local optima do not exist. A number of initiatives to improve mean field VI have been proposed by the physics community and further developed by the machine learning community [190, 150, 143].
An early example of going beyond the mean field theory in a spin glass system is the Thouless-Anderson-Palmer (TAP) equation approach , which introduces perturbative corrections to the variational free energy. A related idea relies on power expansions , which has been extended and applied to machine learning models by various authors [80, 145, 142, 158, 185]. Additionally, information geometry provides an insight into the relation between MFVI and TAP equations [186, 187].  further connects TAP equations with divergence measures. We refer the readers to  for more information. Next, we review the recent advances beyond MFVI based on divergence measures other than the KL divergence.
5.2 VI with Alternative Divergences
The KL divergence often provides a computationally convenient method to measure the distance between two distributions. It leads to analytically tractable expectations for certain model classes. However, traditional Kullback-Leibler variational inference (KLVI) suffers from problems such as underestimating posterior variances . In other cases, it is unable to break symmetry when multiple modes are close , and is a comparably loose bound . Due to these shortcomings, a number of other divergence measures have been proposed which we survey here.
New divergence measures beyond the KL divergence do not only play a role in VI, but also in related approximate inference methods such as EP. Some recent extensions of EP [125, 204, 101, 228] can be viewed as classical EP with alternative divergence measures . While these methods are sophisticated, a practitioner will find them difficult to use due to complex derivations and limited scalability. Recent developments of VI focus mainly on a unified framework in a black box fashion to allow for scalability and accessibility. BBVI rendered the application of other divergence measures, such as the divergence , possible while maintaining the efficiency and simplicity of the method.
In this section, we introduce relevant divergence measures and show how to use them in the context of VI. The KL divergence, as discussed in Section 2.1, is a special form of the -divergence, while the -divergence is a special form of the -divergence. All above divergences can be written in the form of the Stein discrepancy.
The -divergence is a family of divergence measures with interesting properties from an information geometrical and computational perspective [4, 6]. Both the KL divergence and the Hellinger distance are special cases of the -divergence.
where . With this definition of -divergences, a smaller leads to more mass-covering effects, while a larger
results in zero-forcing effects, meaning that the variational distribution avoids areas of low posterior probability. For, we recover standard VI (involving the KL divergence).
For , is a lower bound on the log marginal likelihood. Interestingly, Eq. 18 also admits negative values of , in which case it becomes an upper bound. Note that in this case, is not a divergence. Among various possible definitions of the -divergence, only Renyi’s formulation leads to a bound (Eq. 18) in which the marginal likelihood cancels.
-Divergence and Generalized VI
is a convex function with . For example, the KL divergence is represented by the -divergence with , and the Pearson distance is an -divergence with .
In general, only specific choices of result in a bound that does only trivially depend on the marginal likelihood, and which is therefore useful for VI.
 lower-bounded the marginal likelihood using Jensen’s inequality:
Above, is an arbitrary concave function with . This formulation recovers the true marginal likelihood for , the standard ELBO for , and -VI for . For , the authors propose the following function:
Above, is a free parameter that can be optimized, and which absorbs the bound’s dependence on the marginal likelihood. The authors show that the terms up to linear order in correspond to the KL divergence, whereas higher order polynomials are correction terms which make the bound tighter. This connects to earlier work on TAP equations [190, 150] (see Section 5.1), which generally did not result in a bound.
Stein Discrepancy and VI
Stein’s method  was first proposed as an error bound to measure how well an approximate distribution fits a distribution of interest. The Stein discrepancy has been adapted to modern VI [108, 109, 110, 59]. Here, we introduce the Stein discrepancy and two VI methods that use it: Stein Variational Gradient Descent (SVDG)  and operator VI . These two methods share the same objective but are optimized in different manners.
indicates a set of smooth, real-valued functions. When and are identical, the divergence is zero. More generally, the more similar and are, the smaller is the discrepancy.
The second term in Eq. 20 involves an expectation under the intractable posterior. Therefore, the Stein discrepancy can only be used in VI for classes of functions for which the second term is equal to zero. We can find a suitable class with this property as follows. We define by applying a differential operator on another function , where is only restricted to be smooth:
where . The operator is constructed in such a way that the second expectation in Eq. 20 is zero for arbitrary ; all operators with this property are valid operators . A popular operator that fulfills this requirement is the Stein operator:
Both operator VI  and SVGD  use the Stein discrepancy with the Stein operator to construct the variational objective. The main difference between these two methods lies in the optimization of the variational objective using the Stein discrepancy. Operator VI  uses a minimax (GAN-style) formulation and BBVI to optimize the variational objective directly; while Stein Variational Gradient Descent (SVGD)  uses a kernelized Stein discrepancy. With a particular choice of the kernel and , it can be shown that SVGD determines the optimal perturbation in the direction of the steepest gradient of the KL divergence . SVGD leads to a scheme where samples in the latent space are sequentially transformed to approximate the posterior. As such, the method is reminiscent of, though formally distinct from, a normalizing flow approach  (see Section 6.3).
5.3 Structured Variational Inference
MFVI assumes a fully-factorized variational distribution; as such, it is unable to capture posterior correlations. Fully factorized variational models have limited accuracy, especially when the latent variables are highly dependent such as in models with hierarchical structure. This section examines variational distributions which are not fully factorized, but contain dependencies between the latent variables. These structured variational distributions are more expressive, but often come at higher computational costs.
Allowing a structured variational distribution to capture dependencies between latent variables is a modeling choice; different dependencies may be more or less relevant and depend on the model under consideration. For example, structured variational inference for LDA  shows that maintaining global structure is vital, while structured variational inference for the Beta Bernoulli Process  shows that maintaining local structure is more important. As follows, we review structured inference for hierarchical models, and discuss VI for time series.
For many models, the variational approximation can be made more expressive by maintaining dependencies between latent variables, but these dependencies make it harder to estimate the gradient of the variational bound. Hierarchical variational models (HVM) are a black box VI framework for structured variational distributions which applies to a broad class of models. In order to capture dependencies between latent variables, one starts with a mean-field variational distribution , but instead of estimating the variational parameters , one places a prior over them and marginalizes them out:
The new variational distribution thus captures dependencies through the marginalization procedure. Sampling from this distribution is also possible by simulating the hierarchical process. The resulting ELBO can be made tractable by further lower-bounding the resulting entropy and sampling from the hierarchical model. Notably, this approach is used in the development of the variational Gaussian Process (VGP) , a particular HVM. The VGP applies a Gaussian Process to generate variational estimates, thus forming a Bayesian non-parametric prior. Since GPs can model a rich class of functions, the VGP is able to confidently approximate diverse posterior distributions .
Another method that established dependencies between latent variables is copula VI [199, 60]. Instead of using a fully factorized variational distribution, copula VI assumes the variational family form:
is the copula distribution, which is a joint distribution over the marginal cumulative distribution functions. This copula distribution restores the dependencies among the latent variables.
VI for Time Series
One of the most important model classes in need of structured variational approximations are time series models. Significant examples include Hidden Markov Models (HMM) and Dynamic Topic Models (DTM) . These models have strong dependencies between time steps, leading traditional fully factorized MFVI to produce unsatisfying results. When using VI for time series, one typically employs a structured variational distribution that explicitly captures dependencies between time points, while remaining fully-factorized in the remaining variables [17, 76, 45, 12]. This commonly requires model specific approximations. [76, 45] derive SVI for popular time series models including HMMs, hidden semi-Markov models (HSMM), and hierarchical Dirichlet process-HMMs. Moreover,  derive an accelerated SVI for HSMMs. [12, 11] derive a structured BBVI algorithm for non-conjugate latent diffusion models.
5.4 Other Non-Standard VI Methods
In this section, we cover a number of miscellaneous approaches which fall under the broad umbrella of improving the accuracy of VI, but would not be categorized as alternative divergence measures or structured models.
VI With Mixture Distributions
Mixture distributions form a class of very flexible distributions, and have been used in VI since the 1990s [74, 79]. Due to their flexibility as well as computational difficulties, advancing VI for mixture models has been of continuous interest [51, 170, 58, 124, 8]. To fit a mixture model, we can make use of auxiliary bounds , a fixed point update , or enforce additional assumptions such as using uniform weights . Inspired by boosting methods, recently proposed methods fit mixture components in a successive manner [58, 124]. Here, Boosting VI and variational boosting [58, 124] refine the approximate posterior iteratively by adding one component at a time while keeping previously fitted components fixed. In a different approach, 
utilizes stochastic policy search methods found in the Reinforcement Learning literature for fitting Gaussian mixture models.
VI by Stochastic Gradient Descent
Stochastic gradient descent on the negative log posterior of a probabilistic model can, under certain circumstances, be seen as an implicit VI algorithm. Here we consider SGD with constant learning rates (constant SGD) [113, 114], and early stopping .
Constant SGD can be viewed as a Markov chain that converges to a stationary distribution; as such, it resembles Langevin dynamics . The variance of the stationary distribution is controlled by the learning rate.  shows that the learning rate can be tuned to minimize the KL divergence between the resulting stationary distribution and the Bayesian posterior. Additionally,  derive formulas for this optimal learning rate which resemble AdaGrad  and its relatives. A generalization of SGD that includes momentum and iterative averaging is presented in . In contrast,  interprets SGD as a non-parametric VI scheme. The paper proposes a way to track entropy changes in the implicit variational objective based on estimates of the Hessian. As such, the authors consider sampling from distributions that are not stationary.
Robustness to Outliers and Local Optima
Since the ELBO is a non-convex objective, VI benefits from advanced optimization algorithms that help to escape from poor local optima. Variational tempering  adapts deterministic annealing [164, 136]
to VI, making the cooling schedule adaptive and data-dependent. Temperature can be defined globally or locally, where local temperatures are specific to individual data points. Data points with associated small likelihoods under the model (such as outliers) are automatically assigned a high temperature. This reduces their influence on the global variational parameters, making the inference algorithm more robust to local optima. Variational tempering can also be interpreted as data re-weighting, the weight being the inverse temperature. In this context, lower weights are assigned to outliers. Other means of making VI more robust include the trust-region method , which uses the KL divergence to tune the learning progress and avoids poor local optima, and population VI , which averages the variational posterior over bootstrapped data samples for more robust modeling performance.
6 Amortized variational inference and deep learning
Finally, we review amoritzed variational inference. Consider the setup of Section 3, where each data point is goverend by its latent variable with variational parameter . Traditional VI makes it necessary to optimize a for each data point , which is computationally expensive, in particular when this optimization is embedded a global parameter update loop. The basic idea behind amortized inference is to use a powerful predictor to predict the optimal based on the features of , i.e., . This way, the local variational parameters are replaced by a function of the data whose parameters are shared across all data points, i.e. inference is amortized.
6.1 Amortized Variational Inference
The term amortized inference refers to utilizing inferences from past computations to support future computations [50, 36]. For VI, amortized inference usually refers to inference over local variables. Instead of approximating separate variables for each data point, as shown in Figure 2(a), amortized VI assumes that the local variational parameters can be predicted by a parameterized function of the data. Thus, once this function is estimated, the latent variables can be acquired by passing new data points through the function, as shown in Figure 2(b). Deep neural networks used in this context are also called inference networks. Amortized VI with inference networks thus combines probabilistic modeling with the representational power of deep learning.
As an example, amortized inference has been applied to Deep Gaussian Processes (DGPs) . Inference in these models is intractable, which is why the authors apply MFVI with inducing points (see Section 3.3) . Instead of estimating these latent variables separately, however,  proposes to estimate these latent variables as functions of inference networks, allowing DGPs to scale to bigger datasets and speeding up convergence. Amortization can be also made iterative by feeding back the amortization error into the inference model .
6.2 Variational Autoencoders
Amortized VI has become a popular tool for inference in deep latent Gaussian models (DLGM). This leads to the concept of variational autoencoders (VAEs), which have been proposed independently by two groups [85, 160], and which are discussed in detail below. VAEs apply more generally than to DLGMs, but for simplicity we will restrict our discussion to this model class.
The Generative Model
In this paragraph we introduce the class of deep latent Gaussian models. The corresponding graphical model is depicted in Figure 2(b). The model employs a multivariate normal prior from which we draw a latent variable ,
More generally, this could be an arbitrary prior that depends on additional parameters . The likelihood of the model is:
Most importantly, the likelihood depends on through two non-linear functions and
. These are typically neural networks, which take the latent variables as an input and transform them in a non-linear way. The data are then drawn from a normal distribution centered around the transformed latent variables. The parameter entails the parameters of the networks and .
Deep latent Gaussian models are highly flexible density estimators. There exist many modified versions specific to other types of data. For example, for binary data, the Gaussian likelihood can be replaced by a Bernoulli likelihood. Next, we review how amortized inference is applied to this model class.
Most commonly, VAEs refer to deep latent Gaussian models which are trained using inference networks.
VAEs employ two deep sets of neural networks: a top-down generative model as described above, mapping from the latent variables to the data , and a bottom-up inference model which approximates the posterior . Commonly, the corresponding neural networks are referred to as the generative network and the recognition network, or sometimes as decoder and encoder networks.
In order to approximate the posterior, VAEs employ an amortized mean-field variational distribution:
The conditioning on indicates that the local variational parameters associated with each data point are replaced by a function of the data. This amortized variational distribution is typically chosen as:
Similar to the generative model, the variational distribution employs non-linear mappings and of the data in order to predict the approximate posterior distribution of . The parameter summarizes the corresponding neural network parameters.
The main contribution of [85, 160] was to derive a scalable and efficient training scheme for deep latent variable models. During optimization, both the inference network and the generative network are trained jointly to optimize the ELBO.
The key to training these models is the reparameterization trick (Section 4.3). We focus on the ELBO contribution form a single data point . First, we draw samples from a noise distribution. We also employ a reparameterization function , such that realize samples from the approximate posterior . For Eq. 23, the most common reparametrization function takes the form , where and are parameterized by . One obtains an unbiased Monte Carlo estimator of the VAE’s ELBO by
This stochastic estimate of the ELBO can subsequently be differentiated with respect to and to obtain an estimate of the gradient.
The reparameterization trick also implies that the gradient variance is bounded by a constant . The drawback of this approach however is that we require the approximate posterior to be reparameterizable.
A Probabilistic Encoder-Decoder Perspective
The term variational autoencoder arises from the fact that the joint training of the generative and recognition network resembles the structure of autoencoders, a class of unsupervised, deterministic models. Autoencoders are deep neural networks that are trained to reconstruct their inputs as closely as possible. Importantly, the neural networks involved in autoencoders have an hourglass structure, meaning that there is a small number of units in the inner layers that prevent the neural network from learning the trivial identity mapping. This ’bottleneck’ forces the network to learn a useful and compact representation of the data.
In contrast, VAEs are probabilistic models, but they have a close correspondence to classical autoencoders. It turns out that the hidden variables of the VAE can be thought of as the intermediate representations of the data in the bottleneck of an autoencoder. During VAE training, one injects noise into this intermediate layer, which has a regularizing effect. In addition, the KL divergence term between the prior and the approximate posterior forces the latent representation of the VAE to be close to the prior, leading to a more homogeneous distribution in latent space that generalizes better to unseen data. When the variance of the noise is reduced to zero and the prior term is omitted, the VAE becomes a classical autoencoder.
6.3 Advancements in VAEs
Since the proposal of VAEs, an ever-growing number of extensions have been proposed. While exhaustive coverage of the topic would require a review article in its own right, we summarize a few important extensions. While several model extensions of the VAE have been proposed, this review puts a bigger emphasis on inference procedures. We will report on extensions that modify the variational approximation , the model , and finally discuss the dying units problem when the posterior of some latent units remains close to the prior during the optimization.
Flexible Variational Distributions
Traditional VI, including VAE training, relies on parametric inference models. The approximate posterior can be an explicit parametric distribution, such as a Gaussian or discrete distribution . We can use more flexible distributions, for example by transforming a simple parametric distribution. Here, we review VAE with implicit distributions, normalizing flow, and importance weighted VAE. Using more flexible variational distributions reduces not only the approximation error but also the amortization error, i.e. the error introduced by replacing the local variational parameters by a parametric function .
Implicit distributions can be used in VI since a closed-form density function is not a strict requirement for the inference model; all we need is to be able to sample from it. As detailed below, their reparameterization gradients can still be computed. In addition to the standard reparameterization approach, the entropy contribution to the gradient has to be estimated. Implicit distributions for VI is an active area of research [102, 72, 81, 129, 119, 105, 107, 196, 210].
VI requires the computation of the log density ratio . When is implicit, the standard training procedure faces the problem that log density ratio is intractable. In this case, one can employ a Generative Adversarial Networks (GAN)  style discriminator that discriminates the prior from the variational distribution, [119, 102]. This formulation is very general and can be combined with other ideas, for example a hierarchical structure [202, 217] .
Normalizing flow [40, 159, 41, 29, 84] presents another way to utilize flexible variational distributions. The main idea behind normalizing flow is to transform a simple (e.g. mean field) approximate posterior into a more expressive distribution by a series of successive invertible transformations.
To this end, we first draw a random variable , and transform it using an invertible, smooth function . Let . Then, the new distribution is
It is necessary that we can compute the determinant, since the variational approach requires us to also estimate the entropy of the transformed distribution. By choosing the transformation function such that is easily computable, this normalizing flow constitutes an efficient method to generate multimodal distributions from a simple distribution. As variants, linear time-transformations, Langevin and Hamiltonian flow , as well as inverse autoregressive flow  and autoregressive flow  have been proposed.
Normalizing flow and the previously mentioned implicit distribution share the common idea of using transformations to transform simple distributions into more complicated ones. A key difference is that for normalizing flows, the density of can be estimated due to an invertible transformation function.
One final approach that utilizes flexible variational distributions is the importance weighted variational autoencoder (IWAE) which was originally proposed to tighten the variational bound  and can be reinterpreted to sample from a more flexible distribution . IWAEs require samples from the approximate posterior which are weighted by the ratio
The authors show that the more samples are evaluated, the tighter the variational bound becomes, implying that the true log likelihood is approached in the limit . A reinterpretation of IWAEs, suggests that they are identical to VAEs but sample from a more expressive distribution which converges pointwise to the true posterior as . As IWAEs introduce a biased estimator, additional steps to obtain potentially better variance-bias trade-offs can be taken, such as in[139, 152, 153] .
Modeling Choices of
Modeling choices affect the performance of deep latent Gaussian models. In particular improving the prior model in VAEs can lead to more interpretable fits and better model performance. 
proposed a method to utilize a structured prior for VAEs, combining the advantages of traditional graphical models and inference networks. These hybrid models overcome the intractability of non-conjugate priors and likelihoods by learning variational parameters of conjugate distributions with a recognition model. This allows one to approximate the posterior conditioned on the observations while maintaining conjugacy. As the encoder outputs an estimate of natural parameters, message passing, which relies on conjugacy, is applied to carry out the remaining inference.
Other approaches tackle the drawback of the standard VAE which is the assumption that the likelihood factorizes over dimensions. This can be a poor approximation, e.g., for images, for which a structured output model works better. The Deep Recurrent Attentive Writer  relies on a recurrent structure that gradually constructs the observations while automatically focusing on regions of interest. In comparison, PixelVAE  tackles this problem by modeling dependencies between pixels within an image, using a conditional model that factorizes as , where denotes the th dimension of observation . The dimensions are generated in a sequential fashion, which accounts for local dependencies within the pixels of an image. The expressiveness of the modeling choice comes at a cost. If the decoder is too strong, the inference procedure can fail to learn an informative posterior . This problem, known as the dying units problem, will be discussed in the paragraph below.
The Dying Units Problem
Certain modeling choices and parameter configurations impose problems in VAE training, such that learning a good low-dimensional representation of the data fails. A prominent such problem is known as the dying units problem. As detailed below, two main effects are responsible for this phenomenon: a too powerful decoder, and the KL divergence term.
In some cases, the expressiveness of the decoder can be so strong, that some dimensions of the variables are ignored, i.e. it might model independently of . In this case the true posterior is the prior , and thus the variational posterior tries to match the prior in order to satisfy the KL divergence in Eq. 24. Lossy variational autoencoders  circumvent this problem by conditioning the decoding distribution for each output dimension on partial input information. For example, in the case of images, the likelihood of a given pixel is only conditioned on the values of the immediate surrounding pixels and the global latent state. This forces the distribution to encode global information in the latent variables.
The KL divergence contribution to the VAE loss may exacerbate this problem. To see why, we can rewrite the ELBO as a sum of two KL divergences . If the model is expressive enough, the model is able to render the second term zero (independent of the value of ). In this case, in order to also satisfy the first term, the inference model places its probability mass to match the prior , failing to learn a useful representation of the data. Even if the decoder is not strong, the problem of dying units may arise in the early stages of the optimization where the approximate posterior does not yet carry relevant information about the data . This problem is more severe when the dimension of is high. In this situation, units are regularized towards the prior and might not be reactivated in the later stages of the optimization . To counteract the early influence of the KL constraint, an annealing scheme can be applied to the KL divergence term during training .
We have summarized recent advancements in variational inference. Here we outline some selected active research directions and open questions, including, but not limited to: theory of VI, VI and policy gradients, VI for deep learning (DL), and automatic VI.
Theory of VI
Despite progress in modeling and inference, few authors address theoretical aspects of VI [133, 95, 213]. One important direction is quantifying the approximation errors involved when replacing a true posterior with a simplified variational distribution . A related problem is the predictive error, e.g., when approximating the marginalization involved in a Bayesian predictive distribution using VI.
We also conjecture that VI theory could profit from a deeper connection with information theory. This was already exemplified in [187, 186]. Information theory also inspires the development of new models and inference schemes [13, 193, 2]. For example, the information bottleneck  has recently led to the deep variational information bottleneck . We expect more interesting results to come from fusing these two lines of research.
VI and Deep Learning
Despite its recent successes in various areas, deep learning still suffers from a lack of principled uncertainty estimation, a lack in interpretability of its feature representations, and difficulties in including prior knowledge. Bayesian approaches, such as Bayesian neural networks  and variational autoencoders (as reviewed in Section 6), are improving all these aspects. Recent work aims at using interpretable probabilistic models as priors for VAEs [77, 168, 91, 38]. In such models, VI is an essential component. Making VI computationally efficient and easy to implement in Bayesian deep architectures is becoming an important research direction .
VI and Policy Gradients
Policy gradient estimation is important for reinforcement learning (RL)  and stochastic control. The technical challenges in these applications are similar to VI [98, 99, 173, 110, 211] (See Appendix A.7). As an example, SVGD has been applied in the RL setting as the Stein policy gradient . The application of VI in RL is currently an active area of research.
Probabilistic programming allows practitioners to quickly implement and revise models without having to worry about inference. The user is only required to specify the model, and the inference engine will automatically perform the inference. Popular probabilistic programming tools include but are not limited to: Stan , which covers a large range of advanced VI and MCMC methods, Infer.Net , which is based on variational message passing and EP, Automatic Statistician  and Anglican , which mainly rely on sampling methods, Edward , which supports BBVI as well as Monte Carlo sampling, and Zhusuan , which features VI for Bayesian Deep learning. The longstanding goal of these tools is to change the research methodology in probabilistic modeling, allowing users to quickly revise and improve models and to make them accessible to a broader audience.
Despite current efforts to make VI more accessible to practitioners, its usage is still not straightforward for non-experts. For example, manually identifying posterior symmetries and breaking these symmetries is necessary to work with Infer.Net. Furthermore, variance reduction methods such as control variates can drastically accelerate convergence, but a model specific design of control variates is needed to obtain the best performance. At the time of writing, these problems are not yet addressed in current probabilistic programming toolboxes. We believe these and other directions are important to advance the impact of probabilistic modeling in science and technology.
In this paper, we review the recent major advances in variational inference from four perspectives: scalability, generality, accuracy, and amortized inference. The advancement of variational inference theory and the adoption of approximate inference in new machine learning models are developing rapidly. Although this field has grown in recent years, it remains an open question how to make VI more efficient, more accurate, and easier to use for non-experts. Further development, as discussed in the previous section, is needed.
The authors would like to thank Sebastian Nowozin, Francisco Ruiz, Tianfan Fu, Robert Bamler, and especially Andrew Hartnett and Yingzhen Li for comments and discussions that greatly improved the manuscript.
-  S. Ahn, A. Korattikara, and M. Welling. Bayesian posterior sampling via stochastic gradient fisher scoring. In ICML, 2012.
-  A. Alemi, I. Fischer, J. Dillon, and K. Murphy. Deep variational information bottleneck. In ICLR, 2017.
-  S. M. Ali and S. D. Silvey. A general class of coefficients of divergence of one distribution from another. Journal of the Royal Statistical Society., 1966.
-  S.I. Amari. Differential-geometrical methods in statistics. Springer, 1985.
-  S.I. Amari. Natural gradient works efficiently in learning. Neural computation, 10(2), 1998.
-  S.I. Amari. -divergence is unique, belonging to both -divergence and bregman divergence classes. IEEE Transactions on Information Theory, 55(11), 2009.
-  E. Angelino, M. J. Johnson, and R. P. Adams. Patterns of scalable bayesian inference. Foundations and Trends® in Machine Learning, 9(2-3), 2016.
-  O. Arenz, M.J. Zhong, and G. Neumann. Efficient gradient-free variational inference using policy search. In ICML, 2018.
-  A. Azevedo-Filho and R.D. Shachter. Laplace’s method approximations for probabilistic inferencein belief networks with continuous variables. In UAI, 1994.
-  L. Balles, J. Romero, and P. Hennig. Coupling adaptive batch sizes with learning rates. In UAI, 2017.
-  R. Bamler and S. Mandt. Dynamic word embeddings. In ICML, 2017.
-  R. Bamler and S. Mandt. Structured black box variational inference for latent time series models. In ICML WS, 2017.
-  D. Barber and F. Agakov. The IM algorithm: a variational approach to information maximization. In NIPS, 2003.
-  C. M. Bishop. Pattern recognition and machine learning. springer, 2006.
-  D. M. Blei, A. Kucukelbir, and J. D. McAuliffe. Variational inference: A review for statisticians. 2017.
-  D. M. Blei and J. D. Lafferty. Correlated topic models. In NIPS, volume 18, 2006.
-  D. M. Blei and J. D. Lafferty. Dynamic topic models. In ICML, 2006.
-  L. Bottou. Large-scale machine learning with stochastic gradient descent. Springer, 2010.
-  S. R. Bowman, L. Vilnis, O. Vinyals, A. M. Dai, R. Józefowicz, and S. Bengio. Generating sentences from a continuous space. In CoNLL.
-  P.P. Boyle. Options: A monte carlo approach. Journal of financial economics, 4(3):323–338, 1977.
-  T. Broderick, N. Boyd, A. Wibisono, A. C. Wilson, and M. I. Jordan. Streaming variational bayes. In NIPS, 2013.
-  S. Brooks, A. Gelman, G. Jones, and X.L Meng. Handbook of markov chain monte carlo. CRC press, 2011.
-  A. Buchholz, F. Wenzel, and S. Mandt. Quasi-Monte Carlo Variational Inference. In ICML, 2018.
-  T. D. Bui, J. Yan, and R. E. Turner. A unifying framework for sparse gaussian process approximation using power expectation propagation. arXiv:1605.07066, 2016.
-  Y. Burda, R. Grosse, and R. Salakhutdinov. Importance weighted autoencoders. arXiv:1509.00519, 2015.
-  R.H. Byrd, G.M. Chin, J. Nocedal, and Y.C. Wu. Soamploample size selection in optimization methods for machine learning. Mathematical programming, 134(1), 2012.
-  L. Le Cam. Asymptotic methods in statistical decision theory. Springer Science & Business Media, 2012.
-  B. Carpenter, D. Lee, M. A. Brubaker, A. Riddell, A. Gelman, B. Goodrich, J. Guo, M. Hoffman, M. Betancourt, and P Li. Stan: A probabilistic programming language. Journal of Statistical Software, 2016.
-  X. Chen, D. P. Kingma, T. Salimans, Y. Dua, P. Dhariwal, J. Schulman, I. Sutskever, and P. Abbeel. Variational lossy autoencoder. In ICLR, 2017.
-  C. Cremer, X Li, and D. Duvenaud. Inference suboptimality in variational autoencoders. In ICML, 2018.
-  C. Cremer, Q. Morris, and D. Duvenaud. Reinterpreting importance-weighted autoencoders. In ICLR WS, 2017.
-  D. Csiba and P. Richtárik. Importance sampling for minibatches. arXiv:1602.02283, 2016.
-  I. Csiszár. Eine informationstheoretische ungleichung und ihre anwendung auf beweis der ergodizitaet von markoffschen ketten. Magyer Tud. Akad. Mat. Kutato Int. Koezl., 8, 1964.
-  Z.W. Dai, A. Damianou, J. González, and N. Lawrence. Variational auto-encoded deep Gaussian processes. In ICLR, 2016.
-  A. Damianou and N. Lawrence. Deep gaussian processes. In AISTATS, 2013.
-  P. Dayan, G. E. Hinton, R. M. Neal, and R. S. Zemel. The helmholtz machine. Neural computation, 7(5), 1995.
-  S. De, A. Yadav, D. Jacobs, and T. Goldstein. Automated inference using adaptive batch sizes. In AISTATS, 2017.
-  Z.W. Deng, R. Navarathna, P. Carr, S. Mandt, Y.S. Yue, I. Matthews, and G. Mori. Factorized variational autoencoders for modeling audience reactions to movies. In CVPR, 2017.
-  A. B. Dieng, D. Tran, R. Ranganath, J. Paisley, and D. M. Blei. Variational inference via chi-upper bound minimization. In NIPS, 2017.
-  L. Dinh, D. Krueger, and Y. Bengio. NICE: non-linear independent components estimation. arXiv:1410.8516, 2014.
-  L. Dinh, J. Sohl-Dickstein, and S. Bengio. Density estimation using real nvp. In ICLR, 2017.
-  J. Duchi, E. Hazan, and Y. Singer. Adaptive subgradient methods for online learning and stochastic optimization. JMLR, 12, 2011.
-  D. Duvenaud, D. Maclaurin, and R. Adams. Early stopping as nonparametric variational inference. In AISTATS, 2016.
-  S. R. Eddy. Hidden Markov Models. Current opinion in structural biology, 6(3), 1996.
-  N. Foti, J. Xu, D. Laird, and E. Fox. Stochastic variational inference for hidden Markov models. In NIPS, 2014.
-  M. P. Friedlander and M. Schmidt. Hybrid deterministic-stochastic methods for data fitting. SIAM Journal on Scientific Computing, 34(3), 2012.
-  T.F Fu and Z.H. Zhang. CPSG-MCMC: Clustering-based preprocessing method for stochastic gradient mcmc. In AISTATS, 2017.
-  Y. Gal. Uncertainty in deep learning. PhD thesis, PhD thesis, University of Cambridge, 2016.
-  Y. Gal, M. van der Wilk, and C. Rasmussen. Distributed variational inference in sparse gaussian process regression and latent variable models. In NIPS, 2014.
-  S. Gershman and N. Goodman. Amortized inference in probabilistic reasoning. In CogSci, volume 36, 2014.
-  S. J. Gershman, M. D. Hoffman, and D. M. Blei. Nonparametric variational inference. In ICML, 2012.
Probabilistic Machine Learning and Artificial Intelligence.Nature, 521(7553), 2015.
-  I. Goodfellow, Y. Bengio, and A. Courville. Deep Learning. MIT Press, 2016.
-  I. Goodfellow, J. Pouget-Abadie, M. Mirza, B. Xuand D. Warde-Farley, S. Ozair, A. Courville, and Y. Bengio. Generative adversarial nets. In NIPS, 2014.
-  P. K. Gopalan, S. Gerrish, M. Freedman, D. M. Blei, and D. M. Mimno. Scalable inference of overlapping communities. In NIPS, 2012.
K. Gregor, I. Danihelka, A. Graves, D. J. Rezende, and D. Wierstra.
Draw: A recurrent neural network for image generation.In ICML, 2015.
-  I. Gulrajani, K. Kumar, F. Ahmed, A. A. Taiga, F. V. Visin, D. Vazquez, and A. Courville. PixelVAE: A latent variable model for natural images. 2017.
-  F.J. Guo, X.Y. Wang, K. Fan, T. Broderick, and D. B. Dunson. Boosting variational inference. arXiv:1611.05559, 2016.
-  J. Han and Q. Liu. Stein variational adaptive importance sampling. In UAI, 2017.
-  S. Han, X. Liao, D. Dunson, and L. Carin. Variational gaussian copula inference. 2016.
-  G. Heinrich. Parameter estimation for text analysis, 2008.
-  P. Hennig. Approximate inference in graphical models. PhD thesis, University of Cambridge, 2011.
-  J. Hensman, N. Fusi, and N. D. Lawrence. Gaussian processes for big data. In UAI, 2013.
-  J Hensman, M Rattray, and N. D Lawrence. Fast variational inference in the conjugate exponential family. In NIPS, 2012.
-  M. Hoffman, F. R. Bach, and D. M. Blei. Online learning for Latent Dirichlet Allocation. In NIPS, 2010.
-  M. D. Hoffman and D. M. Blei. Structured stochastic variational inference. In AISTATS, 2015.
-  M D. Hoffman, D. M. Blei, C. Wang, and J. Paisley. Stochastic variational inference. JMLR, 14(1), 2013.
-  R. V. Hogg and A. T. Craig. Introduction to mathematical statistics.(5th edition). Upper Saddle River, New Jersey: Prentice Hall, 1995.
-  A. Honkela, T. Raiko, M. Kuusela, M. Tornio, and J. Karhunen. Approximate riemannian conjugate gradient learning for fixed-form variational bayes. JMLR, 11, 2010.
-  A. Honkela, M. Tornio, T. Raiko, and J. Karhunen. Natural conjugate gradient in variational inference. In ICONIP, 2007.
-  A. Honkela and H. Valpola. Online variational bayesian learning. 2003.
-  F. Huszár. Variational inference using implicit distributions. arXiv:1702.08235, 2017.
-  T. S. Jaakkol, L. K. Saul, and M. I. Jordan. Fast learning by bounding likelihoods in sigmoid type belief networks. NIPS, 1996.
-  T. S. Jaakkola and M. I. Jordan. Improving the mean field approximation via the use of mixture distributions. NATO ASI series D behaviroural and social sciences, 89, 1998.
-  E. Jang, S.X. Gu, and B. Poole. Categorical reparameterization with gumbel-softmax. In ICLR, 2017.
-  M. Johnson and A. Willsky. Stochastic variational inference for bayesian time series models. In ICML, 2014.
-  M. J. Johnson, D. Duvenaud, A. B. Wiltschko, S. R. Datta, and R. P. Adams. Structured VAEs: Composing probabilistic graphical models and variational autoencoders. In NIPS, 2016.
-  R. Johnson and T. Zhang. Accelerating stochastic gradient descent using predictive variance reduction. In NIPS, 2013.
-  M. I. Jordan, Z. Ghahramani, T. S. Jaakkola, and L. K. Saul. An introduction to variational methods for graphical models. Machine learning, 37(2), 1999.
-  H. J. Kappen and W. Wiegerinck. Second order approximations for probability models. In NIPS, 2001.
-  T. Karaletsos. Adversarial message passing for graphical models. In NIPS WS. 2016.
-  M. E. Khan, P. Baqué, F. Fleuret, and P. Fua. Kullback-Leibler Proximal Variational Inference. In NIPS, 2015.
-  N. King and N. Lawrence. Fast variational inference for gaussian process models through KL-correction. In ECML, 2006.
-  D. P. Kingma, T. Salimans, R. Jozefowicz, X. Chen, I. Sutskever, and M. Welling. Improving variational autoencoders with inverse autoregressive flow. In NIPS, 2016.
-  D. P. Kingma and M. Welling. Auto-encoding variational bayes. In ICLR, 2014.
-  D. P. Kingma and M. Welling. Stochastic gradient vb and the variational auto-encoder. In ICLR, 2014.
-  D.P. Kingma and J. Ba. Adam: A method for stochastic optimization. In ICLR, 2015.
-  D. A. Knowles and T. P. Minka. Non-conjugate variational message passing for multinomial and binary regression. In NIPS, 2011.
-  D.A. Knowles. Stochastic gradient variational Bayes for gamma approximating distributions. arXiv, page 1509.01631, 2015.
-  A. Korattikara, V. Rathod, K. Murphy, and M. Welling. Bayesian Dark Knowledge. arXiv:1506.04416, 2015.
R. G. Krishnan, U. Shalit, and D. Sontag.
Deep kalman filters.In NIPS WS, 2015.
-  A. Kucukelbir and D. M. Blei. Population empirical bayes. In UAI, 2015.
-  S. Kullback and R. A. Leibler. On information and sufficiency. The annals of mathematical statistics, 22(1), 1951.
-  K. Kurihara, M. Welling, and Y. W. Teh. Collapsed variational dirichlet process mixture models. In IJCAI, 2007.
-  S. Lacoste-Julien, F. Huszár, and Z. Ghahramani. Approximate inference for the loss-calibrated bayesian. In AISTATS.
-  P. S. Laplace. Memoir on the probability of the causes of events. Statistical Science, 1(3), 1986.
-  M. Lázaro-Gredilla, S. Van Vaerenbergh, and N. D. Lawrence. Overlapping mixtures of gaussian processes for the data association problem. Pattern Recognition, 45(4), 2012.
-  S. Levine. Reinforcement learning and control as probabilistic inference: Tutorial and review. arXiv:1805.00909, 2018.
-  S. Levine and V. Koltun. Variational policy search via trajectory optimization. In NIPS, 2013.
-  Y.Z. Li. Approximate Inference: New Visions. PhD thesis, 2018.
-  Y.Z. Li, J.M. Hernández-Lobato, and R.E. Turner. Stochastic expectation propagation. In NIPS, 2015.
-  Y.Z. Li and Q. Liu. Wild variational approximations. In NIPS WS, 2016.
-  Y.Z Li and r. e. Turner. Rényi divergence variational inference. In NIPS, 2016.
-  Y.Z Li, M. Rowland, T. Bui, D. Hernandez-Lobato, and R. Turner. Black-Box Alpha Divergence Minimization. In ICML, 2016.
-  Y.Z. Li, R. E. Turner, and Q. Liu. Approximate inference with amortised MCMC. arXiv:1702.08343, 2017.
-  Wu. Lin, N. Hubacher, and M.E. Khan. Variational message passing with structured inference networks. arXiv:1803.05589, 2018.
-  Q. Liu and Y. Feng. Two methods for wild variational inference. arXiv:1612.00081, 2016.
-  Q. Liu, J. Lee, and M. Jordan. A kernelized stein discrepancy for goodness-of-fit tests. In ICML, 2016.
-  Q. Liu and D. Wang. Stein variational gradient descent: A general purpose bayesian inference algorithm. In NIPS, 2016.
-  Y. Liu, P. Ramachandran, Q. Liu, and J. Peng. Stein variational policy gradient. In UAI, 2017.
C. J. Maddison, A. Mnih, and Y. W. Teh.
The concrete distribution: A continuous relaxation of discrete random variables, 2017.
-  S. Mandt and D. Blei. Smoothed gradients for stochastic variational inference. In NIPS, 2014.
-  S. Mandt, M. D. Hoffman, and D. M. Blei. A Variational Analysis of Stochastic Gradient Algorithms. In ICML, 2016.
-  S. Mandt, M. D. Hoffman, and D. M. Blei. Stochastic gradient descent as approximate bayesian inference. JMLR, 2017.
-  S. Mandt, J. McInerney, F. Abrol, R. Ranganath, and Blei. Variational Tempering. In AISTATS, 2016.
-  Joseph Marino, Yisong Yue, and Stephan Mandt. Iterative amortized inference. In ICML, 2018.
-  J. McInerney, R. Ranganath, and D. Blei. The population posterior and bayesian modeling on streams. In NIPS, 2015.
-  M.E.¬Khan, D. Nielsen, V. Tangkaratt, W. Lin, Y. Gal, and A. Srivastava. Fast and scalable bayesian deep learning by weight-perturbation in adam. arXiv:1806.04854, 2018.
-  L. Mescheder, S. Nowozin, and A. Geiger. Adversarial variational bayes: Unifying variational autoencoders and generative adversarial networks. In ICML, 2017.
-  L. Mescheder, S. Nowozin, and A. Geiger. The numerics of gans. In NIPS, 2017.
-  M. Mézard, G. Parisi, and M. A. Virasoro. Spin glass theory and beyond. 1990.
-  Y. Miao, L. Yu, and P. Blunsom. Neural variational inference for text processing. In ICML, 2016.
-  A. C. Miller, N. J. Foti, A. D’Amour, and R. P. Adams. Reducing reparameterization gradient variance. In NIPS, 2017.
-  A.C. Miller, N. Foti, and R. P. Adams. Variational boosting: Iteratively refining posterior approximations. In ICML, 2017.
-  T. Minka. Power EP. Technical report, Technical report, Microsoft Research, Cambridge, 2004.
-  T. Minka, J. M. Winn, J. P. Guiver, S. Webster, Y. Zaykov, B. Yangel, A. Spengler, and J. Bronskill. Infer.NET 2.6. http://research.microsoft.com/infernet, 2014.
-  T. P. Minka. Expectation propagation for approximate bayesian inference. In UAI, 2001.
-  T. P. Minka. Divergence measures and message passing. In Microsoft Research Technical Report, 2005.
-  S. Mohamed and B. Lakshminarayanan. Learning in implicit generative models. arXiv:1610.03483, 2016.
-  A. Müller. Integral probability metrics and their generating classes of functions. Advances in Applied Probability, 1997.
-  K. P. Murphy, Y. Weiss, and M. I. Jordan. Loopy belief propagation for approximate inference: An empirical study. In UAI, 1999.
-  C. Naesseth, F. Ruiz, S. Linderman, and D. M. Blei. Reparameterization gradients through acceptance-rejection sampling algorithms. In AISTATS, 2017.
-  S. Nakajima and S. Watanabe. Variational bayes solution of linear neural networks and its generalization performance. Neural Computation, 19(4), 2007.
-  E. Nalisnick and P. Smyth. Stick-breaking variational autoencoders. In ICLR, 2017.
-  R. Nallapati, W. Cohen, and J. Lafferty. Parallelized variational em for latent dirichlet allocation: An experimental evaluation of speed and scalability. In ICDM WS, 2007.
-  R. M. Neal. Probabilistic inference using markov chain monte carlo methods. In Technical Report, 1993.
-  R. M. Neal. Bayesian learning for neural networks, volume 118. Springer Science & Business Media, 2012.
-  W. Neiswanger, C. Wang, and E. Xing. Embarrassingly parallel variational inference in nonconjugate models. arXiv:1510.04163, 2015.
-  S. Nowozin. Debiasing evidence approximations: On importance-weighted autoencoders and jackknife variational inference. In ICLR, 2018.
-  C.J. Oates, M. Girolami, and N. Chopin. Control functionals for monte carlo integration. Journal of the Royal Statistical Society, 79(3):695–718, 2017.
-  M. Opper, M. Fraccaro, U. Paquet, A. Susemihl, and O. Winther. Perturbation theory for variational inference. NIPS WS, 2015.
-  M. Opper, U. Paquet, and O. Winther. Perturbative corrections for approximate inference in gaussian latent variable models. JMLR, 14(1), 2013.
-  M. Opper and D. Saad. Advanced mean field methods: Theory and practice. MIT press, 2001.
M. Opper and O. Winther.
A mean field algorithm for bayes learning in large feed-forward neural networks.In NIPS, 1996.
-  M. Opper and O. Winther. Tractable approximations for probabilistic models: The adaptive thouless-anderson-palmer mean field approach. Physical Review Letters, 86(17), 2001.
-  J. Paisley, D. M. Blei, and M. I. Jordan. Variational bayesian inference with stochastic search. In ICML, 2012.
-  G. Parisi. Statistical field theory. Addison-Wesley, 1988.
-  D Perekrestenko, V Cevher, and M Jaggi. Faster coordinate descent via adaptive importance sampling. In AISTATS, 2017.
-  C. Peterson and J. R. Anderson. A mean field theory learning algorithm for neural networks. Complex systems, 1, 1987.
-  T Plefka. Convergence condition of the TAP equation for the infinite-ranged ising spin glass model. Journal of Physics A: Mathematical and general, 15(6), 1982.
-  I. Porteous, D. Newman, A. Ihler, A. Asuncion, P. Smyth, and M. Welling. Fast collapsed gibbs sampling for latent dirichlet allocation. 2008.
-  T. Rainforth, R. Cornish, H. Yang, A. Warrington, and F. Wood. On Nesting Monte Carlo Estimators. In ICML, 2018.
-  T. Rainforth, A. R. Kosiorek, T. A. Le, C. J. Maddison, M. Igl, F. Wood, and Y. W. Teh. Tighter variational bounds are not necessarily better. In ICML, 2018.
-  R. Ranganath, S. Gerrish, and D. Blei. Black box variational inference. In AISTATS, 2014.
-  R. Ranganath, D. Tran, J. Altosaar, and D. M. Blei. Operator variational inference. In NIPS. 2016.
-  R. Ranganath, D. Tran, and D. M. Blei. Hierarchical variational models. In ICML, 2016.
-  R. Ranganath, C. Wang, D. Blei, and E. Xing. An adaptive learning rate for Stochastic Variational Inference. In ICML, 2013.
-  J. Raymond, A. Manoel, and M. Opper. Expectation propagation. arXiv:1409.6179, 2014.
-  D. Rezende and S. Mohamed. Variational Inference with Normalizing Flows. In ICML, 2015.
D. J. Rezende, S. Mohamed, and D. Wierstra.
Stochastic backpropagation and approximate inference in deep generative models.In ICML, 2014.
-  H. Robbins and S. Monro. A stochastic approximation method. The annals of mathematical statistics, 1951.
-  G. Roeder, Y. Wu, and D. Duvenaud. Sticking the landing: An asymptotically zero-variance gradient estimator for variational inference. In NIPS, 2017.
-  J. T. Rolfe. Discrete variational autoencoders. In ICLR, 2017.
-  K. Rose, E. Gurewitz, and G. Fox. A deterministic annealing approach to clustering. Pattern Recognition Letters, 11(9), 1990.
-  S. M. Ross. Simulation. Elsevier, 2006.
-  F. J.R. Ruiz, M. K. Titsias, and D. M. Blei. The generalized reparameterization gradient. In NIPS. 2016.
-  F. J.R. Ruiz, M. K. Titsias, and D. M. Blei. Overdispersed black-box variational inference. In UAI, 2016.
-  A. Saeedi, M. D. Hoffman, S. J. DiVerdi, A. Ghandeharioun, M. J. Johnson, and R. P. Adams. Multimodal prediction and personalization of photo edits with deep generative models. arXiv:1704.04997, 2017.
-  T. Salimans, D. P. Kingma, and M. Welling. Markov chain monte carlo and variational inference: Bridging the gap. In ICML, 2015.
-  T. Salimans and D. A. Knowles. Fixed-form variational posterior approximation through stochastic linear regression. Bayesian Analysis, 8(4), 2013.
-  M. Sato. Online model selection based on the variational Bayes. Neural Computation, 13(7):1649–1681, 2001.
-  L. K. Saul, T. Jaakkola, and M. I. Jordan. Mean field theory for sigmoid belief networks. Journal of artificial intelligence research, 4(1), 1996.
-  J. Schulman, S. Levine, P. Abbeel, M. Jordan, and P. Moritz. Trust region policy optimization. In ICML, 2015.
-  M. Seeger. Expectation propagation for exponential families. Technical report, 2005.
-  A. Shah, D. A. Knowles, and Z. Ghahramani. An Empirical Study of Stochastic Variational Algorithms for the Beta Bernoulli Process. In ICML, 2015.
-  J. Shi, J. Chen, J. Zhu, S. Sun, Y. Luo, Yihong Gu, and Yuhao Zhou. Zhusuan: A library for bayesian deep learning. arXiv:1709.05870, 2017.
-  E. Snelson and Z. Ghahramani. Sparse gaussian processes using pseudo-inputs. NIPS, 18, 2006.
-  C. K. Sønderby, T. Raiko, L. Maaløe, S. K. Sønderby, and O. Winther. How to train deep variational autoencoders and probabilistic ladder networks. In ICML, 2016.
-  B. K. Sriperumbudur, K. Fukumizu, A. Gretton, B. Schölkopf, and G. RG. Lanckriet. On integral probability metrics,phi-divergences and binary classification. arXiv preprint arXiv:0901.2698, 2009.
-  A. Srivastava and C. Sutton. Autoencoding variational inference for topic models. In ICLR, 2017.
-  C. Stein. A bound for the error in the normal approximation to the distribution of a sum of dependent random variables. In Berkeley Symposium on Mathematical Statistics and Probability, 1972.
-  J. Sung, Z. Ghahramani, and S. Y. Bang. Latent-space variational bayes. TPAMI, 30(12), 2008.
-  R. S. Sutton and A. G. Barto. Reinforcement learning: An introduction, volume 1. MIT press Cambridge, 1998.
-  L. SL. Tan. Stochastic variational inference for large-scale discrete choice models using adaptive batch sizes. Statistics and Computing, 27(1):237–257, 2017.
-  T. Tanaka. Estimation of third-order correlations within mean field approximation. In NIPS, 1998.
-  T. Tanaka. A theory of mean field approximation. In NIPS, 1999.
-  T. Tanaka. Information geometry of mean-field approximation. Neural Computation, 12(8), 2000.
-  Y. W. Teh, D. Newman, and M. Welling. A collapsed variational Bayesian inference algorithm for latent Dirichlet allocation. In NIPS, 2006.
-  L. Theis and M. Hoffman. A Trust-region Method for Stochastic Variational Inference with Applications to Streaming Data. In ICML, 2015.
-  D.J. Thouless, P. W. Anderson, and R. G. Palmer. Solution of ’solvable model of a spin glass’. Philosophical Magazine, 35(3), 1977.
-  T. Tieleman and G. Hinton. Lecture 6.5-rmsprop: Divide the gradient by a running average of its recent magnitude. COURSERA: Neural networks for machine learning, 4(2), 2012.
-  L. Tierney and J.B. Kadane. Accurate approximations for posterior moments and marginal densities. Journal of the american statistical association, 81(393), 1986.
-  N. Tishby, F. C. Pereira, and W. Bialek. The information bottleneck method. physics/0004057, 2000.
-  M. Titsias and M. Lázaro-Gredilla. Local Expectation Gradients for Black Box Variational Inference. In NIPS, 2015.
-  M. K. Titsias. Variational learning of inducing variables in sparse gaussian processes. In AISTATS, 2009.
-  M. K. Titsias. Learning model reparametrizations: Implicit variational inference by fitting mcmc distributions. arXiv:1708.01529, 2017.
M.K. Titsias and M. Lázaro-Gredilla.
Variational heteroscedastic gaussian process regression.In ICML, 2011.
-  D. Tolpin, J. W. van de Meent, H. Yang, and F. Wood. Design and implementation of probabilistic programming language anglican. arXiv:1608.05263, 2016.
-  D. Tran, D. M. Blei, and E. M. Airoldi. Copula variational inference. In NIPS, 2015.
-  D. Tran, A. Kucukelbir, A. B. Dieng, M. Rudolph, D. Liang, and D. M. Blei. Edward: A library for probabilistic modeling, inference, and criticism. arXiv:1610.09787, 2016.
-  D. Tran, R. Ranganath, and D. M. Blei. The Variational Gaussian Process. stat, 1050, 2016.
-  D. Tran, R. Ranganath, and D. M. Blei. Hierarchical implicit models and likelihood-free variational inference. arXiv:1702.08896, 2017.
-  G. Tucker, A. Mnih, C.J. Maddison, J. Lawson, and J. Sohl-Dickstein. REBAR: Low-variance, unbiased gradient estimates for discrete latent variable models. In NIPS, 2017.
-  M. J. Wainwright, T. S. Jaakkola, and A. S. Willsky. A new class of upper bounds on the log partition function. IEEE Transactions on Information Theory, 51(7), 2005.
-  M. J. Wainwright and M. I. Jordan. Graphical models, exponential families, and variational inference. Foundations and Trends® in Machine Learning, 1(1-2), 2008.
-  C. Wang, D. Blei, and L. Fei-Fei. Simultaneous image classification and annotation. In CVPR, 2009.
-  C. Wang and D. M. Blei. Variational inference in non-conjugate models. JMLR, 14(Apr), 2013.
-  C. Wang, X. Chen, A. J. Smola, and E. P. Xing. Variance reduction for stochastic gradient optimization. In NIPS, 2013.
-  C. Wang, J. Paisley, and D. M. Blei. Online variational inference for the hierarchical dirichlet process. In AISTATS, 2011.
-  D. Wang and Q. Liu. Learning to draw samples: With application to amortized MLE for generative adversarial learning. arXiv:1611.01722, 2016.
-  Y. Wang, G. Williams, E. Theodorou, and L. Song. Variational policy for guiding point processes. In ICML, 2017.
-  Y. X. Wang, A. Kucukelbir, and D. M. Blei. Robust Probabilistic Modeling with Bayesian Data Reweighting. In ICML, 2017.
-  Yixin Wang and David M Blei. Frequentist consistency of variational bayes. 2018.
-  M. Welling and Y. W. Teh. Bayesian learning via stochastic gradient langevin dynamics. In ICML, 2011.
-  R. J. Williams. Simple statistical gradient-following algorithms for connectionist reinforcement learning. Machine learning, 8(3-4):229–256, 1992.
-  J. Winn and C. M. Bishop. Variational message passing. JMLR, 6, 2005.
-  M. Yin and M. Zhou. Semi-implicit variational inference. In ICML, 2018.
-  M.D. Zeiler. Adadelta: an adaptive learning rate method. arXiv:1212.5701, 2012.
-  K. Zhai, J. Boyd-Graber, N. Asadi, and M. L. Alkhouja. Mr. lda: A flexible large scale topic modeling package using variational inference in mapreduce. In WWW, 2012.
-  C. Zhang. Structured Representation Using Latent Variable Models. PhD thesis, KTH Royal Institute of Technology, 2016.
-  C. Zhang, R. Bamler, M. Opper, and S. Mandt. Perturbative black box variational inference. In NIPS, 2017.
-  C. Zhang, C. H. Ek, X. Gratal, F. T. Pokorny, and H. Kjellström. Supervised hierarchical Dirichlet processes with variational inference. In ICCV WS, 2013.
-  C. Zhang, H. Kjellstrom, and S. Mandt. Determinantal point processes for mini-batch diversification. In UAI, 2017.
-  C. Zhang, C. Öztireli, S. Mandt, and G. Salvi. Active mini-batch sampling using repulsive point processes. arXiv:1804.02772, 2018.
-  P.L. Zhao and T. Zhang. Accelerating minibatch stochastic gradient descent using stratified sampling. arXiv:1405.3080, 2014.
-  P.L. Zhao and T. Zhang. Stochastic optimization with importance sampling for regularized loss minimization. In ICML, 2015.
-  S. Zhao, J. Song, and S. Ermon. Towards deeper understanding of variational autoencoding models. CoRR, abs/1702.08658, 2017.
-  S.D. Zhe, K.C. Lee, K. Zhang, and J. Neville. Online spike-and-slab inference with stochastic expectation propagation. In NIPS WS, 2016.
-  H. Zhu and R. Rohwer. Information geometric measurements of generalisation. Technical report, 1995.
a.1 ELBO and KL
We show that the difference between the marginal likelihood and the ELBO is the KL divergence between the variational distribution and the target distribution :
With this equivalence, the ELBO can be derived using either Jensen’s inequality as in Eq. 3, or using the KL divergence as .
a.2 Conjugate Exponential family
Many probabilistic models involve exponential family distributions. A random variable
is distributed according to a member of the exponential family if its probability distribution takes the form
is a vector of parameters,is the natural parameter, and are the sufficient statistics. Furthermore, is the base measure and is the log-normalizer. Many distributions fall into this class.
In the context of Bayesian statistics, certain exponential family distributions are conjugate pairs. A likelihood and prior distribution are a conjugate pair if the corresponding posterior distribution is in the same family as the prior. Examples for conjugate pairs include a Gaussian distribution with a Gaussian prior, a Poisson distribution with a gamma prior, or a multinomial distribution with a Dirichlet prior.
a.3 Variational Message Passing
. In a Bayesian network, the update for each node only requires information from the nodes in its Markov blanket, which includes this node’s parents, children, and co-parents of its children,
where indicates the set of parent nodes of , includes the set of the child nodes of , and indicates the th child node. indicates the set of parent nodes of . Hence, the update of one latent variable only depends on its parents, children, and its children’s co-parents.
If we further assume that the model is conjugate-exponential, see Section A.2, a latent variable can be updated by receiving all messages from its parents and children. Here, each child node has already received messages from its co-parents. Thus, to update each node, only nodes in this node’s Markov blanket are involved. Finally, is updated with the following three steps: a) receive messages from all parents and children , ; b) update ’s natural parameter ; c) update the expectation of ’s sufficient statistic .
Variational message passing provides a general message passing formulation for the MFVI. It enjoys all the properties of MFVI, but can be used in large scale Bayesian networks and can be automated easily. Together with EP, it forms the basis for the popular probabilistic programming tool Infer.Net .
a.4 Natural Gradients and SVI
We use the model example as shown in Figure 1 and assume that the true posterior of the global variable is the in exponential family:
We also assume that the variational distribution is in the same family:
Recall that is the variational parameter estimating the global variable . The subscript in and denotes that these are the natural parameter and log-normalizer of the global variable. The natural gradient of a function is given by , where is the Fisher information matrix.
 showed that the ELBO has a closed-form solution in terms of its variational parameters :
The constant contains all those terms that are independent of . The gradient of Equation 29 is given by
Importantly, when is in the exponential family, then it holds that . Thus, the natural gradient simplifies to
Hence, the natural gradient has a simpler form than the regular gradient.
Following the natural gradient has the advantage that we do not optimize in the Euclidean space, which is often not able to represent distances between distributions, but in Riemann space, where distance is defined by the KL divergence, i.e. distance between distributions. More information about the advantages of using natural gradients can be found in .