Laplace's approximation

Laplace's approximation or the quadratic approximation (QUAP) provides an analytical expression for a posterior probability distribution by fitting a Gaussian distribution with a mean equal to the MAP solution and precision equal to the observed Fisher information. The approximation is justified by the Bernstein–von Mises theorem, which states that, under regularity conditions, the error of the approximation tends to 0 as the number of data points tends to infinity.

For example, consider a regression or classification model with data set $\{x_n,y_n{\}}_{n=1,\dots ,N}$ comprising inputs $x$ and outputs $y$ with (unknown) parameter vector $\theta$ of length $D$. The likelihood is denoted $p(𝐲|𝐱,\theta )$ and the parameter prior $p(\theta )$. Suppose one wants to approximate the joint density of outputs and parameters $p(𝐲,\theta |𝐱)$. Bayes' formula reads:

$p(𝐲,\theta |𝐱)=p(𝐲|𝐱,\theta )p(\theta |𝐱)=p(𝐲|𝐱)p(\theta |𝐲,𝐱)\simeq \tilde{q}(\theta )=Zq(\theta ).$

The joint is equal to the product of the likelihood and the prior and by Bayes' rule, equal to the product of the marginal likelihood $p(𝐲|𝐱)$ and posterior $p(\theta |𝐲,𝐱)$. Seen as a function of $\theta$ the joint is an un-normalised density.

In Laplace's approximation, we approximate the joint by an un-normalised Gaussian $\tilde{q}(\theta )=Zq(\theta )$, where we use $q$ to denote approximate density, $\tilde{q}$ for un-normalised density and $Z$ the normalisation constant of $\tilde{q}$ (independent of $\theta$). Since the marginal likelihood $p(𝐲|𝐱)$ doesn't depend on the parameter $\theta$ and the posterior $p(\theta |𝐲,𝐱)$ normalises over $\theta$ we can immediately identify them with $Z$ and $q(\theta )$ of our approximation, respectively.

Laplace's approximation is

$p(𝐲,\theta |𝐱)\simeq p(𝐲,\hat{\theta }|𝐱)\exp (-\frac{1}{2}(\theta -\hat{\theta }{)}^{\mathrm{\top }}{S}^{-1}(\theta -\hat{\theta }))=\tilde{q}(\theta ),$

where we have defined

$\begin{array}{cc}\hat{\theta }& ={\mathrm{argmax}}_{\theta }\log p(𝐲,\theta |𝐱),\\ S^{-1}& =-{{\mathrm{\nabla }}_{\theta }{\mathrm{\nabla }}_{\theta }\log p(𝐲,\theta |𝐱)|}_{\theta =\hat{\theta }},\end{array}$

where $\hat{\theta }$ is the location of a mode of the joint target density, also known as the maximum a posteriori or MAP point and $S^{-1}$ is the $D\times D$ positive definite matrix of second derivatives of the negative log joint target density at the mode $\theta =\hat{\theta }$. Thus, the Gaussian approximation matches the value and the log-curvature of the un-normalised target density at the mode. The value of $\hat{\theta }$ is usually found using a gradient based method.

In summary, we have

$\begin{array}{cc}q(\theta )& =𝒩(\theta |\mu =\hat{\theta },\mathrm{\Sigma }=S),\\ \log Z& =\log p(𝐲,\hat{\theta }|𝐱)+\frac{1}{2}\log |S|+\frac{D}{2}\log (2\pi ),\end{array}$

for the approximate posterior over $\theta$ and the approximate log marginal likelihood respectively.

The main weaknesses of Laplace's approximation are that it is symmetric around the mode and that it is very local: the entire approximation is derived from properties at a single point of the target density. Laplace's method is widely used and was pioneered in the context of neural networks by David MacKay, and for Gaussian processes by Williams and Barber.

Integrated nested Laplace approximation

Template:Vanchor (Template:Vanchor) is a method for approximate Bayesian inference based on Laplace's approximation. It is designed for a class of models called latent Gaussian models (LGMs), for which it can be a fast and accurate alternative for Markov chain Monte Carlo methods to compute posterior marginal distributions. Due to its relative speed even with large data sets for certain problems and models, INLA has been a popular inference method in applied statistics, in particular spatial statistics, ecology, seismology, and epidemiology. It is also possible to combine INLA with a finite element method solution of a stochastic partial differential equation to study e.g. spatial point processes and species distribution models. The INLA method is implemented in the R-INLA R package.

Latent Gaussian models

Let $𝒚=(y_1,\dots ,y_n)$ denote the response variable (that is, the observations) which belongs to an exponential family, with the mean ${\mu }_i$ (of $y_i$) being linked to a linear predictor ${\eta }_i$ via an appropriate link function. The linear predictor can take the form of a (Bayesian) additive model. All latent effects (the linear predictor, the intercept, coefficients of possible covariates, and so on) are collectively denoted by the vector $𝒙$. The hyperparameters of the model are denoted by $𝜽$. As per Bayesian statistics, $𝒙$ and $𝜽$ are random variables with prior distributions.

The observations are assumed to be conditionally independent given $𝒙$ and $𝜽$: $$\pi (𝒚|𝒙,𝜽)=\prod _{i\in ℐ}\pi (y_i|{\eta }_i,𝜽),$$ where $ℐ$ is the set of indices for observed elements of $𝒚$ (some elements may be unobserved, and for these INLA computes a posterior predictive distribution). Note that the linear predictor $𝜼$ is part of $𝒙$.

For the model to be a latent Gaussian model, it is assumed that $𝒙|𝜽$ is a Gaussian Markov Random Field (GMRF) (that is, a multivariate Gaussian with additional conditional independence properties) with probability density $$\pi (𝒙|𝜽)\propto {\left|{𝑸}_{𝜽}\right|}^{1/2}\exp (-\frac{1}{2}{𝒙}^T{𝑸}_{𝜽}𝒙),$$ where ${𝑸}_{𝜽}$ is a $𝜽$-dependent sparse precision matrix and $\left|{𝑸}_{𝜽}\right|$ is its determinant. The precision matrix is sparse due to the GMRF assumption. The prior distribution $\pi (𝜽)$ for the hyperparameters need not be Gaussian. However, the number of hyperparameters, $m=\mathrm{d}\mathrm{i}\mathrm{m}(𝜽)$, is assumed to be small (say, less than 15).

Approximate Bayesian inference with INLA

In Bayesian inference, one wants to solve for the posterior distribution of the latent variables $𝒙$ and $𝜽$. Applying Bayes' theorem $$\pi (𝒙,𝜽|𝒚)=\frac{\pi (𝒚|𝒙,𝜽)\pi (𝒙|𝜽)\pi (𝜽)}{\pi (𝒚)},$$ the joint posterior distribution of $𝒙$ and $𝜽$ is given by $$\begin{array}{cc}\pi (𝒙,𝜽|𝒚)& \propto \pi (𝜽)\pi (𝒙|𝜽)\prod _i\pi (y_i|{\eta }_i,𝜽)\\ & \propto \pi (𝜽){\left|{𝑸}_{𝜽}\right|}^{1/2}\exp (-\frac{1}{2}{𝒙}^T{𝑸}_{𝜽}𝒙+\sum _i\log [\pi (y_i|{\eta }_i,𝜽)]).\end{array}$$ Obtaining the exact posterior is generally a very difficult problem. In INLA, the main aim is to approximate the posterior marginals $$\begin{array}{ccc}\pi (x_i|𝒚)& =& \int \pi (x_i|𝜽,𝒚)\pi (𝜽|𝒚)d𝜽\\ \pi ({\theta }_j|𝒚)& =& \int \pi (𝜽|𝒚)d{𝜽}_{-j},\end{array}$$ where ${𝜽}_{-j}=({\theta }_1,\dots ,{\theta }_{j-1},{\theta }_{j+1},\dots ,{\theta }_m)$.

A key idea of INLA is to construct nested approximations given by $$\begin{array}{ccc}\tilde{\pi }(x_i|𝒚)& =& \int \tilde{\pi }(x_i|𝜽,𝒚)\tilde{\pi }(𝜽|𝒚)d𝜽\\ \tilde{\pi }({\theta }_j|𝒚)& =& \int \tilde{\pi }(𝜽|𝒚)d{𝜽}_{-j},\end{array}$$ where $\tilde{\pi }(\cdot |\cdot )$ is an approximated posterior density. The approximation to the marginal density $\pi (x_i|𝒚)$ is obtained in a nested fashion by first approximating $\pi (𝜽|𝒚)$ and $\pi (x_i|𝜽,𝒚)$, and then numerically integrating out $𝜽$ as $$\begin{array}{c}\tilde{\pi }(x_i|𝒚)=\sum _k\tilde{\pi }(x_i|{𝜽}_k,𝒚)\times \tilde{\pi }({𝜽}_k|𝒚)\times {\mathrm{\Delta }}_k,\end{array}$$ where the summation is over the values of $𝜽$, with integration weights given by ${\mathrm{\Delta }}_k$. The approximation of $\pi ({\theta }_j|𝒚)$ is computed by numerically integrating ${𝜽}_{-j}$ out from $\tilde{\pi }(𝜽|𝒚)$.

To get the approximate distribution $\tilde{\pi }(𝜽|𝒚)$, one can use the relation $$\begin{array}{c}\pi (𝜽|𝒚)=\frac{\pi (𝒙,𝜽,𝒚)}{\pi (𝒙|𝜽,𝒚)\pi (𝒚)},\end{array}$$ as the starting point. Then $\tilde{\pi }(𝜽|𝒚)$ is obtained at a specific value of the hyperparameters $𝜽={𝜽}_k$ with Laplace's approximation $$\begin{array}{cc}\tilde{\pi }({𝜽}_k|𝒚)& \propto {\frac{\pi (𝒙,{𝜽}_k,𝒚)}{{\tilde{\pi }}_G(𝒙|{𝜽}_k,𝒚)}|}_{𝒙={𝒙}^{\ast }({𝜽}_k)},\\ & \propto {\frac{\pi (𝒚|𝒙,{𝜽}_k)\pi (𝒙|{𝜽}_k)\pi ({𝜽}_k)}{{\tilde{\pi }}_G(𝒙|{𝜽}_k,𝒚)}|}_{𝒙={𝒙}^{\ast }({𝜽}_k)},\end{array}$$ where ${\tilde{\pi }}_G(𝒙|{𝜽}_k,𝒚)$ is the Gaussian approximation to $\pi (𝒙|{𝜽}_k,𝒚)$ whose mode at a given ${𝜽}_k$ is ${𝒙}^{\ast }({𝜽}_k)$. The mode can be found numerically for example with the Newton-Raphson method.

The trick in the Laplace approximation above is the fact that the Gaussian approximation is applied on the full conditional of $𝒙$ in the denominator since it is usually close to a Gaussian due to the GMRF property of $𝒙$. Applying the approximation here improves the accuracy of the method, since the posterior $\pi (𝜽|𝒚)$ itself need not be close to a Gaussian, and so the Gaussian approximation is not directly applied on $\pi (𝜽|𝒚)$. The second important property of a GMRF, the sparsity of the precision matrix ${𝑸}_{{𝜽}_k}$, is required for efficient computation of $\tilde{\pi }({𝜽}_k|𝒚)$ for each value ${𝜽}_k$.

Obtaining the approximate distribution $\tilde{\pi }(x_i|{𝜽}_k,𝒚)$ is more involved, and the INLA method provides three options for this: Gaussian approximation, Laplace approximation, or the simplified Laplace approximation. For the numerical integration to obtain $\tilde{\pi }(x_i|𝒚)$, also three options are available: grid search, central composite design, or empirical Bayes.

Further reading

• Template:Cite book
• Template:Cite book
• Template:Cite book