TL;DR
This paper introduces a novel method combining Laplace approximation and variational inference to better capture posterior correlations and improve approximation quality in Bayesian inference.
Contribution
It proposes a mixed variational inference approach that retains applicability to non-conjugate models and preserves posterior correlations, overcoming limitations of traditional variational methods.
Findings
Improved approximation over Laplace method.
Enhanced posterior correlation modeling.
Numerical experiments show superior performance.
Abstract
The Laplace approximation has been one of the workhorses of Bayesian inference. It often delivers good approximations in practice despite the fact that it does not strictly take into account where the volume of posterior density lies. Variational approaches avoid this issue by explicitly minimising the Kullback-Leibler divergence DKL between a postulated posterior and the true (unnormalised) logarithmic posterior. However, they rely on a closed form DKL in order to update the variational parameters. To address this, stochastic versions of variational inference have been devised that approximate the intractable DKL with a Monte Carlo average. This approximation allows calculating gradients with respect to the variational parameters. However, variational methods often postulate a factorised Gaussian approximating posterior. In doing so, they sacrifice a-posteriori correlations. In this…
Click any figure to enlarge with its caption.
Figure 1
Figure 2
Figure 3
Figure 4| MVI | MVIeig | MVIlr | |
| # parameters | D | 2D | 3D |
| mean | |||
| covariance “root” |
| Laplace | MVI | MVIeig | MVIlr | VIdiag | |
|---|---|---|---|---|---|
| LPD | -0.818 | -0.771 | -0.726 | -0.736 | |
| MSE | 0.155 | 0.142 | 0.129 | 0.134 |
| Dataset | Q | Laplace | MVI | MVIeig | MVIlr | VIdiag | ||
|---|---|---|---|---|---|---|---|---|
| Banana | 2 | 400 | 4900 | -1238.76 | -1219.19 | -1221.41 | -1253.36 | |
| Breast cancer | 9 | 200 | 77 | -42.82 | -42.65 | -42.53 | -45.42 | |
| Diabetis | 8 | 468 | 300 | -145.98 | -145.468 | -145.31 | -193.27 | |
| Solar | 9 | 666 | 400 | -232.64 | -232.37 | -232.42 | -234.52 | |
| German | 20 | 700 | 300 | -151.71 | -151.42 | -151.31 | -179.29 | |
| Heart | 13 | 170 | 100 | -39.25 | -38.973 | -38.96 | -48.37 | |
| Image | 18 | 1300 | 1010 | -304.33 | -291.91 | -284.19 | -284.60 | |
| Ringnorm | 20 | 400 | 7000 | -309.87 | -319.67 | -309.80 | -342.627 | |
| Splice | 60 | 1000 | 2175 | -1156.80 | -900.00 | -900.30 | -899.501 | |
| Thyroid | 5 | 140 | 75 | -11.01 | -10.189 | -10.189 | -10.280 | |
| Titanic | 3 | 150 | 2051 | -1019.59 | -1021.7 | -1020.62 | -1023.91 | |
| Twonorm | 20 | 400 | 7000 | -452.57 | -450.716 | -461.10 | -543.115 | |
| Waveform | 21 | 400 | 4600 | -947.66 | -948.62 | -950.31 | -969.61 |
| Dataset | Laplace | MVI | MVIeig | MVIlr | VIdiag |
|---|---|---|---|---|---|
| Banana | 11.76 | 11.49 | 11.53 | 11.74 | |
| Breast cancer | 28.95 | 28.84 | 28.84 | 28.98 | |
| Diabetis | 24.79 | 24.65 | 24.67 | 34.33 | |
| Solar | 35.15 | 35.14 | 35.16 | 35.62 | |
| German | 26.06 | 25.91 | 25.96 | 29.00 | |
| Heart | 18.05 | 17.81 | 17.68 | 23.04 | |
| Image | 13.96 | 13.38 | 12.86 | 13.14 | |
| Ringnorm | 1.95 | 1.90 | 1.91 | 1.93 | |
| Splice | 25.74 | 18.98 | 18.96 | 18.89 | |
| Thyroid | 6.35 | 6.08 | 6.09 | 6.13 | |
| Titanic | 23.31 | 23.33 | 23.32 | 23.62 | |
| Twonorm | 2.92 | 2.86 | 2.92 | 2.96 | |
| Waveform | 10.40 | 10.36 | 10.34 | 10.44 |
| Dataset | K | Q | Laplace | MVI | MVIeig | MVIlr | VIdiag | ||
|---|---|---|---|---|---|---|---|---|---|
| Ecoli | 8 | 7 | 236 | 100 | -50.32 | -48.80 | -49.31 | -51.39 | |
| Crabs | 4 | 5 | 140 | 60 | -64.59 | -64.11 | -64.28 | -68.92 | |
| Iris | 3 | 4 | 105 | 45 | -9.06 | -7.53 | -7.46 | -8.17 | |
| Soybean | 4 | 35 | 33 | 14 | -4.10 | -2.35 | -2.36 | -1.67 | |
| Wine | 3 | 13 | 125 | 53 | -5.72 | -4.01 | -3.94 | -4.66 | |
| Glass | 6 | 9 | 150 | 64 | -61.35 | -60.39 | -60.44 | -76.26 | |
| Vehicle | 4 | 18 | 593 | 293 | -159.783 | -158.60 | -159.15 | -174.725 | |
| Balance | 3 | 4 | 438 | 187 | -23.5734 | -23.2197 | -23.078 | -31.587 |
| Dataset | Laplace | MVI | MVIeig | MVIlr | VIdiag |
|---|---|---|---|---|---|
| Ecoli | 17.22 | 16.75 | 17.23 | 17.87 | |
| Crabs | 55.75 | 55.10 | 55.00 | 58.96 | |
| Iris | 9.24 | 7.86 | 7.58 | 9.76 | |
| Soybean | 13.12 | 4.83 | 4.86 | 9.33 | |
| Wine | 4.97 | 3.31 | 3.14 | 4.72 | |
| Glass | 42.93 | 41.14 | 41.14 | 53.01 | |
| Vehicle | 32.60 | 32.34 | 32.22 | 34.61 | |
| Balance | 6.56 | 6.27 | 6.25 | 8.06 |
Peer Reviews
No public reviews on file for this paper yet. If you reviewed it on a platform where reviews are public (OpenReview, ICLR, NeurIPS, ICML), you can paste yours below so the community can read it here.
Code & Models
Videos
No videos yet. Explain this paper in a talk, walkthrough, or lecture? Add one.
Taxonomy
MethodsDeep Kernel Learning
Mixed Variational Inference
††thanks: The author gratefully acknowledges the generous and invaluable support of the Klaus Tschira Foundation.
Nikolaos Gianniotis
Astroinformatics
*Heidelberg Institute for Theoretical Studies gGmbH
*Heidelberg, Germany
Abstract
The Laplace approximation has been one of the workhorses of Bayesian inference. It often delivers good approximations in practice despite the fact that it does not strictly take into account where the volume of posterior density lies. Variational approaches avoid this issue by explicitly minimising the Kullback-Leibler divergence DKL between a postulated posterior and the true (unnormalised) logarithmic posterior. However, they rely on a closed form DKL in order to update the variational parameters. To address this, stochastic versions of variational inference have been devised that approximate the intractable DKL with a Monte Carlo average. This approximation allows calculating gradients with respect to the variational parameters. However, variational methods often postulate a factorised Gaussian approximating posterior. In doing so, they sacrifice a-posteriori correlations. In this work, we propose a method that combines the Laplace approximation with the variational approach. The advantages are that we maintain: applicability on non-conjugate models, posterior correlations and a reduced number of free variational parameters. Numerical experiments demonstrate improvement over the Laplace approximation and variational inference with factorised Gaussian posteriors.
Index Terms:
Bayesian inference, variational techniques, Laplace approximation, mean field, Gaussian approximation
I Introduction
Bayesian inference provides a way of making use of the complete information available either as data or prior knowledge. It enables us to capture the uncertainty present in the data and model assumptions, and propagate it to further tasks such as prediction and decision making. However, exact Bayesian inference is only possible whenever mathematically convenient priors are combined with particular likelihood functions (e.g. conjugacy). Deviation from such convenience, results in intractable calculations that call for approximations.
The Laplace approximation (LA) has helped the advancement of Bayesian methodology in the machine learning field [1] and has also been an important tool for practitioners [2]. LA produces a Gaussian approximating posterior. In doing so, it operates myopically in the sense that it determines the mean and covariance simply by looking locally around the mode instead of focusing on where the volume of the density actually lies. Despite this shortcoming, it has been found to produce good approximations in a variety of contexts. In a Gaussian process binary classification setting [3], LA is found to be on a par with other approximations in terms of error rate, though it performed poorer on other criteria. The work in [4] employs LA in order to calculate approximate intractable integrals within the Expectation Propagation algorithm [5]. In [6] LA is found to perform well when compared to Expectation Propagation in a bounded regression task. In [7] LA is used to formulate an approximate form of the marginal likelihood that facilitates the update of hyperparameters without the need to update the current Gaussian posterior given by LA.
Variational inference has been put forward (VI) [8, 9] as a solution to calculating posterior distributions in situations where certain expectations are not analytically tractable. Its use has been widespread in eliciting posterior densities in e.g. dimensionality reduction [10], classification [11], regression [9], density estimation [12] and in specialised applications like in astronomy [13]. The applicability of VI depends on choosing a posterior density form that allows the Kullback-Leibler divergence () between the approximating and the true (unnormalised) posterior to be calculated in closed form. However, such a choice may not always exist.
Whenever it is not possible to obtain a closed-form within the VI framework, approximations become necessary. One type of approximation approximates the logarithmic (unnormalised) posterior. In [14] the logarithmic posterior is linearised via a first-order Taylor expansion which allows then calculating the expectation with respect to the approximating posterior in the . In a similar vein, second order Taylor expansions are considered in [15]. Interestingly, [16] considers multiple second order Taylor expansions at different parameter locations of the logarithmic posterior. This results in an approximate posterior density expressed as a mixture of spherical Gaussians that has the potential to capture multiple modes. A second type of approximation [17, 18, 19, 20] approximates the as a Monte Carlo average with samples drawn from the approximating posterior. The resulting expression allows calculating the gradient with respect to the free variational parameters. An update of the variational parameters follows, typically using a small step size in a stochastic gradient descent setting, after which the is approximated with a new Monte Carlo average. In a slightly manner, [21, 22] fix the Monte Carlo average approximation of the throughout the optimisation of the variational parameters. We finally note that often in practice e.g. [17, 19], VI methods choose to work with a factorised Gaussian posterior which has the advantage of reducing the number of free variational parameters that need to be optimised, but also has the inevitable disadvantage of discarding potential parameter correlations in the posterior.
In this work, we propose a method that combines the Laplace approximation with the variational approximation. The method works on non-conjugate models, captures a-posteriori correlations and limits the number of free variational parameters. The main idea is to take the Gaussian posterior obtained from the Laplace approximation, plug it into the variational lower bound and adapt it by optimising the lower bound. The crux of the approach is to allow only a partial update of the Laplace Gaussian posterior.
II Approximate Bayesian Inference
We briefly review methods for approximate inference as a gentle reminder and for the purpose of introducing relevant notation. We write the log-posterior as the sum of the model log-likelihood, log-prior and minus log-evidence:
[TABLE]
where are the data, \mbox{\boldmathw}\in\mathbb{R}^{D} are the model parameters and \mathcal{Z}=\int p(\mathcal{D}|\mbox{\boldmathw},\mbox{\boldmath\theta}_{\ell})p(\mbox{\boldmathw}|\mbox{\boldmath\theta}_{\pi})\mbox{\boldmathdw}. The log-likelihood and log-prior terms have hyperparameters \mbox{\boldmath\theta}_{\ell} and \mbox{\boldmath\theta}_{\pi} which are jointly summarised as in the log-posterior. In the following, the evidence is a constant which we discard. Discarding it, gives us the unnormalised log-posterior \ln\tilde{p}(\mbox{\boldmathw}|\mbox{\boldmath\theta},\mathcal{D}).
II-A Laplace approximation
The Laplace approximation (LA) seeks the mode111Multiple modes may be present. \mbox{\boldmathw}^{*} of the log-posterior density \ln\tilde{p}(\mbox{\boldmathw}|\mbox{\boldmath\theta},\mathcal{D}) where \mbox{\boldmath0}=\nabla_{\tiny\mbox{\boldmathw}}\ln\tilde{p}(\mbox{\boldmathw}|\mbox{\boldmath\theta},\mathcal{D})|_{\scriptstyle\mbox{\boldmath\scriptstyle w}=\mbox{\boldmath\scriptstyle w}^{*}}. This may be carried out with gradient-based optimisation. At the found mode, we calculate the Hessian matrix \mbox{\boldmathH}=\nabla\nabla\ln\tilde{p}(\mbox{\boldmathw}|\mbox{\boldmath\theta},\mathcal{D})|_{\scriptstyle\mbox{\boldmath\scriptstyle w}=\mbox{\boldmath\scriptstyle w}^{*}}. The obtained approximating Gaussian posterior reads:
[TABLE]
We see that the covariance of the approximating posterior q(\mbox{\boldmathw}) is given by the local curvature of the posterior at the found mode. The approximation can be good, if the true posterior concentrates strongly around the mode.
II-B Variational Inference
Variational inference (VI) [8] postulates an approximating posterior q(\mbox{\boldmathw}). VI finds the q(\mbox{\boldmathw}) that maximises the following based objective, also known as the variational lower bound [23, Chapter ] :
[TABLE]
In general, VI does not require that q(\mbox{\boldmathw}) is a Gaussian, but here we choose to work with q(\mbox{\boldmathw})=\mathcal{N}(\mbox{\boldmathw}|\mbox{\boldmath\mu},\mbox{\boldmath\Sigma}). For this choice, the above objective now reads as:
[TABLE]
where the second term is the Gaussian entropy. The free parameters in objective (4) are the variational parameters , and hyperparameters .
II-C Stochastic variational inference
Stochastic variational inference [19] addresses the difficulty that arises when the expectation in the first term of (4) is not tractable. It does so by approximating the expectation by a Monte Carlo average with samples drawn from \mbox{\boldmathw}_{s}\sim q(\mbox{\boldmathw}):
[TABLE]
The variational parameters no longer appear in the approximation, but it is possible to reintroduce them using the reparametrisation \mbox{\boldmathw}_{s}=\mbox{\boldmath\mu}+\mbox{\boldmathC}\mbox{\boldmathz}_{s}, in terms of samples \mbox{\boldmathz}_{s}\sim\mathcal{N}(\mbox{\boldmath0},\mbox{\boldmathI}_{D}):
[TABLE]
where is a matrix222A common choice is the Cholesky decomposition. such that \mbox{\boldmathC}\mbox{\boldmathC}^{T}=\mbox{\boldmath\Sigma}. The free parameters in objective (6) are , and . We note that, typically, one chooses covariance to be a diagonal matrix (e.g. [17, 19]) in order to limit the number of free variational parameters to be optimised. In this case, q(\mbox{\boldmathw}) is a factorised posterior.
III Proposed method
The motivation behind this work is to apply VI on non-conjugate models using an approximating posterior q(\mbox{\boldmathw}) that captures a-posteriori correlations but at the same time limits the number of free variational parameters that need to be optimised. To that end, we make use of the covariance \mbox{\boldmath\Sigma}_{LA} of the approximating posterior q(\mbox{\boldmathw}) obtained via LA and the approximate variational lower bound in (6). Since the proposed method combines LA with the variational lower bound, we name it mixed variational inference (MVI). In the following, we propose three ways that MVI can exploit the correlation structure present in \mbox{\boldmath\Sigma}_{LA}.
III-A Adaptation of mean only - MVI
We perform the following Cholesky decomposition:
[TABLE]
We propose the posterior q(\mbox{\boldmathw})=\mathcal{N}(\mbox{\boldmathw}|\mbox{\boldmath\mu},\mbox{\boldmath\Sigma}_{LA}) and use it in the approximate variational lower bound in (6), which results in the following objective:
[TABLE]
The free parameters in (8) are the mean and hyperparameters . Effectively, the proposed posterior is the Laplace posterior with the added flexibility of shifting its mean while keeping its covariance fixed to \mbox{\boldmath\Sigma}_{LA}. Note, that here the entropy is a constant term that can be discarded during optimisation.
III-B Mean and scaling of covariance - MVIeig
We perform the following eigenvalue decomposition:
[TABLE]
where matrix \mbox{\boldmathQ}_{LA}\in\mathbb{R}^{D\times D} and vector \mbox{\boldmathr}_{LA}\in\mathbb{R}^{D} hold the eigenvectors and square roots of the eigenvalues respectively333Operator creates a diagonal matrix using the vector it is applied to. Notation \mbox{\boldmathr}^{2} implies raising the components of vector to the power of .. We propose q(\mbox{\boldmathw})=\mathcal{N}(\mbox{\boldmathw}|\mbox{\boldmath\mu},\mbox{\boldmathQ}_{LA}\operatorname{diag}(\mbox{\boldmathr}^{2})\mbox{\boldmathQ}_{LA}^{T}) and optimise:
[TABLE]
The free parameters in (10) are , and . Note the simplification in the entropy term due to the orthogonal \mbox{\boldmathQ}_{LA}, i.e. |2\pi e\operatorname{diag}(\mbox{\boldmathr}^{2})\mbox{\boldmathQ}_{LA}^{T}\mbox{\boldmathQ}_{LA}|=|2\pi e\operatorname{diag}(\mbox{\boldmathr}^{2})|. Effectively, the proposed posterior is the Laplace posterior which now has the added flexibility to shift the mean and scale the covariance matrix along its axes by adapting vector .
III-C Mean and low rank update of covariance - MVIlr
We introduce the vectors \mbox{\boldmathU},\mbox{\boldmathV}\in\mathbb{R}^{D}. We use the Cholesky decomposition \mbox{\boldmath\Sigma}_{LA}=\mbox{\boldmathC}_{LA}\mbox{\boldmathC}_{LA}^{T} and form the matrix \mbox{\boldmathL}=\mbox{\boldmathC}_{LA}+\mbox{\boldmathU}\mbox{\boldmathV}^{T}. We propose the posterior q(\mbox{\boldmathw})=\mathcal{N}(\mbox{\boldmathw}|\mbox{\boldmath\mu},\mbox{\boldmathL}\mbox{\boldmathL}^{T}) and the associated objective:
[TABLE]
The free parameters in (11) are , , and . Effectively, the proposed posterior is the Laplace posterior which now has the added flexibility to shift the mean but also modify its covariance matrix via a low-rank update.
The proposed posteriors are summarised in Table I.
III-D Initialisation
We use the mean \mbox{\boldmath\mu}_{LA} and optimised hyperparameters \mbox{\boldmath\theta}_{LA} obtained from the Laplace approximation to initialise the mean in q(\mbox{\boldmathw})=\mathcal{N}(\mbox{\boldmathw}|\mbox{\boldmath\mu}=\mbox{\boldmath\mu}_{LA},\mbox{\boldmath\Sigma}) and hyperparameters \mbox{\boldmath\theta}=\mbox{\boldmath\theta}_{LA} in each of the three proposed objectives. We emphasize that the covariance in MVI is initialised to \mbox{\boldmath\Sigma}_{LA} and fixed. Vector in MVIeig is initialised to the square root of the eigenvalues \mbox{\boldmathr}_{LA}. Vectors , in MVIlr are randomly initialised by drawing them from \mathcal{N}(\mbox{\boldmath0},0.01\mbox{\boldmathI}_{D}).
III-E Optimisation
Following [21, 22] we draw number of samples \mbox{\boldmathz}_{s}\sim\mathcal{N}(\mbox{\boldmath0},\mbox{\boldmathI}_{D}) which we keep fixed throughout the optimisation of the objectives in (8), (10) and (11). This enables the use of scaled-conjugate gradients (SCG) as the optimisation routine[24] in contrast to the typically employed stochastic gradient descent444We note that in [19] the use of stochastic gradient is additionally motivated by the desire to train with “mini-batches”. [19]. We note that the proposed method can in principle also employ the same optimisation scheme as in [19]. The free parameters , and the ones pertaining to the covariance in each proposed posterior are jointly optimised via SCG. In all experiments we fix the number of drawn samples \mbox{\boldmathz}_{s} to .
IV Numerical setup
IV-A Comparisons
The proposed work builds on LA and VI in order to improve the performance (see section IV-B) of the Laplace approximation and do better than VI when employing a factorised Gaussian posterior. Specifically, the proposed posterior for the latter reads \mbox{\boldmathq}(\mbox{\boldmathw})=\mathcal{N}(\mbox{\boldmathw}|\mbox{\boldmath\mu},\operatorname{diag}(\mbox{\boldmath\sigma}^{2})) and has a diagonal covariance matrix whose elements are specified by the vector \mbox{\boldmath\sigma}\in\mathbb{R}^{D}. The associated objective reads:
[TABLE]
The free parameters in (12) are , and . We refer to this method as VIdiag.
In the numerical experiments, we initialise the mean and hyperparameters with \mbox{\boldmath\mu}=\mbox{\boldmath\mu}_{LA} and \mbox{\boldmath\theta}=\mbox{\boldmath\theta}_{LA}. Regarding \mbox{\boldmath\sigma}^{2}, we experimented with two initialisations: either setting the elements of \mbox{\boldmath\sigma}^{2} equal to the diagonal elements of \mbox{\boldmath\Sigma}_{LA}, or all equal to . In the experiments of Section V we report for VIdiag the best performance achieved by either initialisation.
IV-B Measuring performance
In the experiments of Section V, we measure performance in terms of the logarithmic predictive density (LPD) (i.e. marginal log-likelihood) evaluated on test data:
[TABLE]
The LPD is approximated by number of samples drawn from \mbox{\boldmathw}_{s^{\prime}}\sim\mbox{\boldmathq}(\mbox{\boldmathw}), where \mbox{\boldmathq}(\mbox{\boldmathw}) is the respective posterior obtained via LA, MVI or VIdiag. In all numerical experiments we use .
Along LPD, we also report error rates. For the regression problem we report the mean squared error (MSE). For the classification problems, the error rate is the percentage of predicted labels not matching the true labels.
We compare MVI to LA and VIdiag on a number of datasets as detailed in the corresponding sections. Each dataset is split times into a training and testing set. The algorithms are run on each split, hence, we collect samples of the algorithms’ performance in terms of LPD and error rate on the test set. For each dataset, we report the median LPD and median error rate on the test set for each algorithm. The best performance is marked with bold in the tables reporting the results.
Moreover, we attempt to detect whether the observed differences in median, over the collected performances, are statistically significant. In the experiments, we observed that the collected performances are not normally distributed which precludes the use of a paired T-test. The Wilcoxon signed rank test is also precluded as it requires [25, Chapter 4.7] that the distribution of the difference in median of the tested pairs is symmetric. Therefore, we resort to using a sign test [25, Chapter 2.5.2], to check whether a difference in median performance exists (i.e. better or worse, but not by how much), and the confidence intervals constructed by the bootstrap [26].
We test whether the performance of the best algorithm (marked in the tables with bold) is statistically significantly better by checking two conditions: we pair the best performing algorithm with all other algorithms and carry out the sign test. The first condition is satisfied if for each pair, the sign test rejects the null hypothesis that the median of the best algorithm is equal to the median of its respective paired algorithm. The second condition is satisfied if the confidence interval constructed by the bootstrap on the difference of the paired medians does not contain the [math] value. That is, if the confidence interval does not contain [math], then [math] is not a likely value for the difference in the true medians. If both conditions are satisfied, we declare the best performance as statistically significant and mark it with a marker in the respective tables.
V Applications
We first demonstrate the behaviour of the proposed MVI posteriors on two synthetic examples. We then proceed with experiments on benchmark problems.
V-A Illustration with 2D posterior
We illustrate the MVI posteriors on a synthetic 2D example where we specify the true, target posterior as a mixture of two Gaussian components:
[TABLE]
Fig. 1 displays a contour plot of p(\mbox{\boldmathw}). We approximate p(\mbox{\boldmathw}) with LA and the three proposed MVI posteriors and plot them in Fig. 1. The figure also reports the \operatorname{D\textsubscript{KL}}(q(\mbox{\boldmathw})||p(\mbox{\boldmathw})) of each approximating posterior to the true posterior. Here \operatorname{D\textsubscript{KL}}(q(\mbox{\boldmathw})||p(\mbox{\boldmathw})) is calculated numerically as there is no, at least not straightforward, closed-form expression for it. We observe that LA (black), by design, places its mean on the mode of p(\mbox{\boldmathw}). All other approximations place their means on alternative locations, but fairly close to one another. Even MVI (red), whose covariance is constrained to be equal to that of LA, can shift its mean to a better location so that it covers more of the target density. We also note how the axes of MVIeig (green) are parallel by design, but scaled compared to the axes of MVI (red), i.e. we see that the red ellipse of MVI is contained in the green ellipse of MVIeig. By scaling its axes, MVIeig achieves a lower . MVIlr has the flexibility of rotating its covariance and, in this case, achieves the lowest to the true posterior p(\mbox{\boldmathw}).
V-B Robust regression on synthetic task
We experiment with a regression task with input-target pairs , . Inputs are drawn uniformly in . Targets are generated through the expression
[TABLE]
and corrupted with i.i.d. noise drawn from the uniform distribution with support . This is a regression task where adopting a Gaussian likelihood would lead to poor results as it cannot adequately explain the noise. We adopt a Cauchy density instead. The unnormalised log-posterior reads:
[TABLE]
where is the Cauchy density. Additionally, we have calculated a set of radial basis functions on the data inputs:
[TABLE]
where \phi(\mbox{\boldmathx}_{n};r,\mbox{\boldmathc}_{m})=\exp(-\frac{\|\mbox{\boldmathx}_{n}-\mbox{\boldmathc}_{m}\|^{2}}{2r^{2}}). The last element in (17) serves as a bias term. Hence, \mbox{\boldmath\phi}_{n}\in\mathbb{R}^{M+1} and \mbox{\boldmathw}\in\mathbb{R}^{D} with .
In this numerical experiment, we generate datasets with training data items using (15). We also generate test data items in precisely the same way. We report the median log-predictive density (LPD) on test data in Table II. Best performances are marked with bold. We see that the MVIeig approximation performs the best, hence we mark it in bold. When looking at the results, we established that MVIeig performs statistically significantly better than LA, MVI and VIdiag. However, as MVIeig does not outperform MVIlr with statistical significance, we do not mark it additionally with marker. In terms of error rate, we see that MVIlr performs best and marginally better than MVIeig. Finally, in figure 2 we show the true underlying curve, specified in (15), the observed training data as filled circles along with the regressions induced by MVI and MVIeig.
V-C Logistic Regression
We experiment with logistic regression [23, Chapter ]. The data are input-label pairs (\mbox{\boldmathx}_{n},y_{n}) with \mbox{\boldmathx}_{n}\in\mathbb{R}^{Q},\ y_{n}\in\{0,1\}. Just like in Section V-B, equation (17), we calculate a set of radial basis functions \mbox{\boldmath\phi}_{n}\in\mathbb{R}^{M+1} on the data inputs \mbox{\boldmathx}_{n}. The weights are given by \mbox{\boldmathw}\in\mathbb{R}^{D} with . The unnormalised log-posterior reads:
[TABLE]
To avoid inadvertently selecting single datasets on which the proposed algorithm performs well, we experiment with the entire collection of datasets preprocessed by Rätsch et al555http://www.raetschlab.org/Members/raetsch/benchmark. Each dataset has been standardised and split into training and testing instances, except for Image and Splice that have 20 splits. We approximate the log-posterior in (18) with LA, the MVI posteriors and VIdiag. To initialise the hyperparameters \mbox{\boldmath\theta}=(M,r,\alpha,\mbox{\boldmathc}_{1},\dots,\mbox{\boldmathc}_{M}) in Laplace, we proceed as follows: per dataset, we run LA for iterations for each combination of and randomly drawn pairs , , i.e. a total of combinations. The centres \mbox{\boldmathc}_{m} are determined by K-means for each choice of . The combination with the highest lower bound (we are maximising) is declared the winner and used to initialise LA which is then run for a maximum of iterations. The MVI posteriors and hyperparameters in the respective objectives are initialised using the optimised Laplace posterior and hyperparameters, as described in Section III-D.
We report the median log-predictive density (LPD) on test data in Table III, along with details about the datasets. Best performances are marked with bold. Best performances that differ in a statistically significant way to all other performances (see Section IV-B) are additionally marked with a marker. Table III reveals that, in general, the proposed MVI posteriors perform better than LA or VIdiag. In particular, we see that MVIlr scores better on a number of datasets and that the difference in performance is often statistically significant. Table IV displays the results on error rates. We see that all methods achieved more or less the same error rates with no performance being statistically significantly superior. Nonetheless, we do observe a few exceptions, e.g. on datasets Splice and Diabetis LA and VIdiag respectively perform noticeably poorer.
V-D Multiclass Logistic Regression
Similarly to logistic regression, multiclass logistic regression [23, Chapter ] does not allow direct Bayesian inference as the use of the softmax function renders integrals over the likelihood term intractable. The unnormalised log-posterior reads:
[TABLE]
where denotes the total number of classes. The data are input-label pairs (\mbox{\boldmathx}_{n},\mbox{\boldmathy}_{n}) with \mbox{\boldmathx}_{n}\in\mathbb{R}^{Q}. Vectors \mbox{\boldmathy}_{n} are binary vectors encoding class labels using a -of- coding scheme, e.g. encodes class label in a -class problem. The probability p(C_{k}|\mbox{\boldmath\phi}_{n}) of the -th data item belonging to class is modelled via the softmax function:
[TABLE]
where each class is associated with a weight vector \mbox{\boldmathw}_{k}\in\mathbb{R}^{D}, with . The basis functions \mbox{\boldmath\phi}_{n}\in\mathbb{R}^{M+1} are defined in the same way as in Section V-B. We initialise hyperparameters \mbox{\boldmath\theta}=(M,r,\alpha,\mbox{\boldmathc}_{1},\dots,\mbox{\boldmathc}_{M}) in the same way as described in Section V-C.
To avoid inadvertently selecting single datasets on which the proposed algorithm performs well, we experiment with the collection of multiclass datasets used in the work of [27] in a different context. Details of the datasets are shown in Table V. We standardise the data column-wise to zero mean and unit standard deviation. Using random subsampling, we split each dataset times into a training ( of the data) and testing () set. We report the median LPD for each dataset and algorithm in Table V and median error rate in Table VI. Again, best performances are marked in bold. We mark the best performance with a marker if it is found to be statistically significant using the same two conditions described in Section IV-B. The results show an improvement over LA and the use of a factorised posterior in VIdiag. We also see that MVIeig performs well on this set of problems in terms of LPD, though the picture is not as clear in terms of error rate in Table VI. Finally, we note the low performance of all methods on the dataset Crab, evidently in Table VI. This may be perhaps attributed to the particular choice of the RBF kernel made here, though other kernels (cf [27]) may be more appropriate.
VI Discussion and Conclusion
We proposed Mixed Variational Inference (MVI) as a method for approximate Bayesian inference in non-conjugate models. MVI makes use of the posterior obtained via the Laplace approximation and the objective function provided by variational inference. The adoption of the Laplace posterior helps with capturing a-posteriori correlations; the partial adaptation of the Laplace posterior, in the form of the proposed MVI posteriors, helps limit the number of free variational parameters that need to be optimised. The numerical results show that the MVI posteriors have the potential to improve on the performance of the Laplace approximation and on the performance of the commonly adopted factorised Gaussian posterior in variational inference.
Strictly speaking, however, one should be aware of the fact that a posterior q(\mbox{\boldmathw}) that approximates the true posterior better, does not necessarily guarantee improved log-predictive density; vice versa, a “naive” approximating posterior (e.g. factorised) may in principle provide satisfactory predictive performance. This observation has been previously stated in [3] where a variety of approximations are evaluated in the context of Gaussian process binary classification. Therein it is stated that, in principle, even a poor approximation in terms of posterior moments can still provide good predictions. After all, as far as variational approximations are concerned, it is evident in objective (3) that the goal is to find a q(\mbox{\boldmathw}) that is as close as possible to the true posterior; this does not necessarily correlate with improved predictive performance. Nevertheless, one does expect in practice that an approximation that captures posterior correlations in the parameters to be more useful than a factorised approximation that practically draws the parameters independently of one another when making predictions (see (13)). But beyond this expectation, it is admittedly difficult to anticipate what approximation may perform best. Indeed, in the numerical experiments we notice that while MVIlr seems well suited for logistic regression (see Table III), this is not necessarily the case in multiclass logistic regression (see Table V).
In its present form, MVI is limited to Gaussian posteriors. It would be interesting to extend MVI to non-Gaussian posteriors, though at first sight it seems that its dependence on the Laplace approximation considerably limits it. An interesting direction, inspired by [16], would be to postulate a posterior q(\mbox{\boldmathw}) based on a mixture of Gaussians, where the covariance of each Gaussian comes from a Laplace approximation performed at a different mode. Forming an approximating posterior using multiple modes procured by the Laplace approximation has been previously suggested in [28]. However, therein no objective function akin to (6) is guiding the inference of the posterior. This could be potentially addressed by extending MVI so that q(\mbox{\boldmathw}) is now a Gaussian mixture whose covariance matrices are given by the Laplace approximation and partially updated as suggested in Sections III-A, III-B, III-C. We reserve such investigations for future research.
Acknowledgment
The author acknowledges the useful discussions and encouragement of Christoph Schnörr, Kai Polsterer and Ata Kaban.
The reference list from the paper itself. Each links out to its DOI / PubMed record.
- 1[1] D. J. Mac Kay, “A practical bayesian framework for backpropagation networks,” Neural computation , vol. 4, no. 3, pp. 448–472, 1992.
- 2[2] D. Sivia and J. Skilling, Data analysis: a Bayesian tutorial . OUP Oxford, 2006.
- 3[3] H. Nickisch and C. E. Rasmussen, “Approximations for binary gaussian process classification,” Journal of Machine Learning Research , vol. 9, no. Oct, pp. 2035–2078, 2008.
- 4[4] E. Eskin, A. J. Smola, and S. Vishwanathan, “Laplace propagation,” in Advances in neural information processing systems , 2004, pp. 441–448.
- 5[5] T. P. Minka, “Expectation propagation for approximate bayesian inference,” in Proceedings of the Seventeenth conference on Uncertainty in artificial intelligence . Morgan Kaufmann Publishers Inc., 2001, pp. 362–369.
- 6[6] B. S. Jensen, J. B. Nielsen, and J. Larsen, “Bounded gaussian process regression,” in Machine Learning for Signal Processing (MLSP), 2013 IEEE International Workshop on . IEEE, 2013, pp. 1–6.
- 7[7] A. Wu, N. G. Roy, S. Keeley, and J. W. Pillow, “Gaussian process based nonlinear latent structure discovery in multivariate spike train data,” in Advances in Neural Information Processing Systems , 2017, pp. 3496–3505.
- 8[8] M. J. Beal, “Variational Algorithms for Approximate Bayesian Inference,” Ph.D. dissertation, Gatsby Computational Neuroscience Unit, University College London, 2003.
