Fixed-Form Variational Posterior Approximation through Stochastic Linear Regression

Introduction

In Bayesian analysis the form of the posterior distribution is often not analytically tractable. To obtain quantities of interest under such a distribution, such as moments or marginal distributions, we typically need to use Monte Carlo methods or approximate the posterior with a more convenient distribution. A popular method of obtaining such an approximation is structured or fixed-form Variational Bayes, which works by numerically minimizing the Kullback-Leibler divergence of an approximating distribution in the exponential family to the intractable target distribution (Attias, 2000; Beal and Ghahramani, 2006; Jordan et al., 1999; Wainwright and Jordan, 2008). For certain problems, algorithms exist that can solve this optimization problem in much less time than it would take to approximate the posterior using Monte Carlo methods (see e.g. Honkela et al., 2010). However, these methods usually rely on analytic solutions to certain integrals and need conditional conjugacy in the model specification, i.e. the distribution of each variable conditional on its Markov blanket must be an analytically tractable member of the exponential family for these methods to be applicable. As a result this class of methods is limited in the type of approximations and posteriors they can handle.

We show that solving the optimization problem of fixed-form Variational Bayes is equivalent to performing a linear regression with the sufficient statistics of the approximation as explanatory variables and the (unnormalized) log posterior density as the dependent variable. Inspired by this result, we present an efficient stochastic approximation algorithm for solving this optimization problem. In contrast to earlier work, our approach does not require any analytic calculation of integrals, which allows us to extend the fixed-form Variational Bayes approach to problems where it was previously not applicable. Our method can be used to approximate any posterior distribution, provided that it is given in closed form up to the proportionality constant. The type of approximating distribution can be any distribution in the exponential family or any mixture of such distributions, which means that our approximations can in principle be made arbitrarily precise. While our method somewhat resembles performing stochastic gradient descent on the variational objective function in parameter space (Paisley et al., 2012; Nott et al., 2012), the linear regression view gives insights which allow a more computationally efficient approach.

Section 2 introduces fixed-form variational posterior approximation, the optimization problem to be solved, and the notation used in the remainder of the paper. In Section 3 we provide a new way of looking at variational posterior approximation by re-interpreting the underlying optimization as a linear regression problem. We propose a stochastic approximation algorithm to perform the optimization in Section 4. In Section 5 we discuss how to assess the quality of our posterior approximations and how to use the proposed methods to approximate the marginal likelihood of a model. These sections represent the core ideas of the paper.

To make our approach more generally applicable and computationally efficient we provide a number of extensions in two separate sections. Section 6 discusses modifications of our stochastic approximation algorithm to improve efficiency. Section 7 relaxes the assumption that our posterior approximation is in the exponential family, allowing instead mixtures of exponential family distributions. Sections 4, 6, and 7 also contain multiple examples of using our method in practice, and show that despite its generality, the efficiency of our algorithm is highly competitive with more specialized approaches. Code for these examples is available at github.com/TimSalimans/LinRegVB. Finally, Section 8 concludes.

Fixed-form Variational Bayes

Let $x$ be a vector of unknown parameters and/or latent random effects for which we have specified a prior distribution $p(x)$, and let $p(y|x)$ be the likelihood of observing a given set of data, $y$. Upon observing $y$ we can use Bayes’ rule to obtain our updated state of belief, the posterior distribution

An equivalent definition of the posterior distribution is

where the optimization is over all proper probability distributions $q(x)$, and where $D[q(x)|p(x|y)]$ denotes the Kullback-Leibler divergence between $q(x)$ and $p(x|y)$. The KL-divergence is always non-negative and has a unique minimizing solution $q(x)=p(x|y)$ almost everywhere, at which point the KL-divergence is zero. The solution of (2) does not depend on the normalizing constant $p(y)$ of the posterior distribution.

The posterior distribution given in (1) is the exact solution of the variational optimization problem in (2), but except for certain special cases it is not very useful by itself because it does not have an analytically tractable form. This means that we do not have analytic expressions for the posterior moments of $x$, the marginals $p(x_{i}|y)$, or the normalizing constant $p(y)$. One method of solving this problem is to approximate these quantities using Monte Carlo simulation. A different approach is to restrict the optimization problem in (2) to a reduced set of more convenient distributions $Q$. If $p(x,y)$ is of conjugate exponential form, choosing $Q$ to be the set of factorized distributions $q(x)=q(x_{1})q(x_{2})\dots q(x_{k})$ often leads to a tractable optimization problem that can be solved efficiently using an algorithm called Variational Bayes Expectation Maximization (VBEM, Beal and Ghahramani, 2002). Such a factorized solution is attractive because it makes the variational optimization problem easy to solve, but it is also very restrictive: it requires a conjugate exponential model and prior specification and it assumes posterior independence between the different blocks of parameters $x_{i}$. This means that this factorized approach can be used with few models, and that the solution $q(x)$ may be a poor approximation to the exact posterior (see e.g. Turner et al., 2008).

An alternative choice for $Q$ is the set of distributions of a certain parametric form $q_{\eta}(x)$, where $\eta$ denotes the vector of parameters governing the shape of the posterior approximation. This approach is known as structured or fixed-form Variational Bayes (Honkela et al., 2010; Storkey, 2000; Saul and Jordan, 1996). Usually, the posterior approximation is chosen to be a specific member of the exponential family of distributions:

where $T(x)$ is a $1\times k$ vector of sufficient statistics, $U(\eta)$ takes care of normalization, and $\nu(x)$ is a base measure. The $k\times 1$ vector $\eta$ is often called the set of natural parameters of the exponential family distribution $q_{\eta}(x)$. Using this approach, the variational optimization problem in (2) reduces to a parametric optimization problem in $\eta$:

Variational Bayes as linear regression

For notational convenience we will write our posterior approximation in the adjusted form,

Setting this expression to zero in order to find the minimum gives

For notational simplicity, we will assume a constant base measure $\nu(x)=1$ in the remaining discussion, but the linear regression analogy continues to hold if the base measure $\nu(x)$ is non-constant in $x$. In that case, the fixed point condition (9) simply becomes

i.e. we perform the linear regression on the residual of the base model $\log\nu(x)$.

A stochastic approximation algorithm

If the algorithm is run for additional iterations after the true posterior is recovered, the approximation will not change. This is to be contrasted with other stochastic gradient descent algorithms which have non-vanishing variance for a finite number of samples, and is due to the fact that our regression in itself is noise free: only its support points are stochastic. This exact convergence will not hold for cases of actual interest, where $p$ and $q$ will not be of the exact same functional form, but we generally still observe a dramatic improvement when using Algorithm 1 instead of more conventional stochastic gradient descent algorithms. A deeper analysis of the variance of our stochastic approximation is given in Appendix D.

Contrary to most applications in the literature, Algorithm 1 uses a fixed step size $w=1/\sqrt{N}$ rather than a declining one in updating our statistics. The analyses of Robbins and Monro (1951) and Amari (1997) show that a sequence of learning rates $w_{t}=ct^{-1}$ is asymptotically efficient in stochastic gradient descent as the number of iterations $N$ goes to infinity, but this conclusion rests on strong assumptions on the functional form of the objective function (e.g. strong convexity) that are not satisfied for the problems we are interested in. Moreover, with a finite number of iterations, the effectiveness of a sequence of learning rates that decays this fast is highly dependent on the proportionality constant $c$. If we choose $c$ either too low or too high, it may take a very long time to reach the efficient asymptotic regime of this learning rate sequence.

Nemirovski et al. (2009) show that a more robust approach is to use a constant learning rate $w=1/\sqrt{N}$ and that this is optimal for finite $N$ without putting stringent requirements on the objective function. In order to reduce the variance of the last iterate with this non-vanishing learning rate, they propose to use an average of the last $L$ iterates as the final output of the optimization algorithm. The value of $L$ should grow with the total number of iterations, and is usually chosen to be equal to $N/2$. Remarkably, they show that such an averaging procedure can match the asymptotic efficiency of the optimal learning sequence $w_{t}=ct^{-1}$.

Like other optimization algorithms for Variational Bayes, Algorithm 1 will only find a local minimum of the KL-divergence. This is generally not a problem when approximating unimodal posterior distributions, such as with the examples in this paper, since the optimization problem then often only has a single optimum (depending on the type of approximation, see Bishop, 2006, Ch. 10). If the true posterior distribution is multimodal and the approximation is unimodal, however, the variational approximation will tend to pick one of the posterior modes and ignore the others (Minka, 2005). Although this is often desirable (see e.g. Stern et al., 2009), there is no guarantee that the recovered local minimum of the KL-divergence is then also a global minimum.

These seemingly similar alternatives perform dramatically worse than Algorithm 1. We set the true $\lambda:=2$, and initialize $\eta:=1$ and $C:=I_{2}$, the identity matrix. Figure 1 shows the mean and variance of the estimates of $\log(\eta)$ across $100$ repeat runs of each method with varying number of iterations $N$. We see it takes option i (“different randomness”) and ii (“analytic”) well over $1000$ iterations to give a reasonable answer, and even with $N=10^{4}$ samples, option i) estimates $\hat{\eta}=2.04\pm 0.15$ and option ii) $2.01\pm 0.11$.

Marginal likelihood and approximation quality

The stochastic approximation algorithm presented in the last section serves to minimize the Kullback-Leibler divergence between $q_{\eta}(x)$ and $p(x|y)$, given by

As discussed before, we do not need to know $p(y)$ (the marginal likelihood) in order to minimize $D(q_{\eta}|p)$ as $p(y)$ does not depend on $\eta$, but we do need to know it if we want to determine the quality of the approximation, as measured by the final Kullback-Leibler divergence. In addition, the constant $p(y)$ is also essential for performing Bayesian model comparison or model averaging. This section presents a method for approximating the marginal likelihood and final Kullback-Leibler divergence.

When our algorithm has converged, we have the following identity

where $r(x)$ is the ‘residual’ or ‘error term’ in the linear regression of $\log p(x,y)$ on the sufficient statistics of $q_{\eta}(x)$, and where $U(\eta)$ is the normalizer of $q_{\eta}(x)$. The intercept of the regression is

which we need to integrate with respect to $x$ in order to find the marginal likelihood $p(y)$. Doing so gives

so that $\hat{\eta}_{0}+U(\eta)$ is indeed a lower bound on the log marginal likelihood. If our approximation is perfect, the KL-divergence is zero and $r(x)$ is zero almost everywhere. In that case the residual term vanishes and the lower bound will be tight, otherwise it will underestimate the true marginal likelihood. The lower bound $\hat{\eta}_{0}+U(\eta)$ is often used in model comparison, which works well if the KL-divergence between the approximate and true posterior distribution is of approximately the same size for all models that are being compared. However, if we compare two very different models this will often not be the case, and the model comparison will be biased as a result. In addition, as opposed to the exact marginal likelihood, the lower bound gives us no information on the quality of our posterior approximation. It would therefore be useful to obtain a better estimate of the marginal likelihood.

One approach to doing this would be to evaluate the expectation in (18) using Monte Carlo sampling. Some analysis shows that this corresponds to approximating $p(y)$ using importance sampling, with $q_{\eta}(x)$ as the candidate distribution. It is well known that this estimator of the marginal likelihood may have infinite variance, unless $r(x)$ is bounded from above (see e.g. Geweke, 2005, p. 114). In general, we cannot guarantee the boundedness of $r(x)$ for our approach, so we will instead approximate the expectation in (18) using something that is easier to calculate.

At convergence, we know that the mean of $r(x)$ is zero when sampling from $q_{\eta}(x)$. The variance of $r(x)$ can be estimated using the mean squared error of the regressions we perform during the optimization, with relatively low variance. We denote our estimate of this variance by $s^{2}$. The assumption we will then make in order to approximate $\log p(y)$ is that $r(x)$ is approximately distributed as a normal random variable with these two moments. This leads to the following simple estimate of the log marginal likelihood

That is, our estimate of the marginal likelihood is equal to its lower bound plus a correction term that captures the error in our posterior approximation $q_{\eta}(x)$. Similarly, we can approximate the KL-divergence of our posterior approximation as

The KL-divergence is approximately equal to half the mean squared error in the regression of $\log p(x,y)$ on the sufficient statistics of the approximation. This relationship should not come as a surprise: this mean squared error is exactly what we minimize when we perform linear regression.

The scale of the KL-divergence is highly dependent on the amount of curvature in $\log p(x|y)$ and is therefore not easily comparable across different problems. If we scale the approximate KL-divergence to account for this curvature, this naturally leads to the R-squared measure of fit for regression modeling:

The R-squared measure corrects for the amount of curvature in the posterior distribution and is therefore comparable across different models and data sets. In addition it is a well-known measure and easily interpretable. We therefore propose to use the R-squared as the measure of approximation quality for our variational posterior approximations. Although we find the R-squared to be a useful measure for the majority of applications, it is important to realize that it mostly contains information about the mass of the posterior distribution and its approximation, and not directly about their moments. It is therefore possible to construct pathological examples in which the R-squared is relatively high, yet the (higher) moments of the posterior and its approximation are quite different. This may for example occur if the posterior distribution has very fat tails.

Section 7.2.1 provides an application of the methods developed here. In that section, Figure 6 shows that the approximation of the KL-divergence is quite accurate, especially when the approximation $q_{\eta}(x)$ is reasonably good. The same figure also shows that the approximation of the marginal likelihood proposed here (19) is much more accurate than the usual lower bound. In Sections 6 and 7, we also calculate the R-squared measure of approximation quality for a number of different posterior approximations, and we conclude that it corresponds well to visual assessments of the approximation accuracy.

The discussion up to this point represents the core ideas of this paper. To make our approach more general and computationally efficient we now provide a number of extensions in two separate sections. Section 6 discusses modifications of our stochastic approximation algorithm to improve efficiency, and Section 7 generalizes the exponential family approximations $q(x)$ used so far to include mixtures of exponential family distributions. Examples are given throughout. Finally, Section 8 concludes.

Extensions I: Improving algorithmic efficiency

We use data simulated from the model, with $N=100$ and $M=5$, to be able to show the performance averaged over $500$ datasets and many different settings of the algorithm. We compare our algorithm to the VBEM algorithm of Ormerod and Wand (2010) which makes use of the fact that the expectations required for this model can be calculated analytically. We choose not to do this for our method to investigate how effective our MC estimation strategy can be. For completeness we also compare to variational message passing (VMP, Winn and Bishop, 2006), a message passing implementation of VBEM, and expectation propagation (EP, Minka, 2001), which is known to have excellent performance on binary classification problems (Nickisch and Rasmussen, 2008). These last two alternatives are both implemented in Infer.NET (Minka et al., 2010) a library for probabilistic inference in graphical models, whereas we implement VBEM and our approximation algorithm ourselves in MATLAB. VMP and VBEM use a different variational approximation to our methods, introducing auxiliary variables $z_{i}\sim N(x^{\prime}v_{i},1)$, with $z_{i}$ constrained to be positive if $y_{i}=1$ and negative otherwise. A factorized variational posterior $q(x)\prod_{i}q(z_{i})$ is used, where $q(x)$ is multivariate normal and each $q(z_{i})$ can be thought of as a truncated univariate Gaussian.

For all implementations of our algorithm, we initialize the posterior approximation to the prior. All algorithms then use a single random draw to update the parameters during each iteration. This is often not the best implementation in terms of computational efficiency, since the contributions of multiple draws can often be calculated in parallel at little extra cost, and using antithetic sampling (i.e. sampling of negatively correlated draws) can reduce the variance of our approximations. By using the most basic implementation, however, we can more clearly compare the different stochastic approximations proposed in this section. Since the time required to run the different algorithms is strongly dependent on their precise implementation (e.g. the chosen programming language), we choose to perform this comparison by looking at statistical efficiency, as measured by the accuracy as a function of the number of likelihood evaluations, rather than the running time of the algorithms.

Since this experiment is on synthetic data we are able to assess performance in terms of the method’s ability to recover the known regression coefficients $x$, which we quantify as the root mean squared error (RMSE) between the variational mean and the true regression weights, and the “log score”: the log density of the true weights under the approximate variational posterior. The log score is useful because it rewards a method for finding good estimates of the posterior variance as well as the mean, which should of course be central to any approximate Bayesian method.

Figure 2 shows the performance of the different versions of our algorithm as presented in the following discussion, as well as the performance of the VBEM algorithm of Ormerod and Wand (2010). As can be seen from this graph, our approximation method achieves a lower RMSE than the VBEM algorithm. This is because of the extra factorization assumptions made by VBEM when introducing the $z_{i}$ variables. Where the improvement in the RMSE is noticeable, the difference in log score is dramatic: $0.193$ versus $-4.46$ (not shown), indicating that our approximation gives significantly better estimates of the variance than VBEM. The average R-squared obtained by our variational approximation was $0.97$, indicating a close fit to the exact posterior distribution. In terms of accuracy, our results are very similar to those of EP, which obtained an RMSE and log score identical to those of our approximation (up to 3 significant digits). As expected, VMP gave consistent results with VBEM: an RMSE of $0.265$ and a log score of $-4.56$.

As can be seen from Figure 2, our basic algorithm is considerably slower than VBEM in terms of the number of likelihood evaluations that are required to achieve convergence. In terms of wall clock time, our basic algorithm ran about an order of magnitude slower than VBEM, although it could easily be sped up by using multiple random draws in parallel. The basic algorithm was about as fast as EP and VMP, needing about 15 milliseconds to converge on this small data set, but note that the system set ups were not completely comparable: EP and VMP were run on a laptop rather than a desktop, and Infer.NET is implemented in C# rather than MATLAB.

The remainder of this section introduces the other implementations of our variational approximation, presented in Figure 2, some of which are much faster and more computationally efficient than both our basic algorithm and VBEM.

1 Making use of factor structure

For most statistical problems, including our probit regression model, the log posterior can be decomposed into a number of additive factors, i.e. $\log p(x,y)=\sum_{j=1}^{N}\log\phi_{j}(x,y)$. The optimality condition in (9) can then also be written as a sum:

This means that rather than performing one single linear regression we can equivalently perform $N$ separate regressions.

One benefit of this is that some of the factors $\phi_{j}(x,y)$ may be conjugate to the posterior approximation, such as the prior $p(x)$ in our probit regression example. The regression coefficients $\hat{\eta}^{j}$ for these conjugate factors are known analytically and do not need to be approximated.

Our probit regression model provides a straightforward example, for which the log joint density of $x$ and $y$ has the following factor structure

Here, each likelihood factor $p(y_{i}|v_{i},x)$ depends on all the parameters $x$, but only through the univariate product $f_{i}=x^{\prime}v_{i}$. We can emphasize this by writing our model as

and regressing these against the likelihood factors $\log p(y_{i}|v_{i},x_{i})$. At each iteration of Algorithm 1, we can then update the natural parameters of each approximate likelihood term using

Rather than performing a single regression of dimension $1+M(M+3)/2$, we may thus equivalently perform $N$ regressions of dimension 3. Performing these lower dimensional regressions is computationally more efficient as long as $N$ is not very large, and it is also statistically more efficient. Figure 2 shows that this factorized regression implementation of our approximation indeed needs far fewer random draws to achieve convergence. All $N$ regressions can be performed in parallel, which offers further opportunities for computational gain on multicore machines or computer clusters.

So far, we have assumed that we sample $x^{*}$ and then form the $f_{i}$ by multiplying with the $v_{i}$, but note that we can equivalently sample the $f_{i}$ directly and separately from their univariate Gaussian approximate posteriors $q_{\eta}(f_{i})=N[\mu_{i}(\eta,v_{i}),\sigma^{2}_{i}(\eta,v_{i})]$. For the current example we find that both implementations are about equally efficient.

2 Using the gradient of the log posterior

Furthermore, using the properties of the exponential family of distributions, we know that

By using the same random number seed $z^{*}$ in both Monte Carlo approximations we once again get the beneficial variance reduction effect described in Section 4.

Performing a single iteration using (27) provides about the same information as doing $2\times\dim(x)$ iterations with the basic algorithm, making it more computationally efficient if the gradients can be obtained analytically.

We can also do updates of this form while still making use of the factor structure of the posterior distribution, as proposed above for the probit regression example. Using this example, and assuming we sample the $f_{i}$ separately (see last paragraph of Section 6.1), this gives the following regression statistics for each of the $N$ low dimensional regressions:

Figure 2 shows the performance of this approximation on our probit example, showing again a large gain in efficiency with respect to the approximations introduced earlier. Empirically, we find that using gradients also leads to more efficient stochastic optimization algorithms for many other applications. For some problems the posterior distribution will not be differentiable in some of the elements of $x$, for example when $x$ is discrete. In that case the stochastic approximations presented here may be combined with the basic approximation of Section 4.

In addition, for many samplers $\nabla_{\eta}s(\eta,z^{*})$ may be not defined, e.g. rejection samplers. However, for the gradient approximations it does not matter what type of sampler is actually used to draw $x^{*}$, only that it is from the correct distribution. A correct strategy is therefore to draw $x^{*}$ using any desired sampling algorithm, and then proceeding as if we had used a different sampling algorithm for which $\nabla_{\eta}s(\eta,z^{*})$ is defined. For example, we might use a nondifferentiable rejection sampler to draw a univariate $x^{*}$, and then calculate (27) as if we had used an inverse-transform sampler, for which we have

for all natural parameters $\eta_{i}$, with $Q_{\eta}(x)$ the cdf and $q_{\eta}(x)$ the pdf of $x$. Similarly, it does not matter for the probit example whether we sample the $f_{i}$ jointly by sampling $x$, or whether we sample them directly and independently. After sampling the $f_{i}$, we can use $s_{i}(\eta,z_{i}^{*})=\mu_{i}+\sigma_{i}z_{i}^{*}$ as proposed above, but we might equivalently proceed using (42), or something else entirely. Finding the most efficient strategy we mostly leave for future work, although Sections 6.3 and 6.4 offer some further insights into what is possible.

3 Using the Hessian of the log posterior

When we have both first and second order gradient information for $\log p(x,y)$ and if we choose our approximation to be multivariate Gaussian, i.e. $q_{\eta}(x)=N(m(\eta),V(\eta))$, we have a third option for approximating the statistics used in the regression. For Gaussian $q(x)$ and twice differentiable $\log p(x,y)$, Minka (2001) and Opper and Archambeau (2009) show that

where $\nabla_{x}\nabla_{x}\log p(x,y)$ denotes the Hessian matrix of $\log p(x,y)$ in $x$.

For the multivariate Gaussian distribution we know that the natural parameters are given as $\eta_{1}=V^{-1}m$ and $\eta_{2}=V^{-1}$. Using this relationship, we can derive Monte Carlo estimators $\hat{g}$ and $\hat{C}$ using the identities (25, 26). We find that these stochastic approximations are often even more efficient than the ones in Section 6.2, provided that the Hessian matrix of $\log p(x,y)$ can be calculated cheaply. This type of approximation is especially powerful when combined with the extension presented in the next section.

4 Linear transformations of the regression problem

It is well known that classical linear least squares regression is invariant to invertible linear transformations of the explanatory variables. We can use the same principle in our stochastic approximation algorithm to allow us to work with alternative parameterizations of the approximate posterior $q(x)$. These alternative forms can be easier to implement or lead to more efficient algorithms, as we show in this section.

Using the same principle, we can rewrite the optimality condition of (9) as

for any invertible matrix $K$, which may depend on the variational parameters $\eta$. Instead of solving our original least squares regression problem, we may thus equivalently solve this transformed version. When we perform the linear regression in (45) for a fixed set of parameters $\eta$, the result will be identical to that of the original regression with $K(\eta)=\operatorname{I}$, as long as we use the same random numbers for both regressions. However, when the Monte Carlo samples (‘data points’ in our regression) are generated using different values of $\eta$, as is the case with the proposed stochastic approximation algorithm, the two regressions will not necessarily give the same solution for a finite number of samples. If the true posterior $p(x|y)$ is of the same functional form as the approximation $q_{\eta}$, the exact convergence result of Section 4 holds for any invertible $K(\eta)$, so it is not immediately obvious which $K(\eta)$ is best for general applications.

We hypothesize that certain choices of $K(\eta)$ may lead to statistically more efficient stochastic approximation algorithms for certain specific problems, but we will not pursue this idea here. What we will discuss is the observation that the stochastic approximation algorithm may be easier to implement for some choices of $K(\eta)$ than for others, and that the computational costs are not identical for all $K(\eta)$. In particular, the transformation $K(\eta)$ allows us to use different parameterizations of the variational approximation. Let $q_{\phi}$ be such a reparameterization of the approximation, let the new parameter vector $\phi(\eta)$ be an invertible and differentiable transformation of the original parameters $\eta$, and set $K(\eta)$ equal to the inverse transposed Jacobian of this transformation, i.e. $K(\eta)=[\nabla_{\eta}\phi(\eta)]^{-1}$. Using the properties of the exponential family of distributions, we can then show that

for any differentiable function $h(x)$. Using this result, the stochastic approximations of Section 6.2 for the transformed regression problem are

These new expressions for $\hat{g}$ and $\hat{C}$ may be easier to calculate than the original ones (27), and the resulting $\hat{C}$ may have a structure making it easier to invert in some cases. An example of this occurs when we use a Gaussian approximation in combination with the stochastic approximations of Section 6.3, using the gradient and Hessian of $\log p(x,y)$. In this case we may work in the usual natural parameterization, but doing so gives a dense matrix $\hat{C}$ with dimensions proportional to $M^{2}$, where $M$ is the dimension of $x$. For large $M$, such a stochastic approximation is expensive to store and invert. However, using the stochastic approximations above, we may alternatively parameterize our approximation in terms of the mean $m$ and variance $V$. Working in this parameterization, we can express the update equations for the natural parameters in terms of the gradient and Hessian of $\log p(x,y)$ and the average sampled $x$ value, instead of the (higher dimensional) $g$ and $C$ statistics. The resulting algorithm, as derived in Appendix B, is therefore more efficient in terms of both computation and storage. Pseudocode for the new algorithm is given below.

Instead of storing and inverting the full $C$ matrix, this algorithm uses the sparsity induced by the transformation $K(\eta)$ to work with the precision matrix $P$ instead. The dimensions of this matrix are equal to the dimension of $x$, rather than its square, providing great savings. Moreover, while the $C$ matrix in the original parameterization is always dense, $P$ will have the same sparsity pattern as the Hessian of $\log p(x,y)$, which may reduce the costs of storing and inverting it even further for many applications.

Figure 2 shows the performance of Algorithm 2 as applied to our probit regression example. As is typical for this version of the algorithm, it performs even better than the algorithm using only the gradient and factor structure of the posterior distribution. Since this type of approximation is also very easy to implement efficiently in a matrix programming language like MATLAB, it also runs significantly faster than the VBEM algorithm for this example. Moreover, the algorithm is now again completely general and does not make any assumptions as to the structure of the posterior distribution (other than it being twice differentiable). This means it can easily be used for Gaussian variational approximation of almost any posterior distribution.

5 Subsampling the data: double stochastic approximation

The stochastic approximations derived above are all linear functions of $\log p(x,y)$ and its first and second derivatives. This means that these estimates are still unbiased even if we take $\log p(x,y)$ to be a noisy unbiased estimate of the true log posterior, rather than the exact log posterior. For most statistical applications $\log p(x,y)$ itself is a separable additive function of a number of independent factors, i.e. $\log p(x,y)=\sum_{j=1}^{N}\log\phi_{j}(x,y)$ as explained in Section 6.1. Using this fact we can construct an unbiased stochastic approximation of $\log p(x,y)$ as

For our probit regression example we implement subsampling by dividing the sample into 10 equally sized ‘minibatches’ of data. During each iteration of the algorithm, these minibatches are processed in random order, using Algorithm 2 combined with (50) to update the variational parameters after each minibatch. As can be seen in Figure 2 this approach allows us to get a good approximation to the posterior very quickly: reaching the accuracy of converged VBEM now only requires three passes over the training data, although final convergence is not much faster than when using the full sample.

Extensions II: Using mixtures of exponential family distributions

So far, we have assumed that the approximating distribution $q_{\eta}(x)$ is a member of the exponential family. Here we will relax that assumption. If we choose a non-standard approximation, certain moments or marginals of $q_{\eta}(x)$ are typically no longer available analytically, which should be taken into account when choosing the type of approximation. However, if we can at least sample directly from $q_{\eta}(x)$, it is often still much cheaper to approximate these moments using Monte Carlo than it would be to approximate the corresponding moments of the posterior using MCMC or other indirect sampling methods. We have identified two general strategies for constructing useful non-standard posterior approximations which are discussed in the following two sections.

If we split our vector of unknown parameters $x$ into $p$ non-overlapping blocks, our approximating posterior may be decomposed as

If we then choose every conditional posterior $q(x_{i}|x_{1},\ldots,x_{i-1})$ to be an analytically tractable member of the exponential family, we can easily sample from the joint $q(x)$, while still having much more freedom in capturing the dependence between the different blocks of $x$. In practice, such a conditionally tractable approximation can be achieved by specifying the sufficient statistics of each exponential family block $q(x_{i}|x_{1},\ldots,x_{i-1})$ to be a function of the preceding elements $x_{1},x_{2},\ldots,x_{i-1}$. This leads to a natural type of approximation for hierarchical Bayesian models, where the hierarchical structure of the prior often suggests a good hierarchical structure for the posterior approximation.

If every conditional $q(x_{i}|x_{1},\ldots,x_{i-1})$ is in the exponential family, the joint may not be if the normalizing constant of any of those conditionals is a non-separable function of the preceding elements $x_{1},x_{2},\ldots,x_{i-1}$ and the variational parameters. However, because the conditionals are still in the exponential family, our optimality condition still holds separately for the variational parameters of each conditional with only slight modification. Taking again the derivative of the KL-divergence and setting it to zero yields:

where $T_{i}(x_{i})$ and $\eta_{i}$ denote the sufficient statistics and corresponding natural parameters of the $i$-th conditional approximation $q(x_{i}|x_{1},\ldots,x_{i-1})$, and where $r_{-i}(x)$ can be seen as the residual of the approximation with the $i$-th block left out. Note that we cannot rewrite this expression as a linear regression any further, like we did in Section 2, since the intercept of such a regression is related to the normalizing constant of $q(x_{i}|x_{1},\ldots,x_{i-1})$ which may now vary in $x_{1},\ldots,x_{i-1}$. However, $C_{i}$ and $g_{i}$ can still be approximated straightforwardly using Monte Carlo, and Algorithm 1 can still be used with these approximations, performing separate ‘regressions’ for all conditionals during each iteration like we proposed for factorized $p(x,y)$ in Section 6.1. Alternatively, Algorithm 2 or any of the extensions in Section 6 may be used to fit the different blocks of $q_{\eta}(x)$.

Using this type of approximation, the marginals $q(x_{i})$ will generally be mixtures of exponential family distributions, which is where the added flexibility of this method comes from. By allowing the marginals $q(x_{i})$ to be mixtures with dependency on the preceding elements of $x$, we can achieve much better approximation quality than by forcing them to be a single exponential family distribution. A similar idea was used in the context of importance sampling by Hoogerheide et al. (2012). A practical example of this is given below.

Stochastic volatility models for signals with time varying variances are considered extremely important in finance. Here we apply our methodology to the model and prior specified in Girolami and Calderhead (2011). The data we will use, from Kim et al. (1998), is the percentage change $y_{t}$ in GB Pound vs. US Dollar exchange rate, modeled as:

The relative volatilities, $v_{t}$ are governed by the autoregressive AR(1) process

The distributions of the error terms are given by $\epsilon_{t}\sim N(0,1)$ and $\xi_{t}\sim N(0,\sigma^{2})$. The prior specification is as in Girolami and Calderhead (2011):

Following the strategy outlined above, we use the hierarchical structure of the prior to suggest a hierarchical structure for the approximate posterior:

The prior of $\phi$ is in the exponential family, so we choose the posterior approximation $q_{\eta}(\phi)$ to be of the same form:

The prior for $\sigma^{2}$ is inverse-Gamma, which is also in the exponential family. We again choose the same functional form for the posterior approximation, but with a slight modification in order to capture the posterior dependency between $\phi$ and $\sigma^{2}$:

where the extra term $\eta_{5}\phi^{2}$ was chosen by examining the functional form of the exact full conditional $p(\sigma^{2}|\phi,v)$.

Using the notation $f=(\log(\beta),v^{\prime})^{\prime}$, the conditional prior $p(f|\phi,\sigma^{2})$ can be seen as the diffuse limit of a multivariate normal distribution. We therefore also use a multivariate normal conditional approximate posterior:

with $p(f|\phi,\sigma^{2})$ the Gaussian prior, $q_{\eta}(y|f)$ a Gaussian approximate likelihood of the form

with $\eta_{6}$ a $T\times T$ positive-definite matrix and $\eta_{7}$ a $T\times 1$ vector, and where

is the normalizing constant of our posterior approximation $q_{\eta}(f|\phi,\sigma^{2})$.

Now that we have defined the functional form of the approximate posterior, we can fit its parameters by applying (51) to each of the blocks $q_{\eta}(\phi)$, $q_{\eta}(\sigma^{2}|\phi)$, and $q_{\eta}(f|\phi,\sigma^{2})$. We approximate the statistics of the first two blocks using gradients as proposed in Section 6.2. The last (multivariate Gaussian) block is updated using both the gradient and the Hessian of $p(y|f)$ via the optimized expressions of Algorithm 2.

For the first block $q_{\eta}(\phi)$ this gives us the following stochastic approximations:

where $T_{1}(\phi^{*})$ are the sufficient statistics of $q_{\eta}(\phi)$, and where we make use of the fact that

Cancelling the prior term $p(f|\phi,\sigma^{2})$ in $p()$ and $q()$ then allows us to go from (7.1.1) to (7.1.1). The approximate marginal likelihood $q_{\eta}(y|\phi,\sigma^{2})$ and the expectations with respect to $q_{\eta}(f|\phi,\sigma^{2})$ can be evaluated analytically using the Kalman filter and smoother (e.g. Durbin and Koopman, 2001), which means we do not have to sample $f$ for this problem. Note that (7.1.1) includes both the direct effect of $\phi$, as well as its indirect effects through $q_{\eta}(\sigma^{2}|\phi)$ and $q_{\eta}(f|\phi,\sigma^{2})$. If the functional form of $q()$ is close to that of $p()$, the relative importance of these indirect effects is low. In most cases we can therefore ignore these indirect effects with little to no loss of accuracy. For the current application we find that using (57) instead of (7.1.1) gives virtually identical results.

The stochastic approximations for the second block $q_{\eta}(\sigma^{2}|\phi)$ are given by

where $T_{2}(\sigma^{2*})$ are the sufficient statistics of $q_{\eta}(\sigma^{2}|\phi)$.

Finally, the updates for the likelihood approximation (using Algorithm 2) are given by

Here again, the expectations with respect to the approximate posterior $q_{\eta}(f|\phi,\sigma^{2})$ can be calculated analytically using the Kalman filter/smoother and do not have to be approximated by sampling. Furthermore we know that the Hessian of the log likelihood is sparse, which means that only a relatively small number of the parameters in $\eta_{6}$ will be non-zero: all elements on the diagonal and all elements in the column and row belonging to $\log(\beta)$. This sparsity is also what makes fitting this posterior approximation feasible, since inverting a dense $T\times T$ precision matrix would be much too expensive. Even with this sparsity, our optimization problem is still fairly high dimensional with about 2000 free parameters. Nevertheless, we find that our approximation converges very quickly using 250 iterations of our algorithm, with a single $(\phi,\sigma^{2})$ sample per iteration, which takes our single-threaded MATLAB implementation half a second to complete on a 3GHz processor. This is more than two orders of magnitude faster than the running time required by advanced MCMC algorithms for this problem.

We compare the results of our posterior approximation against the “true” posterior, provided by a very long run of the MCMC algorithm of Girolami and Calderhead (2011). As can be seen from Figures 3, 4 and 5, the posterior approximations for the model parameters are nearly exact. Similarly, the posterior approximations for the latent volatilities $v$ (not shown) are also indistinguishable from the exact posterior.

Our approach to doing inference in the stochastic volatility model shares some characteristics with the approach of Liesenfeld and Richard (2008). They fit a Gaussian approximation to the posterior of the volatilities for given $\phi,\sigma^{2},\beta$ parameters, using the importance sampling algorithm of Richard and Zhang (2007), which is based on auxiliary regressions somewhat similar to those in Algorithm 1. They then infer the model parameters using MCMC methods. The advantage of our method is that we are able to leverage the information in the gradient and Hessian of the posterior, and that our stochastic approximation algorithm allows us to fit the posterior approximation very quickly for all volatilities simultaneously, while their approach requires optimizing the approximation one volatility at a time. Unique to our approach is also the ability to concurrently fit a posterior approximation for the model parameters $\phi,\sigma^{2},\beta$ and have the approximate posterior of the volatilities depend on these parameters, while Liesenfeld and Richard (2008) need to re-construct their approximation every time a new set of model parameters is considered. As a result, our approach is significantly faster for this problem.

2 Using auxiliary variables

Another approach to constructing flexible posterior approximations is using the conditional exponential family approximation of Section 7.1, but letting the first block of variables be a vector of auxiliary variables $u$, that are not part of the original set of model parameters and latent variables, $x$. The posterior approximation then has the form

The factors $q(u)$ and $q(x|u)$ should both be analytically tractable members of the exponential family, which allows the marginal approximation $q(x)$ to be a general mixture of exponential family distributions, like a mixture of normals for example. If we use enough mixture components, the approximation $q(x)$ could then in principle be made arbitrarily close to $p(x|y)$. Note, however, that if $p(x|y)$ is multimodal our optimization problem might suffer from multiple local minima, which means that we are generally not guaranteed to find the optimal approximation.

The mixture approximation $q(x)$ can be fitted by performing the standard KL-divergence minimization:

From (59) it becomes clear that an additional requirement of this type of approximation is that we can integrate out the auxiliary variables $u$ from the joint $q(x,u)$ in order to evaluate the marginal density $q(x)$ at a given point $x$. Fortunately this is easy to do for many interesting approximations, such as discrete mixtures of normals or continuous mixtures like Student’s t distributions. Also apparent from (59) is that we cannot use this approximation directly with the stochastic approximation algorithms proposed in the last sections since $q(x)$ is itself not part of the exponential family of distributions. However, we can rewrite (59) as

Equation 60 now once again has the usual form of a KL-divergence minimization where the approximation, $q_{\eta}(x,u)$, consists of exponential family blocks $q_{\eta}(u)$ and $q_{\eta}(x|u)$. By including the auxiliary variables $u$ in the ‘true’ posterior density, we can thus once again make use of our efficient stochastic optimization algorithms. Including $u$ in the posterior does not change the marginal posterior $p(x|y)$ which is what we are interested in. We now describe a practical example of this approach using an approximation consisting of a mixture of normals.

Albert (2009, Section 5.4) considers the problem of estimating the rates of death from stomach cancer for the largest cities in Missouri. This cancer mortality data is available from the R package LearnBayes, and consists of 20 pairs $(n_{j},y_{j})$ where $n_{j}$ contains the number of individuals that were at risk in city $j$, and $y_{j}$ is the number of cancer deaths that occurred in that city. The counts $y_{j}$ are overdispersed compared to what one could expect under a binomial model with constant probability, so Albert (2009) assumes the following beta-binomial model with mean $m$ and precision $K$:

where $B(\cdot,\cdot)$ denotes the Beta-function. The parameters $m$ and $K$ are given the following improper prior:

The resulting posterior distribution is non-standard and extremely skewed. To ameliorate this, Albert (2009) proposes the reparameterization

The form of the posterior distribution $p(x|y)$ still does not resemble any standard distribution, so we will approximate it using a finite mixture of $L$ bivariate Gaussians. In order to do this, we first introduce an auxiliary variable $u$, to which we assign a categorical approximate posterior distribution with $L$ possible outcomes:

where $\delta(.)$ is the indicator function and $U(\eta)$ is the normalizer.

Conditional on $u$, we assign $x$ a Gaussian approximate posterior

By adapting the true posterior to include $u$ as described above, we can fit this approximate posterior to $p(x|y)$. Here, the auxiliary variable $u$ is discrete, and hence our posterior approximation is not differentiable with respect to this variable. We must therefore use the basic stochastic approximation of Section 4 to fit $q_{\eta}(u)$. In order to reduce the variance of the resulting stochastic approximations, we Rao-Blackwellize them by taking expectations with respect to $q_{\eta}(u|x)$. If we then also take advantage of the sparsity in the covariance matrix of the sufficient statistics, this leads to the following update equations:

Conditional on $u$, the approximate posterior for $x$ is Gaussian, and we can therefore once again use the optimized expressions from Algorithm 2 to update $q_{\eta}(x|u)$:

for each mixture component $i$. Here we have once again Rao-Blackwellized the stochastic approximations with respect to $q_{\eta}(u|x)$, which introduced the extra variable $\hat{C}_{t,i}$ compared to Algorithm 2. Also note the presence of the $\log q_{\eta_{t}}(u=i|x^{*})$ term, which enters our equations as a result of expanding the posterior to include $u$. This term has the effect of pushing apart the different mixture components of the approximation.

We fit the approximation $q_{\eta}(x)$ using a varying number of mixture components and examine the resulting KL-divergence to the true posterior density. Since this is a low dimensional problem, we can obtain this divergence very precisely using quadrature methods. Figures 6 and 7 show that we can indeed approximate this skewed and heavy-tailed density very well using a large enough number of Gaussians. The R-squared of the mixture approximation with 8 components is 0.997.

Also apparent is the inadequacy of an approximation consisting of a single Gaussian for this problem, with an R-squared of only 0.82. This clearly illustrates the advantages of our approach which allows us to use much richer approximations than was previously possible. Furthermore, Figure 6 shows that the KL-divergence of the approximation to the true posterior can be approximated quite accurately using the measure developed in Section 5, especially if the posterior approximation is reasonably good.

The variational optimization problem for this approximation has multiple solutions, since all Gaussian mixture components are interchangeable. Since $p(x|y)$ is unimodal, however, we find that all local optima (that we find) are equally good, and are presumably also global optima. In this case, we find that we can therefore indeed approximate $p(x|y)$ arbitrarily well by using a large enough number of mixture components.

Conclusion and future work

We have introduced a stochastic optimization scheme for variational inference inspired by a novel interpretation of fixed-form variational approximation as linear regression of the target log density against the sufficient statistics of the approximating family. Our scheme allows very generic implementation for a wide class of models since in its most basic form only the unnormalized density of the target distribution is required, although we have shown how gradient or even Hessian information can be used if available. The generic nature of our methodology would lend itself naturally to a software package for Bayesian inference along the lines of Infer.NET (Minka et al., 2010) or WinBUGS (Gilks et al., 1994), and would allow inference in a considerably wider range of models. Incorporating automatic differentiation in such a package could clearly be beneficial. Automatic selection of the approximating family would be very appealing from a user perspective, but could be challenging in general.

Despite its general applicability, the performance of our approach was demonstrated to be very competitive for problems where we can either decompose the posterior distribution into low dimensional factors (Section 6.1), or where we can make use of the gradient and Hessian of the log posterior (Section 6.3). For those rare cases where this is not the case (e.g. high dimensional discrete distributions without factor structure) we cannot presently recommend the optimization algorithm presented in this paper. The extension of our approach to this class of problems is an important direction for future work.

We have shown it is straightforward to extend our methodology to use hierarchical structured approximations and more flexible approximating families such as mixtures. This closes the gap considerably relative to MCMC methods. Perhaps the biggest selling point of MCMC methods is that they are asymptotically exact: in practice this means simply running the MCMC chain for longer can give greater accuracy, an option not available to a researcher using variational methods. However, if we use a mixture approximating family then we can tune the computation time vs. accuracy trade off simply by varying the number of mixture components used. Another interesting direction of research along this line would be to use low rank approximating families such as factor analysis models.

Variational inference usually requires that we use conditionally conjugate models: since our method removes this restriction several possible avenues of research are opened. For example, for MCMC methods collapsed versions of models (i.e. with certain parameters or latent variables integrated out) sometimes permit much more efficient inference (Porteous et al., 2008) but adapting variational methods to work with collapsed models is complex and requires custom per model methodology (Teh et al., 2006). However, our method is indifferent to whether the model is collapsed or not, so it would be straightforward to experiment with different representations of the same model.

It is also possible to mix our method with VBEM, for example using our method for any non-conjugate parts of the model and VBEM for variables that happen to be conditionally conjugate. This is closely related to the non-conjugate variational message passing (NCVMP) algorithm of Knowles and Minka (2011) implemented in Infer.NET, which aims to fit non-conjugate models while maintaining the convenient message passing formalism. NCVMP only specifies how to perform the variational optimization, not how to approximate required integrals: in Infer.NET where analytic expectations are not available quadrature or secondary variational bounds are used, unlike the Monte Carlo approach proposed here. It is still an open question how these different methods could best be combined into a joint framework.

Acknowledgements

Tim Salimans wishes to acknowledge his advisors Richard Paap and Dennis Fok, as well as the anonymous referees, for their substantial help in improving the paper. He thanks The Netherlands Organization for Scientific Research (NWO) for financially supporting this project. DAK thanks Wolfson College, Cambridge, Microsoft Research Cambridge, and the Stanford Univeristy Center for Cancer Systems Biology for funding.

Appendix A Unnormalized to normalized optimality condition

The unnormalized optimality condition in (8) is

where $Y:=\log p(x,y)$. Rearranging gives

Note that (78), combined with $\mathop{\rm Cov}[T(x),\log q_{\eta}(x))]=\mathop{\rm Cov}[T(x),T(x)]\eta$ also implies that $\mathop{\rm Cov}[T(x),\log p(x,y)-\log q_{\eta}(x)]=0$ at a solution of the KL-divergence minimization. This is the same fixed point condition used in other applications of stochastic approximation variational Bayes such as Paisley et al. (2012).

Appendix B Derivation of Gaussian variational approximation

For notational simplicity we will derive our stochastic approximation algorithm for Gaussian variational approximation (Algorithm 2) under the assumption that $x$ is univariate. The extension to multivariate $x$ is conceptually straightforward but much more tedious in terms of notation.

Let $p(x,y)$ be the unnormalized posterior distribution of a univariate random variable $x$, and let $q(x)=N(m,V)$ be its Gaussian approximation with sufficient statistics, $T(x)=(x,-0.5x^{2})$. In order to find the mean $m$ and variance $V$ that minimize the KL-divergence between $q(x)$ and $p(x|y)$ we solve the transformed regression problem defined in (45), i.e.

with $\phi=(\phi_{1},\phi_{2})=(m,V)$ the usual mean-variance parameterization and where the natural parameters are given by $\eta=(V^{-1}m,V^{-1})$. Recall identity (43) which states that

with $\phi_{1}=m$ the first element of the parameter vector $\phi$, and $g(x)$ any differentiable function. Similarly, identity (44) reads

with $\phi_{2}=V$ the second element of the parameter vector. Using these identities we find that the regression statistics for this optimization problem are given by

where $Pm$ and $P=V^{-1}$ are the natural parameters (mean times precision and precision) of the approximation. Thus the quantities we need to stochastically approximate are

Appendix C Connection to Efficient Importance Sampling

It is worth pointing out the connection between fixed-form variational Bayes and Richard and Zhang’s (2007) Efficient Importance Sampling (EIS) algorithm. Although these authors take a different perspective (that of importance sampling) their goal of approximating the intractable posterior distribution with a more convenient distribution is shared with variational Bayes. Specifically, Richard and Zhang (2007) choose their posterior approximation to minimize the variance of the log-weights of the resulting importance sampler. This leads to an optimization problem obeying a similar fixed-point condition as in (9), but with the expectation taken over $p(x|y)$ instead of $q(x)$. Since sampling from $p(x|y)$ directly is not possible, they evaluate this expectation by sampling from $q(x)$ and weighting the samples using importance sampling. In practice however, these ‘weights’ are often kept fixed to one during the optimization process in order to improve the stability of the algorithm. When all weights are fixed to one, Richard and Zhang’s (2007) fixed-point condition becomes identical to that of (9) and the algorithm is in fact fitting a variational posterior approximation.

The connection between EIS and variational Bayes seems to have gone unnoticed until now, but it has some important consequences. It is for example well known (e.g. Minka, 2005; Nickisch and Rasmussen, 2008; Turner et al., 2008) that the tails of variational posterior approximations tend to be thinner than those of the actual posterior unless the approximation is extremely close, which means that using EIS with the importance-weights fixed to one is not to be recommended for general applications: In the case that the posterior approximation is nearly exact, one might as well use it directly instead of using it to form another approximation using importance sampling. In cases where the approximation is not very close, the resulting importance sampling algorithm is likely to suffer from infinite variance problems. The literature on variational Bayes offers some help with these problems. Specifically, de Freitas et al. (2001) propose a number of ways in which variational approximations can be combined with Monte Carlo methods, while guarding for the aforementioned problems.

Much of the recent literature (e.g. Teh et al., 2006; Honkela et al., 2010) has focused on the computational and algorithmic aspects of fitting variational posterior approximations, and this work might also be useful in the context of importance sampling. Algorithmically, the ‘sequential EIS’ approach of Richard and Zhang (2007) is most similar to the non-conjugate VMP algorithm of Knowles and Minka (2011). As these authors discuss, such an algorithm is not guaranteed to converge, and they present some tricks that might be used to improve convergence in some difficult cases.

The algorithm presented in this paper for fitting variational approximations is provably convergent, as discussed in Section 4. Furthermore, Sections 5 and 6 present multiple new strategies for variance reduction and computational speed-up that might also be useful for importance sampling. In this paper we will not pursue the application of importance sampling any further, but exploring these connections more fully is a promising direction for future work.

Appendix D Choosing an estimator

As discussed in Section 4, the particular estimator used in our stochastic approximation is not the most obvious choice, but it seems to provide a lower variance approximation than other choices. In this section we consider three different MC estimators for approximating (9) to see why this might be the case.

The first separately approximates the two integrals and then calculates the ratio:

with $S$ the number of Monte Carlo samples. The second approximates both integrals using the same samples from $q$:

Only this estimator is directly analogous to the linear regression estimator. The third estimator is available only when the first expectation is available analytically:

We wish to understand the bias/variance tradeoff inherent in each of these estimators. To keep notation manageable consider the case with only $k=1$ sufficient statisticThese results extend in a straightforward manner to the case where $k>1$. and let

We can now write the three estimators of $\eta$ more concisely as

From this we can derive expressions for the expectation of $f(y)$:

so that $\eta_{2}=f(y)=\frac{y_{2}}{y_{1}}$. Note that $\mathop{\rm Cov}(y)=\frac{1}{S}\mathop{\rm Cov}([a,b]^{\prime})$ and

We now turn to the variances. The analytic estimator is a standard MC estimator with variance

Consider only the linear terms of the Taylor expansion:

Substituting this into the formula for variance gives

We will calculate the variance of the second estimator and derive the variance of the first estimator from this. Again let $y$ be as in (94). Note that $\mathop{\rm Var}(y)=\mathop{\rm Cov}(a,b)/S$. We find

The final term is equal to that for the analytic estimator. The second term is not present in the variance of the first estimator, since then $a$ and $b$ have no covariance under the sampling distribution, i.e.

The first term is always positive, suggesting that $\hat{\eta}_{1}$ is dominated by the analytic estimator.

Note that the first term is shared, but the first estimator does not have the covariance term as a result of the independent sampling in approximating the numerator and denominator. In contrast $\hat{\eta}_{a}$ is unbiased. Now consider the variances

All three estimators have the same final term (the variance of the “analytic” estimator). Again the second estimator has an additional term resulting from the covariance between $a$ and $b$ which we find is typically beneficial in that it results in the variance of $\hat{\eta}$ being significantly smaller. It is worth recalling that the mean squared error (MSE) of an estimator is given by

We see that in this case for $\hat{\eta}_{2}$ the positive and negative contributions to both the bias and variance cancel. While this result will not hold exactly for cases of interest, it suggests that for exponential families which are capable of approximating $p$ reasonably well, $\hat{\eta}_{2}$ should perform significantly better than $\hat{\eta}_{1}$ or even $\hat{\eta}_{a}$. If $q$ and $p$ are of the same exponential family, it is actually possible to see that $\hat{\eta}_{2}$ will in fact give the exact solution in $k+1$ samples (with $k$ the number of sufficient statistics), while the other estimators have non-vanishing variance for a finite number of samples. This means that the approximate equality in (113) can be replaced by exact equality. Using $k+1$ samples $x_{i},i=1,...,k+1$, assumed to be unique (which holds almost surely for continuous distributions $q$), we have

That is, the algorithm has recovered $p(x,y)$ exactly with probability one. If we assume we know how to normalize $q$, this means we also have $p(x|y)$ exactly in this case. Note that we recover the exact answer here because the $p(x,y)$ function evaluations are in themselves noise free, so the regression analogy really corresponds to a noise free regression.

We test the three estimators in (79), (80) and (81) on the trivial exponential example of Section 4 when the true exponential rate is $\lambda=1.5$, and sampling from the optimal $q$ distribution with $\eta=1.5$. The results confirm that $\hat{\eta}_{2}$ finds the exact rate using just $S=2$ MC samples, as predicted by (115). We would expect $\hat{\eta}_{a}$ to be unbiased, and this is borne out by the results shown in Figure 8. The estimator $\hat{\eta}_{1}$ has both poor bias and such large variance that it often gives an invalid negative rate if fewer than 10 MC samples are used. While this is clearly a very simple example it hopefully emphasizes the potential benefit to be gained from using estimators related to $\hat{\eta}_{2}$.

References