Laplace Approximation

Background : Taylor Expansion

Consider a function $f: \mathcal X \rightarrow \reals$ and $\mathbf x \in \mathcal X.$ Let's say we would to approximate $f(x)$ around a local maximum value of $f$. We start by finding the optimal value using calculus, i.e. $f'(x_0) = 0$; note also that at $x_0$ we have $f''(x_0) < 0$. We then do the Taylor expansion to the 2nd order around $x_0$:

Let's try to approximate $f(x) = e^{\frac{\sin x}{x}}$ . For this function, $x_0 = 1$ and $f''(x_0) = -1/3$. From derivation, we have

Fig 1: Approximating a function with Taylor expansion.

Laplace Approximation

In brief, Laplace (sometimes also called Gaussian) approximation applies the Taylor expansion on$f(w)$ to compute

such that $\int p(w) dw = 1$. In other words, $Z= \int_{-\infty}^\infty e^{f(w)} dw$ and $e^{f(w)}$ is a likelihood.

Let's assume that we only care about the behaviour of $p(w)$ around the mode $w_\text{MAP}$. Define $\gamma := - f''(x_0)$. For the previous section, we know that

Thus, we get

In other word, $p(w) \approx \hat p(w)= \mathcal N(w_\text{MAP}, \sigma^2= (1/\gamma^2))$.

Application 1: Approximating Posterior Distributions

Consider a dataset $\mathcal D = \{ ( \mathbf x_i, y_i)\}_{i=1}^n$ and $\mathbf x_i, \mathbf w \in \reals^d$. We want to estimate estimate the posterior

Using the identity $e^{\log a} =a$, we can write

where $f(\mathbf w) =- (\log p(\mathcal D | \mathbf w) + \log p(\mathbf w) )$ and $Z = p(\mathcal D)$. For the sake of completeness, we repeat the same derivation again but for the $d$ dimensional case. We perform the Taylor explanation on $f (\mathbf w)$ at the mode $\mathbf w_\text{MAP}$ , yielding

where the first order is zero by construction and $\mathbf H_{\mathbf w_\text{MAP}} := \nabla^2 f(\mathbf w_\text{MAP})$. Therefore, the posterior is

which is in the form of a Gaussian distribution. If we set

we get

Fig. 2: Two posterior approximations using Laplace approximation.

Proton Counter Problem (Mackay (2013), Ex. 27.1)

Define $r \in \mathbb N$ the number of protons measured and we want to infer their arrival rate $\lambda \in \reals_+$. We assume that $r$ follows a Poisson distribution, hence

We further assume that we have an improper prior, $p(\lambda) = 1/\lambda.$ Define

Taking derivate $f'(\lambda) \overset{!}{=} 0$ yields

We also have $f''(\lambda) = \frac{r-1}{\lambda^2}$, thus

Therefore, the approximated posterior is

where $\mu := \lambda_\text{MAP} = r-1$ and $\sigma^2 := 1/f''(\lambda_\text{MAP}) = \lambda_\text{MAP}^2$.

Application 2: Stirling’s formula

Consider the Euler Gamma function

$\forall t > 0$. We can see that $\Gamma(t+1) = t\Gamma(t)$ and $\Gamma(1) = 1$; in other words, this function is the factorial function

Although we know the formula for the factorial, our goal is to approximate the function with a closed form without explicitly computing the factorial. We start by defining a new variable

thus $x = e^y$ , $dx = e^y dy$, and $y \in \reals$. Using integration by substitution yields

where $g(y) = (t-1)y + y - e^y = ty - e^y$. Computing the first and second derivatives, we get

• $g'(y) = t - e^y$; setting it to zero yields the maximum value $y_0 := \ln t$.
• $g''(y) = e^y$ ; therefore, the second derivative at $y_0$ is $g''(y_0) = t$.

Therefore, using the Laplace approximation, we arrive at

With this approximation, we can compute

Fig. 3: Relative error of t! approximation

Conclusion

One obvious limitation of the Laplace approximation is inherit from the Taylor expansion that it captures only the local behaviour of the function near the root point: in the case of approximating the posterior distribution, this root point is the mode of the distribution.

I consulted the following materials while writing this article:

• Murphy (2012), "Machine Learning: A Probabilistic Perspective", Section 8.4
• Manfred Opper's Probabilistic Bayesian Modeling course, Lecture Laplace Approximation (TU Berlin, Summer 2020).
• Gundersen (2019), "Laplace's Method"

Figures in this article are made using Google Colab.