在机器学习问题中，很多时候无法确定一个概率分布的具体密度函数，因而在对这种分布进行后续操作（例如，贝叶斯学派求后验概率时）时难度很大，无法进行。为了简化问题经常需要对这种复杂分布进行近似，从而方便计算或操作。目前常用的近似算法主要有三种：拉普拉斯近似、变分近似、Gibbs采样。拉普拉斯近似便是一种简单且广泛应用的近似方法，并且是很多采样方法的基础思想，本文介绍拉普拉斯近似（Laplace Approximation）。
拉普拉斯近似的基本思想是使用一个高斯分布来近似复杂分布，求解过程即为求正态分布的期望$\mu$和方差$\sigma^2$(或精度$\frac{1}{\sigma^2}$)的过程。注意：该方法是用于对连续变量的概率密度函数进行近似的。
一维变量的情形
假设有一个变量$z$（标量），其分布为$p(z)$，该分布的具体形式未知，定义如下：
$$
p(z)=\frac{1}{Z}f(z) \\
Z=\int f(z) dz
$$
其中$Z$为归一化项。
使用拉普拉斯近似寻找一个正态分布$q(z)$来近似这个分布$p(z)$。只需要求出$q(z)$的期望$\mu$和方差$\sigma^2$(或精度$\frac{1}{\sigma^2}$)。
怎么在数学上定义这个近似呢？
在拉普拉斯近似中分布$q(z)$与分布$p(z)$近似定义为：$q(z)$的峰值与$p(z)$的峰值相等。即$q(z)$的期望等于$p(z)$最大值。
先观察一下单变量正态分布公式的形式：
$$
\mathcal{N}(x|\mu,\sigma^{2})=\frac{1}{(2\pi \sigma ^{2})^{1/2}} \exp \{-\frac{1}{2\sigma ^{2}}(x-\mu)^{2}\} \tag1
$$
指数项中包含$(x-\mu)^{2}$，在正态分布$q(z)$中的均值$\mu$应该为$p(z)$的最大值点。
假设$z_0$是$p(z)$的最大值点，则在$z_0$处的一阶导数为0：
$$
\frac{dp(z)}{dz} \bigg|_{z=z_0} = 0 \Leftrightarrow \frac{df(z)}{dz} \bigg|_{z = z_0} =0
$$
借用$\ln f(z)$在$z_0$处的二阶Taylor展开来构建指数中的平方项，如下：
$$
\begin {aligned}
\ln f(z) &\simeq \ln f(z_{0}) -\frac{1}{2} A (z-z_{0})^{2} \\
where: & A=-\frac{d^{2}}{dz^2} \ln f(z)\bigg|_{z=z_0}
\end {aligned} \tag 2
$$
note：Taylor展开的一阶项为0。
对式2两边去掉$\ln$， 则：
$$
f(z) \simeq f(z_{0})\exp \{-\frac{A}{2}(z-z_{0})^{2}\} \tag 3
$$
note：上式是未归一化的。
对比式3中的指数项与式1中的指数项可以发现，如果$A= \frac{1}{\sigma^2}$ ，则可得归一化的高斯函数$q(z)$为下式：
$$
q(z)=(\frac{A}{2\pi })^{1/2} \exp \{-\frac{A}{2}(z-z_{0})^{2\}}
$$
总结：
$q(z)$的
期望$\mu=z_0$ ，即，$f(z)$的最大值点；
精度$\frac{1}{\sigma^2}=A$，即，$\ln f(z)$在$z_0$处的二阶导数，由于$z_0$是极大值点，所以要求$A>0$。
多维变量的情形
多维的情形推导过程和一维的情形类似，简单推导如下。
多维变量的拉普拉斯近似的目的是寻找一个多维正态分布$q(\mathbf z)$来近似这个分布$p(\mathbf z)$。
首先看一下多维变量$\mathbf x$的高斯分布公式的形式:
$$
\mathcal{N}(\mathbf{x}|\mathbf{\mu },\mathbf{\Sigma })=\frac{1}{(2\pi )^{D/2}}\frac{1}{|\mathbf{\Sigma }|^{1/2}} \exp \{-\frac{1}{2}(\mathbf{x}-\mathbf{\mu })^{\mathrm{T}}\mathbf{\Sigma }^{-1}(\mathbf{x}-\mathbf{\mu })\} \tag 4
$$
其中$D$为$\mathbf x$的维度。
假设有一个变量$\mathbf z$（$M$维向量，粗体），其分布为$p(\mathbf z)$，该分布的具体形式未知，定义如下：
$$
p(\mathbf z)=\frac{1}{\mathbf Z}f(\mathbf z) \\
\mathbf Z=\int f(\mathbf z) dz
$$
其中$\mathbf Z$为归一化项。
令$\mathbf z_0$为$f(\mathbf z) $的驻点，
$$
\ln f(\mathbf{z}) \simeq \ln f(\mathbf{z}_{0})-\frac{1}{2}(\mathbf{z}- \mathbf{z}_{0})^{\mathrm{T}}\mathbf{A}(\mathbf{z}-\mathbf{z}_{0}) \tag 5
$$
其中$\mathbf A$为$M \times M$的矩阵，是$\ln f(\mathbf z)$在驻点$\mathbf z_0$处的海森矩阵(Hessian matrix)。
$$
\mathbf{A}= \nabla\nabla \ln f(\mathbf{z}) \bigg|_{\mathbf{z}= \mathbf{z_0}}
$$
对于公式5两边去掉$\ln$，则：
$$
f(\mathbf{z})\simeq f(\mathbf{z}_{0})\exp\{-\frac{1}{2}(\mathbf{z}-\mathbf{z}_{0})^{\mathrm T}\mathbf{A}(\mathbf{z}-\mathbf{z}_{0})\} \tag6
$$
对比式4中的指数项与式6中的指数项可以发现，如果$A= {\Sigma }^{-1}$ ，则可得归一化的高斯函数$q(\mathbf z)$为下式：
$$
q(\mathbf z)=\frac{A^{1/2}}{(2\pi )^{M/2}} \exp \{-\frac{1}{2}(\mathbf{z}-\mathbf{z}_0)^{\mathrm{T}}\mathbf{\Sigma }^{-1}(\mathbf{z}-\mathbf{z}_0)\}
$$
其中$|\mathbf A|$为矩阵$\mathbf A$的行列式。
总结
求一个未知分布的拉普拉斯近似的过程分为两个步骤：
- 求其分布函数的最大值(mode)，即，驻点。通常是通过一些数值优化算法求的。实际情况下会存在多峰情况的分布，那么可以对不同的波峰进行拉普拉斯近似。
- 求海森矩阵。
根据中心极限定理，随着数据量的增加，一个模型的后验分布会越来越与高斯分布相似，所以数据集越大，拉普拉斯近似的作用约理想。
参考资料
Pattern Recognition and Machine Learning