机器学习之变分推断

一、前言

概率模型的中心任务就是求潜变量的后验概率分布$p(Z|X)$，以及关于该后验的期望。前者是推断，后者是决策。
为什么决策是关于后验的期望？ 见如下推导：
$$\begin{array}{rl}p(\tilde{x}|X)& ={\int }_{Z}p(\tilde{x}|Z)p(Z|X)\mathrm{d}Z\\ & ={\mathbb{E}}_{p(Z|X)}[p(\tilde{x}|Z)]\end{array}$$由此可见，对于后验概率分布的求解是十分重要的。但是，实际上后验分布很难求，原因主要有两点：潜在空间维度过高；后验分布形式复杂。具体地，对于连续型随机变量，积分可能没有解析解，同时由于空间维度过高和被积函数复杂导致没有数值解；对于离散型随机变量，理论上是可以计算的，但是隐含状态的数量可能有指数多个，从而导致精确计算的计算代价过大。
因此精确推断一般是不可行的，需要使用近似推断。大体分为随机近似和确定近似：
（1）随机近似 ：这种近似方法的特点是假如我们能给无限的计算资源，他可以生成精确的结果。马尔科夫链蒙特卡洛方法就是一种随机近似。
（2）确定近似：这种方法的特点是永远无法生成精确的解。该方法一般假设后验分布可以通过一种特定的方式分解，或者有一个具体的参数形式（例如，高斯分布）。本节介绍的变分推断就是一种确定近似。

二、变分推断

在进行推导之前，需要明确两个概念：泛函和变分。
泛函： 简单来说，泛函就是一个映射，该映射接受一个函数作为输入，返回一个值作为输出。机器学习中一个典型的例子就是熵$H[p]$。
$$H[P]=-\int p(x)\ln p(x)\mathrm{d}x$$变分法: 传统的微积分中，我们讨论的是$x$值的微小变化对函数$y(x)$的影响。在变分法中，我们关心的是函数$y(x)$的微小变化对泛函$F(y)$的影响。从而，在变分法中，我们可以寻找一个$y(x)$来最大化或者最小化泛函$F(y)$。

---

现在，我们利用变分法的思想，将后验概率分布看做是一类函数，我们的目的在于寻找一个后验概率分布，使得其对应的泛函最大或最小。从而我们关注，这个泛函具体是什么？ 很自然地，在概率角度，我们一般都是希望对数似然能够最大，所以我们先将对数似然作为泛函，然后看看是否能够进一步化简。
我们定义$X$为观测变量，$Z$为潜在变量和参数，推导如下：
$$\begin{array}{rl}\log p(X)& =\log p(X,Z)-\log p(Z|X)\\ & =\log \frac{p(X,Z)}{q(Z)}-\log \frac{p(Z|X)}{q(Z)}\end{array}$$等式两边同时对$q(Z)$求期望
$$左边={\int }_{Z}q(Z)\log p(X)\mathrm{d}Z=\log p(X){\int }_{Z}q(Z)\mathrm{d}Z=\log p(X)$$$$\begin{array}{c}右边=\underset{ELBO}{\underbrace{{\int }_{Z}q(Z)\log \frac{p(X,Z)}{q(Z)}\mathrm{d}Z}}\underset{KL(q||p)}{\underbrace{-{\int }_{Z}q(Z)\log \frac{p(Z|X)}{q(Z)}\mathrm{d}Z}}\end{array}$$因此
$$\log p(X)=ELBO+KL(q||p)$$
由KL Divergence的性质可知，其大于等于0。当且仅当两个分布相等时，KL Divergence为0。所以，当KL Divergence最小时，$ELBO$最大。那么换句话说，只要我们取一个$q(Z)$使得$ELBO$最大，那么此时的KL Divergence最小，即可以认为$q(Z)$接近于$p(Z|X)$。从而，我们的泛函由对数似然变为了$ELBO$。
即，此时我们的目标为
$$q(Z)=\underset{q}{\mathrm{argmax}}ELBO$$从而，我们需要限制$q(Z)$的类别范围，使得其能够被我们计算，并使得这个范围足够大，能够对后验分布进行较好的近似。限制范围的方法有很多，本文只介绍经典的变分推断。

---

我们假设随机变量$Z$可以划分为若干个互不相关、彼此独立的组，即 Z = { Z i ∣ i = 1 , 2 , . . . , M } Z=\{Z_i|i=1,2,...,M\} Z={Zi​∣i=1,2,...,M}。该假设的这种分解形式在统计物理学中对应于平均场理论。
则$q(Z)$可以表示为
$$q(Z)=\prod _{i=1}^M{q}_i$$式中，我们将$q_i(Z_i)$记为$q_i$。
因此，我们现在目的是找一系列的$q_i$使得$EBLO$达到最大。在变分推断中，我们对每个$q_i$进行最优化，从而完成整体的最优化。因此，我们在一次优化过程中，只需要关注一个因子就可以，假设当前我们关注的因子是$q_j$，我们将$ELBO$进行如下化简
$$\begin{array}{rl}ELBO& ={\int }_Z\prod _{i=1}^M{q}_i[\log p(X,Z)-\log \prod _{k=1}^M{q}_k]\mathrm{d}Z\\ & \begin{array}{c}=\underset{Left}{\underbrace{{\int }_Z{q}_j[\prod _{i\ne j}{q}_i\log p(X,Z)]\mathrm{d}Z}}-\underset{Right}{\underbrace{{\int }_Z\prod _{i=0}^M{q}_i\sum _{k=1}^M\log q_k\mathrm{d}Z}}\end{array}\end{array}$$我们分别对$Left$和$Right$进行计算
$$\begin{array}{rl}Left& ={\int }_{Z_j}{q}_j[{\int }_{Z_{i\ne j}}\prod _{i\ne j}{q}_i\log p(X,Z)\mathrm{d}{Z}_{i\ne j}]\mathrm{d}{Z}_j\\ & ={\int }_{Z_j}{q}_j{\mathbb{E}}_{i\ne j}[\log p(X,Z)]\mathrm{d}{Z}_j\\ Right& ={\int }_Z\prod _{i=0}^M{q}_i\log q_j\mathrm{d}Z+C_1\\ & ={\int }_{Z_j}{q}_j\log q_j\mathrm{d}{Z}_j\prod _{i\ne j}{\int }_{Z_i}{q}_i\mathrm{d}{Z}_i+C_1\\ & ={\int }_{Z_j}{q}_j\log q_j\mathrm{d}{Z}_j+C_1\end{array}$$其中，$C_1=-{\int }_Z{\prod }_{i=0}^M{q}_i{\sum }_{i\ne j}\log q_j\mathrm{d}Z$
同时，我们令$\log \tilde{p}(X,Z_j)={\mathbb{E}}_{i\ne j}[\log p(X,Z)]+C_2$（注意这里的$C_2$是用于归一化的，因为等式左边是对一个概率取对数，所以需要保证概率的性质，后面也会展示$C_2$的值），则$ELBO$最终化为
$$\begin{array}{rl}ELBO& ={\int }_{Z_j}{q}_j{\mathbb{E}}_{i\ne j}[\log p(X,Z)]\mathrm{d}{Z}_j-{\int }_{Z_j}{q}_j\log q_j\mathrm{d}{Z}_j-C_1\\ & ={\int }_{Z_j}{q}_j\log \frac{\tilde{p}(X,Z_j)}{q_j}\mathrm{d}{Z}_j+C_3\\ & =-KL(q_j||\tilde{p}(X,Z_j))+C_3\end{array}$$显然，当$q_j$逼近于$\tilde{p}(X,Z_j)$时，$ELBO$能够取得最大值，即
$$\log q_j^{\ast }(Z_j)={\mathbb{E}}_{i\ne j}[\log p(X,Z)]+C_2$$式中，$q_j^{\ast }(Z_j)$为当前时刻$q_j$的最优解。
注意到，上式中有一个归一化常数，我们推导一下他的表达式：
q j ∗ ( Z j ) = e x p { E i ≠ j [ l o g   p ( X , Z ) ] + C 2 } ∫ Z j q j ∗ ( Z j ) d Z j = e x p { C 2 } ∫ Z j e x p { E i ≠ j [ l o g   p ( X , Z ) ] } d Z j C 2 = l o g 1 ∫ Z j e x p { E i ≠ j [ l o g   p ( X , Z ) ] } d Z j \begin{align*} q_j^*(Z_j) &=exp\{ \mathbb{E}_{i \neq j}[{\rm log}\,p(X,Z)]+C_2\} \\ \int_{Z_j} q_j^*(Z_j) {\rm d}Z_j &= exp\{ C_2 \} \int_{Z_j} exp\{ \mathbb{E}_{i \neq j}[{\rm log}\,p(X,Z)]\}{\rm d}Z_j \\ C_2&= {\rm log} \frac{1}{ \int_{Z_j} exp\{ \mathbb{E}_{i \neq j}[{\rm log}\,p(X,Z)]\}{\rm d}Z_j} \\ \end{align*} qj∗​(Zj​)∫Zj​​qj∗​(Zj​)dZj​C2​​=exp{Ei=j​[logp(X,Z)]+C2​}=exp{C2​}∫Zj​​exp{Ei=j​[logp(X,Z)]}dZj​=log∫Zj​​exp{Ei=j​[logp(X,Z)]}dZj​1​​至此$q_j$求解完毕。

---

最后，我们使用坐标上升（Coordinate Ascent）法更新所有的因子。即

1. 初始化所有的$q_i$，令$t=0$
2. 第$t$次迭代：
   固定$q_2^{(t-1)},q_3^{(t-1)},...,q_M^{(t-1)}$，计算$q_1^{(t)}$
   固定$q_1^{(t)},q_3^{(t-1)},...,q_M^{(t-1)}$，计算$q_2^{(t)}$
   固定$q_1^{(t)},q_2^{(t)},q_4^{(t-1)}...,q_M^{(t-1)}$，计算$q_3^{(t)}$
   $⋮$
   固定$q_1^{(t)},q_2^{(t)},...,q_{M-1}^{(t)}$，计算$q_M^{(t)}$
3. 回到第2步，进行第$t+1$次迭代，直至收敛
   该算法能够保证收敛，因为下界关于每个$q_i$都是一个凸函数。

---

但是经典的基于平均场理论的变分推断仍然存在问题：

1. 平均场理论假设过强
2. 期望的计算在某些情况下仍然是不可行的

参考资料

《Pattern Recognition and Machine Learning》
机器学习（白板推导系列）

---

如有错误，感谢指正！！！