变分贝叶斯近似-编程知识

马尔可夫链蒙特卡洛方法（MCMC）是一个非常有用和重要的工具，但在用于估计大型数据集的复杂后验分布或模型时可能会遇到困难。变分近似（variational approximations）或变分推断（variational inference）可用于拟合贝叶斯模型（Jordan et al. 1999）。对于完全贝叶斯推断，变分近似通常比MCMC快得多，在某些情况下，它有助于估计其他方法无法估计的模型。作为一种确定性后验近似方法，变分近似保证了收敛性，很容易对收敛性进行评估。变分近似法最适合于其他方法无法解决的完全贝叶斯推断问题。当然，变分近似确实也有一些限制。它可能低估后验分布的变异性。这种被低估的可变性是变分近似的一个缺点，然而，它是直接可控的，取决于一组透明且易于修改的假设。

变分近似为贝叶斯模型的估计提供了一种不同的方法。与EM算法一样，变分近似也是确定性优化算法，他保证了收敛性，通过检查标量中的变化可以轻松地进行评估。与MCMC算法一样，变分近似估计完全后验分布，不需要额外的步骤来执行判断。下面描述变分近似的基础知识。

变分近似的目标是用某一个称为近似分布（approximating distribution）的工作分布 $q(\theta;\phi)$ 来近似后验分布 $p(\theta|X)$ （Bishop, 2006）。为了使这一近似尽可能接近，我们在近似分布的空间上搜索，找 $q$ 和 $p$ 之间最小的相对熵（Kullback-Leibler divergence, KL-divergence），也称KL散度（KL散度常被用来衡量两个概率分布之间的相似性或差异性。当KL散度为0时，表示两个分布完全相同；当KL散度大于0时，表示两个分布之间存在差异。）。形式上，我们在近似分布集合中搜索的分布 $q(\theta;\phi)$ 使下式最小（下面积分运算对离散 $\theta$ 情况应该为和运算）：

$KL(q(\theta;\phi)||p(\theta|X))=E_{q(\theta;\phi)}[\log q(\theta;\phi)-\log p(\theta|X)]\\ \\ =-\int q(\theta)\log\{\frac{p(\theta|X)}{q(\theta;\phi))}\}d\theta \quad......\ (1)$

或者更确切地说寻求 $\phi^*$ ，使得

$\phi^*=\mathop{\arg \min}\limits_{\phi\in\Phi}KL(q(\theta;\phi)||p(\theta|X))$

直接最小化KL是困难的（这里 $\Phi$ 为 $\phi$ 的空间）。因此，等价地，人们试图最大化证据下限（evidence lower bound, ELBO）。ELBO定义为：

$ELBO(\phi)=E_{q(\theta;\phi)}[\log p(X|\theta)+\log p(\theta)-\log q(\theta;\phi)\\\\=E_{q(\theta;\phi)}[\log p(X,\theta)-\log q(\theta;\phi)] \quad...... \ (2)$

因此，问题归结为寻求 $\phi^*$ ，使得

$\phi^*=\mathop{\arg \min}\limits_{\phi\in\Phi}ELBO(\phi)$

显然，当 $q(\theta;\phi)=p(\theta|X)$ 时，（1）式和（2）式被最大化。当然，这并不是特别有帮助，但可以作为选择 $q$ 的一种参考。我们在近似分布中引入了额外的假定以使得推断可以操作。这里使用类似于Jordan et al. (1999) 和 Bishop (2006) 所使用的分布形式。这里假定所谓的近似分布在参数之间的独立性，这种独立性对于真实后验分布可能并不存在。但除此之外没有任何其他假定。选择这种近似分布，也是因为它在应用于大量模型时已经被证明表现良好。Wang and Titterington (2004) 证明，在足够多的观测结果下，这种近似分布家族将正确地描述后验均值，而这对基于抽样的推断方法一般是不成立的。

我们首先将 $\theta$ 划分为 $K$ 块： $\theta=(\theta_1,\theta_2,...,\theta_n)$ ，近似分布在参数之间独立意味着

$q(\theta)=\prod_{k=1}^Kq(\theta_k)$

变分算法将识别（而不是假定）构成分解分布的每个组成部分的特定参数族。类似于EM算法的迭代，我们的步骤为下面两步的迭代（直到收敛）：

（1）假定有了第 $i$ 步的估计： $\{q^{(i)}(\theta_k)\}_{k=1}^K$

（2）为更新第 $k$ 个因子（ $k=1,2,...,K$ ），定义：

$E_{j\neq k}[\log p(\theta|X)]=\int \prod_{j\neq k}\log p(\theta,X)q^{(i)}(\theta_j)d\theta_j$

对 $k=1,2,...,K$ ，令

$q^{(i+1)}(\theta_k)=\frac{\exp\{E_{j\neq k}[\log p(\theta|X)]\}}{\int\exp\{E_{j\neq k}[\log p(\theta|X)]\}d\theta_k}$

变分推断可能非常耗时，并且在某些情况下是不可行的。关于这部分的细节请参看Blei et al.

最大化ELBO的另一种方法是自动微分变分推断（automatic differentation variational inference, ADVI）。这是一种基于梯度的方法，仍然使用迭代优化程序来获得 $\phi^*$ 。人们使用类似（随机）梯度下降的理念。这需要根据参数计算ELBO关于参数 $\phi$ 的导数。在复杂的模型中，这可能非常繁琐。在软件存在自动微分（AD）库的情况下，可以依赖库来准确计算导数。AD依赖成为“对偶数（dual numbers）”的变量表示来有效地计算梯度。

首先假设所有模型参数都是连续的。在ADVI中，ELBO被重写为

$ELBO(\phi)=E_{q(\zeta;\phi)}\{\log p[X,T^{-1}(\zeta)]+\log|\det J_{T^{-1}}(\zeta)|-\log q(\zeta;\phi)\}$

这里的 $T$ 是一个把 $\theta$ 转换到 $\zeta$ 的函数，其值域 $\zeta \in R^{\dim(\theta)}$ 。因此，需要将数据和模型参数的对数联合密度与雅可比行列式的对数绝对值相加。因为变换后的模型参数 $\zeta$ 的支撑在实轴，所以多元正态分布为一个合适的变分分布。

将多元正态作为变分分布允许我们使用蒙特卡洛估计在ELBO中计算期望及其梯度。也就是说，为估计ELBO，可以从变分正态分布中抽样，并评估上述期望值内的表达式。

为了最大化ELBO，需要ELBO相对于变分参数的梯度。即

$\nabla_{\phi}ELBO(\phi)\\\\=\nabla_{\phi}E_{q(\zeta;\phi)}\{\log p[X,T^{-1}(\zeta)]+\log|\det J_{T^{-1}}(\zeta)|-\log q(\zeta;\phi)\}\\\\\approx \log p[X,T^{-1}(\hat \zeta)]+\log|\det J_{T^{-1}}(\hat \zeta)|-\log q(\hat \zeta;\phi)$