21 变分推断¶

21.1 导论¶

本书已经探讨了多种计算后验分布及其函数的算法。例如：对于离散型变量的概率图，可以使用链接树算法（JTA: junction tree algorithm）来进行精确推断（见第 20 章）。但该算法时间复杂度是链接树宽度的指数函数，导致确切推断常常是不可用。而对于连续型变量的高斯图模型，确切推断的时间复杂度是树宽的三次函数。此外，该算法对多维变量会变得奇慢无比，而且对非高斯的连续型变量以及离散 / 连续混合型变量束手无策。

对于某些简单共用 $x->D$ 形式的的双节点图模型，可在先验分布 $p(x)$ 和似然函数共轭的情况下，计算其后验分布 $p(x|D)$ 的确切闭式解（也称解析解），更多例子见第 5 章。（注：本章中 $x$ 代表未知变量，而在第五章中用的是 $\theta$ 表示未知模型参数）。

更一般的情况下，必须使用近似的推断方法。在第 8.4.1 节，我们曾讨论过高斯近似在非共轭双节点概率图中的应用，在 8.4.3 节中也讨论了其在逻辑回归上的应用。高斯近似方法很简单，但其无法很好地模拟一些复杂的后验分布。例如：在推断多项分布的参数时，狄利克雷分布会是更好的选择；推断离散图模型的状态变量时，类别分布会是一个更好的选择。

在本章节中，我们将考虑一个更加一般的确定性近似推断算法 — 变分推断法。其基本理念是：从易处理、可参数化的某个分布族中，选取一个分布 $q (x)$，然后设法让 $q (x)$ 尽可能地接近真正的后验 $p^*(x)\triangleq p (x|D)$ ，使得原本复杂后验的推断问题简化成了该分布族的参数最优化问题。通过放宽限制或对目标函数作近似，可以在精度和速度之间寻找平衡。因此，变分推断可以用最大后验估计算法的速度实现贝叶斯方法的统计学优越性。

21.2 变分推断法¶

假设 $p^*(x)$ 是真实却难以处理的分布，而 $q (x)$ 是某个便于处理的近似分布（例如：一个可被因子分解的分布族，或多维高斯分布）。假设 $q$ 被可由某些参数控制，并且可通过优化这些参数使得 $q$ 更加像 $p^*$。显然可以最小化损失函数 KL 散度：

\[ \mathbb{KL}(p^*||q) = \sum_{x} p^*（x） \log \frac {p^*(x)}{q (x)} ~~\tag {21.1}\label {eqn:21.1} \]

该公式非常难算，因为在分布 $p^*$ 上求期望回归本身就是我们要求解的问题。一个替代方案是最小化 逆 KL 散度``：

\[ \mathbb{KL}(q||p^*) = \sum_{x} q (x) \log \frac {q (x)}{p^*(x)} ~~\tag {21.2}\label {eqn:21.2} \]

该目标函数优势是：由于 $q$ 的形式可以提前假设，因此在其上的期望是可计算的。有关两个散度的区别，会在第 21.2.2 节讨论。

不幸的是，公式 21.2 仍然没有那么好算，因为按照贝叶斯公式，即便逐点计算后验分布 $p (x|D)$ 也很困难，因为归一化常数 $Z=p (D)$ 的计算非常棘手。好在未归一化的后验分布 $\tilde {p}(x)\triangleq p (x,D)=p^*(x) Z$ 很好计算。所以将目标函数改为如下：

\[ J(q) \triangleq \mathbb{KL}(q \| \tilde{p}) ~~\tag {21.3}\label {eqn:21.3} \]

当然该写法有点滥用记号的意思，因为 $\tilde {p}$ 严格意义上并非一个概率分布。不过无所谓，让我们带入 KL 散度的定义：

\[\begin{split} \begin{aligned} J(q) &=\sum_{\mathbf{x}} q(\mathbf{x}) \log \frac{q(\mathbf{x})}{\tilde{p}(\mathbf{x})} \\ &=\sum_{\mathbf{x}} q(\mathbf{x}) \log \frac{q(\mathbf{x})}{Z p^{*}(\mathbf{x})} \\ &=\sum_{\mathbf{x}} q(\mathbf{x}) \log \frac{q(\mathbf{x})}{p^{*}(\mathbf{x})}-\log Z \\ &=\mathbb{K} \mathbb{L}\left(q \| p^{*}\right)-\log Z \end{aligned} \end{split}\]

因为对于确定的 $p^*$， $Z$ 是常数，所以最小化 $J (q)$ 也就达到了使 $q$ 趋近 $p^*$ 的目的。因为 KL 散度总是非负的，可以看出 $J (q)$ 是负对数似然的上界：

\[ J (q) = \mathbb{KL} (q||p^*) - \log Z \ge -\log Z = -\log p（D） ~~\tag {21.8}\label {eqn:21.8} \]

换句话说，可以尝试最大化如下能量泛函的量。它同时也是对数边缘似然 $\log p(D)$ 的下界： $$ L (q) \triangleq -J (q) = - \mathbb{KL} (q||p^*) + \log Z \le \log Z = \log p (D) ~~\tag {21.9}\label {eqn:21.9} $$

该界在 $q=p^*$ 时是紧致的，由此可看出变分推断法和 EM 算法联系之紧密（见第 11.4.7 节）。

21.2.1 变分目标函数的其他意义¶

前述目标函数还有几种同样深刻的写法。其中一种如下：

\[ J (q) = \mathbb{E}_q [\log q (x) ] + \mathbb{E} [-\log \tilde {p}(x)] = - \mathbb {H}(q) + \mathbb{E}_q [E (x)] ~~\tag {21.10}\label {eqn:21.10} \]

也就是能量的期望（因为 $E (x)=-\log \tilde { p} (x)$）减去系统的熵。在统计物理里，$J (q)$ 被称为变分自由能。

另一写法如下：

\[\begin{split} \begin {align*} J (q) &= \mathbb{E}_q { [\log q (x) - \log p (x) p (D|x) ] } ~~\tag {21.11}\label {eqn:21.11} \\ &= \mathbb{E}_q { [\log q (x) - \log p (x) - \log p (D|x) ] } ~~\tag {21.12}\label {eqn:21.12} \\ &= \mathbb{E}_q [ -\log p (D|x) ] + \mathbb{KL} (q (x)||p (x)) ~~\tag {21.13}\label {eqn:21.13} \end {align*} \end{split}\]

也就是负对数似然的期望，加上一个表示后验分布到确切先验距离的惩罚项。

21.2.2 前向还是后向 KL 散度？¶

因为 KL 散度是非对称的，最小化 $\mathbb{KL} (q||p)$ 和 $\mathbb{KL} (p||q)$ 会给出不同的结果。

首先考虑后向 KL， $\mathbb{KL} (q||p)$，也称 I-投影 或 信息投影。根据定义有

\[ \mathbb{KL} (q||p) = \sum_x q (x) \ln {q (x) \over p (x)} ~~\tag {21.14}\label {eqn:21.14} \]

该量在 $p (x)=0$ 且 $q (x)>0$ 是无穷的。因此如果 $p (x)=0$ 就必须要有 $q (x) = 0$ ，后向 KL 被称作是 迫零的，而且近似分布 $q$ 常常会欠估计 $p$ 的支撑集。

其次考虑前向 KL, 也称 M-投影 或者 矩投影

\[ \mathbb{KL} (p||q) = \sum_x p (x) \ln {p (x) \over q (x)} ~~\tag {21.15}\label {eqn:21.15} \]

该量在 $q (x)=0$ 且 $p (x)>0$ 是无穷的。因此如果 $p (x)>0$ 就必须要有 $q (x) > 0$。前向 KL 被称作 避零的，且近似分布 $p$ 常常会过估计 $p$ 的支撑集。

两者的区别请见图 21.1。可以发现当真实分布 $p$ 是多模态的时候，前向 KL 是一个很差的选择（此处假设使用了一个单模态的 $q$）, 因为它给出的后验众数和平均数会落在一个低密度区域，恰恰落在两个模态的峰值之间。而后向 KL 不仅更便于计算，也具有更好的统计性质。

图 21.1 说明了双峰分布的正向KL与反向KL。蓝色曲线是真实分布 $p$ 的轮廓线，红色曲线是单峰近似 $q$ 的轮廓线。(a) 最小化前向 KL：$q$ 趋向于“覆盖” $p$ 。(b-c) 最小化反向 KL：$q$ 锁定在两种模式之一。基于（Bishop 2006b）的图10.3。

图 21.2 说明了对称高斯上的正向 KL 和反向 KL。蓝色曲线是真实分布 $p$ 的轮廓，红色曲线是因子分解的近似 $q$ 。(a) 最小化后向 $KL(q||p)$ 的轮廓。(b) 最小化前向 $KL(p||q)$。基于（Bishop 2006b）的图10.2。

另一个区别显示在图 21.2 中。此处真实分布是一个拉长的二维高斯分布而近似分布是两个一维高斯的乘积。换句话说 $p (x) = \mathcal {N}(x|\mu,\Lambda^{-1})$，且 $$ \mu = \begin {pmatrix} \mu_1 \\ \mu_2 \end {pmatrix} , \Lambda = \begin {pmatrix} \Lambda_{11} & \Lambda_{12} \\ \Lambda_{21} & \Lambda_{22} \end {pmatrix} $$

在图 21.2a 中给出了最小化后向 KL（ $\mathbb{KL} (q||p)$）的结果。在该例中，可以证明有如下解：

\[\begin{split} \begin {align*} q (x) &= \mathcal {N}(x_1|m_1,\Lambda_{11}^{-1}) \mathcal {N}(x_2|m_2,\Lambda_{22}^{-1}) ~~\tag {21.17}\label {eqn:21.17} \\ m1 &= \mu_1 - \Lambda_{11}^{-1} \Lambda_{12} (m_2 - \mu_2) ~~\tag {21.18}\label {eqn:21.18} \\ m2 &= \mu_2 - \Lambda_{22}^{-1} \Lambda_{21} (m_1 - \mu_1) ~~\tag {21.19}\label {eqn:21.19} \end {align*} \end{split}\]

图 21.2a 的结果显示我们正确地估计了平均值，但是该近似太紧凑了：其方差是由 $p$ 的方差中最小的那个方向决定的。事实上，通常情况下，在 $q$ 被因子化后，最小化 $\mathbb{KL} (q||p)$ 会给出一个过于置信的近似。

图 21.2b 给出了最小化前向 $\mathbb{KL} (p||q)$ 的结果。在习题 21.7 中，我们证明对一个因子化的近似，最小化正向 KL 给出的最优解正好是其各个因子边缘分布的乘积，也就是说有：

\[ q (x) = \mathcal {N}(x_1 | \mu_1,\Lambda_{11}^{-1}) \mathcal {N}(x_2 | \mu_2,\Lambda_{22}^{-1}) ~~\tag {21.20}\label {eqn:21.20} \]

图 21.2b 显示出该估计太宽泛了，因为它过估计了 $p$ 的支撑集。

在本章剩余部分，以及本书接下来大部分章节，我们会专注于最小化 后向 KL （ $\mathbb{KL} (q||p)$ ）。在第 22.5 节对期望传播算法的介绍中，会探讨 前向 KL （$\mathbb{KL} {p||q}$ ）在局部优化中的作用。

21.2.3 其他相关的度量¶

通过引入参数 $\alpha \in \mathbb {R}$ ，我们还可以定义如下 $\alpha$ 散度族：

\[ D_{\alpha}( p || q ) \triangleq \frac {4}{1-\alpha^2} \left ( 1 - \int p (x) ^{(1+\alpha)/2} q (x) ^{(1-\alpha)/2} dx \right) ~~\tag {21.21}\label {eqn:21.21} \]

该量满足 $D_\alpha (p || q) =0 \iff p=q$，但是显然也是不对称的，因而不是同一个度量。 $\mathbb{KL} (p||q)$ 对应极限 $\alpha \rightarrow 1$ ，而 $\mathbb{KL} (q||p)$ 对应极限 $\alpha \rightarrow -1$。当 $\alpha=0$，我们可以取得一个和海灵格距离线性相关的对称散度，定义如下

\[ D_H (p||q) \triangleq \int \left ( p (x)^{1\over2} - q (x) ^{1\over 2}\right)^2 dx ~~\tag {21.22}\label {eqn:21.22} \]

注意到 $\sqrt {D_H (p||q)}$ 是一个有效的距离度量。也就是说，它对称、非负且满足三角不等式。

21.3 平均场方法¶

21.3.1 基本思想¶

最流行的变分推断形式之一被称为平均场近似（Opper, Saad 2001）。在这种方法中，假设后验在形式上可以被完全分解近似： $$ q(x) = \prod \limits_i q_i(x_i) ~~\tag {21.23}\label {eqn:21.23} $$ 而变分推断的目标是解决如下优化问题： $$ \min _{q_{1}, \ldots, q_{D}} \mathbb{K} \mathbb{L}(q \| p) ~~\tag {21.24}\label {eqn:21.24} $$

因此平均场近似方法本质上是：将原优化问题转化为对近似分布的每个因子（维度）上的边缘分布的某些参数进行优化。在第 21.3.2 节中，将推导出一种坐标下降法，其中每一步都进行如下更新：

\[ \log q_{j}\left(\mathbf{x}_{j}\right)=\mathbb{E}_{-q_{j}}[\log \tilde{p}(\mathbf{x})]+\text { const } ~~\tag {21.25}\label {eqn:21.25} \]

其中： $\tilde{p}(x)= p(x,D)$ 是未归一化的后验，符号 $E_{-q_j}[f(x)]$ 表示除 $x_j$ 之外的所有变量对 $f(x)$ 的期望。例如：如果有三个变量，则： $$ \mathbb{E}_{-q_{2}}[f(\mathbf{x})]=\sum_{x_{1}} \sum_{x_{3}} q_1\left(x_{1}\right) q_{3}\left(x_{3}\right) f\left(x_{1}, x_{2}, x_{3}\right) ~~\tag {21.26}\label {eqn:21.26} $$

其中 $\sum$ 在必要时可采用积分代替。

当更新 $q_j$ 时，只需要推断与 $x_j$ 共享因子的变量，即 $j$ 的马尔可夫毯中的项（见第 10.5.3 节）；其他项被吸收到常项$const$ 中。由于用平均值替换相邻的值，所以该方法又被称为均值场。这与吉布斯采样（第 24. 2节）非常相似，但不是在相邻节点之间发送采样值作为消息，而是发送平均值。这往往更有效，因为平均值可以用作大量样本的代理。另一方面，平均场消息是密集的，而样本是稀疏的，这使采样更容易扩展到非常大的模型。

当然，每次更新一个分布 $q_j(\mathbf{x}_j)$ 会很慢。已经提出了几种方法来加速这种基本方法，包括：使用模式搜索、基于参数扩展等技术。需要注意：平均场方法可以使用各种参数形式的 $q_i$ 来推断离散或连续型隐变量。与其他变分方法相比，平均场方法适用性更受限制。表 21.1 列出了本书涉及的一些平均场应用案例。

表 21.1 本书中的一些模型，书中提供了平均场推断算法的详细推导。

21.3.2 平均场更新过程¶

回想一下，变分推断的目标是最小化上界 $J(q)≥-\log p(D)$。等效地，可以尝试最大化下限： $$ L(q) \triangleq-J(q)=\sum_{\mathbf{x}} q(\mathbf{x}) \log \frac{\tilde{p}(\mathbf{x})}{q(\mathbf{x})} \leq \log p(\mathcal{D}) ~~\tag {21.27}\label {eqn:21.27} $$ 如果写出目标，挑出涉及 $q_j$ 的项，并将所有其他项视为常量，就会得到： $$ \begin{array}{l}

\[\begin{align*} L\left(q_{j}\right) &=\sum_{\mathbf{x}} \prod_{i} q_{i}\left(\mathbf{x}_{i}\right)\left[\log \tilde{p}(\mathbf{x})-\sum_{k} \log q_{k}\left(\mathbf{x}_{k}\right)\right] ~~\tag {21.28}\label {eqn:21.28}\\ =& \sum_{\mathbf{x}_{j}} \sum_{\mathbf{x}-j} q_{j}\left(\mathbf{x}_{j}\right) \prod_{i \neq j} q_{i}\left(\mathbf{x}_{i}\right)\left[\log \tilde{p}(\mathbf{x})-\sum_{k} \log q_{k}\left(\mathbf{x}_{k}\right)\right] ~~\tag {21.29}\label {eqn:21.29}\\ =& \sum_{\mathbf{x}_{j}} q_{j}\left(\mathbf{x}_{j}\right) \sum_{\mathbf{x}_{-j}} \prod_{i \neq j} q_{i}\left(\mathbf{x}_{i}\right) \log \tilde{p}(\mathbf{x}) \\ &-\sum_{\mathbf{x}_{j}} q_{j}\left(\mathbf{x}_{j}\right) \sum_{\mathbf{x}_{-j}} \prod_{i \neq j} q_{i}\left(\mathbf{x}_{i}\right)\left[\sum_{k \neq j} \log q_{k}\left(\mathbf{x}_{k}\right)+q_{j}\left(\mathbf{x}_{j}\right)\right] ~~\tag {21.30}\label {eqn:21.30}\\ =& \sum_{\mathbf{x}_{j}} q_{j}\left(\mathbf{x}_{j}\right) \log f_{j}\left(\mathbf{x}_{j}\right)-\sum_{\mathbf{x}_{j}} q_{j}\left(\mathbf{x}_{j}\right) \log q_{j}\left(\mathbf{x}_{j}\right)+\mathrm{const} ~~\tag {21.31}\label {eqn:21.31}\\\\ \text { where }\\ \log f_{j}\left(\mathbf{x}_{j}\right) &\triangleq \sum_{\mathbf{x}_{-j}} \prod_{i \neq j} q_{i}\left(\mathbf{x}_{i}\right) \log \tilde{p}(\mathbf{x})=\mathbb{E}_{-q_{j}}[\log \tilde{p}(\mathbf{x})] ~~\tag {21.32}\label {eqn:21.32} \end{align*}\]

\end{array} $$ 所以取除 $\mathbf{x}{j}$ 外的所有隐变量的平均值。进而可以重写 $ L\left(q{j}\right)$ 如下： $$ L\left(q_{j}\right)=-\mathbb{K} \mathbb{L}\left(q_{j} \| f_{j}\right) ~~\tag {21.33}\label {eqn:21.33} $$ 可通过最小化该 KL 来最大化 $L$，这可以通过设置 $q_j=f_j$ 来实现，如下所示： $$ q_{j}\left(\mathbf{x}_{j}\right)=\frac{1}{Z_{j}} \exp \left(\mathbb{E}_{-q_{j}}[\log \tilde{p}(\mathbf{x})]\right) ~~\tag {21.34}\label {eqn:21.34} $$ 通常可以忽略局部归一化常数 $Z_j$，因为知道 $q_j$ 一定是归一化分布。因此，通常使用形式： $$ \log q_{j}\left(\mathbf{x}_{j}\right)=\mathbb{E}_{-q_{j}}[\log \tilde{p}(\mathbf{x})]+\text { const } ~~\tag {21.35}\label {eqn:21.35} $$ { }

$q_j$ 分布的函数形式由变量 $\mathbf{x}_j$ 的类型以及模型的形式来确定（这有时被称为自由形式优化)。如果 $\mathbf{x}_j$ 是离散型随机变量，则 $q_j$ 将是离散分布；如果 $x_j$ 是连续型随机变量，则 $q_j$ 将是某种 $pdf$ 。

21.3.3 示例：因子化的伊辛模型¶

考虑 19.4.1 节中的图像去噪示例，其中 $x_i \in \{−1，+1\}$ 是“净”图像的隐像素值。有如下形式的联合模型： $$ p(\mathbf{x},\mathbf{y})=p(\mathbf{x})p(\mathbf{y}|\mathbf{x}) ~~\tag {21.36}\label {eqn:21.36} $$ 其中，先验具有如下形式： $$ \begin{aligned} p(\mathbf{x}) &=\frac{1}{Z_{0}} \exp \left(-E_{0}(\mathbf{x})\right) \\ E_{0}(\mathbf{x}) &=-\sum_{i=1}^{D} \sum_{j \in \mathrm{nbr}_{i}} W_{i j} x_{i} x_{j} \end{aligned} $$ 似然具有如下形式： $$ p(\mathbf{y} \mid \mathbf{x})=\prod_{i} p\left(\mathbf{y}_{i} \mid x_{i}\right)=\sum_{i} \exp \left(-L_{i}\left(x_{i}\right)\right) $$ 后验具有如下形式： $$ \begin{aligned} p(\mathbf{x} \mid \mathbf{y}) &=\frac{1}{Z} \exp (-E(\mathbf{x})) \\ E(\mathbf{x}) &=E_{0}(\mathbf{x})-\sum_{i} L_{i}\left(x_{i}\right) \end{aligned} $$ 现在可以用完全因子化来近似后验： $$ q(\mathbf{x})=\prod_{i} q\left(x_{i}, \mu_{i}\right) $$ 其中 $ \mu_{i} $ 表示第 $i$ 个节点的平均值。对 $q_i$ 的更新，其实就是对变分参数 $μ_i$ 的更新，首先写出 $\log \tilde{p}(\mathbf{x})=-E(\mathbf{x})$，丢弃不涉及 $x_i$ 的项： $$ \log \tilde{p}(\mathbf{x})=x_{i} \sum_{j \in \mathrm{nbr}_{i}} W_{i j} x_{j}+L_{i}\left(x_{i}\right)+\text { const } $$ 这仅取决于相邻节点的状态。现在期望相对于 $\prod_{j \neq i} q_{j}\left(x_{j}\right)$，可以获得： $$ q_{i}\left(x_{i}\right) \propto \exp \left(x_{i} \sum_{j \in \mathrm{nbr}_{i}} W_{i j} \mu_{j}+L_{i}\left(x_{i}\right)\right) $$ 因此，将邻居的状态替换为它们的平均值。让 $$ m_{i}=\sum_{j \in \mathrm{nbr}_{i}} W_{i j} \mu_{j} $$

是节点 $i$ 上的平均场影响。还有，设 $ L_{i}^{+} \triangleq L_{i}(+1)$ 并且 $L_{i}^{-} \triangleq L_{i}(-1) $。则近似后验的边缘分布由下式给出：

\[\begin{split} \begin{aligned} q_{i}\left(x_{i}=1\right) &=\frac{e^{m_{i}+L_{i}^{+}}}{e^{m_{i}+L_{i}^{+}}+e^{-m_{i}+L_{i}^{-}}}=\frac{1}{1+e^{-2 m_{i}+L_{i}^{-}-L_{i}^{+}}}=\operatorname{sigm}\left(2 a_{i}\right) \\ a_{i} & \triangleq m_{i}+0.5\left(L_{i}^{+}-L_{i}^{-}\right) \end{aligned} \end{split}\]

类似地，我们有。由此，可以计算节点 $i$ 的新平均值： $$ \begin{aligned} \mu_{i} &=\mathbb{E}_{q_{i}}\left[x_{i}\right]=q_{i}\left(x_{i}=+1\right) \cdot(+1)+q_{i}\left(x_{i}=-1\right) \cdot(-1) \\ &=\frac{1}{1+e^{-2 a_{i}}}-\frac{1}{1+e^{2 a_{i}}}=\frac{e^{a_{i}}}{e^{a_{i}}+e^{-a_{i}}}-\frac{e^{-a_{i}}}{e^{-a_{i}}+e^{a_{i}}}=\tanh \left(a_{i}\right) \end{aligned} $$ 因此更新方程变为： $$ \mu_{i}=\tanh \left(\sum_{j \in \mathrm{nbr}_{i}} W_{i j} \mu_{j}+0.5\left(L_{i}^{+}-L_{i}^{-}\right)\right) $$

可以通过以下方式将上述方程转化为迭代算法： $$ \mu_{i}^{t}=\tanh \left(\sum_{j \in \mathrm{nbr}_{i}} W_{i j} \mu_{j}^{t-1}+0.5\left(L_{i}^{+}-L_{i}^{-}\right)\right) $$

通常情况下，最好使用抑制更新的形式： $$ \mu_{i}^{t}=(1-\lambda) \mu_{i}^{t-1}+\lambda \tanh \left(\sum_{j \in \mathrm{nbr}_{i}} W_{i j} \mu_{j}^{t-1}+0.5\left(L_{i}^{+}-L_{i}^{-}\right)\right) $$

对于 $0<λ<1$ ，可以并行更新所有节点，也可以异步更新。

图 21.3 显示了该方法的实际应用，应用于具有齐次引力势的二维伊辛模型，$W_{ij}=1$。我们使用阻尼因子为 $λ=0.5$ 时的并行更新。

图 21.3 使用平均场进行图像去噪的示例（并行更新，阻尼因子为0.5）。我们使用伊辛先验（$W_{ij}=1$ ）和高斯噪声模型（ $σ=2$ ）。我们显示了图像经过1、3和15次迭代后的结果。请参阅图 24.1。

21.4 结构化平均场¶

21.4.1 基本思路¶

假设所有变量在后验都是独立的，这是一个非常强有力的假设，可能会导致糟糕的结果。有时可以在问题中利用易处理的子结构，从而可以有效地处理某些类型的依赖关系。这被称为结构化平均场方法(Saul和Jordan，1995)。方法与以前类似，只是将变量集进行了分组，并按照组同时更新（通过简单地将第 $i$ 组中的所有变量视为单个“巨型变量”，然后重复21.3.2 节中的推导来实现）。只要能在每个 $q_j$ 中进行有效的推断，该方法总体上是容易掌握的。下面给出一个例子。见(Bouchard-Cote and Jordan，2009) 以了解这一领域的最新工作。

图 21.4 （a）具有3个链的因子化隐马尔可夫模型。（b）全因子分解近似。（c）链乘积近似。基于 （Ghahramani And Jordan 1997) 图 2。

21.4.2 示例：因子化的HMM¶

考虑 17.6.5 节中介绍的因子化HMM模型 (Ghahramani和Jordan，1997)。假设有 $M$ 个链，每个链的长度为 $T$，并且假设每个隐变量有 $K$ 个状态。模型定义如下： $$ p(\mathbf{x}, \mathbf{y})=\prod_{m} \prod_{t} p\left(x_{t m} \mid x_{t-1, m}\right) p\left(\mathbf{y}_{t} \mid x_{t m}\right) $$

其中 $p\left(x_{t m}=k \mid x_{t-1, m}=j\right)=A_{m j k}$ 是链 $m$ 的转移矩阵中的条目，$p\left(x_{1 m}=k \mid x_{0 m}\right)=p\left(x_{1 m}=k\right)=\pi_{m k}$ 是链 $m$ 的初始状态分布，并且： $$ p\left(\mathbf{y}_{t} \mid \mathbf{x}_{t}\right)=\mathcal{N}\left(\mathbf{y}_{t} \mid \sum_{m=1}^{M} \mathbf{W}_{m} \mathbf{x}_{t m}, \mathbf{\Sigma}\right) $$ 是观测模型。其中 $\mathbf{x}{t m} $ 是 $x{t m} $ 的一个 $1 \text{ of } K$ 编码，$\mathbf{W}{m}$ 是 $D \times K$ 阶矩阵（假设 $ \mathbf{y}{t} \in \mathbb{R}^{D} $）。图 21.4(a) 说明了 $M=3$ 的情况下的模型。即使每条链都是先验独立的，但由于有可观测的共同孩子 $ \mathbf{y}_{t} $ ，它们在后验中被耦合了。应用于该概率图的连接树算法所需时间为 $O\left(T M K^{M+1}\right) $ 。下面将推导出一个结构化的平均场算法，它需要 $O\left(T M K^{2} I\right)$ 时间，其中 $I$ 是平均场迭代的次数（通常 $I∼10$ 就足以获得良好的性能）。

可以用下面的形式写出确切的后验：

确切的后验概率可以写成如下形式： $$ \begin{aligned} p(\mathbf{x} \mid \mathbf{y})=& \frac{1}{Z} \exp (-E(\mathbf{x}, \mathbf{y})) \\ E(\mathbf{x}, \mathbf{y})=& \frac{1}{2} \sum_{t=1}^{T}\left(\mathbf{y}_{t}-\sum_{m} \mathbf{W}_{m} \mathbf{x}_{t m}\right)^{T} \boldsymbol{\Sigma}^{-1}\left(\mathbf{y}_{t}-\sum_{m} \mathbf{W}_{m} \mathbf{x}_{t m}\right) \\ &-\sum_{m} \mathbf{x}_{1 m}^{T} \tilde{\boldsymbol{\pi}}_{m}-\sum_{t=2}^{T} \sum_{m} \mathbf{x}_{t m}^{T} \tilde{\mathbf{A}}_{m} \mathbf{x}_{t-1, m} \end{aligned} $$

其中 $ \tilde{\mathbf{A}}_{m} \triangleq \log \mathbf{A}_{m}$ 并且 $\tilde{\boldsymbol{\pi}}_{m} \triangleq \log \boldsymbol{\pi}_{m} $ （两者都逐元素解释）。

可以将后验近似为边缘的乘积，如图 21.4(b) 所示，但更好的近似是使用链的乘积，如图 21.4(c) 所示。每条链都可以使用向前-向后算法分别进行跟踪更新。更准确地说，假设

\[\begin{split} \begin{aligned} q(\mathbf{x} \mid \mathbf{y}) &=\frac{1}{Z_{q}} \prod_{m=1}^{M} q\left(x_{1 m} \mid \boldsymbol{\xi}_{1 m}\right) \prod_{t=2}^{T} q\left(x_{t m} \mid x_{t-1, m}, \boldsymbol{\xi}_{t m}\right) \\ q\left(x_{1 m} \mid \boldsymbol{\xi}_{1 m}\right) &=\prod_{k=1}^{K}\left(\xi_{1 m k} \pi_{m k}\right)^{x_{1 m k}} \\ q\left(x_{t m} \mid x_{t-1, m}, \boldsymbol{\xi}_{t m}\right) &=\prod_{k=1}^{K}\left(\xi_{t m k} \prod_{j=1}^{K}\left(A_{m j k}\right)^{x_{t-1, m, j}}\right)^{x_{t m k}} \end{aligned} \end{split}\]

可以看到 $ξ_{tmk}$ 参数起到近似局部证据的作用，平均了其他链的影响。这与将所有的链连接在一起的精确局部证据形成了鲜明对比。

可将近似后验公式重写为 $q(\mathbf{x})=\frac{1}{Z_{q}} \exp \left(-E_{q}(\mathbf{x})\right)$ ，其中： $$ E_{q}(\mathbf{x})=-\sum_{t=1}^{T} \sum_{m=1}^{M} \mathbf{x}_{t m}^{T} \tilde{\boldsymbol{\xi}}_{t m}-\sum_{m=1}^{M} \mathbf{x}_{1 m}^{T} \tilde{\boldsymbol{\pi}}_{m}-\sum_{t=2}^{T} \sum_{m=1}^{M} \mathbf{x}_{t m}^{T} \tilde{\mathbf{A}}_{m} \mathbf{x}_{t-1, m} $$ 其中 $\tilde{\boldsymbol{\xi}}{t m}=\log \boldsymbol{\xi}{t m}$。可以看到，这与精确后验具有相同的时间要素，但局部证据项有所不同。目标函数由 $$ \mathbb{K} \mathbb{L}(q \| p)=\mathbb{E}[E]-\mathbb{E}\left[E_{q}\right]-\log Z_{q}+\log Z $$

给出，其中期望取相对于 $q$ 。更新具有以下形式： $$ \begin{aligned} \boldsymbol{\xi}_{t m} &=\exp \left(\mathbf{W}_{m}^{T} \boldsymbol{\Sigma}^{-1} \tilde{\mathbf{y}}_{t m}-\frac{1}{2} \boldsymbol{\delta}_{m}\right) \\ \boldsymbol{\delta}_{m} & \triangleq \operatorname{diag}\left(\mathbf{W}_{m}^{T} \mathbf{\Sigma}^{-1} \mathbf{W}_{m}\right) \\ \tilde{\mathbf{y}}_{t m} & \triangleq \mathbf{y}_{t}-\sum_{\ell \neq m}^{M} \mathbf{W}_{\ell} \mathbb{E}\left[\mathbf{x}_{t, \ell}\right] \end{aligned} $ξ_{tm}$ 参数起到局部证据的作用，在相邻链上取平均。在为每个链计算了这个之后，可以并行地执行`向前-向后`，使用这些近似的局部证据项来为每个 $m$ 和 $t$ 计算 $q(x_t,m|y_{1:T})$。

对所有变分参数完成“扫描”所需更新成本为 $0(TMK^2)$，因为必须对每个链独立地向前-向后运行 $M$ 次。这与完全因子分解近似法的成本相同，但要精确得多。

21.5 变分贝叶斯（VB）¶

目前为止，我们主要专注于在模型参数 $\theta$ 已知的情况下推断隐变量 $z_i$ 的分布。现在起开始考虑推断模型参数其本身。做一个完全因子分解近似（即平均场） $p (\theta | D) \approx \prod_k q (\theta_k)$，实际上就是在做 变分贝叶斯（VB: Variational Bayes） 。接下来会给出一些在假设不存在隐变量情况下，应用变分贝叶斯的例子，如果想要同时推断隐变量和模型参数，并作如下近似 $p (\theta,z_{1:N}|D) \approx q (\theta) \prod_i q_i (z_i)$，则实际上就在使用变分贝叶斯期望最大化算法 (VBEM)，详见 21.6 节。

21.5.1 从一维变量的变分贝叶斯开始¶

根据 (MacKay 2003,p429)，考虑推断一维高斯分布的参数分布 $p (\mu,\lambda | D)$, 此处 $\lambda=1/\sigma^2$ 指高斯精度。为简便起间，我们使用一个共轭先验分布：

\[ p (\mu,\lambda) = \mathcal {N}(\mu|\mu_0,(\kappa \lambda)^{-1})\text {Ga}(\lambda|a_0,b_0)\tag {21.65}\label {eqn:21.65} \]

而近似后验是如下做了因子化的分布：

\[ q (\mu,\lambda) = q_\mu (\mu) q_\lambda (\lambda) \tag {21.66}\label {eqn:21.66} \]

此处不需要指定分布 $q_\mu,\,q_\lambda$ 的确切形式，因为其最优形式会在推导时自动 “出现”（实际上它们恰巧是相应的高斯分布和伽马分布）。

你可能想知道为什么要这样做，因为知道如何计算该模型的精确后验（第4.6.3.7节）。原因有两条：首先这是一个有教学意义的练习，可以乘机通过与确切分布的对比，来分析近似分布的质量；其次，该方法可以很自然地加以改造以适用一个半共轭的先验 $p (\mu,\lambda)=\mathcal {N}(\mu|\mu_0,\tau_0)\text {Ga}(\lambda|a_0,b_0)$，而这种先验下是无法进行确切推断的。

21.5.1.1 目标分布¶

未归一化的对数后验具有形式：

\[\begin{split} \begin {align*} \log \tilde {p}(\mu,\lambda) = \log p (\mu,\lambda,D) &= \log p (D|\mu,\lambda) + \log p (\mu|\lambda) + \log p (\lambda) \tag {21.67}\label {eqn:21.67} \\ &=\frac {N}{2}\log\lambda - {\lambda \over 2 } \sum_{i=1}^N (x_i - \mu)^2 - \frac {\kappa_0 \lambda}{2}(\mu-\mu_0)^2 \notag\\ & ~~ +{1\over 2}\log (\kappa_0 \lambda) + (a_0 - 1)\log \lambda - b_0 \lambda + \text {const} \tag {21.68}\label {eqn:21.68} \end {align*} \end{split}\]

21.5.1.2 更新 $q_\mu（\mu）$¶

通过对 $\lambda$ 求期望可以得到 $q_\mu (\mu)$ 的最优形式：

\[\begin{split} \begin {align*} \log q_\mu (\mu) &=\mathbb{E}_{q_\lambda}[\log p (D|\mu,\lambda) + \log p (\mu|\lambda) ] + \text {const} \tag {21.69}\label {eqn:21.69} \\ &= -\frac {\mathbb{E} [\lambda]}{2} \left\{ \kappa_0 (\mu - \mu_0)^2 + \sum_{i=1}^N (x_i - \mu)^2 \right \} + \text {const} \tag {21.70}\label {eqn:21.70} \end {align*} \end{split}\]

通过展开平方项（？疑）可以证明 $q_\mu (\mu) = \mathcal {N}(\mu | \mu_N , \kappa^{-1}_N)$, 且有

\[ \mu_N={ \kappa_0 \mu_0 + N \bar {x} \over \kappa_0 + N }, \kappa_N = (\kappa_0 + N )\mathbb{E}_{q_\lambda}[\lambda] \tag {21.71}\label {eqn:21.71} \]

此时我们还不知道 $q_\lambda (\lambda)$ 的确切分布形式，所以无法计算 $\mathbb{E} [\lambda]$，但会解决它。

21.5.1.3 更新 $q_\lambda (\lambda)$¶

$q_\lambda (\lambda)$ 的最优分布有如下表示

\[\begin{split} \begin {align*} \log q_\lambda (\lambda) &=\mathbb{E}_{q_\mu}[ \log p (D|\mu,\lambda) + \log p (\mu | \lambda ) + \log p (\lambda) ] + \text{ const} \tag {21.72}\label {eqn:21.72} \\ &= (a_0 - 1 )\log \lambda -　b_0 \lambda + {1\over2}\log\lambda + {N\over2} \log\lambda \\ &~~-{\lambda\over2}\mathbb{E}_{q_\mu}\left [ \kappa_0 (\mu - \mu_0)^2 + \sum_{i=1}^N (x_i-\mu)^2 \right] + \text{ const} \tag {21.73}\label {eqn:21.73} \end {align*} \end{split}\]

可以发现这对应一个伽马分布，因而有　$q_\lambda (\lambda)=Ga (\lambda|a_N,b_N)$, 且有

\[\begin{split} \begin {align*} a_N &= a_0 + \frac {N+1}{2} \tag {21.74}\label {eqn:21.74} \\ b_N &= b_0 + \frac {1}{2}\mathbb{E}_{q_\mu}\left [ \kappa_0 (\mu - \mu_0)^2 + \sum_{i=1}^N (x_i-\mu)^2 \right] \tag {21.75}\label {eqn:21.75} \end {align*} \end{split}\]

21.5.1.4 计算期望值¶

为实现上述更新，必须确定如何计算诸多期望值。鉴于 $q (\mu) = \mathcal (\mu | \mu_N,\kappa_N^{-1})$，可以得出：

\[\begin{split} \begin {align*} \mathbb{E}_{q_(\mu)}[\mu] &= \mu_N \tag {21.76}\label {eqn:21.76} \\ \mathbb{E}_{q (\mu)}[\mu^2] &= {1\over \kappa_N} + \mu_N^2 \tag {21.77}\label {eqn:21.77} \end {align*} \end{split}\]

鉴于 $q (\lambda) = \text {Ga}(\lambda \mid a_N, b_N)$，有：

\[ \mathbb{E}_{q (\lambda)}[\lambda] = \frac {a_N}{b_N} \tag {21.78}\label {eqn:21.78} \]

如此可以显式给出更新方程。对 $q (\mu)$ 有

\[\begin{split} \begin {align*} \mu_N &= \frac {\kappa_0 \mu_0 + N \bar {x}}{\kappa_0 + N} \tag {21.79}\label {eqn:21.79} \\ \kappa_N &= (\kappa_0 + N ) \frac {a_N}{b_N} \tag {21.80}\label {eqn:21.80} \end {align*} \end{split}\]

相应地对 $q (\lambda)$ 有

\[\begin{split} \begin {align*} a_N &= a_0 + \frac {N+1}{2} \tag {21.81}\label {eqn:21.81} \\ b_N &= b_0 + \kappa_0 (\mathbb{E} [\mu^2] +\mu_0^2 -2\mathbb{E} [\mu]\mu_0 ) + \frac {1}{2}\sum_{i=1}^N (x_i^2 +\mathbb{E} [\mu^2] -2\mathbb{E} [\mu ] x_i) \tag {21.82}\label {eqn:21.82} \end {align*} \end{split}\]

可以看到 $\mu_N$ 和 $a_N$ 实际上是固定的常数，且只有 $\kappa_N\,,b_N$ 需要被迭代更新。（事实上，可以利用迭代方程，解析地求出 $\kappa_N\,\,b_N$ 的固定点，但是在此不作展示，只介绍迭代更新法）

21.5.1.5 展示¶

图 21.5 给出了一个工作示意。绿色的轮廓表示确切的高斯 - 伽马后验分布，虚线的轮廓表示的是变分近似在不同的迭代步的结果。可以看到最终的近似是和确切解较为相似的。然而，近似解比真实分布更为 “紧致”。事实上平均场近似推断常常欠估计后验的不确定性，更多有关讨论请见 21.2.2 节。

图 21.5 高斯-伽马分布(绿色)的因子化变分近似（红色）。（a）初始猜测。（b）更新 $q_μ$ 后。（c）更新 $q_λ$ 后。（d）收敛时（5次迭代后）。基于（Bishop 2006b）第10.4节。

21.5.1.6 下界¶

在变分贝叶斯中，最大化的 $L (q)$ 是数据边缘似然值（证据）的下界。

\[ L (q) \le \log p (D) = \log \int \int p (D\mu) p (\mu, \lambda) d\mu d\lambda \tag {21.83}\label {eqn:21.83} \]

出于多种原因，需要计算该下界本身。（1）下界可用来检查算法是否收敛了；（2）下界可用来检查算法的正确性。与期望最大算法 (EM) 的情况一样，如果下界不是单调增加的话，那一定是哪里出了问题；（3）下界可用作对数据边缘似然（证据）的近似，因而可以用于贝叶斯模型的选择（注：模型选择是贝叶斯方法的一个重要组成部分）。

不幸的是，计算这个下界需要大量繁琐的代数运算。本节会为例子计算出细节，但是对于其他模型，将只陈述结果而不做证明，或者完全省略对边界的讨论。

对该模型，可以如下计算似然函数 $L (q)$ ：

\[\begin{split} \begin {align*} L (q) &= \int \int q (\mu, \lambda) \log \frac { p (D,\mu,\lambda) } {q (\mu,\lambda) } d\mu d\lambda \tag {21.84}\label {eqn:21.84} \\\\ &=\mathbb{E} [\log p (D\ gvn \mu,\lambda) ] +\mathbb{E} [ \log p (\mu \mid \lambda) ] +\mathbb{E} [ \log p (\lambda) ] \notag \\\\ & ~~ -\mathbb{E} [\log q (u)] -\mathbb{E} [ \log q (\lambda) ] \tag {21.85}\label {eqn:21.85} \end {align*} \end{split}\]

上式中所有期望值都是关于 $q (\mu,\lambda)$ 的。可看出最后两项不过是高斯分布和伽马分布的熵，因此可以直接写出：

\[\begin{split} \begin {align*} \mathbb{H} ( \mathcal {N}(\mu_N,\kappa_N^{-1} ) &= -\frac {1}{2} \log \kappa_N + \frac {1}{2} (1 + \log (2\pi)) \tag {21.86}\label {eqn:21.86}\\\\ \mathbb{H} (\text {Ga}(a_N,b_N)) &= \log \Gamma (a_N) - (a_N - 1) \psi (a_N) - \log (b_N) + a_N \tag {21.87}\label {eqn:21.87} \end {align*} \end{split}\]

此处 $\psi ()$ 是双伽马函数。

为计算其它项，要用到如下结果：

\[\begin{split} \begin {align*} \mathbb{E} [ \log x \mid x \sim \text {Ga}(a,b)] &= \psi (a) - \log (b) \tag {21.88}\label {eqn:21.88} \\ \mathbb{E} [ x \mid x \sim \text {Ga}(a,b) ] &= {a \over b} \tag {21.89}\label {eqn:21.89} \\ \mathbb{E} [ x \mid x \sim \mathcal {N}(\mu, \sigma^2) ] &= \mu \tag {21.90}\label {eqn:21.90} \\ \mathbb{E} [ x^2 \mid x \sim \mathcal {N} (\mu, \sigma^2)] &= \mu^2 + \sigma^2 \tag {21.91}\label {eqn:21.91} \end {align*} \end{split}\]

可以证明对数似然的期望符合下式：

\[\begin{split} \begin {align*} \mathbb{E}_{q (\mu,\lambda)}&[ \log p (D \mid \mu \lambda) ] \tag {21.92}\label {eqn:21.92} \\ &= -{N\over 2 }\log (2\pi) + {N\over 2}\mathbb{E}_{q_\lambda}[\log \lambda] - \frac {\mathbb{E}_{q_\lambda}[ \lambda]}{2} \sum_{i=1}^N \mathbb{E}_{q (\mu)} [ (x_i - \mu)^2 ] \tag {21.93}\label {eqn:21.93} \\ &= -{N\over 2 }\log (2\pi) + {N\over 2} (\psi (a_N) - \log b_N) \\ & ~~ - {N a_N \over 2 b_N} \left ( \hat {\sigma}^2 + \bar {x}^2 - 2\mu_N \bar {x} + \mu^2_N + \frac {1}{\kappa_N} \right) \tag {21.94}\label {eqn:21.94} \end {align*} \end{split}\]

此处的 $\bar {x},\,\hat {\sigma}^2$ 对应经验均值和经验方差。

对 $\lambda$ 的对数先验的期望可做如下处理：

\[\begin{split} \begin {align*} \mathbb{E}_{q (\lambda)}[ \log p (\lambda) ] ~~&=~~ (a_0 - 1)\mathbb{E} [ \log \lambda] -b_0\mathbb{E} [\lambda ] + a_0\log b_0 -\log \Gamma ( a_0 ) \tag {21.95}\label {eqn:21.95} \\ ~~&=~~ (a_0 - 1)(\psi (a_N) - \log b_N) -b_0 \frac {a_N}{b_N} + a_0\log b_0 -\log \Gamma ( a_0 ) \tag {21.96}\label {eqn:21.96} \end {align*} \end{split}\]

对 $\mu$ 的对数先验的期望可做如下处理：

\[\begin{split} \begin {align*} \mathbb{E}_{q (\mu, \lambda)}[ \log p (\mu \mid \lambda) ] ~~&=~~ \frac {1}{2} \log \frac {\kappa_0}{2\pi} + \frac {1}{2} \mathbb{E} [ \log \lambda] q (\lambda) - \frac {1}{2}\mathbb{E}_{q (\mu,\lambda)} [ (\mu - \mu_0)^2 \kappa_0 \lambda ] \\ ~~&=~~ \frac {1}{2} \log \frac {\kappa_0}{2\pi} + \frac {1}{2} (\psi (a_N) - \log b_N) \\ & ~~ - \frac {\kappa_0}{2} \frac {a_N}{b_N} \left [ (\mu_N - \mu_0)^2 + \frac {1}{\kappa_N} \right] \tag {21.97}\label {eqn:21.97} \end {align*} \end{split}\]

将以上结果结合，可得到对数似然的下界

\[ L (q) = {1 \over 2 }\log {1 \over \kappa_N} + \log \Gamma (a_N) - a_N \log b_N + \text{ const} \tag {21.98}\label {eqn:21.98} \]

该量随着变分贝叶斯的迭代更新而单调递增。

21.5.2 示例：线性回归的变分贝叶斯¶

在第 7.6.4 节中，曾经讨论了一种经验贝叶斯方法来设置岭回归的超参数。特别地，我们假设了似然的形式为 $p(y|X，θ) = \mathcal{N}(X\mathbf{w}，λ^{-1})$ ，先验的形式为 $p(\mathbf{w}) = \mathcal{N}(\mathbf{w}|0，α^{-1}I)$ 。然后计算了 $α$ 和 $λ$ 的 II 型估计。同样的方法在第 13.7 节中得到扩展，以处理形式为 $N(\mathbf{w}|0，\text{diag}(α)^{-1})$ 的先验，它允许每个特征有一个超参数，该技术被称为自动相关性确定。

本节将为该模型导出了一个VB算法。首先，用以下先验知识： $$ p(\mathbf{w}, \lambda, \alpha)=\mathcal{N}\left(\mathbf{w} \mid \mathbf{0},(\lambda \alpha)^{-1} \mathbf{I}\right) \mathrm{Ga}\left(\lambda \mid a_{0}^{\lambda}, b_{0}^{\lambda}\right) \mathrm{Ga}\left(\alpha \mid a_{0}^{\alpha}, b_{0}^{\alpha}\right) $$

选择如下的后验因子分解近似： $$ q(\mathbf{w}, \alpha, \lambda)=q(\mathbf{w}, \lambda) q(\alpha) $$

给定这些假设，可以证明（略）后验的最优形式是： $$ q(\mathbf{w}, \alpha, \lambda)=\mathcal{N}\left(\mathbf{w} \mid \mathbf{w}_{N}, \lambda^{-1} \mathbf{V}_{N}\right) \operatorname{Ga}\left(\lambda \mid a_{N}^{\lambda}, b_{N}^{\lambda}\right) \mathrm{Ga}\left(\alpha \mid a_{N}^{\alpha}, b_{N}^{\alpha}\right) $$

其中： $$

\begin{aligned} \mathbf{V}{N}^{-1} &=\overline{\mathbf{A}}+\mathbf{X}^{\mathbf{X}} \ \mathbf{w}{N} &=\mathbf{V}{N} \mathbf{X}^{T} \mathbf{y} \ a{N}^{\lambda} &=a_{0}^{\lambda}+\frac{N}{2} \ b_{N}^{\lambda} &=b_{0}^{\lambda}+\frac{1}{2}\left(|\mathbf{y}-\mathbf{X} \mathbf{w}|^{2}+\mathbf{w}{N}^{T} \overline{\mathbf{A}} \mathbf{w}{N}\right) \ a_{N}^{\alpha} &=a_{0}^{\alpha}+\frac{D}{2} \ b_{N}^{\alpha} &=b_{0}^{\alpha}+\frac{1}{2}\left(\frac{a_{N}^{\lambda}}{b_{N}^{\lambda}} \mathbf{w}{N}^{T} \mathbf{w}{N}+\operatorname{tr}\left(\mathbf{V}{N}\right)\right) \ \overline{\mathbf{A}} &=\langle\alpha\rangle \mathbf{I}=\frac{a{N}^{\alpha}}{b_{N}^{\alpha}} \mathbf{I} \end{aligned}

\[ 这种方法可以通过使用以下先验以直接的方式扩展到 $ARD$ 案例： \]

\begin{aligned} p(\mathbf{w}) &=\mathcal{N}\left(\mathbf{0}, \operatorname{diag}(\alpha)^{-1}\right) \ p(\boldsymbol{\alpha}) &=\prod_{j=1}^{D} \mathrm{Ga}\left(\alpha_{j} \mid a_{0}^{\alpha}, b_{0}^{\alpha}\right) \end{aligned} $$

$\mathbf{w}$ 和 $λ$ 的后验像以前一样计算，除了用 $\bar A = \text{diag} (a^α_N/b^α_{N_j})$ 代替 $a^α_N/b^α_NI$。$α$ 的后验具有以下形式： $$ \begin{aligned} q(\boldsymbol{\alpha}) &=\prod_{j} \mathrm{Ga}\left(\alpha_{j} \mid a_{N}^{\alpha}, b_{N_{j}}^{\alpha}\right) \\ a_{N}^{\alpha} &=a_{0}^{\alpha}+\frac{1}{2} \\ b_{N_{j}}^{\alpha} &=b_{0}^{\alpha}+\frac{1}{2}\left(\frac{a_{N}^{\lambda}}{b_{N}^{\lambda}} w_{N, j}^{2}+\left(\mathbf{V}_{N}\right)_{j j}\right) \end{aligned} $$ 该算法在更新 $q(\mathbf{w}, \lambda)$ 和 $q(\boldsymbol{\alpha})$ 之间交替。一旦推断出 $\mathbf{w}$ 和 $\lambda$ ，后验预测就是学生分布，如等式 7.76 所示。具体来说，对于单个数据案例有： $$ p(y \mid \mathbf{x}, \mathcal{D})=\mathcal{T}\left(y \mid \mathbf{w}_{N}^{T} \mathbf{x}, \frac{b_{N}^{\lambda}}{a_{N}^{\lambda}}\left(1+\mathbf{x}^{T} \mathbf{V}_{N} \mathbf{x}\right), 2 a_{N}^{\lambda}\right) $$ 用于模型选择的精确边缘似然由下式给出： $$ p(\mathcal{D})=\iiint p(\mathbf{y} \mid \mathbf{X}, \mathbf{w}, \lambda) p(\mathbf{w} \mid \alpha) p(\lambda) d \mathbf{w} d \alpha d \lambda $$ 可以计算 $\log p(\mathcal{D})$ 上的下界： $$ \begin{aligned} L(q)=&-\frac{N}{2} \log (2 \pi)-\frac{1}{2} \sum_{i=1}^{N}\left(\frac{a_{N}^{\lambda}}{b_{N}^{\lambda}}\left(y_{i}-\mathbf{w}_{N}^{T} \mathbf{x}_{i}\right)^{2}+\mathbf{x}_{i}^{T} \mathbf{V}_{N} \mathbf{x}_{i}\right) \\ &+\frac{1}{2} \log \left|\mathbf{V}_{N}\right|+\frac{D}{2} \\ &-\log \Gamma\left(a_{0}^{\lambda}\right)+a_{0}^{\lambda} \log b_{0}^{\lambda}-b_{0}^{\lambda} \frac{a_{N}^{\lambda}}{b_{N}^{\lambda}}+\log \Gamma\left(a_{N}^{\lambda}\right)-a_{N}^{\lambda} \log b_{N}^{\lambda}+a_{N}^{\lambda} \\ &-\log \Gamma\left(a_{0}^{\alpha}\right)+a_{0}^{\alpha} \log b_{0}^{\alpha}+\log \Gamma\left(a_{N}^{\alpha}\right)-a_{N}^{\alpha} \log b_{N}^{\alpha} \end{aligned} $$ 在 ARD 案例中，最后一行变为： $$ \sum_{j=1}^{D}\left[-\log \Gamma\left(a_{0}^{\alpha}\right)+a_{0}^{\alpha} \log b_{0}^{\alpha}+\log \Gamma\left(a_{N}^{\alpha}\right)-a_{N}^{\alpha} \log b_{N_{j}}^{\alpha}\right] $$ 图 21.6 在多项式回归的模型选择问题上比较了 `VB` 和 `EB` 。可以看到 `VB` 给出了与 `EB` 类似的结果，但精确行为取决于样本大小。当 $N = 5$ 时，`VB` 对模型后验的估计比 `EB` 更发散，因为 `VB` 要对超参数的不确定性建模。当 `N = 3 0` 时，超参数的后验估计变得更加确定。事实上，如果有一个非信息先验 $a^α_0=b^α_0=0$ 时，计算 $\mathbb{E}[α|D]$，会得到： $$ \bar \alpha = \frac{a^\alpha_N}{b^\alpha_N}=\frac{D/2}{\frac{1}{2}(\frac{a^\lambda_N}{a^\lambda_N}\mathbf{w}^T_N\mathbf{w}_N+tr(\mathbf{V}_N)}将其与EB的公式13.167进行比较: $$ 将其与 `EB` 的公式 13.167 进行比较： $$ \hat \alpha=\frac{D}{\mathbb{E}[\mathbf{w}^T\mathbf{w}]}=\frac{D}{\mathbf{w}^T_N\mathbf{w}_N+tr(\mathbf{V}_N)} $$ 除 $a^λ_N$ 和 $b^λ_N$ 项外，两者基本相同。这并不奇怪，因为 `EB` 试图最大化 $\log p(D)$，而`VB`试图最大化 $\log p(D)$ 的下界。

图 21.6 我们绘制了均匀先验假设（ $p(m) ∝ 1$）下的多项式回归模型的后验图（分别为1，2和3次多项式）。用 VB （a, c）和 EB（b，d）来近似边缘似然。在 (a-b) 中，使用 $N = 5$ 个数据点。在 (c-d) 中，使用 $N = 3 0$ 个数据点。

21.6 变分贝叶斯期望最大算法¶

21.6.1 基本思路¶

现在让考虑具有形式 $z_i \rightarrow x_i \leftarrow \theta$ 的隐变量模型。这囊括了混合模型、主成分分析、隐马尔可夫模型等。未知变量现在有两种：模型参数 $\theta$，以及隐变量 $z_i$。在章节 11.4 讨论过这类模型而且它们通常是用 “期望最大化（EM）” 算法来拟合的。

EM 算法 的 “期望” 步会推断一个关于隐变量的后验分布 $p (z_i \mid x_i, \theta)$；在 “最大化” 步，计算一个关于模型参数 $\theta$ 的点估计。作此处理的原因有两个：（1）这会给出一个简单的算法；（2）模型参数 $\theta$ 的后验不确定性通常小于隐变量 $z_i$ 的，鉴于 $\theta$ 从所有 $N$ 个样本中得到信息，而 $z_i$ 只从 $x_i$ 获得信息。这就使得对 $\theta$ 而不是 $z_i$ 作点估计显得更为合理。

不过，变分贝叶斯通过同时对模型参数 $\theta$ 和隐变量 $z_i$ 的不定性建模，给出了一个更加 “贝叶斯化” 的算法，同时又具备与 EM 算法差不多的计算速度。该算法被称为 变分贝叶斯 EM (VBEM)。

其基本理念是用平均场来近似后验分布： $$ p (\theta, z_{1:N} \mid D)\approx q (\theta) q (z) = q (\theta) \prod_i q (z_i) \tag {21.120}\label {eqn:21.120} $$

第一个关于 $\theta,\, z_i$ 的因子分解式是简化算法的关键假设。第二个因子分解式可以从在给定模型参数 $\theta$ 后，模型对隐变量独立同分布 的假设中得出。

在 VBEM 中，交替更新 $q (z_i \mid D)$ （变分“E-step”）和 $q (\theta \mid D)$ （变分“M-step”）。要从VBEM 回到 EM 算法，只需用狄拉克函数来表示对 $\theta$ 的点估计即可，即 $q (\theta \mid D) \approx \delta_{\hat {\theta}}(\theta)$

变分 E-step 和标准 E-step 差不多，不过变分 E-step 不像标准 E-step 那样插入参数的最大后验估计和 $p (z_i \mid D,\hat {\theta})$ 计算，而是对参数求期望（均值）。简单地说，可以用模型参数的后验平均值代替最大后验估计，然后计算隐变量的后验分布 $p (z_i \mid D, \hat {\theta})$ 。不幸的是，虽然思想很简单，但事情并不总是那么简单，更多细节常需要根据模型形式决定。

变分 M-setp 和标准 M-step 相似，不过此时更新的是模型参数分布的超参数而不是它的点估计，这是借由充分统计量的期望实现的。该步骤和正常的点估计很相似，但是细节还是要等模型定下来才能讨论。

VBEM 相对于经典 EM 的主要优势在于：通过边缘化模型参数，可以计算数据边缘似然（证据），而该值模型选择非常重要。我们会在 21.6.2.6 看到一个实例。同时，VBEM 奉行 “平等主义”，因为它并不把模型参数视为 “二等公民”，而是和其他未知量平等地联立；相比之下，EM 算法人为区分了模型参数和隐变量。

21.6.2 案例：高斯混合模型的 VBEM¶

现在让我们考虑用变分贝叶斯期最来 “拟合” 一个高斯混合模型（此处的引号是因为我们并不是在估计模型参数，而是在推断它们的后验分布）。我们会基于 (@Bishop2006b,Sec 10.2) 阐述。不幸的是，这些细节都十分的繁杂。不过万幸的是，同期最算法一样，算多了你就自然而然知道该怎么搞了 (?-_-),（和大多数数学一样，光读公式是没啥大用处滴，你一定要尝试自己推导出这些结果（或者试一试课后练习也好啊）不然是没法深入理解这些玩意儿的）

21.6.2.1 高斯混合模型的变分后验¶

高斯混合的似然函数就是还是熟悉的那个

\[ p (z,X \mid \theta ) = \prod_i \prod_k \pi_k^{z_{ik}} \mathcal {N}(x_i \mid \mu_k, \Lambda_k^{-1})^{z_{ik}} \tag {21.121}\label {eqn:21.121} \]

此处有指示函数 $z_{ik} = 1$ 如若第 $i$ 个样本属于类别 $k$, 反之则有 $z_{ik}=0$。

我们假设有如下的因子化的共轭先验

\[ p (\theta) = \text {Dir}(\pi \mid \alpha_0) \pi_k \mathcal {N} (\mu_k \mid m_0, (\beta_0 \Lambda_k)^{-1}) \text {Wi}(\Lambda_k \mid L_0, \nu_0) \tag {21.122}\label {eqn:21.122} \]

此处 $\Lambda_k$ 是类别 $k$ 的精度矩阵。下标 0 对应着先验参数，并假设所有的类别共享同样的先验参数。对于混合权重，则使用对称先验 $\pmb {\alpha}_0=\alpha_0\pmb {1}$。

该模型后验分布 $p (z,\theta \mid {D})$ 的确切形式是由 $K^N$ 个对应所有可能分类 / 标注方法 $\{z\}$ 的分布混合而成的。此处仅在其中某个标注模态附近做近似。考虑基本的变分贝叶斯近似后验： $$ p (\theta,z_{1:N} \mid D) \approx q (\theta)\prod_i q (z_i) \tag {21.123}\label {eqn:21.123} $$

此处我们还没有选定函数 $q$ 的具体形式，它们会由似然和先验的形式共同决定。接下来我们会证明有如下的最优形式

\[\begin{split} \begin {align*} q (z,\theta) &= q (z \mid \theta) q (\theta) = \left [ \prod_i \text {Cat}(z_i \mid r_i ) \right] \cdot \tag {21.124}\label {eqn:21.124} \\ &~~\left [ \text {Dir}(\pi \mid \alpha) \prod_k \mathcal {N} ( \mu_k \mid m_k, (\beta_k \Lambda_k)^{-1}) \text {Wi}(\Lambda_k \mid L_k,\nu_k) \right] \tag {21.125}\label {eqn:21.125} \end {align*} \end{split}\]

（注意上式中不含下标 0 因而都是后验而非先验的参数）。接下来会给出的这些变分参数的更新公式。

21.6.2.2 变分 E-step¶

后验期望 $q (z)$ 的形式可以通过考虑数据似然函数：忽略那些不含 $z$ 的那些项目，并将剩下来的项目对除去 $z$ 的所有隐变量求期望，亦即

\[\begin{split} \begin {align*} \log q (z) &=\mathbb{E}_{q (\theta)} [ \log p (x,z,\theta) ] + \text{ const} \tag {21.126}\label {eqn:21.126} \\ &= \sum_k \sum_i z_{ik} +\text{ const} \tag {21.127}\label {eqn:21.127} \end {align*} \end{split}\]

并定义

\[\begin{split} \begin {align*} \log \rho_{ik} &\triangleq\mathbb{E}_{q (\theta)}[ \log \pi_k ] + \frac {1}{2}\mathbb{E}_{q (\theta)}[ \log |\Lambda_k| ] -\frac {D}{2}\log (2\pi) \\ &~~ - \frac {1}{2}\mathbb{E}_{q (\theta)} [(x_i-\mu_k)^T \Lambda_k (x_i - \mu_k)] \tag {21.128}\label {eqn:21.128} \end {align*} \end{split}\]

由于已经有 $q (\pi) = \text {Dir}(\pi)$，可以得出

\[ \log \tilde {\pi} \triangleq\mathbb{E} [\log \pi_k] = \psi (\alpha_k) - \psi ( \sum_{k'} \alpha_{k'} ) \tag {21.129}\label {eqn:21.129} \]

此处 $\psi ()$ 是双伽马函数。（详细推导见 Exer21.5\ref {exer:21.5}）接下来利用如下事实

\[ q (\mu_k, \Lambda_k) = \mathcal {N} ( \mu_k \mid m_k, (\beta_k \Lambda_k)^{-1}) \text {Wi}(\Lambda_k \mid L_k, \nu_k) \tag {21.130}\label {eqn:21.130} \]

然后得出

\[ \log \tilde {\Lambda_k} \triangleq \mathbb{E} [\log|\Lambda_k|] = \sum_{j=1}^D \psi\left ({v_k + 1 -j \over 2}\right) + D\log 2 + \log |\Lambda_k| \tag {21.131}\label {eqn:21.131} \]

最后，对二次项求期望可以得出

\[ \mathbb{E} [ (x_i - \mu_k)^T \Lambda_k (x_i - \mu_k) ] = D\beta_k^{-1} + \nu_k (x_i - m_k)^T \Lambda_k (x_i - m_k) \tag {21.132}\label {eqn:21.132} \]

综合考虑以上种种，得出

\[ r_{ik} \propto \tilde {\pi_k}\tilde {\Lambda}_k^{\frac {1}{2}} \exp\left ( -{D\over 2\beta_k} - {\nu_k \over 2}(x_i - m_k)^T \Lambda_k (x_i - m_k) \right) \tag {21.133}\label {eqn:21.133} \]

将以上形式与常用的期最算法作对比

\[ r_{ik}^{EM} \propto \hat {\pi}_k |\hat {\Lambda}|_k^{1/2} \exp\left (-{1\over 2}(x_i - \hat {\mu}_k)^T \hat {\Lambda}_k (x_i - \hat {\mu}_k) \right) \tag {21.134}\label {eqn:21.134} \]

两者的区别会晚些在 Sec21.6.1.7\ref {sec:21.6.1.7} 中探讨。

21.6.2.3 变分 M-step¶

按照平均场的标准菜谱，你可以炒出

\[\begin{split} \begin {align*} \log q (\theta) &= \log p (\pi) + \sum_k \log p (\mu_k, \Lambda_k) + \sum_i\mathbb{E}_{q (z)} [ \log p (z_i \mid \pi)] \\ &~~+ \sum_k \sum_i\mathbb{E}_{q (z)}[ z_{ik} ] \log \mathcal {N}(x_i \mid \mu_k, \Lambda_k^{-1}) + \text{ const} \tag {21.135}\label {eqn:21.135} \end {align*} \end{split}\]

这玩意可以因子分解成

\[ q (\theta) = q (\pi) \prod_k q (\mu_k, \Lambda_k) \tag {21.136}\label {eqn:21.136} \]

收集所有含 $\pi$ 的项可以得出

\[ \log q (\pi) = (\alpha_0 - 1 ) \sum_k \log \pi_k + \sum_k \sum_i r_{ik} \log \pi_k + \text{ const} \tag {21.137}\label {eqn:21.137} \]

两边取指数，可以看出这是一个狄利克雷分布 (?-_-)

\[\begin{split} \begin {align*} q (\pi) &= \text {Dir} ( \pi \mid \alpha) &\tag {21.138}\label {eqn:21.138} \\ \alpha_k &= \alpha_0 + N_k &\tag {21.139}\label {eqn:21.139} \\ N_k &= \sum_i r_{ik} &\tag {21.138}\label {eqn:21.140} \end {align*} \end{split}\]

收集所有含 $\mu_k,\,\Lambda_k$ 的项目，可以得出

\[\begin{split} \begin {align*} q (\mu_k, \Lambda_k) &= \mathcal {N}(\mu_k \mid m_k, (\beta_k \Lambda_k)^{-1} ) \text {Wi}(\Lambda_k \mid L_k,\nu_k) &\tag {21.141}\label {eqn:21.141} \\ \beta_k &= \beta_0 + N_k &\tag {21.142}\label {eqn:21.142} \\ m_k &= (\beta_0 m_0 + N_k \bar {x_k})/\beta_k &\tag {21.143}\label {eqn:21.143} \\ L_k^{-1} &= L_0^{-1} + N_k S_k + \frac {\beta_0 N_k}{\beta_0 + N_k} (\bar {x}_k -m_0)(\bar {x}_k - m_0)^T &\tag {21.144}\label {eqn:21.144} \\ \nu_k &= \nu_0 + N_k + 1 &\tag {21.145}\label {eqn:21.145} \\ \bar {x}_k &= \frac {1}{N_k} \sum_i r_{ik} x_i &\tag {21.146}\label {eqn:21.146} \\ S_k &= \frac {1}{N_k} \sum_i r_{ik} (x_i - \bar {x}_k)(x_i - \bar {x}_k)^T &\tag {21.147}\label {eqn:21.147} \end {align*} \end{split}\]

以上结果和最大似然后验点估计（见 Sec11.4.2.8\ref {sec:11.4.2.8}) 的最大化步是很相似的，不过这里我们计算的是模型参数 $\theta$ 之后验分布的超参数，而不是最大似然点估计。

21.6.2.4 边缘似然（证据）下界¶

该算法设法最大化的是如下的似然函数下界

\[ \mathcal {L} = \sum_z \int q (z,\theta) \log \frac { p (x,z,\theta) }{ q (z, \theta) } d\theta \le \log p (\mathcal {D}) \tag {21.148}\label {eqn:21.148} \]

该量应该随着每次迭代单调增加，如图 Fig21.7\ref {fig:21.7} 所示。不幸的是，该界限的推导过程有些杂乱，因为需要同时取未正规化的对数后验的期望并计算 $q$ 的熵。该函数的推导细节（其实和 Sec21.5.1.6\ref {sec:21.5.1.6} 很相似）留作习题 Exer21.4\ref {exer:21.4}。

21.6.2.5 后验预测性分布¶

我们已经证得近似后验有如下形式

\[ q (\theta) = \text {Dir}(\pi\mid\alpha)\prod_k \mathcal {N}( \mu_k \mid m_k, (\beta_k\Lambda_k)^{-1}) \text {Wi}(\Lambda_k \mid L_k,\nu_k) \tag {21.149}\label {eqn:21.149} \]

因此后验预测分布可以参照 Sec4.6.3.6\ref {sec:4.6.3.6} 的结果作如下近似

\[\begin{split} \begin {align*} p (x\mid \mathcal {D}) &\approx \sum_z \int p (x\mid z,\theta) p (z\mid \theta) q (\theta) d \theta &\tag {21.150}\label {eqn:21.150} \\ &= \sum_k \int \pi_k \mathcal {N}(x\mid\mu_k, \Lambda_k^{-1}) q (\theta) d\theta &\tag {21.151}\label {eqn:21.151} \\ &= \sum_k \frac {\alpha_k}{\sum_{k'}\alpha_{k'}} \mathcal {T}(x\mid m_k, M_k,\nu_k + 1 -D) &\tag {21.152}\label {eqn:21.152} \\ M_k &= {(v_k + 1 -D)\beta_k\over 1+\beta_k} L_k &\tag {21.153}\label {eqn:21.153} \end {align*} \end{split}\]

这实际上是一个学生 t - 分布的加权混合。如果我们现在改作一个 Plug-In 近似，则会得出高斯分布的加权和。

21.6.2.6 使用 VBEM 做模型选择¶

为变分贝叶斯模型选择分类数 $K$ 的最简单的方法是多拟合几个模型，然后用对数边际似然的变分下界 $\mathcal {L}(K)\le \log p (\mathcal {D}\mid K)$ 来近似 $p (K\mid \mathcal {D})$

\[ p (K\mid \mathcal {D}) \approx \frac {\exp^{\mathcal {L}(K)}}{\sum_{K'}\exp^{\mathcal {L}(K')}} \tag {21.154}\label {eqn:21.154} \]

实际上，该下界还要考虑到参数的不可分辨性而做一些修改 (Sec11.3.1\ref {sec:11.3.1})。特别地，尽管分贝能够近似模型参数后验附近的概率密度 $K$, 它事实上只能对付其中一个局域模态。对于 $K$ 组分的混合分布，有 $K!$ 个由于标注的可置换性而等价的模态。因此我们应该使用如下修正： $\log p (\mathcal {D}\mid K) \approx \mathcal {L} + \log (K!)$。

21.6.1.7 VBEM 的自动稀疏诱导效应¶

尽管变分贝叶斯提供了一个还凑活的关于边际似然的近似（好过 BIC,@Beal&Ghahramani2006)，它仍然需要为多个分类数 $K$ 各自拟合一个模型。另外一个选择是直接拟合仅仅一个 $K$ 数巨大而 $\alpha_0$ 巨小的模型 $\alpha_0 \ll 1$。从图 2.14d\ref {fig:2.14d] 可以看出，此时混合权重 $\pi$ 的先验在单纯形的顶点附近有 “尖峰”，因而会偏爱一个稀疏的混合权重。

在常规期最算法中，混合权重的最大后验点估计具有形式 $\hat {\pi}^k \propto (\alpha_k - 1 )$, 此处 $\alpha_k=\alpha_0 + N_k$。不幸的是，这玩意在 $\alpha_0=0$ 且 $N_k=0$ 会取到负值（@Figueiredo&Jain2002）。不过呢，在变分贝叶斯期最中相应的后验估计是

\[ \tilde {\pi}_k = \frac {\exp [\Psi (\alpha_k)]}{\exp [ \Psi ( \sum_{k'}\alpha_{k'} ) ]} \tag {21.155}\label {eqn:21.155} \]

对于 $x>1$ 有近似 $\exp {\Psi (x)} \approx x-0.5$。因而如果 $\alpha_k=0$, 那么就相当于我们在计算 $\tilde {\pi}_k$ 的时候从后验个数中减去了 0.5。该效果对于没几个有分量成员的的小分类来说是更加严重的（就像累退税那样）。其最终结果就是，随着不断地迭代，小分类变得越来越空而成员众多的分类变得越来越大。这也被称作是马太效应，这在 Sec25.2\ref {sec:25.2} 讨论狄利克雷过程混合模型的时候会被再次提起。

该自动剪枝方法在图 21.8\ref {fig:21.8} 中进行了展示。一个 6 分类的高斯分布被用来拟合 OldFaithful 数据集，但是数据实际上只 “需要” 两个分类，因而多出来的那些分类就被 “扼杀” 了。在该例子里，我们采用了 $\alpha_0 = 0.001$。如果我们用一个更大的 $\alpha_0$，我们不一定会观察到稀疏化现象。在图 21.9\ref {fig:21.9} 中我们展示了 $q (\alpha\mid \mathcal {D})$ 在不同迭代步的形状。可以见得多余的成分逐渐灭绝了。这相较于对分类数 $K$ 的离散搜索是一个更有效率的替代方案。

21.7 变分消息传递和 VIBES¶

我们已经看到，平均场方法，至少是完全因子分解的方法，都非常相似:只需计算每个节点的完全条件，并平均出相邻节点。这与吉布斯采样(第24.2节)非常相似，只是方程的推导通常要做更多的工作。幸运的是，可以推导出一组通用的更新方程，适用于所有CPD都在指数族中并且所有父节点都具有共轭分布的任何DGM(Ghahramani和Beal 2001)。(参见(Wand等人，2011)最近的扩展，以处理非共轭先验。)然后可以扫描图，以类似于吉布斯抽样的方式一次更新一个节点。这被称为变分消息传递或VMP (Winn和Bishop 2005)，并且已经在开源程序VIBES5中实现。这是一个模拟BUGS的VB程序，BUGS是Gibbs抽样的通用程序，在24.2.6节中讨论过。

VMP/均值场最适合于一个或多个隐藏节点是连续的推断(例如，当执行“贝叶斯学习”时)。对于所有隐藏节点都是离散的模型，可以使用更精确的近似推断算法，正如我们在第22章中讨论的那样

21.8 局部变分界¶

到目前为止，我们一直关注平均场推断，这是一种基于最小化 $\mathbb{K}(q||\tilde p)$ ，其中 $q$ 为近似后验，通常是因子化的，而 $\tilde p$ 是精确后验（通常是未归一化的）。然而，还存在另一种变分推断方法，用更简单的函数代替联合分布中的特定项，以简化后验的计算。这种方法有时被称为局部变分近似，因为只修改了模型的一部分，而不像平均场那样做全局近似。本节中将研究这种方法的几个例子。

21.8.1 动机¶

在解释如何导出局部变分界之前，先给出一些说明局部变分推断用处的例子。

（1）变分逻辑回归¶

对于高斯先验下的多类别逻辑回归模型，考虑其参数后验的近似问题。高斯（拉普拉斯）近似是其中一种选择，如第 8.4.3 节所述。但变分法可以产生更精确的后验近似，因为它具有可调参数。而且变分方法单调地优化了数据边缘似然（证据）下界。

要明白为什么需要一个界，注意模型的似然函数（注意不是边缘似然）可以写成: $$ p(\mathbf{y} \mid \mathbf{X}, \mathbf{w})=\prod_{i=1}^{N} \exp \left[\mathbf{y}_{i}^{T} \boldsymbol{\eta}_{i}-\operatorname{lse}\left(\boldsymbol{\eta}_{i}\right)\right] $$ 其中 $\boldsymbol{\eta}{i}=\mathbf{X}{i} \mathbf{w}{i}=\left[\mathbf{x}{i}^{T} \mathbf{w}{1}, \ldots, \mathbf{x}{i}^{T} \mathbf{w}{M}\right]$, $M=C-1$ （为独立性设 $\mathbf{w}{C}=\mathbf{0}$ ）,并且： $$ \operatorname{lse}\left(\boldsymbol{\eta}_{i}\right) \triangleq \log \left(1+\sum_{m=1}^{M} e^{\eta_{i m}}\right) $$ 此模型最大的问题是：似然函数与高斯先验不共轭。如果能够计算一个与该似然相似的“类高斯”下限，则似乎有可能解析地生成近似的高斯后验。

（2）多任务学习¶

贝叶斯统计在逻辑回归中的一个重要作用是：有多个相关的分类器要拟合，而我们希望在每个分类器的模型参数之间共享信息，这要求维持参数的后验分布，因此我们同时有一个置信度的度量和一个值的估计。我们可以将上述变分方法嵌入到一个更大的层次模型中，以便执行这种多任务学习，例如，如（Braun and McAuliffe ，2010）所述（即多任务的共享参数学习）。

（3）离散因子分析¶

当将因子分析模型拟合到离散型数据时，另一种有用的情况出现了。这个模型就像多项式逻辑回归一样，只是输入变量是隐变量。我们需要对隐变量和回归权重进行推断。为简单起见，可能会对权重进行点估计，然后积分出隐变量。可以用变分 EM 来做，在 E-step 中使用变分界。详见第 12.4 节。

（4）相关主题模型¶

主题模型是文本文档和其他形式离散数据的隐变量模型；详见第 27.3 节。通常假设主题上的分布具有狄利克雷先验，但一个更强大的模型（称为相关主题模型）使用高斯先验，可以更容易地建模相关性（详见第 27.4.1 节）。不幸的是，这会涉及到 lse 函数。不过，我们可以在变分 EM 算法的上下文中使用变分界来解决这个问题。

21.8.2 log-sum-exp 函数的Bohning二次界¶

上述所有例子都需要处理高斯先验与多项式似然相乘的问题；而由于存在 log-sum-exp 项(lse 项)，处理起来很困难。本节导出了获得此类似然“类高斯”下界的方法。

机器学习：概率视角

21 变分推断

Contents