Variational Bayesian methods

Variational Bayesian methods are a family of techniques for approximating intractable integrals arising in Bayesian inference and machine learning. They are typically used in complex statistical models consisting of observed variables as well as unknown parameters and latent variables, with various sorts of relationships among the three types of random variables, as might be described by a graphical model. As typical in Bayesian inference, the parameters and latent variables are grouped together as "unobserved variables". Variational Bayesian methods are primarily used for two purposes:

To provide an analytical approximation to the posterior probability of the unobserved variables, in order to do statistical inference over these variables.
To derive a lower bound for the marginal likelihood of the observed data. This is typically used for performing model selection, the general idea being that a higher marginal likelihood for a given model indicates a better fit of the data by that model and hence a greater probability that the model in question was the one that generated the data.

In the former purpose, variational Bayes is an alternative to Monte Carlo sampling methods—particularly, Markov chain Monte Carlo methods such as Gibbs sampling—for taking a fully Bayesian approach to statistical inference over complex distributions that are difficult to evaluate directly or sample. In particular, whereas Monte Carlo techniques provide a numerical approximation to the exact posterior using a set of samples, variational Bayes provides a locally-optimal, exact analytical solution to an approximation of the posterior.
Variational Bayes can be seen as an extension of the expectation–maximization algorithm from maximum likelihood or maximum a posteriori estimation of the single most probable value of each parameter to fully Bayesian estimation which computes the entire posterior distribution of the parameters and latent variables. As in EM, it finds a set of optimal parameter values, and it has the same alternating structure as does EM, based on a set of interlocked equations that cannot be solved analytically.
For many applications, variational Bayes produces solutions of comparable accuracy to Gibbs sampling at greater speed. However, deriving the set of equations used to update the parameters iteratively often requires a large amount of work compared with deriving the comparable Gibbs sampling equations. This is the case even for many models that are conceptually quite simple, as is demonstrated below in the case of a basic non-hierarchical model with only two parameters and no latent variables.

Mathematical derivation

Problem

In variational inference, the posterior distribution over a set of unobserved variables given some data is approximated by a so-called variational distribution,
The distribution is restricted to belong to a family of distributions of simpler form than , selected with the intention of making similar to the true posterior,.
The similarity is measured in terms of a dissimilarity function and hence inference is performed by selecting the distribution that minimizes.

KL divergence

The most common type of variational Bayes uses the Kullback–Leibler divergence of Q from P as the choice of dissimilarity function. This choice makes this minimization tractable. The KL-divergence is defined as
Note that Q and P are reversed from what one might expect. This use of reversed KL-divergence is conceptually similar to the expectation–maximization algorithm.

Intractability

Variational techniques are typically used to form an approximation for:
The marginalization over to calculate in the denominator is typically intractable, because, for example, the search space of is combinatorially large. Therefore, we seek an approximation, using.

Evidence lower bound

Given that, the KL-divergence above can also be written as
Because is a constant with respect to and because is a distribution, we have
which, according to the definition of expected value, can be written as follows
which can be rearranged to become
As the log-evidence is fixed with respect to, maximizing the final term minimizes the KL divergence of from. By appropriate choice of, becomes tractable to compute and to maximize. Hence we have both an analytical approximation for the posterior, and a lower bound for the log-evidence .
The lower bound is known as the variational free energy in analogy with thermodynamic free energy because it can also be expressed as a negative energy plus the entropy of. The term is also known as Evidence Lower Bound, abbreviated as ELBO, to emphasize that it is a lower bound on the log-evidence of the data.

Proofs

By the generalized Pythagorean theorem of Bregman divergence, of which KL-divergence is a special case, it can be shown that:
where is a convex set and the equality holds if:
In this case, the global minimizer with can be found as follows:
in which the normalizing constant is:
The term is often called the evidence lower bound in practice, since, as shown above.
By interchanging the roles of and we can iteratively compute the approximated and of the true model's marginals and respectively. Although this iterative scheme is guaranteed to converge monotonically, the converged is only a local minimizer of.
If the constrained space is confined within independent space, i.e. the above iterative scheme will become the so-called mean field approximation as shown below.

Mean field approximation

The variational distribution is usually assumed to factorize over some partition of the latent variables, i.e. for some partition of the latent variables into,
It can be shown using the calculus of variations that the "best" distribution for each of the factors satisfies:
where is the expectation of the logarithm of the joint probability of the data and latent variables, taken with respect to over all variables not in the partition: refer to Lemma 4.1 of for a derivation of the distribution.
In practice, we usually work in terms of logarithms, i.e.:
The constant in the above expression is related to the normalizing constant and is usually reinstated by inspection, as the rest of the expression can usually be recognized as being a known type of distribution.
Using the properties of expectations, the expression can usually be simplified into a function of the fixed hyperparameters of the prior distributions over the latent variables and of expectations of latent variables not in the current partition. This creates circular dependencies between the parameters of the distributions over variables in one partition and the expectations of variables in the other partitions. This naturally suggests an iterative algorithm, much like EM, in which the expectations of the latent variables are initialized in some fashion, and then the parameters of each distribution are computed in turn using the current values of the expectations, after which the expectation of the newly computed distribution is set appropriately according to the computed parameters. An algorithm of this sort is guaranteed to converge.
In other words, for each of the partitions of variables, by simplifying the expression for the distribution over the partition's variables and examining the distribution's functional dependency on the variables in question, the family of the distribution can usually be determined. The formula for the distribution's parameters will be expressed in terms of the prior distributions' hyperparameters, but also in terms of expectations of functions of variables in other partitions. Usually these expectations can be simplified into functions of expectations of the variables themselves ; sometimes expectations of squared variables, or expectations of higher powers also appear. In most cases, the other variables' distributions will be from known families, and the formulas for the relevant expectations can be looked up. However, those formulas depend on those distributions' parameters, which depend in turn on the expectations about other variables. The result is that the formulas for the parameters of each variable's distributions can be expressed as a series of equations with mutual, nonlinear dependencies among the variables. Usually, it is not possible to solve this system of equations directly. However, as described above, the dependencies suggest a simple iterative algorithm, which in most cases is guaranteed to converge. An example will make this process clearer.

A duality formula for variational inference

The following theorem is referred to as a duality formula for variational inference. It explains some important properties of the variational distributions used in variational Bayes methods.
Consider two probability spaces and with. Assume that there is a common dominating probability measure such that and. Let denote any real-valued random variable on that satisfies. Then the following equality holds
Further, the supremum on the right-hand side is attained if and only if it holds
almost surely with respect to probability measure, where and denote the Radon–Nikodym derivatives of the probability measures and with respect to, respectively.

A basic example

Consider a simple non-hierarchical Bayesian model consisting of a set of i.i.d. observations from a Gaussian distribution, with unknown mean and variance. In the following, we work through this model in great detail to illustrate the workings of the variational Bayes method.
For mathematical convenience, in the following example we work in terms of the precision — i.e. the reciprocal of the variance — rather than the variance itself.

The mathematical model

We place conjugate prior distributions on the unknown mean and precision, i.e. the mean also follows a Gaussian distribution while the precision follows a gamma distribution. In other words:
The hyperparameters and in the prior distributions are fixed, given values. They can be set to small positive numbers to give broad prior distributions indicating ignorance about the prior distributions of and.
We are given data points and our goal is to infer the posterior distribution of the parameters and