Variational Inference

Why Variational Inference?

Bayesian inference, in its idealized form, is a single-step procedure. We write down a prior \(p(\boldsymbol{z})\) over latent variables and a likelihood \(p_{\boldsymbol{\theta}}(\boldsymbol{x} \mid \boldsymbol{z})\) generating observations, so that the joint factorizes as \[ p_{\boldsymbol{\theta}}(\boldsymbol{x}, \boldsymbol{z}) = p_{\boldsymbol{\theta}}(\boldsymbol{x} \mid \boldsymbol{z}) p(\boldsymbol{z}). \] Here the prior is taken to be independent of \(\boldsymbol{\theta}\), in line with the notation of the variational inference literature. Bayes' rule then gives the posterior and the evidence, \[ \begin{align*} p_{\boldsymbol{\theta}}(\boldsymbol{z} \mid \boldsymbol{x}) &= \frac{p_{\boldsymbol{\theta}}(\boldsymbol{x}, \boldsymbol{z})} {p_{\boldsymbol{\theta}}(\boldsymbol{x})} \\\\ p_{\boldsymbol{\theta}}(\boldsymbol{x}) &= \int p_{\boldsymbol{\theta}}(\boldsymbol{x}, \boldsymbol{z}) \, d\boldsymbol{z}. \end{align*} \] Here and throughout, \(\boldsymbol{x}\) is an observation with \(0 \lt p_{\boldsymbol{\theta}}(\boldsymbol{x}) \lt \infty\), so that the ratio is defined. This density form of Bayes' rule is the analogue of Bayes' theorem for events, and we take it for granted here.

We use \(\int \cdot \, d\boldsymbol{z}\) throughout to denote integration against the appropriate reference measure on the latent space (Lebesgue measure for continuous \(\boldsymbol{z}\), counting measure for discrete \(\boldsymbol{z}\), product for mixed cases). All densities below are understood as densities with respect to this reference measure. The difficulty hidden in the formula is the evidence \(p_{\boldsymbol{\theta}}(\boldsymbol{x})\). When, for instance, \(\boldsymbol{z}\) is the parameter of an exponential-family likelihood and the prior is conjugate to that likelihood, the integral admits a closed form. Outside such regimes, as with neural-network likelihoods, non-conjugate priors, or discrete-continuous interactions, there is in general no analytic solution, and numerical quadrature scales prohibitively in the dimension of \(\boldsymbol{z}\). The posterior itself is a conditional distribution. Each of its probabilities is a conditional expectation of an indicator, and selecting these versions coherently across all sets of latent values is the question of regular conditional distributions, discussed but not developed on the linked page. What variational inference addresses is the intractability of computing with the posterior, not its existence.

Variational inference (VI) replaces the intractable computation with a tractable one of a different shape. We restrict attention to a family \(\mathcal{Q}\) of distributions over \(\boldsymbol{z}\) that we have decided in advance we can compute with. Typically \(\mathcal{Q}\) is a parametric family \(\{q_{\boldsymbol{\psi}}(\boldsymbol{z}) : \boldsymbol{\psi} \in \Psi\}\) indexed by variational parameters \(\boldsymbol{\psi}\). Inside this family we look for the member closest to the true posterior in the Kullback-Leibler divergence: \[ \boldsymbol{\psi}^* = \arg\min_{\boldsymbol{\psi} \in \Psi} D_{\mathrm{KL}}\!\left( q_{\boldsymbol{\psi}}(\boldsymbol{z}) \,\|\, p_{\boldsymbol{\theta}}(\boldsymbol{z} \mid \boldsymbol{x}) \right). \] The linked definition is stated for finite sample spaces. Below we use its analogue for densities, \(D_{\mathrm{KL}}(q \,\|\, p) = \int q \log (q/p) \, d\boldsymbol{z}\), with the same conventions for zeros.

The geometric picture is clean. The family \(\mathcal{Q}\) is a set of distributions, the true posterior is a single distribution generally lying outside \(\mathcal{Q}\), and \(\boldsymbol{\psi}^*\) selects the member of \(\mathcal{Q}\) closest to the posterior in the sense of KL. Because the KL divergence is not symmetric, this notion of "closest" depends on the order of the arguments. The choice \(D_{\mathrm{KL}}(q_{\boldsymbol{\psi}} \,\|\, p_{\boldsymbol{\theta}}(\cdot \mid \boldsymbol{x}))\) above, rather than the reverse, is itself a modeling decision, with consequences we return to in the closing section.

The substitution carries an immediate complication. The objective \(D_{\mathrm{KL}}(q_{\boldsymbol{\psi}} \,\|\, p_{\boldsymbol{\theta}}(\cdot \mid \boldsymbol{x}))\) is itself defined in terms of the intractable posterior, so it cannot be evaluated directly. The headline payoff of the next section is that an algebraic rearrangement converts KL minimization into an equivalent objective, the evidence lower bound (ELBO), which depends only on the joint \(p_{\boldsymbol{\theta}}(\boldsymbol{x}, \boldsymbol{z})\) and on \(q_{\boldsymbol{\psi}}(\boldsymbol{z})\). Both are tractable by construction. Maximizing the ELBO is exactly the same problem as minimizing the KL, but stated in a form one can actually compute.

Where needed, we adopt the absolute-continuity hypothesis \(q_{\boldsymbol{\psi}} \ll p_{\boldsymbol{\theta}}(\cdot \mid \boldsymbol{x})\), under which the KL is well-defined as an element of \([0, \infty]\), together with the integrability hypothesis, namely that \(\log p_{\boldsymbol{\theta}}(\boldsymbol{x}, \boldsymbol{z})\) and \(\log q_{\boldsymbol{\psi}}(\boldsymbol{z})\) are \(q_{\boldsymbol{\psi}}\)-integrable. Whether they hold depends on the model as well as on the family. A Gaussian \(q_{\boldsymbol{\psi}}\), whose density is positive on the whole latent space, violates the first against a posterior confined to a bounded region. The theorems below state which of these hypotheses, if any, they use.

With the ELBO substitution in hand, what was an integration problem has become an optimization problem over \(\boldsymbol{\psi}\). This trade is favorable because optimization in high dimensions has far more developed machinery than integration: stochastic gradient methods, automatic differentiation, and the entire toolkit of modern deep learning are available to attack it. Variational inference trades the asymptotic exactness of MCMC-style sampling for tractable, fast, gradient-based optimization in a chosen surrogate family, at the price that the result is only as good as the family allows.

The reader familiar with the variational autoencoder literature will recognize this framework as the foundation of the VAE. A VAE combines the amortized specialization of the variational program below, in which a neural network (the inference network) reads \(\boldsymbol{x}_n\) and produces the variational parameters \(\boldsymbol{\psi}_n = f_{\boldsymbol{\phi}}(\boldsymbol{x}_n)\), with the learning of the model parameters \(\boldsymbol{\theta}\) by maximizing the same ELBO. We treat the non-amortized case throughout. The amortized reduction is developed in Variational Autoencoders and builds on this framework rather than replacing it.

The notation is consistent with the earlier measure-theoretic probability pages and with the advanced machine learning literature. The probability space is \((\Omega, \mathcal{F}, \mathbb{P})\), with \(\boldsymbol{x}\) and \(\boldsymbol{z}\) measurable functions on \(\Omega\) taking values in observation- and latent-spaces respectively. Lowercase letters \(p_{\boldsymbol{\theta}}, q_{\boldsymbol{\psi}}\) denote densities with respect to the reference measures on the latent and observation spaces, extending to the observation space the convention adopted above for the latent space. When the reference measure is clear from context, we follow the variational-inference literature and sometimes call these densities "distributions".

The joint density \(p_{\boldsymbol{\theta}}(\boldsymbol{x}, \boldsymbol{z})\) is parameterized by the model parameters \(\boldsymbol{\theta}\), regarded as fixed throughout the inference problem (their estimation is a separate task that this page does not address). The variational density \(q_{\boldsymbol{\psi}}(\boldsymbol{z})\) is parameterized by the variational parameters \(\boldsymbol{\psi}\), and the variational family \(\mathcal{Q} = \{q_{\boldsymbol{\psi}} : \boldsymbol{\psi} \in \Psi\}\) is the set of all such densities as \(\boldsymbol{\psi}\) ranges over its parameter space \(\Psi\).

The ELBO is denoted \(\mathcal{L}(\boldsymbol{\theta}, \boldsymbol{\psi} \mid \boldsymbol{x})\) throughout. The symbol \(\mathcal{L}\) refers to the ELBO itself, which is maximized, and not to its negation, the variational free energy, which is minimized. The two conventions differ only by a sign and yield equivalent optimization problems, and this page uses the maximization convention. Where convenient, we follow a standard abuse of notation in the variational-inference literature and write \(q_{\boldsymbol{\psi}} \ll p_{\boldsymbol{\theta}}(\cdot \mid \boldsymbol{x})\) and \(D_{\mathrm{KL}}(q_{\boldsymbol{\psi}} \,\|\, p_{\boldsymbol{\theta}}(\cdot \mid \boldsymbol{x}))\) with the densities themselves as arguments. Such expressions are statements about the induced probability measures \(Q_{\boldsymbol{\psi}}\) and \(P_{\boldsymbol{\theta}}(\cdot \mid \boldsymbol{x})\) on the latent space.

The ELBO and the KL Decomposition

We define the evidence lower bound, prove that it is indeed a lower bound on the log-evidence by a direct application of Jensen's inequality, and then prove the central exact identity \[ \log p_{\boldsymbol{\theta}}(\boldsymbol{x}) = \mathcal{L}(\boldsymbol{\theta}, \boldsymbol{\psi} \mid \boldsymbol{x}) + D_{\mathrm{KL}}(q_{\boldsymbol{\psi}} \,\|\, p_{\boldsymbol{\theta}}(\cdot \mid \boldsymbol{x})). \] The Jensen bound and the exact decomposition are two views of the same fact, one an inequality and the other an equality. We present the inequality first because it is the standard entry point. The equality absorbs the inequality as a corollary and is the working tool in everything that follows.

The ELBO Definition

The two forms of the ELBO are connected by the identity \[ \log p_{\boldsymbol{\theta}}(\boldsymbol{x}, \boldsymbol{z}) = \log p_{\boldsymbol{\theta}}(\boldsymbol{x} \mid \boldsymbol{z}) + \log p(\boldsymbol{z}) \] together with the definition of the KL divergence. The first form is what we use to prove the lower-bound property and the exact decomposition. The second form drives variational autoencoders, where the expected log-likelihood is the "reconstruction" term and the KL to the prior is the "regularization" term. We work with the first form from here on.

Three remarks before proceeding. First, the ELBO depends on \(\boldsymbol{x}\) through the joint and on \(\boldsymbol{\psi}\) through the variational distribution. It depends on \(\boldsymbol{\theta}\) only via the model's joint density. Throughout this section \(\boldsymbol{\theta}\) is fixed, and we occasionally drop it from the notation when the dependence is not in play. Second, the expectation \(\mathbb{E}_{q_{\boldsymbol{\psi}}(\boldsymbol{z})}[\,\cdot\,]\) is against the variational distribution, not against the intractable posterior. That choice is the computational point of the construction, since \(q_{\boldsymbol{\psi}}\) has been chosen precisely so that such expectations can be computed, analytically for Gaussian \(q_{\boldsymbol{\psi}}\) with simple likelihoods and by Monte Carlo otherwise.

Third, the integrand \(\log p_{\boldsymbol{\theta}}(\boldsymbol{x}, \boldsymbol{z}) - \log q_{\boldsymbol{\psi}}(\boldsymbol{z})\) is the logarithm of the density ratio \(p_{\boldsymbol{\theta}}(\boldsymbol{x}, \boldsymbol{z}) / q_{\boldsymbol{\psi}}(\boldsymbol{z})\), which one should think of as the unnormalized "importance weight" of sampling \(\boldsymbol{z}\) from \(q_{\boldsymbol{\psi}}\) when the target is the joint \(p_{\boldsymbol{\theta}}(\boldsymbol{x}, \boldsymbol{z})\). The connection to importance sampling is more than notational. Tighter bounds built on it, such as IWAE, sharpen the ELBO by averaging several importance weights before taking the logarithm.

The Jensen Lower Bound

The ELBO is not just any quantity associated with \(q_{\boldsymbol{\psi}}\). It is a lower bound on the log-evidence \(\log p_{\boldsymbol{\theta}}(\boldsymbol{x})\), and the inequality holds for every \(q_{\boldsymbol{\psi}}\) in the variational family under the integrability of the two logarithms alone. The proof is a one-step application of Jensen's inequality.

The Jensen step deserves a moment of reflection. The inequality \(\mathbb{E}[\log Y] \leq \log \mathbb{E}[Y]\) is a one-direction relation, and beyond its equality case Jensen tells us nothing in general about how tight it is. The bound \(\mathcal{L} \leq \log p_{\boldsymbol{\theta}}(\boldsymbol{x})\) it produces is sharp in a precise sense. Because \(-\log\) is strictly convex, equality holds in the Jensen step if and only if the ratio \[ p_{\boldsymbol{\theta}}(\boldsymbol{x}, \boldsymbol{z}) / q_{\boldsymbol{\psi}}(\boldsymbol{z}) \] is \(q_{\boldsymbol{\psi}}\)-almost-everywhere equal to a constant \(c\), by the equality case of the Jensen statement linked above. The constant is then \(c = \int_{\{q_{\boldsymbol{\psi}} \gt 0\}} p_{\boldsymbol{\theta}}(\boldsymbol{x}, \boldsymbol{z}) \, d\boldsymbol{z}\), and the step governed by the covering property is an equality as well exactly when \(c = p_{\boldsymbol{\theta}}(\boldsymbol{x})\). Both steps are equalities precisely when, for almost every \(\boldsymbol{z}\), \[ q_{\boldsymbol{\psi}}(\boldsymbol{z}) = p_{\boldsymbol{\theta}}(\boldsymbol{z} \mid \boldsymbol{x}). \] In that case the ratio equals \(p_{\boldsymbol{\theta}}(\boldsymbol{x})\) wherever \(q_{\boldsymbol{\psi}}\) is positive, and the Jensen bound recovers the log-evidence exactly. The variational family \(\mathcal{Q}\) may or may not be rich enough to meet this condition.

The Exact KL Decomposition

The Jensen bound is suggestive but, on its own, gives no information about the gap between the ELBO and the log-evidence. The exact decomposition fills this in by identifying the gap as precisely the KL divergence from the variational distribution to the true posterior. This identity is the central result of variational inference.

We add a few words on the rigorous backing of this manipulation. The ratio \(q_{\boldsymbol{\psi}}(\boldsymbol{z}) / p_{\boldsymbol{\theta}}(\boldsymbol{z} \mid \boldsymbol{x})\) is, at the level of measures, the Radon-Nikodym derivative \(d q_{\boldsymbol{\psi}} / d p_{\boldsymbol{\theta}}(\cdot \mid \boldsymbol{x})\). The absolute-continuity hypothesis of the theorem guarantees its existence, and the integrability hypothesis guarantees finiteness of the expectations involved. The remaining steps split the logarithm by linearity of expectation and pull out the constant \(\log p_{\boldsymbol{\theta}}(\boldsymbol{x})\). Both are standard manipulations of Lebesgue integrals, justified term by term once each integrand is \(q_{\boldsymbol{\psi}}\)-integrable.

Two Consequences

The decomposition has two consequences that drive everything that follows on this page.

The Jensen bound is recovered as a corollary, under the decomposition's hypotheses.
Under the absolute-continuity hypothesis \(q_{\boldsymbol{\psi}} \ll p_{\boldsymbol{\theta}}(\cdot \mid \boldsymbol{x})\) of the decomposition theorem, the KL term is non-negative by Gibbs' inequality, which the linked page proves for finite sample spaces and which we use here in its density form without proof. The decomposition then gives \[ \log p_{\boldsymbol{\theta}}(\boldsymbol{x}) \geq \mathcal{L}(\boldsymbol{\theta}, \boldsymbol{\psi} \mid \boldsymbol{x}), \] with equality if and only if \(q_{\boldsymbol{\psi}} = p_{\boldsymbol{\theta}}(\cdot \mid \boldsymbol{x})\) almost everywhere under \(q_{\boldsymbol{\psi}}\). The corollary costs nothing beyond the hypotheses of the Jensen-based proof. The integrability hypothesis of that proof already makes \(p_{\boldsymbol{\theta}}(\boldsymbol{x}, \boldsymbol{z})\) positive for \(q_{\boldsymbol{\psi}}\)-almost every \(\boldsymbol{z}\), so every set of posterior probability zero has \(q_{\boldsymbol{\psi}}\)-probability zero, which is the absolute-continuity hypothesis. What the corollary adds is the exact gap, and with it the fact that the bound is tight exactly when \(q_{\boldsymbol{\psi}}\) is the posterior.

ELBO maximization is KL minimization.
Fix \(\boldsymbol{x}\) and \(\boldsymbol{\theta}\). Then \(\log p_{\boldsymbol{\theta}}(\boldsymbol{x})\) is a constant in \(\boldsymbol{\psi}\), and the decomposition becomes \[ \mathcal{L}(\boldsymbol{\theta}, \boldsymbol{\psi} \mid \boldsymbol{x}) = \log p_{\boldsymbol{\theta}}(\boldsymbol{x}) - D_{\mathrm{KL}}(q_{\boldsymbol{\psi}} \,\|\, p_{\boldsymbol{\theta}}(\cdot \mid \boldsymbol{x})). \] Maximizing \(\mathcal{L}\) over \(\boldsymbol{\psi}\) is therefore exactly the same optimization problem as minimizing the KL divergence over \(\boldsymbol{\psi}\). This equivalence licenses the entire variational program. The program converts the original problem \[ \boldsymbol{\psi}^* = \arg\min_{\boldsymbol{\psi}} D_{\mathrm{KL}}(q_{\boldsymbol{\psi}} \,\|\, p_{\boldsymbol{\theta}}(\cdot \mid \boldsymbol{x})), \] which involves the intractable posterior, into the equivalent problem of maximizing the ELBO, and the conversion involves no approximation. The ELBO involves only the joint \(p_{\boldsymbol{\theta}}(\boldsymbol{x}, \boldsymbol{z})\) and the variational distribution \(q_{\boldsymbol{\psi}}(\boldsymbol{z})\), and we can compute with both by construction.

The remaining work, namely the choice of variational family, gradient computation, and coordinate-ascent updates, occupies the rest of this page.

Choice of Variational Family

The variational family \(\mathcal{Q}\) is a modeling choice that the practitioner makes before running any algorithm, and it determines both the quality of the approximation and the cost of finding it. Two families are sufficiently common that they organize much of the variational literature. We introduce them here informally, and each receives a fuller treatment in the algorithmic sections that follow.

The first is the fixed-form family. One chooses a parametric distribution and optimizes its parameters by gradient methods on the ELBO. The most common choice is a multivariate Gaussian \[ q_{\boldsymbol{\psi}}(\boldsymbol{z}) = \mathcal{N}(\boldsymbol{z} \mid \boldsymbol{\mu}, \boldsymbol{\Sigma}) \] with \(\boldsymbol{\psi} = (\boldsymbol{\mu}, \boldsymbol{\Sigma})\). The covariance structure of \(\boldsymbol{\Sigma}\) (full, diagonal, or low-rank-plus-diagonal) trades expressivity against the cost of storage and inversion. A diagonal \(\boldsymbol{\Sigma}\) gives the mean-field Gaussian that reappears below. The fixed-form approach is the path taken by variational autoencoders and is the natural target of the gradient-based machinery developed in the next section.

The second is the mean-field family, in which the variational distribution factorizes across coordinate groups, \[ q_{\boldsymbol{\psi}}(\boldsymbol{z}) = \prod_{j=1}^{J} q_j(\boldsymbol{z}_j), \] with \(\boldsymbol{z}\) partitioned into blocks \(\boldsymbol{z}_1, \ldots, \boldsymbol{z}_J\) and each \(q_j\) a distribution on the \(j\)-th block. The factors \(q_j\) are not constrained to lie in any particular parametric family. Instead, their functional forms are derived from the model and the factorization by the coordinate-ascent procedure of a later section, which maximizes the ELBO with respect to one factor at a time while holding the others fixed.

The mean-field family is historically the original variational construction and remains the workhorse of classical Bayesian VI, particularly for models with conjugate updates between blocks. Its principal weakness is that it tends to underestimate posterior variance when the true posterior has inter-block dependencies. The weakness is visible already in two dimensions, where a mean-field Gaussian cannot represent any posterior correlation between the two coordinates.

The two families are not mutually exclusive. A mean-field factorization can have Gaussian factors, and mean-field Gaussian VI is a standard baseline. A fixed-form parametric family can also be made richer by adopting structured rather than diagonal covariances. Richer alternatives generalize both directions and are an active research frontier. They include normalizing flows, structured mean-field with learned dependency graphs, and implicit posteriors specified by sampling rather than by a density. We mention some of them in the closing section but do not develop them here.

Computing the ELBO Gradient

The decomposition theorem of the previous section reduces variational inference to the optimization problem \[ \arg\max_{\boldsymbol{\psi}} \mathcal{L}(\boldsymbol{\theta}, \boldsymbol{\psi} \mid \boldsymbol{x}). \] For the parametric variational families used in practice, such as Gaussians, mean-field exponential families, and neural-network-parameterized distributions, this is a continuous optimization in \(\boldsymbol{\psi}\), and the natural attack is gradient ascent on the ELBO. The technical obstacle is that the ELBO is an expectation against a distribution whose parameters are themselves what we are differentiating with respect to. Naive interchange of the gradient and the expectation is illegitimate, and recovering a valid gradient estimator is the work of this section.

Two standard estimators address the obstacle. Both are unbiased, and they differ in variance and in what they require of \(q_{\boldsymbol{\psi}}\). The reparameterization gradient applies when \(q_{\boldsymbol{\psi}}\) admits a sampling representation \(\boldsymbol{z} = g_{\boldsymbol{\psi}}(\boldsymbol{\epsilon})\) with \(\boldsymbol{\epsilon}\) drawn from a fixed base distribution. It pushes the \(\boldsymbol{\psi}\)-dependence inside the expectation, where ordinary differentiation under the integral sign applies, and typically yields a low-variance estimator at the cost of requiring \(g_{\boldsymbol{\psi}}\) and the log-joint to be differentiable. The score-function estimator, also called REINFORCE, applies much more widely, including when \(\boldsymbol{z}\) is discrete or \(g_{\boldsymbol{\psi}}\) is non-differentiable, provided \(\log q_{\boldsymbol{\psi}}\) is differentiable in \(\boldsymbol{\psi}\) and its support does not move with \(\boldsymbol{\psi}\). The estimator exploits an algebraic identity that converts the gradient of the expectation into a different expectation that can be estimated by Monte Carlo. The price is high variance and a corresponding need for variance-reduction techniques. We develop both estimators in turn, with the rigorous justification (in each case, an application of dominated convergence) explicit.

The Gradient Interchange Problem

Write the ELBO as an explicit integral: \[ \begin{align*} \mathcal{L}(\boldsymbol{\theta}, \boldsymbol{\psi} \mid \boldsymbol{x}) &= \int q_{\boldsymbol{\psi}}(\boldsymbol{z}) \, f_{\boldsymbol{\theta}}(\boldsymbol{x}, \boldsymbol{z}; \boldsymbol{\psi}) \, d\boldsymbol{z}, \\\\ f_{\boldsymbol{\theta}}(\boldsymbol{x}, \boldsymbol{z}; \boldsymbol{\psi}) &= \log p_{\boldsymbol{\theta}}(\boldsymbol{x}, \boldsymbol{z}) - \log q_{\boldsymbol{\psi}}(\boldsymbol{z}). \end{align*} \] The \(\boldsymbol{\psi}\)-dependence enters in two places. The integrand \(f_{\boldsymbol{\theta}}\) contains \(\log q_{\boldsymbol{\psi}}\), and the integrating measure is itself \(q_{\boldsymbol{\psi}}(\boldsymbol{z}) \, d\boldsymbol{z}\). The naive identity \[ \nabla_{\boldsymbol{\psi}} \mathbb{E}_{q_{\boldsymbol{\psi}}(\boldsymbol{z})}\!\left[ h(\boldsymbol{z}) \right] \stackrel{?}{=} \mathbb{E}_{q_{\boldsymbol{\psi}}(\boldsymbol{z})}\!\left[ \nabla_{\boldsymbol{\psi}} h(\boldsymbol{z}) \right] \] fails in general because the right-hand side ignores the contribution of the gradient acting on the measure \(q_{\boldsymbol{\psi}}\) itself. For an integrand that does not depend on \(\boldsymbol{\psi}\) (so that \(\nabla_{\boldsymbol{\psi}} h = 0\)), the right-hand side would be identically zero, while the left-hand side is generally non-zero whenever \(q_{\boldsymbol{\psi}}\) actually depends on \(\boldsymbol{\psi}\).

The legitimate operation underlying both estimators we develop below is differentiation under the integral sign, which the dominated convergence theorem licenses provided one can find a single integrable function dominating the \(\boldsymbol{\psi}\)-difference quotients uniformly in a neighborhood of the parameter value. The two estimators differ in where they apply this principle. The reparameterization gradient first rewrites the expectation against a fixed \(\boldsymbol{\psi}\)-independent measure (eliminating the second source of \(\boldsymbol{\psi}\)-dependence) and then differentiates under the integral. By contrast, the score-function estimator differentiates under the integral directly and pays for the \(\boldsymbol{\psi}\)-dependent measure with an algebraic identity that produces an additional term.

The Reparameterization Gradient

Suppose \(q_{\boldsymbol{\psi}}\) admits a reparameterization: a fixed base distribution \(p_0\) on some space \(\mathcal{E}\), typically the standard multivariate Gaussian \(\mathcal{N}(\mathbf{0}, \mathbf{I})\), and a deterministic map \(g_{\boldsymbol{\psi}} : \mathcal{E} \to \mathcal{Z}\) such that the law of \(g_{\boldsymbol{\psi}}(\boldsymbol{\epsilon})\) under \(\boldsymbol{\epsilon} \sim p_0\) coincides with \(q_{\boldsymbol{\psi}}\). The canonical example is the location-scale family for a Gaussian variational distribution. If \(q_{\boldsymbol{\psi}}(\boldsymbol{z}) = \mathcal{N}(\boldsymbol{z} \mid \boldsymbol{\mu}, \operatorname{diag}(\boldsymbol{\sigma}^2))\) with \(\boldsymbol{\psi} = (\boldsymbol{\mu}, \boldsymbol{\sigma})\), where \(\boldsymbol{\sigma}\) has positive entries, one takes \(p_0 = \mathcal{N}(\mathbf{0}, \mathbf{I})\) and \[ g_{\boldsymbol{\psi}}(\boldsymbol{\epsilon}) = \boldsymbol{\mu} + \boldsymbol{\sigma} \odot \boldsymbol{\epsilon}, \] where \(\odot\) denotes element-wise multiplication. This construction is the reparameterization trick of variational autoencoders, and the present treatment isolates the gradient identity from any specific parametric form.

Applied to the ELBO with \[ h_{\boldsymbol{\psi}}(\boldsymbol{z}) = \log p_{\boldsymbol{\theta}}(\boldsymbol{x}, \boldsymbol{z}) - \log q_{\boldsymbol{\psi}}(\boldsymbol{z}), \] the theorem gives, whenever its hypotheses hold for this choice of \(h_{\boldsymbol{\psi}}\), \[ \nabla_{\boldsymbol{\psi}} \mathcal{L}(\boldsymbol{\theta}, \boldsymbol{\psi} \mid \boldsymbol{x}) = \mathbb{E}_{p_0(\boldsymbol{\epsilon})}\!\left[ \nabla_{\boldsymbol{\psi}} \!\left( \log p_{\boldsymbol{\theta}}(\boldsymbol{x}, g_{\boldsymbol{\psi}}(\boldsymbol{\epsilon})) - \log q_{\boldsymbol{\psi}}(g_{\boldsymbol{\psi}}(\boldsymbol{\epsilon})) \right) \right], \] and the right-hand side is amenable to Monte Carlo estimation by drawing samples \(\boldsymbol{\epsilon}^{(1)}, \ldots, \boldsymbol{\epsilon}^{(K)} \sim p_0\) and averaging the bracketed gradient.

Evaluating the integrand at a sample \(\boldsymbol{\epsilon}\) requires computing \(\log q_{\boldsymbol{\psi}}(g_{\boldsymbol{\psi}}(\boldsymbol{\epsilon}))\). When the variational family has a known closed-form density, as for Gaussian, Gamma, Beta, or exponential-family distributions, one simply substitutes \(\boldsymbol{z} = g_{\boldsymbol{\psi}}(\boldsymbol{\epsilon})\) into that closed form.

When the family is defined only through the transform \(g_{\boldsymbol{\psi}}\), as in normalizing flows where the variational density is induced by composing Jacobian-tractable diffeomorphisms with the base distribution, the closed form is unavailable and the change-of-variables formula is the only route. If \(g_{\boldsymbol{\psi}} : \mathcal{E} \to \mathcal{Z}\) is a diffeomorphism onto its image with Jacobian \(\partial g_{\boldsymbol{\psi}} / \partial \boldsymbol{\epsilon}\), the densities are related by the multivariable change-of-variables formula, which we use without proof, \[ \log q_{\boldsymbol{\psi}}(\boldsymbol{z}) = \log p_0(\boldsymbol{\epsilon}) - \log \!\left| \det \frac{\partial g_{\boldsymbol{\psi}}}{\partial \boldsymbol{\epsilon}}(\boldsymbol{\epsilon}) \right|, \quad \boldsymbol{z} = g_{\boldsymbol{\psi}}(\boldsymbol{\epsilon}). \] The two routes agree where they overlap. For the location-scale Gaussian \[ g_{\boldsymbol{\psi}}(\boldsymbol{\epsilon}) = \boldsymbol{\mu} + \boldsymbol{\sigma} \odot \boldsymbol{\epsilon}, \] the Jacobian is \(\operatorname{diag}(\boldsymbol{\sigma})\) and the log-determinant is \(\sum_k \log \sigma_k\), and the change-of-variables expression reduces to the closed-form Gaussian log-density that one would have written down directly.

The practical relevance of the formula is therefore that it lifts the reparameterization trick to families without closed-form densities. Designing \(g_{\boldsymbol{\psi}}\) so that its Jacobian determinant is tractable is the central constraint on the variational family in such cases. Full-covariance Gaussians use a triangular Cholesky factor \(\mathbf{L}\) with positive diagonal, whose log-determinant is \(\sum_k \log L_{kk}\), and richer families such as normalizing flows compose Jacobian-tractable maps to extend expressiveness while preserving the closed-form log-determinant. For the closed-form Gaussian case, the reparameterization trick is developed on the variational autoencoder page, where the ELBO decomposition labels the KL to the prior as the regularization term.

A complementary technique, automatic differentiation variational inference (ADVI), addresses the case where the latent space \(\mathcal{Z}\) is a constrained subset of \(\mathbb{R}^D\) (positivity, simplex, bounded interval) by composing a fixed bijection \(T : \mathcal{Z} \to \mathbb{R}^D\) with a Gaussian variational distribution on the unconstrained space. The density on \(\mathcal{Z}\) is recovered by the same change-of-variables formula, and the resulting ELBO is differentiable in the variational parameters by automatic differentiation through \(T^{-1}\) and the model. ADVI is implemented in probabilistic programming systems such as Stan and PyMC, and it is a special case of the reparameterization gradient with the bijection chosen to handle parameter constraints rather than to expand expressiveness.

The crucial feature of all of the above is that the entire chain \[ \boldsymbol{\epsilon} \mapsto g_{\boldsymbol{\psi}}(\boldsymbol{\epsilon}) \mapsto \log p_{\boldsymbol{\theta}}(\boldsymbol{x}, \cdot) - \log q_{\boldsymbol{\psi}}(\cdot) \] is differentiable in \(\boldsymbol{\psi}\), so the gradient can be computed by ordinary automatic differentiation through the model and the variational density. This is what makes the reparameterization estimator the workhorse of variational autoencoders and, more broadly, of deep generative models trained with variational objectives whose latent variables admit a location-scale or related differentiable parameterization. In the amortized setting of the variational autoencoder the variational parameters \(\boldsymbol{\psi}_n = f_{\boldsymbol{\phi}}(\boldsymbol{x}_n)\) themselves depend on the data through an inference network, and the reparameterization map becomes \(g_{\boldsymbol{\phi}}(\boldsymbol{x}, \boldsymbol{\epsilon})\). This is a special case of the present theorem in which the gradient with respect to the encoder parameters \(\boldsymbol{\phi}\) flows through both the data-dependence of \(\boldsymbol{\psi}_n\) and the deterministic transform \(g\).

The Score-Function (REINFORCE) Estimator

When \(q_{\boldsymbol{\psi}}\) does not admit a differentiable reparameterization, the reparameterization theorem is unavailable. Examples include discrete latents and mixture distributions whose component assignment is itself a random variable. Implicit distributions, specified only by their sampler, lie outside both estimators as developed here, since each needs the log-density \(\log q_{\boldsymbol{\psi}}\). For the remaining cases one needs an estimator that requires no differentiability in \(\boldsymbol{z}\), only samples from \(q_{\boldsymbol{\psi}}\) and evaluations of the log-joint, of \(\log q_{\boldsymbol{\psi}}(\boldsymbol{z})\), and of its \(\boldsymbol{\psi}\)-gradient. The score-function estimator delivers this by exchanging the gradient of an expectation for an expectation of a different gradient.

The quantity \(\nabla_{\boldsymbol{\psi}} \log q_{\boldsymbol{\psi}}(\boldsymbol{z})\) is the score function of the variational distribution at \(\boldsymbol{z}\). Monte Carlo sampling of the right-hand side, drawing \(\boldsymbol{z}^{(k)} \sim q_{\boldsymbol{\psi}}\) and averaging the bracketed integrand, gives the score-function or REINFORCE estimator. It is the same construction that appears in policy-gradient reinforcement learning under the policy \(q_{\boldsymbol{\psi}}\) and reward \(h\).

Applied to the ELBO with \[ h_{\boldsymbol{\psi}}(\boldsymbol{z}) = \log p_{\boldsymbol{\theta}}(\boldsymbol{x}, \boldsymbol{z}) - \log q_{\boldsymbol{\psi}}(\boldsymbol{z}), \] the theorem gives, whenever its hypotheses hold for this choice of \(h_{\boldsymbol{\psi}}\), \[ \nabla_{\boldsymbol{\psi}} \mathcal{L} = \mathbb{E}_{q_{\boldsymbol{\psi}}(\boldsymbol{z})}\!\left[ \nabla_{\boldsymbol{\psi}} h_{\boldsymbol{\psi}}(\boldsymbol{z}) + h_{\boldsymbol{\psi}}(\boldsymbol{z}) \, \nabla_{\boldsymbol{\psi}} \log q_{\boldsymbol{\psi}}(\boldsymbol{z}) \right]. \]

The first term simplifies. Since \[ \nabla_{\boldsymbol{\psi}} h_{\boldsymbol{\psi}}(\boldsymbol{z}) = -\nabla_{\boldsymbol{\psi}} \log q_{\boldsymbol{\psi}}(\boldsymbol{z}) \] (only the \(\log q_{\boldsymbol{\psi}}\) term in \(h_{\boldsymbol{\psi}}\) carries \(\boldsymbol{\psi}\)-dependence), and the score has zero mean under its own distribution, \[ \mathbb{E}_{q_{\boldsymbol{\psi}}}[\nabla_{\boldsymbol{\psi}} \log q_{\boldsymbol{\psi}}] = 0 \] (the theorem applied with \(h \equiv 1\), whose domination hypothesis asks for \(q_{\boldsymbol{\psi}}(\boldsymbol{z})\) to be Lipschitz in \(\boldsymbol{\psi}\) near the parameter value with an integrable constant), the first term contributes zero. The ELBO gradient is therefore \[ \nabla_{\boldsymbol{\psi}} \mathcal{L} = \mathbb{E}_{q_{\boldsymbol{\psi}}(\boldsymbol{z})}\!\left[ \big( \log p_{\boldsymbol{\theta}}(\boldsymbol{x}, \boldsymbol{z}) - \log q_{\boldsymbol{\psi}}(\boldsymbol{z}) \big) \nabla_{\boldsymbol{\psi}} \log q_{\boldsymbol{\psi}}(\boldsymbol{z}) \right], \] and the score-function estimator of the ELBO gradient is the Monte Carlo average of this expression over samples drawn from \(q_{\boldsymbol{\psi}}\). Crucially, this formula requires only that we can sample from \(q_{\boldsymbol{\psi}}\), evaluate the log-joint, and evaluate the log-density of \(q_{\boldsymbol{\psi}}\) and its score. Neither \(q_{\boldsymbol{\psi}}\) nor the integrand needs to be differentiable in \(\boldsymbol{z}\), which is exactly the regime where reparameterization fails.

Bias, Variance, and the Choice of Estimator

As noted at the start of this section, both estimators are unbiased whenever the hypotheses of their theorems hold, with expectation equal to the exact ELBO gradient. They differ sharply in variance, and this difference dictates the choice in practice.

The reparameterization estimator typically achieves low variance, when the integrand is smooth with moderate derivatives, because the \(\boldsymbol{\psi}\)-dependence has been pushed into a single deterministic transform \(g_{\boldsymbol{\psi}}\), and the gradient is taken pathwise. Each Monte Carlo sample evaluates the derivative of a smooth composition, and no factor of the score \(\nabla_{\boldsymbol{\psi}} \log q_{\boldsymbol{\psi}}\) multiplies the integrand, which is the main source of variance in the score-function estimator. For Gaussian \(q_{\boldsymbol{\psi}}\) and smooth log-joint \(\log p_{\boldsymbol{\theta}}\), estimates from one or two samples per gradient step are typically sufficient to drive stochastic gradient optimization to convergence, and standard variational autoencoder training relies on this.

The score-function estimator, by contrast, suffers from variance that scales with the magnitude of \(h(\boldsymbol{z})\), which the control variates below can reduce toward the scale of \(h(\boldsymbol{z}) - \mathbb{E}_{q_{\boldsymbol{\psi}}}[h]\), and with the variance of the score function itself. In high-dimensional latent spaces, or when the integrand has wide range, single-sample estimates can have variance orders of magnitude larger than the reparameterization counterpart. Both situations are common in modern probabilistic models.

The standard remedy is variance reduction by control variates. We replace each component \[ \tilde{g}_i(\boldsymbol{z}) = h(\boldsymbol{z}) \, \partial_{\psi_i} \log q_{\boldsymbol{\psi}}(\boldsymbol{z}) \] of the naive estimator by \[ \tilde{g}_i^{\mathrm{cv}}(\boldsymbol{z}) = \tilde{g}_i(\boldsymbol{z}) - c_i \, b_i(\boldsymbol{z}), \] where \(b_i(\boldsymbol{z})\) is a baseline with \[ \mathbb{E}_{q_{\boldsymbol{\psi}}}[b_i(\boldsymbol{z})] = 0 \] and \(c_i \in \mathbb{R}\) is a coefficient to be chosen. The zero-mean condition guarantees that the modified estimator remains unbiased. A natural baseline is \[ b_i(\boldsymbol{z}) = \partial_{\psi_i} \log q_{\boldsymbol{\psi}}(\boldsymbol{z}) \] itself, which is zero-mean by the same identity that simplified the ELBO gradient derivation.

The coefficient \(c_i\) is then chosen to minimize the variance of \(\tilde{g}_i^{\mathrm{cv}}\). Since \[ \operatorname{Var}(\tilde{g}_i - c_i b_i) = \operatorname{Var}(\tilde{g}_i) - 2 c_i \operatorname{Cov}(\tilde{g}_i, b_i) + c_i^2 \operatorname{Var}(b_i) \] is a quadratic in \(c_i\), for \(\operatorname{Var}(b_i) \gt 0\) the optimum is \[ c_i^* = \frac{\operatorname{Cov}(\tilde{g}_i(\boldsymbol{z}), \, b_i(\boldsymbol{z}))}{\operatorname{Var}(b_i(\boldsymbol{z}))}, \] which can be estimated from earlier iterations or from a separate batch of samples, so that the coefficient does not depend on the samples it multiplies. Choosing \(b_i\) and \(c_i^*\) in this way can reduce variance by orders of magnitude while preserving unbiasedness. A further refinement, common in policy-gradient reinforcement learning, subtracts from \(h\) a learned scalar \(b(\boldsymbol{x})\) that depends on the observation but not on the sampled \(\boldsymbol{z}\) (a function of the state rather than of the action). This amounts to the control variate \(b(\boldsymbol{x}) \, \partial_{\psi_i} \log q_{\boldsymbol{\psi}}(\boldsymbol{z})\), which has zero mean because the score does. Effective control variates are an extensive research topic that we do not pursue further here.

The choice of estimator is therefore largely dictated by the structure of \(q_{\boldsymbol{\psi}}\). When a reparameterization is available, it is usually preferred. The literature on "black-box" variational inference, in which the score-function estimator is paired with control variates and applied to models for which no reparameterization is known, was developed to extend the variational program to models without model-specific derivations, including the discrete and combinatorial models that the reparameterization cannot reach. The two estimators cover complementary regimes, and a practitioner whose model contains both reparameterizable and non-reparameterizable latents can apply each estimator to the appropriate component.

Coordinate Ascent Variational Inference

The gradient-based estimators of the previous section operate within a fixed-form variational family. The practitioner chooses a parametric distribution \(q_{\boldsymbol{\psi}}\), such as a Gaussian, a Gaussian mixture, or a normalizing-flow-induced density, and optimizes its parameters by stochastic gradient ascent on Monte Carlo estimates of \(\nabla_{\boldsymbol{\psi}} \mathcal{L}\). The present section develops a complementary approach in which the variational family is the mean-field family \(q(\boldsymbol{z}) = \prod_{j=1}^{J} q_j(z_j)\) introduced in the earlier discussion of the variational family, with \(z_j\) now denoting the \(j\)-th block, and the functional form of each factor \(q_j\) is not chosen in advance but emerges from the optimization itself. This is sometimes called free-form variational inference, in contrast to the fixed-form setting of the previous section.

The free-form approach exchanges flexibility in the variational family for flexibility in the functional form of each factor. Under the mean-field assumption, the ELBO admits an explicit decomposition that makes coordinate-wise optimization tractable. With each \(q_j\) treated as the variational degree of freedom (rather than a finite-dimensional parameter \(\boldsymbol{\psi}_j\)), the optimal form of \(q_j^*\) holding the other factors fixed is determined by the model's log-joint via a closed-form expression. Cycling through the factors and updating each in turn yields the coordinate ascent variational inference (CAVI) algorithm. It is one of the most classical and historically important variational algorithms, and it remains the standard algorithm for mean-field VI in conjugate-exponential-family models.

The Mean-Field ELBO and the Coordinate Update

Substituting the mean-field factorization \(q(\boldsymbol{z}) = \prod_j q_j(z_j)\) into the ELBO and using \(\log q(\boldsymbol{z}) = \sum_j \log q_j(z_j)\) and the fact that each factor integrates to one, the ELBO decomposes as \[ \mathcal{L}(q) = \int q(\boldsymbol{z}) \log p_{\boldsymbol{\theta}}(\boldsymbol{x}, \boldsymbol{z}) \, d\boldsymbol{z} + \sum_{j=1}^{J} \mathbb{H}(q_j), \] where \(\mathbb{H}(q_j) = -\int q_j(z_j) \log q_j(z_j) \, dz_j\) is the (differential) entropy of the \(j\)-th factor. The first term is the expected log-joint under the factorized \(q\). The second is a sum of entropies that decouples completely across factors. This decomposition is what makes coordinate-wise optimization feasible. Holding all factors except \(q_j\) fixed, we find the dependence of the ELBO on \(q_j\) isolated to the \(j\)-th term of the entropy sum and to the contribution of \(q_j\) to the expectation in the first term.

The optimal coordinate update is given by the following theorem, the central result of mean-field variational inference.

Two structural features of this update deserve emphasis. First, the optimal \(q_j^*\) depends on the other factors \(q_{-j}\) only through the expected log-joint. In graphical-model language, the dependence is mediated entirely by the Markov blanket of node \(j\), since \[ \log p_{\boldsymbol{\theta}}(\boldsymbol{x}, \boldsymbol{z}) = \log p_{\boldsymbol{\theta}}(z_j \mid \boldsymbol{x}, \boldsymbol{z}_{-j}) + \log p_{\boldsymbol{\theta}}(\boldsymbol{x}, \boldsymbol{z}_{-j}), \] where the first term depends on the other latents only through that blanket and the second is constant in \(z_j\) and drops out of the proportionality.

Second, the functional form of \(q_j^*\) is dictated by the model's log-joint. When each complete conditional \(p_{\boldsymbol{\theta}}(z_j \mid \boldsymbol{x}, \boldsymbol{z}_{-j})\) is a member of an exponential family, as in models with conjugate prior structure, the expectation \(\mathbb{E}_{q_{-j}}[\log p]\), as a function of \(z_j\), is a linear function of the sufficient statistics of \(z_j\) plus the log base measure and a constant, and \(q_j^*\) inherits the same exponential-family form as the corresponding conditional. This is the source of the closed-form update equations that make CAVI practical. In effect, the model rather than the practitioner dictates the variational family.

The CAVI Algorithm

The optimal-update theorem yields an iterative algorithm. We initialize the factors \(\{q_j\}\) arbitrarily, cycle through \(j = 1, \ldots, J\), updating each \(q_j\) to its optimal form given the current values of the other factors, and repeat until convergence. As long as the hypotheses of the theorem hold at each update, each individual update increases the ELBO (or leaves it unchanged at a fixed point), so the sequence of ELBO values is monotonically non-decreasing. Provided the initial ELBO is finite, the sequence converges to a finite limit, since the ELBO is bounded above by the log-evidence.

We emphasize that this is convergence of the objective values, not of the iterates themselves. Monotone boundedness of \(\mathcal{L}(q^{(t)})\) does not by itself guarantee that the factors \((q_1^{(t)}, \ldots, q_J^{(t)})\) converge in the space of distributions. Convergence of the iterates to a stationary point requires additional regularity conditions, such as compactness of the variational parameter space, continuity of the update map, or the standard hypotheses for block coordinate ascent. They are not automatic, even for the conjugate-exponential models on which CAVI is typically run. In practice the iterates are observed to settle at a local optimum that depends on the initialization, since the ELBO is concave with respect to each \(q_j\) individually but generally non-concave jointly.

Algorithmically, after \(T\) sweeps:

The expectation in the inner loop is, in conjugate-exponential models, available in closed form as a function of the parameters of the other factors \(q_{-j}\). The update therefore reduces to a finite-dimensional parameter update on each iteration, and the apparent functional optimization collapses to ordinary numerical computation. In non-conjugate models the expectation is typically intractable and CAVI is replaced by gradient-based methods of the previous section, or by hybrid schemes that Monte-Carlo-estimate the inner expectation while preserving the coordinate ascent structure of the outer loop.

Bayesian Gaussian Mixture Models and Automatic Sparsity

The canonical application of CAVI is Bayesian inference in conjugate latent-variable models, of which the Bayesian Gaussian mixture model is the most pedagogically important example. With a Dirichlet prior on the mixing weights, conjugate Normal-Wishart priors on the component parameters, and the standard latent cluster-assignment variables, the mean-field factorization \[ q(\boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Lambda}, z_{1:N}) = q(\boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Lambda}) \prod_{n=1}^{N} q_n(z_n) \] treats the mixing weights \(\boldsymbol{\pi}\) and the component parameters \((\boldsymbol{\mu}_k, \boldsymbol{\Lambda}_k)\) as latent variables, part of \(\boldsymbol{z}\) in the notation of this page, rather than as fixed model parameters. It yields CAVI updates in which the factor over \((\boldsymbol{\pi}, \boldsymbol{\mu}, \boldsymbol{\Lambda})\) splits further and every factor's optimal form is itself a tractable member of the exponential family: Categorical for the cluster-assignment factors \(q_n(z_n)\), Dirichlet for the mixing-weight factor \(q(\boldsymbol{\pi})\), and Normal-Wishart for the component-parameter factors \(q(\boldsymbol{\mu}_k, \boldsymbol{\Lambda}_k)\). The update equations involve expected sufficient statistics under the other factors. We do not derive them on this page and instead highlight a single emergent phenomenon that is characteristic of the Bayesian treatment of the mixing weights.

Limitations of Mean-Field Inference

The mean-field assumption is structurally restrictive in a way that cannot be remedied by clever choice of factor functional form. The family \(q(\boldsymbol{z}) = \prod_j q_j(z_j)\) is, by construction, incapable of representing posterior dependence between latent variables. When the true posterior \(p_{\boldsymbol{\theta}}(\boldsymbol{z} \mid \boldsymbol{x})\) exhibits strong correlation between the \(z_j\), mean-field VI returns an approximation whose marginal variances typically underestimate the true posterior uncertainty. Such correlation is common, since the data couples the latents. The direction of the KL divergence, combined with the factorization, is responsible. Minimizing \(D_{\mathrm{KL}}(q \,\|\, p(\cdot \mid \boldsymbol{x}))\) penalizes configurations where \(q\) places mass in regions of low posterior probability much more strongly than configurations where \(q\) misses regions of high posterior probability, which produces the well-documented mode-seeking, variance-shrinking behavior of mean-field approximations.

These limitations motivate several directions of research. Structured mean-field approximations retain dependencies between blocks of latent variables and recover much of the correlation structure that fully-factorized mean-field discards. A hidden Markov sequence, for instance, can be treated as a single block while sequences are assumed independent. Richer fixed-form families for use with the gradient-based methods of the previous section include Gaussians with full or low-rank-plus-diagonal covariance, Gaussian mixtures, and the family of normalizing flows, in which a sequence of Jacobian-tractable diffeomorphisms transforms a simple base distribution into a highly expressive variational posterior.

Tighter bounds than the ELBO are available when one is willing to draw multiple samples per gradient step. The importance-weighted autoencoder bound (IWAE) and its descendants exploit \(K\)-sample importance averaging to reduce the gap to the log-evidence. Expectation propagation works with the reverse direction \(D_{\mathrm{KL}}(p(\cdot \mid \boldsymbol{x}) \,\|\, q)\), which it minimizes locally, one approximating factor at a time, rather than globally, exchanging the mode-seeking behavior of standard VI for a moment-matching, mass-covering behavior at the cost of a more delicate algorithmic implementation. We do not develop these directions here, since each is a substantial research area in its own right.

The broader role of variational inference in the contemporary toolkit is best understood by contrast with the alternatives. Markov chain Monte Carlo, treated on the Monte Carlo page, produces samples whose distribution converges to the true posterior, under suitable conditions on the chain, in the limit of infinite computation. It sacrifices exactness in finite time and pays in autocorrelated samples and slow mixing on multimodal posteriors. Variational inference, conversely, returns an approximation, typically much faster to compute, whose quality is bounded by the expressiveness of the variational family rather than by sampling diagnostics. Its price is asymptotic exactness, which it trades for tractability on problems where MCMC is prohibitive. The two are complementary tools, not competitors, and a working Bayesian practitioner draws on both.

Among gradient-based methods, the variational autoencoder page builds on the reparameterization trick and develops the amortized, deep-network specialization that is central to much of modern deep generative modeling, and the natural gradient page develops the geometry of parameter space that underlies many advanced variational optimizers. The variational foundation laid on this page, namely the ELBO, its KL decomposition, and the gradient and coordinate-ascent algorithms that operate on it, remains the structural backbone underneath these specializations.