Approximate Variational Inference Based on a Finite Sample of Gaussian Latent Variables

TL;DR

This paper introduces a new variational inference scheme that simplifies approximation for complex integrals involving Gaussian latent variables, especially in cases where traditional bounds are hard to derive.

Contribution

It proposes a novel approach to variational inference that overcomes difficulties in deriving lower bounds for certain integrals, broadening applicability.

Findings

Effective on synthetic and real datasets

Applicable to geophysical models where standard methods fail

Simplifies variational inference in complex scenarios

Abstract

Variational methods are employed in situations where exact Bayesian inference becomes intractable due to the difficulty in performing certain integrals. Typically, variational methods postulate a tractable posterior and formulate a lower bound on the desired integral to be approximated, e.g. marginal likelihood. The lower bound is then optimised with respect to its free parameters, the so called variational parameters. However, this is not always possible as for certain integrals it is very challenging (or tedious) to come up with a suitable lower bound. Here we propose a simple scheme that overcomes some of the awkward cases where the usual variational treatment becomes difficult. The scheme relies on a rewriting of the lower bound on the model log-likelihood. We demonstrate the proposed scheme on a number of synthetic and real examples, as well as on a real geophysical model for which the standard variational approaches are inapplicable.

Tables

Table 1. Table 1: Classification performance for Bayesian logistic regression (higher is better).

Dataset	Bound in [11]	Proposed
Banana	0.8893 $\pm$ 0.0055	0.8893 $\pm$ 0.0054
Cancer	0.7119 $\pm$ 0.0456	0.7116 $\pm$ 0.0454
Heart	0.5528 $\pm$ 0.0445	0.5500 $\pm$ 0.0476
Solar	0.6449 $\pm$ 0.0172	0.6445 $\pm$ 0.0167

Table 2. Table 2: Classification performance for Bayesian multiclass classification (higher is better).

Dataset	C	Suggested Kernel	${\text{mRVM}}_2$	Proposed
Ecoli	8	Gaussian	$0.855\pm 0.047$	$0.852\pm 0.060$
Glass	6	Polynomial	$0.581\pm 0.143$	$0.667\pm 0.112$
Iris	3	Gaussian	$0.913\pm 0.071$	$0.947\pm 0.053$
Wine	3	Linear	$0.959\pm 0.062$	$0.976\pm 0.030$
Soybean (small)	4	Linear	$1.000\pm 0.000$	$1.000\pm 0.000$
Vehicle	4	Polynomial	$0.481\pm 0.054$	$0.539\pm 0.103$
Balance	3	Polynomial	$0.931\pm 0.044$	$0.947\pm 0.029$
Crabs	4	Linear	$0.915\pm 0.071$	$0.950\pm 0.033$

Table 3. Table 3: Prediction error for ground motion problem (lower is better).

$N$	Laplace	Proposed
100	$0.8550\pm 0.4008$	$0.6040\pm 0.0307$
200	$0.7411\pm 0.4925$	$0.5776\pm 0.0403$
500	$0.6926\pm 0.4790$	$0.5496\pm 0.0451$
1000	$0.5395\pm 0.0230$	$0.5323\pm 0.0273$

Equations

$$\displaystyle p(\mbox{\boldmath$Y$}|\mbox{\boldmath$X$},\mbox{\boldmath$w$},\beta)$$

$$p(\mbox{\boldmath$w$}|\alpha)={\mathcal{N}}(\mbox{\boldmath$w$}|\mbox{\boldmath$0$},\alpha^{-1}\mbox{\boldmath$I$}_{M})\ .$$

$$q(\mbox{\boldmath$w$})={\mathcal{N}}(\mbox{\boldmath$w$}|\mbox{\boldmath$\mu$},\mbox{\boldmath$\Sigma$})\ .$$

$$\displaystyle q(\mbox{\boldmath$w$})={\mathcal{N}}(\mbox{\boldmath$w$}|\mbox{\boldmath$\mu$},\mbox{\boldmath$L$}\mbox{\boldmath$L$}^{T})\ .$$

$$\displaystyle\log p(\mbox{\boldmath$Y$}|\mbox{\boldmath$X$},\alpha,\beta)$$

$$\displaystyle\frac{1}{S}\sum_{s=1}^{S}\log p(\mbox{\boldmath$Y$}|\mbox{\boldmath$X$},\mbox{\boldmath$w$}_{(s)},\beta)\ ,$$

$$\displaystyle\mbox{\boldmath$w$}_{(s)}=\mbox{\boldmath$\mu$}+\mbox{\boldmath$L$}\mbox{\boldmath$z$}_{(s)}\ ,$$

$$\displaystyle\frac{1}{S}\sum_{s=1}^{S}\log p(\mbox{\boldmath$Y$}|\mbox{\boldmath$X$},\mbox{\boldmath$\mu$}+\mbox{\boldmath$L$}\mbox{\boldmath$z$}_{(s)},\beta)\ .$$

$$\displaystyle\frac{1}{S}\sum_{s=1}^{S}\log p(\mbox{\boldmath$Y$}|\mbox{\boldmath$X$},\mbox{\boldmath$\mu$}+\mbox{\boldmath$L$}\mbox{\boldmath$z$}_{(s)},\beta)=$$

$$\displaystyle\frac{1}{S}\sum_{s=1}^{S}\log\mathcal{N}(\mbox{\boldmath$Y$}|f(\mbox{\boldmath$X$};\mbox{\boldmath$\mu$}+\mbox{\boldmath$L$}\mbox{\boldmath$z$}_{(s)}),\beta^{-1}\mbox{\boldmath$I$}_{N})=$$

$$\displaystyle\frac{1}{S}\sum_{s=1}^{S}\frac{N}{2}\log(\beta)-\frac{\beta}{2}\|\mbox{\boldmath$Y$}-f(\mbox{\boldmath$X$};\mbox{\boldmath$\mu$}+\mbox{\boldmath$L$}\mbox{\boldmath$z$}_{(s)})\|^{2}$$

$$+const.$$

$$\displaystyle\frac{1}{2}\bigg{(}tr(\alpha\mbox{\boldmath$L$}^{T}\mbox{\boldmath$L$})+\alpha\mbox{\boldmath$\mu$}^{T}\mbox{\boldmath$\mu$}-M-\ln|\alpha\mbox{\boldmath$L$}\mbox{\boldmath$L$}^{T}|\bigg{)}\ .$$

$$\displaystyle{\mathcal{L}}_{(FS)}(\mbox{\boldmath$\mu$},\mbox{\boldmath$L$},\alpha,\beta,Z)=$$

$$\displaystyle\frac{1}{S}\sum_{s=1}^{S}\frac{N}{2}\log(\beta)-\frac{\beta}{2}\|\mbox{\boldmath$Y$}-f(\mbox{\boldmath$X$};\mbox{\boldmath$\mu$}+\mbox{\boldmath$L$}\mbox{\boldmath$z$}_{(s)})\|^{2}$$

$$\displaystyle-\frac{1}{2}\bigg{(}tr(\alpha\mbox{\boldmath$L$}\mbox{\boldmath$L$}^{T})+\alpha\mbox{\boldmath$\mu$}^{T}\mbox{\boldmath$\mu$}-\ln|\alpha\mbox{\boldmath$L$}\mbox{\boldmath$L$}^{T}|\bigg{)}\ .$$

$$\displaystyle\nabla_{\mbox{\boldmath$\mu$}}\mathcal{L}_{(FS)}(\mbox{\boldmath$\mu$},\mbox{\boldmath$L$},\alpha,\beta,Z)=$$

$$\displaystyle\frac{1}{S}\sum_{s=1}^{S}\beta(\nabla_{\mbox{\boldmath$w$}}f)_{(s)}^{T}(\mbox{\boldmath$Y$}-f(\mbox{\boldmath$X$};\mbox{\boldmath$\mu$}+\mbox{\boldmath$L$}\mbox{\boldmath$z$}_{(s)}))-\alpha\mbox{\boldmath$\mu$}\ ,$$

$$\displaystyle\nabla_{\mbox{\boldmath$L$}}\mathcal{L}_{(FS)}(\mbox{\boldmath$\mu$},\mbox{\boldmath$L$},\alpha,\beta,Z)=$$

$$\displaystyle\frac{1}{S}\sum_{s=1}^{S}\beta(\nabla_{\mbox{\boldmath$w$}}f)_{(s)}^{T}(\mbox{\boldmath$Y$}-f(\mbox{\boldmath$X$};\mbox{\boldmath$\mu$}+\mbox{\boldmath$L$}\mbox{\boldmath$z$}_{(s)}))\mbox{\boldmath$z$}_{(s)}^{T}$$

$$\displaystyle-\alpha\mbox{\boldmath$L$}+{\mbox{\boldmath$L$}^{+}}^{T}\ ,$$

$$\alpha=\frac{M}{\mbox{\boldmath$\mu$}^{T}\mbox{\boldmath$\mu$}+\mathrm{tr}(\mbox{\boldmath$L$}\mbox{\boldmath$L$}^{T})}\ ,$$

$$\beta=\frac{SN}{\sum_{s=1}^{S}\|\mbox{\boldmath$Y$}-f(\mbox{\boldmath$X$};\mbox{\boldmath$\mu$}+\mbox{\boldmath$L$}\mbox{\boldmath$z$}_{(s)})\|^{2}}\ .$$

$$f(w_{1},w_{2})=2{\mathcal{N}}(w_{1},w_{2}|0,\mbox{\boldmath$I$}_{2})\Phi(h(w_{1},w_{2}))\ ,$$

$$y=2\cos (x)\sin (x)-0.1x^2,$$

$$\mbox{\boldmath$\phi$}_{n}=[\phi(x_{n};r,c_{1})\dots\phi(x_{n};r,c_{M-1})\ 1]^{T}$$

$$\mbox{\boldmath$y$}=\mbox{\boldmath$W$}\mbox{\boldmath$x$}+\mbox{\boldmath$\xi$}+\mbox{\boldmath$\epsilon$}\ .$$

$${\left(\pi \gamma \left[1+{\left(\frac{x-x_0}{\gamma }\right)}^2\right]\right)}^{-1}.$$

Full text

1 Introduction

Bayesian inference is becoming the standard mode of inference as computational resources increase, algorithms advance and scientists across fields become aware of the importance of uncertainty. However, exact Bayesian inference is hardly ever possible whenever the model likelihood function deviates from mathematically convenient forms (i.e. conjugacy). Deterministic approximations are constantly gaining ground on the ubiquitous and computationally intensive Monte Carlo sampling methods that are capable of producing high quality approximations to otherwise intractable quantities such as posterior densities or marginal likelihoods.

Often, however, deterministic schemes are tailored to a particular model, or family of models, and hence previously derived methods might not be transferable to a new setting (e.g. in [4] a specialised algorithm for Bayesian inference in neural networks is considered). This introduces practical difficulties, mostly in fields beyond machine learning, whenever an implementation of Bayesian inference is required particularly in early explorative stages where the model formulation is likely to keep changing. In this work we introduce a scheme for approximate Bayesian inference that aims to be general enough so that it can accommodate a variety of models. The proposed scheme is conceptually simple and requires only that the gradient of the model log-likelihood with respect to its parameters be available.

Specifically, we consider the task of computing a global Gaussian approximation q(\mbox{\boldmathw})=\mathcal{N}(\mbox{\boldmathw}|\boldsymbol{\mu},\boldsymbol{\Sigma}) to a given intractable posterior distribution representing a probabilistic model by maximizing a standard variational lower bound [14]. This lower bound involves the expectation of the log-likelihood of the model distribution with respect to the approximating distribution q(\mbox{\boldmathw}). To enable the computation of these expectations either in closed form or through computationally tractable numerical approximations, the likelihoods are typically restricted to conditionally factorized forms. This step of the approach specifically depends on the model at hand.

By contrast, our conceptually simple method presented below is more generally applicable to various models in the same way. We demonstrate this by working out examples for a variety of models. In particular, we demonstrate empirically that our approach results in Gaussian approximations that are superior to the basic Laplace approximation [6], which is the typical objective of variational approximation schemes. A formal comparison to related state-of-the-art methods for computing improved Gaussian approximations, e.g. through the nested Laplace approximation [18] or by expectation propagation [9], is beyond the scope of this paper, however.

Our paper is organized as follows: in Section 2 we introduce in general terms the proposed scheme. Starting from a standard formulation, that typically variational methods use, we postulate a Gaussian posterior q(\mbox{\boldmathw}) and show how to form an approximation to the lower bound of the marginal log-likelihood. The obtained approximation allows us to optimise the variational parameters $\boldsymbol{\mu}$ and $\boldsymbol{\Sigma}$ of q(\mbox{\boldmathw}) by gradient optimisation. For the reader that wishes to refresh her/his memory or obtain a more detailed explanation of the equations presented in Section 2, we refer to [5, 21].

In Section 3 we demonstrate the proposed scheme on a number of applications and compare it against exact inference, Laplace and variational approximations. Specifically, in Subsection 3.1 we show a visual comparison of approximating flexible bivariate densities using the Laplace approximation and the proposed scheme. In Subsection 3.2 we apply our approach on the problem of Bayesian linear regression which actually does admit an analytical and exact solution. This is useful as it allows us to empirically verify the correctness of our scheme against the posterior obtained by exact inference. Subsequently, in Subsections 3.3 and 3.4 we compare the proposed scheme with approaches that take into account the functional forms of classification problems. We show that despite its general formulation, the proposed scheme performs up to par in this setting without exploiting any such problem specific knowledge. In Subsection 3.5 we show how a change in the model likelihood of probabilistic principal component analysis [20], that renders inference problematic, can easily be accommodated by the proposed scheme. This demonstrates the versatility of our approach in handling such cases in a direct and simple manner. Finally, in Subsection 3.6 we show how the proposed scheme can be applied beyond the usual statistical models, namely on a real geophysical model [7]. We believe that the proposed method raises a range of interesting questions and directions of research; we briefly discuss them in Section 4.

2 Proposed Scheme for Approximate Variational Inference

Assume an observed dataset of inputs \mbox{\boldmathX}=(x_{1},\dots,x_{N})^{T} and outputs \mbox{\boldmathY}=(y_{1},\dots,y_{N})^{T} modelled by a model $f$ parametrised by \mbox{\boldmathw}\in\mathbb{R}^{M}. For observed outputs corrupted by Gaussian noise of precision $\beta$, the following likelihood111Besides the likelihood based on the Gaussian density, others can also be accommodated as shown in Sections 3.3-3.5. The choice of the Gaussian is made for ease of exposition. arises:

[TABLE]

where f(\mbox{\boldmathX};\mbox{\boldmathw})=(f(x_{1};\mbox{\boldmathw}),\dots,f(x_{N};\mbox{\boldmathw})) is the vector of model outputs calculated on all inputs $x_{n}$. Furthermore, assume a Gaussian prior on the parameters $w$:

[TABLE]

Our wish is to approximate the true posterior of the parameters p(\mbox{\boldmathw}|\mbox{\boldmathY},\mbox{\boldmathX},\alpha,\beta). We do not make any assumptions about the model having conjugate priors for the parameters $w$. Model $f$ may have a complex functional form that hinders exact Bayesian inference, or even the application of an approximate Bayesian scheme such as VBEM [5] with a factorised prior. However, we do have to make an assumption on the form of the posterior. We choose to postulate an approximate Gaussian posterior for the parameters $w$:

[TABLE]

Parameters \mbox{\boldmath\mu}\in\mathbb{R}^{M} and \mbox{\boldmath\Sigma}\in\mathbb{R}^{M\times M} of the posterior are called variational parameters. For reasons that will become obvious shortly, we choose to parametrise covariance matrix as \mbox{\boldmath\Sigma}=\mbox{\boldmathL}\mbox{\boldmathL}^{T} with \mbox{\boldmathL}\in\mathbb{R}^{M\times M}. The postulated posterior now reads:

[TABLE]

Hence, the actual variational parameters are $\mu$ and $L$.

2.1 Approximate Lower Bound

The first step in introducing the proposed scheme, is writing the marginal log-likelihood and lower-bounding it in the standard way using Jensens’ inequality [5, Eq. (2.46)]:

[TABLE]

An alternative motivation of the lower bound is provided in [21, Eq. (15)] Maximising the lower bound $\mathcal{L}$ in Eq. (5) with respect to the free variational parameters $\mu$ and $L$ of q(\mbox{\boldmathw}), results in the best Gaussian approximation to the true posterior. Term $(1)$, the integrated likelihood in Eq. (5), is a potentially intractable integral. We approximate term $(1)$ using Monte Carlo sampling:

[TABLE]

where we draw $S$ samples \mbox{\boldmathw}_{(s)} from the postulated posterior q(\mbox{\boldmathw}). Due to the sampling, however, the variational parameters no longer appear in the approximation Eq. (6) . Nevertheless, it is possible to re-introduce them by rewriting the sampled weights \mbox{\boldmathw}_{(s)} as222This is where the parametrisation \mbox{\boldmath\Sigma}=\mbox{\boldmathL}\mbox{\boldmathL}^{T} becomes useful.:

[TABLE]

where variables \mbox{\boldmathz}_{(s)} are sampled from the standard normal \mbox{\boldmathz}\sim\mathcal{N}(\mbox{\boldmath0},\mbox{\boldmathI}_{M}). We summarise all samples \mbox{\boldmathz}_{(s)} by Z=\{\mbox{\boldmathz}_{(1)}\dots,\mbox{\boldmathz}_{(S)}\}. Hence, we can rewrite Eq. (6) as:

[TABLE]

Hence, the variational parameters $\mu$ and $L$ are now made explicit in this approximation. We expand the approximation of term $(1)$ further:

[TABLE]

Term $(2)$ in Eq. (5) is simply the Kullback-Leibler divergence (KLD) between the two Gaussian densities q(\mbox{\boldmathw}) and p(\mbox{\boldmathw}|\alpha), and can be calculated in closed form:

[TABLE]

We can now put together the approximated term $(1)$ in Eq. (9) and the KLD term $(2)$ in Eq. (10), to formulate the following objective function333The subscript $FS$ stands for finite sample. ${\mathcal{L}}_{(FS)}$. Discarding constants, the proposed approximate lower bound reads:

[TABLE]

Objective ${\mathcal{L}}_{(FS)}$ is an approximation to the intractable lower bound ${\mathcal{L}}$ in Eq. (5). It consists of two parts, the approximation to the integrated likelihood (term $(1)$) and the exact KLD (term $(2)$). The proposed lower bound ${\mathcal{L}}_{(FS)}$ becomes more accurate when the number $S$ of samples \mbox{\boldmathz}_{(s)} is large.

2.2 Optimisation of Approximate Lower Bound

Gradients of ${\mathcal{L}}_{(FS)}$ can be calculated with respect to the variational parameters $\mu$ and $L$ in order to find the best approximating posterior q(\mbox{\boldmathw}):

[TABLE]

[TABLE]

where \nabla_{\mbox{\boldmathw}}f denotes the Jacobian matrix of $f$ and \mbox{\boldmathL}^{+} is the pseudo-inverse of $L$ due to the derivation of the log-determinant444See derivative rule $55$ in [15]. in Eq. (11). Analogous equations for the case of exact variational inference can be found at [5, Eq. (2.64)]. Given the current posterior q(\mbox{\boldmathw}), hyperparameters $\alpha$ and $\beta$ have analytical updates:

[TABLE]

[TABLE]

Again, analogous equations for the above hyperparameter updates can be found in [21, Eqs. (38), (39)].

The proposed scheme is summarised with the pseudocode in Algorithm 1. Convergence was established in the experiments by checking whether the difference between the objective function values ${\mathcal{L}}_{(FS)}$ between two successive iterations is less than $10^{-4}$. Gradient optimisation of $\mu$ and $L$ was carried out using the scaled conjugate gradient algorithm [13]. The outcome of the above scheme is the approximation ${\mathcal{L}}_{(FS)}$ to the marginal log-likelihood \log p(\mbox{\boldmathY}|\mbox{\boldmathX},\alpha,\beta) (also called log-evidence). The scheme imparts us with the approximate Gaussian posterior q(\mbox{\boldmathw})={\mathcal{N}}(\mbox{\boldmathw}|\mbox{\boldmath\mu},\mbox{\boldmathL}\mbox{\boldmathL}^{T}).

2.3 Monitoring Generalisation Performance

For large values of $S$ the proposed lower bound $\mathcal{L}_{(FS)}$ approximates the true bound $\mathcal{L}$ in Eq. (5) closely. Therefore, we expect that optimising $\mathcal{L}_{(FS)}$ will yield approximately the same optimal variational parameters \mbox{\boldmath\mu},\mbox{\boldmathL} as the optimisation of the intractable true lower bound $\mathcal{L}$ would.

The proposed scheme exhibits some fluctuation as \mathcal{L}_{(FS)}(\mbox{\boldmath\mu},\mbox{\boldmathL},\alpha,\beta,Z) is a function of the random set of samples $Z$. Hence, if the algorithm, as summarised in Algorithm 1, is run again, a new set $Z^{(new)}$ will be drawn and a different function \mathcal{L}_{(FS)}(\mbox{\boldmath\mu},\mbox{\boldmathL},\alpha,\beta,Z^{(new)}) will be optimised. However, for large enough $S$ the fluctuation due to $Z$ will be innocuous and approximately the same variational parameters will be found for any drawn $Z$555Discounting other sources of randomness like initialisation, etc..

However, if on the other hand we choose a small $S$, then the variational parameters will overly depend on the small set of samples $Z$ that happened to be drawn at the beginning of the algorithm. As a consequence, $\mathcal{L}_{(FS)}$ will approximate $\mathcal{L}$ poorly, and the resulting posterior q(\mbox{\boldmathw}) will also be a poor approximation to the true posterior. Hence, the variational parameters will be overfitted to the small set of samples $Z$ that happened to be drawn.

Naturally, the question arises of how to choose a large enough $S$ in order avoid overfitting the variational parameters on $Z$. A practical answer to this question is the following: at the beginning of the algorithm we draw a second independent set of samples Z^{\prime}=\{\mbox{\boldmathz}^{\prime}_{(1)},\dots,\mbox{\boldmathz}^{\prime}_{(S^{\prime})}\} where $S^{\prime}$ is preferably a number larger than $S$. At each (or every few) iteration(s) we monitor the quantity \mathcal{L}_{(FS)}(\mbox{\boldmath\mu},\mbox{\boldmathL},\alpha,\beta,Z^{\prime}) on the independent666We stress that $Z^{\prime}$ is not used in training. sample set $Z^{\prime}$. If the variational parameters are not overfitting the drawn $Z$, then we should see that as the lower bound \mathcal{L}_{(FS)}(\mbox{\boldmath\mu},\mbox{\boldmathL},\alpha,\beta,Z) increases, the quantity \mathcal{L}_{(FS)}(\mbox{\boldmath\mu},\mbox{\boldmathL},\alpha,\beta,Z^{\prime}) should also display a tendency to increase. If on the other hand the variational parameters are overfitting the drawn $Z$, then though \mathcal{L}_{(FS)}(\mbox{\boldmath\mu},\mbox{\boldmathL},\alpha,\beta,Z) is increasing, we will notice that \mathcal{L}_{(FS)}(\mbox{\boldmath\mu},\mbox{\boldmathL},\alpha,\beta,Z^{\prime}) is actually deteriorating. This is a clear sign that a larger $S$ is required.

The described procedure is reminiscent of monitoring the generalisation performance of a learning algorithm on a validation set during training. A significant difference, however, is that while validation sets are typically of limited size, here we can set $S^{\prime}$ arbitrarily large. For practical purposes, we found that $S^{\prime}=5S$ was good enough to detect overfitting. An illustration of overfitting the variational parameters in provided in Sec. 3.2.

3 Applications

In this section we apply the proposed approach on a variety of applications, namely regression, classification, denoising and geophysical modelling. In particular, the geophysical example shows how the method can be applied beyond standard statistical models.

3.1 Fitting Bivariate Posteriors

We test our proposed scheme on some artificially constructed posteriors by using a flexible parametric class of densities, due to [3], which reads:

[TABLE]

where $\Phi$ is the cumulative distribution function of the standard normal distribution, and in general $h$ is a real-value function such that $h(-w)=-h(w)$. Here, we take $h$ to be the dot product of a row and column vector, h(w_{1},w_{2})=(w_{1},w_{2},w_{1}w_{2}^{2},w_{1}^{2}w_{2},w_{1}^{3},w_{2}^{3})\mbox{\boldmatha} as in [3]. The goal is to find the best Gaussian approximation to instances of Eq. (16) for different column vectors $a$. To that end, we tried to find the best Gaussian using the Laplace approximation and the proposed scheme. We used $S=50$. The results are shown in Fig. 1. The Gaussian approximations are drawn as black dashed curves with their mean marked as a red cross. The goodness of each approximation has been measured as the Kullback-Leibler divergence, and it is noted in the respective captions. We note that the proposed scheme fares better than the Laplace approximation as the latter evidently focuses on the mode of the target density instead on where the volume of the density lies. The KLD in these examples was calculated numerically as there is no closed form between a Gaussian and a member of the densities in Eq. (16)

3.2 Bayesian Linear Regression

Bayesian linear regression constitutes a useful example for corroborating that the proposed scheme works correctly as we can compare the obtained posterior q(\mbox{\boldmathw}) to the posterior obtained by exact Bayesian inference [6].

We consider data targets $y_{n}$ generated by the equation

[TABLE]

with inputs $x_{n}$ uniformly drawn in the range $[-6,6]$. We add white Gaussian noise to the data targets with a standard deviation of $\sigma=0.2$. We calculate a set of radial basis functions on the data inputs

[TABLE]

where $\phi(x_{n};r,c_{m})=\exp(-\frac{\|x_{n}-c_{m}\|^{2}}{2r^{2}})$. The last element $1$ in \mbox{\boldmath\phi}_{n} serves as a bias term. We set $r=1$, and adopt the linear model y=\mbox{\boldmath\phi}^{T}\mbox{\boldmathw} where \mbox{\boldmathw}\in\mathbb{R}^{M}. We complete the model by choosing the following densities:

• •

Likelihood: p(\mbox{\boldmathY}|\mbox{\boldmathX},\mbox{\boldmathw},\beta)=\prod_{n=1}^{N}\mathcal{N}(y_{n}|\mbox{\boldmath\phi}_{n}^{T}\mbox{\boldmathw},\beta^{-1}).

• •

Prior: p(\mbox{\boldmathw})=\mathcal{N}(\mbox{\boldmathw}|\mbox{\boldmath0},\alpha^{-1}\mbox{\boldmathI}_{M}) .

• •

Postulated posterior: q(\mbox{\boldmathw})=\mathcal{N}(\mbox{\boldmathw}|\mbox{\boldmath\mu},\mbox{\boldmathL}\mbox{\boldmathL}^{T}), where \mbox{\boldmath\mu}\in\mathbb{R}^{M} and \mbox{\boldmathL}\in\mathbb{R}^{M\times M}.

We set the number of samples of variables $z$ to $S=100$. We inferred the Gaussian posterior of the weights $w$ using both exact Bayesian inference [6] and the proposed scheme. In Fig. 2 we plot the mean predictions as obtained by the two schemes and note that they are very similar, especially in the areas where enough data items are present. Similarly, in Fig. 3 we plot the covariance matrices found by the two schemes and note their close similarity. Hence, we conclude that the proposed scheme stands in close agreement to the exact Bayesian solution.

Finally, we demonstrate on the same dataset the effect of overfitting the variational parameters when $S$ is set too low. In Fig. 4 we are monitoring the lower bounds \mathcal{L}_{(FS)}(\mbox{\boldmath\mu},\mbox{\boldmathL},\alpha,\beta,Z) and \mathcal{L}_{(FS)}(\mbox{\boldmath\mu},\mbox{\boldmathL},\alpha,\beta,Z^{\prime}), see Sec. 2.3. In the left, we run the algorithm for $S=100$ and $S^{\prime}=500$: we see that as \mathcal{L}_{(FS)}(\mbox{\boldmath\mu},\mbox{\boldmathL},\alpha,\beta,Z) increases with each iteration, so does \mathcal{L}_{(FS)}(\mbox{\boldmath\mu},\mbox{\boldmathL},\alpha,\beta,Z^{\prime}). This means that the fitted variational parameters \mbox{\boldmath\mu},\mbox{\boldmathL} generalise well. On the right hand side, we run the algorithm for $S=10$ but kept $S^{\prime}=500$. Here we clearly see that while \mathcal{L}_{(FS)}(\mbox{\boldmath\mu},\mbox{\boldmathL},\alpha,\beta,Z) increases, the lower bound \mathcal{L}_{(FS)}(\mbox{\boldmath\mu},\mbox{\boldmathL},\alpha,\beta,Z^{\prime}) is deteriorating. This a clear sign that a larger $S$ is required and that the variational parameters are overfitted.

3.3 Bayesian Logistic Regression

In this section we apply the proposed scheme to Bayesian logistic regression and compare with the variational approach presented in [11]. The data are input-label pairs (\mbox{\boldmathx},y) with $y\in\{0,1\}$. Again, like in Sec. 3.2, we calculate basis functions \mbox{\boldmath\phi}_{n} on the input data \mbox{\boldmathx}_{n} and take $r=0.5$. We set $S=200$. We complete the model by choosing the following densities:

• •

\fleqn[0pt]\begin{aligned} \mbox{Likelihood:\ }&p(\mbox{\boldmathY}|\mbox{\boldmathX},\mbox{\boldmathw})=\\ &\prod_{n=1}^{N}\sigma(\mbox{\boldmath\phi}_{n}^{T}\mbox{\boldmathw})^{y_{n}}(1-\sigma(\mbox{\boldmath\phi}_{n}^{T}\mbox{\boldmathw}))^{1-y_{n}}\ .\end{aligned}

• •

Prior: p(\mbox{\boldmathw})=\mathcal{N}(\mbox{\boldmathw}|\mbox{\boldmath0},\alpha^{-1}\mbox{\boldmathI}_{M}) .

• •

Postulated posterior: q(\mbox{\boldmathw})=\mathcal{N}(\mbox{\boldmathw}|\mbox{\boldmath\mu},\mbox{\boldmathL}\mbox{\boldmathL}^{T}) .

We evaluated both schemes on datasets preprocessed by Rätsch et al777http://www.raetschlab.org/Members/raetsch/benchmark. Each preprocessed dataset has been randomly partitioned into a 100 non-overlapping training and testing sets. Hence the performance of both schemes was evaluated as the accuracy on the test set, that is the ratio of correctly classified test samples over all test samples. The predictive distribution for the proposed scheme was approximated using a Monte Carlo estimate. We drew 200 parameter samples from the fitted Gaussian posterior $q$ and measured performance on the testing set as the average accuracy under each sample of parameters. The results reported in Table 1, mean squared error and standard deviation, show that both the proposed schemes and the variational bound in [11] perform virtually the same. We note that, as opposed to [11] which exploits the functional form of logistic regression in order to design a bespoke lower bound, the proposed method does not take into account any such knowledge and still is capable of delivering comparable performance. Hence, we find the results in this section encouraging.

3.4 Bayesian Multiclass Classification

In this section we apply the proposed scheme on Bayesian multiclass classification. The data are input-label pairs (\mbox{\boldmathx},\mbox{\boldmathy}). Vectors $y$ are binary vectors encoding class labels using $1$-of-$K$ coding scheme, e.g. $[0\ 1\ 0]$ encodes class label $2$ in a $3$-class problem. A typical way of formulating multiclass classification is multiclass logistic regression (MLR), see [6, Chapter $4$] for more details. MLR models the probability p(C_{k}|\mbox{\boldmath\phi}_{n}) of the $n$-th data item belonging to class $C_{k}$ via the softmax function p(C_{k}|\mbox{\boldmath\phi}_{n})=\frac{\exp(\mbox{\boldmath\phi}_{n}^{T}\mbox{\boldmathw}_{k})}{\sum_{\ell=1}^{K}\exp(\mbox{\boldmath\phi}_{n}^{T}\mbox{\boldmathw}_{\ell})}. $K$ denotes the total number of classes, and each class $C_{k}$ is associated with a weight vector \mbox{\boldmathw}_{k}. Similarly to logistic regression, MLR does not allow direct Bayesian inference as the use of the softmax function renders integrals over the likelihood term intractable. Thus, Bayesian MLR is a good candidate problem for the proposed approach. We specify the following model:

• •

\fleqn[0pt]\begin{aligned} \mbox{Likelihood:\ }&p(\mbox{\boldmathY}|\mbox{\boldmathX},\mbox{\boldmathw}_{1},\dots,\mbox{\boldmathw}_{K})=\\ &\prod_{n=1}^{N}\prod_{k=1}^{K}p(C_{k}|\mbox{\boldmath\phi}_{n})^{y_{nk}}\ .\end{aligned}

• •

Prior: \prod_{k=1}^{K}p(\mbox{\boldmathw}_{k})=\mathcal{N}(\mbox{\boldmathw}_{k}|\mbox{\boldmath0},\alpha^{-1}\mbox{\boldmathI}_{M}) .

• •

Postulated posterior: q(\mbox{\boldmathw}_{1},\dots,\mbox{\boldmathw}_{K})=\prod_{k=1}^{K}q(\mbox{\boldmathw}_{k}), with q(\mbox{\boldmathw}_{k})=\mathcal{N}(\mbox{\boldmathw}_{k}|\mbox{\boldmath\mu}_{k},\mbox{\boldmathL}_{k}\mbox{\boldmathL}_{k}^{T}) .

As a corroboration of the usefulness of our approximation to Bayesian MLR, we compare with the multiclass generalisation of the relevance vector machine (RVM) [19] presented in [16]. We use the multiclass UCI datasets suggested therein. Amongst the two generalisations of the RVM suggested in [16], we use the $\mbox{mRVM}_{2}$ version. We also follow the suggestion of [16] concerning the choice of kernels for the different datasets. We set $S=200$. Table 2 summarises the results of our numerical simulations along with details of the datasets. While the $\mbox{mRVM}_{2}$ designs a refined probabilistic model in order to make probabilistic multiclass classification possible, the proposed scheme does not take into account any kind of such knowledge and is still able to deliver competitive performance, in terms of predictive accuracy, as seen in Table 2. The good performance demonstrates both the usefulness and versatility of the proposed method.

3.5 Probabilistic Image Denoising

In this section we further demonstrate how the proposed method can take in its stride a change in the model-likelihood that complicates computations. The model considered here is the ubiquitous probabilistic principal component analysis (PPCA) introduced in [20]. PPCA assumes that the observed high-dimensional data \mbox{\boldmathy}\in\mathbb{R}^{d} are manifestations of low-dimensional latent variables \mbox{\boldmathx}\in\mathbb{R}^{q}, under a linear mapping expressed by a matrix \mbox{\boldmathW}\in\mathbb{R}^{d\times q} and an offset \mbox{\boldmath\xi}\in\mathbb{R}^{d}, corrupted by Gaussian noise $\epsilon$:

[TABLE]

PPCA formulates a computationally amenable linear-Gaussian model which allows integrating out the latent variables $x$ and obtaining the marginal log-likelihood. Estimating $W$ and $\mu$ follows by maximising the marginal log-likelihood [20]. Various works extend PPCA by replacing the noise model with other choices, e.g. [2] uses the Student-t distribution, in order to deal with different types of noise. A recent interesting suggestion is the choice of the Cauchy density as the noise model [22], albeit in a non-probabilistic formulation. The Cauchy density with location $x_{0}$ and scale $\gamma>0$ parameters reads:

[TABLE]

Choosing the Cauchy density as the noise model leads to a version of PPCA where the marginal log-likelihood is no longer tractable and so the latent variables $x$ cannot be integrated out. This is simply because the prior on $x$ is not conjugate to the Cauchy likelihood. However, the proposed method can be used to approximate this intractable marginal log-likelihood. Formally, we specify the following Cauchy-PPCA model:

• •

\fleqn[0pt]\begin{aligned} &\mbox{Likelihood:\ }p(\mbox{\boldmathY}|\mbox{\boldmathX},\mbox{\boldmathW},\mbox{\boldmath\xi},\gamma)=\\ &\prod_{n=1}^{N}\left(\pi\gamma\left[1+\left(\frac{\mbox{\boldmathy}_{n}-\mbox{\boldmathW}\mbox{\boldmathx}_{n}-\mbox{\boldmath\xi}}{\gamma}\right)^{2}\right]\right)^{-1}\ .\end{aligned}

• •

Prior: p(\mbox{\boldmathX})=\mathcal{N}(\mbox{\boldmathX}|\mbox{\boldmath0},\mbox{\boldmathI}_{N}) .

• •

Postulated posterior: q(\mbox{\boldmathX})=\mathcal{N}(\mbox{\boldmathX}|\mbox{\boldmath\mu},\mbox{\boldmathL}\mbox{\boldmathL}^{T}) .

Parameters $W$, $\xi$ and $\gamma$ are obtained by gradient optimisation of the proposed lower bound $\mathcal{L}_{(FS)}$.

We applied the original PPCA [20] and the proposed Cauchy-PPCA on a task concerning the denoising of images that have undergone pixel corruption. The aim here is to show the ease with which the proposed method can accommodate a change in the specification (i.e. change from Normal to Cauchy likelihood) and deliver a well-performing model.

The data $Y$ are $2,414$ face images of $38$ individuals obtained from the Extended Yale B Database [10]. There are $64$ images per individual under $9$ poses and $64$ illumination conditions. The images are grayscale images whose pixels have values between [math] and $255$. We rescale the images to $96\times 84$ pixels. Hence $d=96\times 84=8064$ and we set $q=2$, i.e. both PCA schemes project the images to a latent space of dimension equal to $2$. We corrupt $33.33\%$ of the pixels in each image by drawing a new value uniformly in the range $[0,\dots,255]$. For each individual, we use half of the corrupted images as the training set and the other half as the test set.

Fig. 5 presents results obtained on test images from $4$ individuals. The figure shows from left to right the original and corrupted test image followed by the two reconstructions obtained by PPCA and Cauchy-PPCA respectively. In order to quantify the quality of reconstruction, we use the following measure between the original and reconstructed images: \|\mbox{\boldmathy}_{\mbox{\tiny orig}}-\mbox{\boldmathy}_{\mbox{\tiny rec}}\|^{2}/\|\mbox{\boldmathy}_{\mbox{\tiny orig}}\|^{2}. This measure is quoted below each image. The results in Fig. 5 evidently show that Cauchy-PPCA achieves better denoising levels than PPCA. In actual fact, in our numerical experiments we found that Cauchy-PPCA outperformed PPCA on all $38$ individuals.

The present numerical experiment demonstrates the versatility of the proposed method in how it can easily extend PPCA to incorporate a Cauchy likelihood. This is achieved without exploiting any particular knowledge pertaining to the probabilistic specification of the model.

3.6 Bayesian Inference for the Stochastic Model by Boore

In this section we apply the proposed scheme on a geophysical model called the stochastic model. The stochastic model, due to Boore [7], is used to predict ground motion at a given site of interest caused by an earthquake. Ground motion is simply the shaking of the earth and it is a fundamental quantity in estimating the seismic hazard of structures. From a physical point of view, the stochastic model describes the radiated source spectrum and its amplitude changes in the frequency domain due to wave propagation from the source to the site of interest. The inputs to the stochastic model are the distance $R$ of the site of interest to the seismic source, the magnitude $M_{w}$ of the earthquake, and the frequency $f$ of ground motion. The stochastic model, in its simple form, has a parameter associated with the seismic source known in the literature as stress parameter ($\Delta\sigma$), two parameters associated with the path attenuation called geometrical spreading ($\eta$) and quality factor ($Q$), and one more parameter associated with the site called near-surface attenuation ($\kappa_{0}$). All aforementioned parameters are bounded within a physical range. In the case of multiple seismic sources, each source is associated with its own distinct stress parameter. The scalar output of the model $y$ is the mean Fourier amplitude of the ground motion. The type of ground motion we consider here is acceleration. We denote the stochastic model as a function y=g(M_{w},R,f;\mbox{\boldmathw}), where \mbox{\boldmathw}=[\Delta\sigma_{1},\dots,\Delta\sigma_{E},\eta,Q,\kappa_{0}], where $E$ is the number of seismic sources. This situation is depicted in Fig. 6. We refer the interested reader to [7] for more details. Estimating the posterior uncertainty of the model parameters is important in seismic hazard analysis as the propagation of uncertainty in the parameters can have an impact on the estimated hazard curve [17]. A discussion of how these posteriors can be utilised in later stages of seismic hazard analysis is beyond the scope of this work.

We specify the model by choosing the following densities:

• •

\fleqn[0pt]\begin{aligned} \mbox{Likelihood:\ }&p(\mbox{\boldmathY}|\mathcal{D},\mbox{\boldmathw})=\\ &\prod_{n=1}^{N}\mathcal{N}(y_{n}|g({M_{w}}_{n},R_{n},f_{n};\mbox{\boldmathw}),\sigma^{2})\ .\end{aligned}

• •

Flat prior: p(\mbox{\boldmathw})\propto 1 .

• •

Postulated posterior: q(\mbox{\boldmathw})=\mathcal{N}(\mbox{\boldmathw}|\mbox{\boldmath\mu},\mbox{\boldmathL}\mbox{\boldmathL}^{T}) .

In contrast to the previous applications, here we choose a very flat prior. Ground motion data inputs are denoted by $\mathcal{D}$ and targets by $Y$.

We performed experiments on a subset of the recently compiled RESORCE database [1] which is a pan-European database of strong-motion recordings. In particular we focused on data records originating from a station in the region of L’Aquilla for $E=8$ seismic sources. Hence, the total number of free model parameters in $w$ is $11$. We experimented with varying number of data records, $N\in\{100,200,500,1000\}$, in order to test the robustness of the Laplace and the proposed approximation in scenarios of limited data. Such situations arise in geophysical studies when data recordings are incomplete due to distortions in frequencies caused by instrumentation errors. The performance of Laplace and the proposed scheme was evaluated as the prediction error on test sets. Both schemes were run $10$ times, and each run involved a new random realisation of the training and testing set. Parameter $S$ was set to $1000$ for all experiments in this section. The predictive distribution for the proposed scheme was approximated using a Monte Carlo estimate. We drew $200$ parameter samples from the Gaussian fitted posterior $q$ and estimated performance on the testing set as the average of the mean squared error under each parameter sample. The results are reported in Table 3.

The results show that the proposed approximation fares better than Laplace, although at $N=1000$ the performances are virtually identical. For lower $N$, however, the Laplace approximation exhibits much higher variance than the proposed scheme.

4 Discussion and Conclusion

We have presented a scheme for Bayesian variational inference that is applicable in cases where the likelihood function renders more standard approaches difficult. The scheme is conceptually simple as it relies on a simple Monte Carlo average of the intractable part of the variational lower bound, see Eq. (5), and the re-introduction of the variational parameters resulting in the objective of Eq. (11). The scheme can thus be generally applied to other models where variational inference is difficult requiring only the gradients of the log-likelihood function with respect to the parameters.

In the numerical experiments we have shown that (a) the proposed scheme stands in close agreement with exact inference in Bayesian linear regression, (b) it performs up to par in classification tasks against methods that design bespoke model formulations, (c) it fares better than the Laplace approximation in a number of cases, and (d) it is very versatile and can be applied to a variety of problems.

Future work will address the relationship of our approach to variational approaches [18, 9] that provide alternative ways to compute improved Gaussian approximations to intractable posteriors relative to the Laplace approximation. Another aspect concerns ways to cope with very large problems that would require a large number of samples $S$ to obtain a sufficiently accurate approximation in Eq. (8). A natural choice would be to turn the scheme in Algorithm 1 into a recursive stochastic optimisation scheme [12] that employs small sample sets computed at each iteration, akin to stochastic gradient-based large-scale empirical risk minimisation [8]. These two approaches should not be confused, however. The latter employs subsampling of the data (\mbox{\boldmathX},\mbox{\boldmathY}) whereas our scheme generates samples based on current parameter estimates of the approximate posterior. Clearly, our scheme could incorporate sequential subsampling of the data as well. The problem of proving convergence of such an overall stochastic approximation approach in a suitable sense [12] seems to be open.

Acknowledgement

The RESORCE database [1] was used in this work with the kind permission of the SIGMA project888http://www.projet-sigma.com. N. Gianniotis was partially funded by the BMBF project “Potsdam Research Cluster for Georisk Analysis, Environmental Change and Sustainability”. C. Molkenthin and S. S. Bora were funded by the graduate research school GeoSim of the Geo.X initiative999http://www.geo-x.net.