PRML学习总结（10）——Approximate Inference

在有些情况下，精确推断是难以实现，只能寻求近似推断，根据近似⽅法依赖于随机近似还是确定近似，⽅法⼤体分为两⼤类。随机⽅法，例如第11章介绍的马尔科夫链蒙特卡罗⽅法，使得贝叶斯⽅法能够在许多领域中⼴泛使⽤。这些⽅法通常具有这样的性质：给定⽆限多的计算资源，它们可以⽣成精确的结果，近似的来源是使⽤了有限的处理时间。在实际应⽤中，取样⽅法需要的计算量会相当⼤，经常将这些⽅法的应⽤限制在了⼩规模的问题中。并且，判断⼀种取样⽅法是否⽣成了服从所需的概率分布的独⽴样本是很困难的。
本章中，我们介绍了⼀系列的确定性近似⽅法，有些⽅法对于⼤规模的数据很适⽤。这些⽅法基于对后验概率分布的解析近似，例如通过假设后验概率分布可以通过⼀种特定的⽅式分解，或者假设后验概率分布有⼀个具体的参数形式，例如⾼斯分布。对于这种情况，这些⽅法永远⽆法⽣成精确的解，因此这些⽅法的优点和缺点与取样⽅法是互补的。

10.1 Variational Inference

假设我们有⼀个纯粹的贝叶斯模型，其中每个参数都有⼀个先验概率分布。这个模型也可以有潜在变量以及参数，我们会把所有潜在变量和参数组成的集合记作$\mathbf{Z}$。类似地，我们会把所有观测变量的集合记作$\mathbf{X}$。例如，我们可能有$N$个独⽴同分布的数据，其中 X = { x 1 , … , x N } \mathbf{X}=\left\{\mathbf{x}_{1}, \dots, \mathbf{x}_{N}\right\} X={x1​,…,xN​}且 Z = { z 1 , … , z N } \mathbf{Z}=\left\{\mathbf{z}_{1}, \dots, \mathbf{z}_{N}\right\} Z={z1​,…,zN​}。我们的概率模型确定了联合分布$p(\mathbf{X},\mathbf{Z})$，我们的目标是寻找后验分布$p(\mathbf{Z}|\mathbf{X})$和模型evidence$p(\mathbf{X})$的近似分布，之前在EM算法中介绍了$$\ln p(\mathbf{X})=\mathcal{L}(q)+\mathrm{K}\mathrm{L}(q\Vert p)$$其中$$\begin{array}{rl}\mathcal{L}(q)& =\int q(\mathbf{Z})\ln \left\{\frac{p(\mathbf{X},\mathbf{Z})}{q(\mathbf{Z})}\right\}\mathrm{d}\mathbf{Z}\\ \mathrm{K}\mathrm{L}(q\Vert p)& =-\int q(\mathbf{Z})\ln \left\{\frac{p(\mathbf{Z}|\mathbf{X})}{q(\mathbf{Z})}\right\}\mathrm{d}\mathbf{Z}\end{array}$$与EM不同的是，参数$\mathbf{\theta }$一起包括在$\mathbf{Z}$中。和之前的分析一样，当最大化lower bound$\mathcal{L}(q)$时，就是在最小化KL部分，也就是说，当$q(\mathbf{Z})$足够自由时，KL部分就会消失，从而$q(\mathbf{Z})=p(\mathbf{Z}|\mathbf{X})$。但是现在我们考虑的模型十分复杂，真是后验是没法计算的。于是，我们转⽽考虑概率分布$q(\mathbf{Z})$的⼀个受限制的类别，然后寻找这个类别中使得KL散度达到最⼩值的概率分布。

10.1.1 Factorized distributions

下面介绍一种平均场理论，也就是将$q(\mathbf{Z})$分解为$$q(\mathbf{Z})=\prod _{i=1}^M{q}_i\left({\mathbf{Z}}_i\right)$$在这种分解的限制下，优化ELBO，只考虑$q_j$$$\begin{array}{rl}\mathcal{L}(q)& =\int \prod _i{q}_i\left\{\ln p(\mathbf{X},\mathbf{Z})-\sum _i\ln q_i\right\}\mathrm{d}\mathbf{Z}\\ & =\int q_j\left\{\int \ln p(\mathbf{X},\mathbf{Z})\prod _{i̸=j}{q}_i\mathrm{d}{\mathbf{Z}}_i\right\}\mathrm{d}{\mathbf{Z}}_j-\int q_j\ln q_j\mathrm{d}{\mathbf{Z}}_j+\text{ const }\\ & =\int q_j\ln \tilde{p}\left(\mathbf{X},{\mathbf{Z}}_j\right)\mathrm{d}{\mathbf{Z}}_j-\int q_j\ln q_j\mathrm{d}{\mathbf{Z}}_j+\text{ const }\end{array}$$当我们固定 { q i ≠ j } \left\{q_{i \neq j}\right\} {qi̸​=j​}时，可以发现最大ELBO就是在$q_j\left({\mathbf{Z}}_j\right)=\tilde{p}\left(\mathbf{X},{\mathbf{Z}}_j\right)$。$$\ln q_j^{\star }\left({\mathbf{Z}}_j\right)={\mathbb{E}}_{i̸=j}[\ln p(\mathbf{X},\mathbf{Z})]+\text{ const. }$$注意的是，const就是一个配比因子，归一化。$$q_j^{\star }\left({\mathbf{Z}}_j\right)=\frac{\exp \left({\mathbb{E}}_{i̸=j}[\ln p(\mathbf{X},\mathbf{Z})]\right)}{\int \exp \left({\mathbb{E}}_{i̸=j}[\ln p(\mathbf{X},\mathbf{Z})]\right)\mathrm{d}{\mathbf{Z}}_j}$$得到这样的结果，我们就可以采用迭代的方式得到最终结果。

10.1.2 Properties of factorized approximations

我们的变分推断的⽅法基于的是真实后验概率分布的分解近似。让我们现在考虑⼀下使⽤分解概率分布的⽅式近似⼀个⼀般的概率分布的问题。⾸先，我们讨论使⽤分解的⾼斯分布近似⼀个⾼斯分布的问题，这会让我们认识到在使⽤分解近似时会引⼊的不准确性有哪些类型。需要近似的二维高斯分布为$p(\mathbf{z})=\mathcal{N}\left(\mathbf{z}|\mathbf{\mu },{\mathbf{\Lambda }}^{-1}\right)$，利用$q(\mathbf{z})=q_1\left(z_1\right)q_2\left(z_2\right)$近似它。$$\begin{array}{rl}\ln q_1^{\star }\left(z_1\right)& ={\mathbb{E}}_{z_2}[\ln p(\mathbf{z})]+\text{ const }\\ & ={\mathbb{E}}_{z_2}\left[-\frac{1}{2}{\left(z_1-{\mu }_1\right)}^2{\Lambda }_{11}-\left(z_1-{\mu }_1\right){\Lambda }_{12}\left(z_2-{\mu }_2\right)\right]+\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\\ & =-\frac{1}{2}{z}_1^2{\Lambda }_{11}+z_1{\mu }_1{\Lambda }_{11}-z_1{\Lambda }_{12}\left(\mathbb{E}\left[z_2\right]-{\mu }_2\right)+\text{ const. }\end{array}$$上式的右边刚好关于$z_1$的二次函数，则说明$q^{\star }\left(z_1\right)$刚好也是高斯分布，需要说明的是，平均场理论，我们并没有规定各个分解后的分布的类型，这儿就能推出分布是个什么类型。$$q^{\star }\left(z_1\right)=\mathcal{N}\left(z_1|m_1,{\Lambda }_{11}^{-1}\right)$$其中$$m_1={\mu }_1-{\Lambda }_{11}^{-1}{\Lambda }_{12}\left(\mathbb{E}\left[z_2\right]-{\mu }_2\right)$$同理可得$$q_2^{\star }\left(z_2\right)=\mathcal{N}\left(z_2|m_2,{\Lambda }_{22}^{-1}\right)$$$$m_2={\mu }_2-{\Lambda }_{22}^{-1}{\Lambda }_{21}\left(\mathbb{E}\left[z_1\right]-{\mu }_1\right)$$可以发现$q^{\star }\left(z_1\right)$和$q_2^{\star }\left(z_2\right)$相互耦合，即互相包括了期望值，因此可以采用迭代的方式进行求解。通过分析可以发现，最终均值就是$\mathbb{E}\left[z_1\right]={\mu }_1$，$\mathbb{E}\left[z_2\right]={\mu }_2$。结果如图10.2(a)，我们看到，均值被正确地描述了，但是$q(\mathbf{z})$的⽅差由$p(\mathbf{z})$的最⼩⽅差的⽅向所确定，沿着垂直⽅向的⽅差被强烈地低估了。这是⼀个⼀般的结果，即分解变分近似对后验概率分布的近似倾向于过于紧凑。
作为⽐较，假设我们最⼩化相反的Kullback-Leibler散度$\mathrm{K}\mathrm{L}(p\Vert q)$。正如我们将看到的那样，这种形式的KL散度被⽤于另⼀种近似推断的框架中，这种框架被称为期望传播（expectation propagation）。于是，我们考虑⼀般的最⼩化$\mathrm{K}\mathrm{L}(p\Vert q)$的问题，其中$q(\mathbf{Z})$是形式为$q(\mathbf{Z})={\prod }_{i=1}^M{q}_i\left({\mathbf{Z}}_i\right)$的分解近似。这样，KL散度可以写成$$\mathrm{K}\mathrm{L}(p\Vert q)=-\int p(\mathbf{Z})\left[\sum _{i=1}^M\ln q_i\left({\mathbf{Z}}_i\right)\right]\mathrm{d}\mathbf{Z}+\text{ const }$$$$q_j^{\star }\left({\mathbf{Z}}_j\right)=\int p(\mathbf{Z})\prod _{i̸=j}\mathrm{d}{\mathbf{Z}}_i=p\left({\mathbf{Z}}_j\right)$$以上结果为一个解析解，不需要迭代，结果如图10.2(b)。我们再⼀次看到，对均值的近似是正确的，但是它把相当多的概率质量放到了实际上具有很低的概率的变量空间区域中。之所以会造成这样的不同$$\mathrm{K}\mathrm{L}(q\Vert p)=-\int q(\mathbf{Z})\ln \left\{\frac{p(\mathbf{Z})}{q(\mathbf{Z})}\right\}\mathrm{d}\mathbf{Z}$$当$p(\mathbf{Z})$为0的区域，$q(\mathbf{Z})$必须也为0，否则会让上式引入一个很大的正数贡献。而反KL散度，当$p(\mathbf{Z})$为0的区域，当$q(\mathbf{Z})$不为0时，不会引入正数贡献，并且$p(\mathbf{Z})$不为0的区域，$q(\mathbf{Z})$必须不为0。
当用单高斯近似多元高斯时，基于最⼩化$\mathrm{K}\mathrm{L}(q\Vert p)$的变分⽅法倾向于找到这些峰值中的⼀个。相反，如果我们最⼩化$\mathrm{K}\mathrm{L}(p\Vert q)$，那么得到的近似会在所有的均值上取平均。在混合模型问题中，这种⽅法会给出较差的预测分布（因为两个较好的参数值的平均值通常不是⼀个较好的参数值）。可以使⽤$\mathrm{K}\mathrm{L}(p\Vert q)$定义⼀个有⽤的推断步骤，但是这需要⼀种与这⾥讨论的内容相当不同的⽅法。当我们讨论期望传播的时候，我们会仔细讨论这⼀点。
以上两个KL散度其实为alpha家族$${\mathrm{D}}_{\alpha }(p\Vert q)=\frac{4}{1-{\alpha }^2}\left(1-\int p(x{)}^{(1+\alpha )/2}q(x{)}^{(1-\alpha )/2}\mathrm{d}x\right)$$

10.1.3 Example: The univariate Gaussian

考虑一个高斯分布，我们的目标是去近似均值和精度的后验分布。似然为$$p(\mathcal{D}|\mu ,\tau )={\left(\frac{\tau }{2\pi }\right)}^{N/2}\exp \left\{-\frac{\tau }{2}\sum _{n=1}^N{\left(x_n-\mu \right)}^2\right\}$$引入共轭先验为$$\begin{array}{rl}p(\mu |\tau )& =\mathcal{N}\left(\mu |{\mu }_0,{\left({\lambda }_0\tau \right)}^{-1}\right)\\ p(\tau )& =\mathrm{Gam}\left(\tau |a_0,b_0\right)\end{array}$$对于这个问题我们能够计算出精确的后验分布，这儿我们采用平均场变分推理$$q(\mu ,\tau )=q_{\mu }(\mu )q_{\tau }(\tau )$$可得$$\begin{array}{rl}\ln q_{\mu }^{\star }(\mu )& ={\mathbb{E}}_{\tau }[\ln p(\mathcal{D}|\mu ,\tau )+\ln p(\mu |\tau )]+\text{ const }\\ & =-\frac{\mathbb{E}[\tau ]}{2}\left\{{\lambda }_0{\left(\mu -{\mu }_0\right)}^2+\sum _{n=1}^N{\left(x_n-\mu \right)}^2\right\}+\text{ const. }\end{array}$$可以发现该分布为一个高斯分布，$$\mathcal{N}\left(\mu |{\mu }_N,{\lambda }_N^{-1}\right)$$其中$$\begin{array}{rl}{\mu }_N& =\frac{{\lambda }_0{\mu }_0+N\overline{x}}{{\lambda }_0+N}\\ {\lambda }_N& =\left({\lambda }_0+N\right)\mathbb{E}[\tau ]\end{array}$$注意当$N\to \infty$时，就变为了MLE。同样地$$\begin{array}{rl}\ln q_{\tau }^{\star }(\tau )=& {\mathbb{E}}_{\mu }[\ln p(\mathcal{D}|\mu ,\tau )+\ln p(\mu |\tau )]+\ln p(\tau )+\text{ const }\\ =& \left(a_0-1\right)\ln \tau -b_0\tau +\frac{N}{2}\ln \tau \\ & -\frac{\tau }{2}{\mathbb{E}}_{\mu }\left[\sum _{n=1}^N{\left(x_n-\mu \right)}^2+{\lambda }_0{\left(\mu -{\mu }_0\right)}^2\right]+\text{ const }\end{array}$$此时为一个gamma分布$\mathrm{Gam}\left(\tau |a_N,b_N\right)$$$\begin{array}{rl}{a}_N& =a_0+\frac{N}{2}\\ b_N& =b_0+\frac{1}{2}{\mathbb{E}}_{\mu }\left[\sum _{n=1}^N{\left(x_n-\mu \right)}^2+{\lambda }_0{\left(\mu -{\mu }_0\right)}^2\right]\end{array}$$需要说明的是，我们并没有事先规定各个分布的类型，而是很自然地推导出。由于以上两个分布相互耦合，同样可以采取迭代的方式

10.1.4 Model comparison

10.2. Illustration: Variational Mixture of Gaussians

这儿讨论的是全贝叶斯的GMM

$$p(\mathbf{Z}|\mathbf{\pi })=\prod _{n=1}^N\prod _{k=1}^K{\pi }_k^{z_{nk}}$$$$p(\mathbf{X}|\mathbf{Z},\mathbf{\mu },\mathbf{\Lambda })=\prod _{n=1}^N\prod _{k=1}^K\mathcal{N}{\left({\mathbf{x}}_n|{\mathbf{\mu }}_k,{\mathbf{\Lambda }}_k^{-1}\right)}^{z_{nk}}$$$$p(\mathbf{\pi })=\mathrm{Dir}\left(\mathbf{\pi }|{\mathbf{\alpha }}_0\right)=C\left({\mathbf{\alpha }}_0\right)\prod _{k=1}^K{\pi }_k^{{\alpha }_0-1}$$$$\begin{array}{rl}p(\mathbf{\mu },\mathbf{\Lambda })& =p(\mathbf{\mu }|\mathbf{\Lambda })p(\mathbf{\Lambda })\\ & =\prod _{k=1}^K\mathcal{N}\left({\mathbf{\mu }}_k|{\mathbf{m}}_0,{\left({\beta }_0{\mathbf{\Lambda }}_k\right)}^{-1}\right)\mathcal{W}\left({\mathbf{\Lambda }}_k|{\mathbf{W}}_0,{\nu }_0\right)\end{array}$$

10.2.1 Variational distribution

首先写成该模型的联合分布$$p(\mathbf{X},\mathbf{Z},\mathbf{\pi },\mathbf{\mu },\mathbf{\Lambda })=p(\mathbf{X}|\mathbf{Z},\mathbf{\mu },\mathbf{\Lambda })p(\mathbf{Z}|\mathbf{\pi })p(\mathbf{\pi })p(\mathbf{\mu }|\mathbf{\Lambda })p(\mathbf{\Lambda })$$利用平均场变分$$q(\mathbf{Z},\mathbf{\pi },\mathbf{\mu },\mathbf{\Lambda })=q(\mathbf{Z})q(\mathbf{\pi },\mathbf{\mu },\mathbf{\Lambda })$$首先更新第一个因子$$\ln q^{\star }(\mathbf{Z})={\mathbb{E}}_{\mathbf{\pi },\mathbf{\mu },\mathbf{\Lambda }}[\ln p(\mathbf{X},\mathbf{Z},\mathbf{\pi },\mathbf{\mu },\mathbf{\Lambda })]+\text{ const. }$$则$$\ln q^{\star }(\mathbf{Z})={\mathbb{E}}_{\pi }[\ln p(\mathbf{Z}|\mathbf{\pi })]+{\mathbb{E}}_{\mathbf{\mu },\mathbf{\Lambda }}[\ln p(\mathbf{X}|\mathbf{Z},\mathbf{\mu },\mathbf{\Lambda })]+\text{ const. }$$进一步$$\ln q^{\star }(\mathbf{Z})=\sum _{n=1}^N\sum _{k=1}^K{z}_{nk}\ln {\rho }_{nk}+\text{ const }$$其中$$\begin{array}{rl}\ln {\rho }_{nk}=& \mathbb{E}\left[\ln {\pi }_k\right]+\frac{1}{2}\mathbb{E}\left[\ln \left|{\mathbf{\Lambda }}_k\right|\right]-\frac{D}{2}\ln (2\pi )\\ & -\frac{1}{2}{\mathbb{E}}_{{\mathbf{\mu }}_k,{\mathbf{\Lambda }}_k}\left[{\left({\mathbf{x}}_n-{\mathbf{\mu }}_k\right)}^{\mathrm{T}}{\mathbf{\Lambda }}_k\left({\mathbf{x}}_n-{\mathbf{\mu }}_k\right)\right]\end{array}$$从而$$q^{\star }(\mathbf{Z})\propto \prod _{n=1}^N\prod _{k=1}^K{\rho }_{nk}^{z_{nk}}$$归一化后$$q^{\star }(\mathbf{Z})=\prod _{n=1}^N\prod _{k=1}^K{r}_{nk}^{z_{nk}}$$其中$$r_{nk}=\frac{{\rho }_{nk}}{{\sum }_{j=1}^K{\rho }_{nj}}$$有$$\mathbb{E}\left[z_{nk}\right]=r_{nk}$$预先定义$$\begin{array}{rl}{N}_k& =\sum _{n=1}^N{r}_{nk}\\ {\overline{\mathbf{x}}}_k& =\frac{1}{N_k}\sum _{n=1}^N{r}_{nk}{\mathbf{x}}_n\\ {\mathbf{S}}_k& =\frac{1}{N_k}\sum _{n=1}^N{r}_{nk}\left({\mathbf{x}}_n-{\overline{\mathbf{x}}}_k\right){\left({\mathbf{x}}_n-{\overline{\mathbf{x}}}_k\right)}^{\mathrm{T}}\end{array}$$考虑第二个因子$q(\mathbf{\pi },\mathbf{\mu },\mathbf{\Lambda })$，$$\begin{array}{c}\ln q^{\star }(\mathbf{\pi },\mathbf{\mu },\mathbf{\Lambda })=\ln p(\mathbf{\pi })+{\sum }_{k=1}^K\ln p\left({\mathbf{\mu }}_k,{\mathbf{\Lambda }}_k\right)+{\mathbb{E}}_{\mathbf{Z}}[\ln p(\mathbf{Z}|\mathbf{\pi })]+{\sum }_{k=1}^K{\sum }_{n=1}^N\mathbb{E}\left[z_{nk}\right]\ln \mathcal{N}\left({\mathbf{x}}_n|{\mathbf{\mu }}_k,{\mathbf{\Lambda }}_k^{-1}\right)+\text{ const. }\end{array}$$可以看到一些项只与$\mathbf{\pi }$有关，有些与$\mathbf{\mu },\mathbf{\Lambda }$有关，说明可以将该因子进一步分解为$q(\mathbf{\pi })q(\mathbf{\mu },\mathbf{\Lambda })$，$$q(\mathbf{\pi },\mathbf{\mu },\mathbf{\Lambda })=q(\mathbf{\pi })\prod _{k=1}^{K}q\left({\mathbf{\mu }}_k,{\mathbf{\Lambda }}_k\right)$$从而可得$$\ln q^{\star }(\mathbf{\pi })=\left({\alpha }_0-1\right)\sum _{k=1}^K\ln {\pi }_k+\sum _{k=1}^K\sum _{n=1}^N{r}_{nk}\ln {\pi }_k+\text{ const }$$整理后$$q^{\star }(\mathbf{\pi })=\mathrm{Dir}(\mathbf{\pi }|\mathbf{\alpha })$$其中$${\alpha }_k={\alpha }_0+N_k$$最后，变分后验概率分布$q^{\star }\left({\mathbf{\mu }}_k,{\mathbf{\Lambda }}_k\right)$⽆法分解成边缘概率分布的乘积，但是我们总可以使⽤概率的乘积规则，将其写成$q^{\star }\left({\mathbf{\mu }}_k,{\mathbf{\Lambda }}_k\right)=q^{\star }\left({\mathbf{\mu }}_k|{\mathbf{\Lambda }}_k\right)q^{\star }\left({\mathbf{\Lambda }}_k\right)$。跟预期地一样$$q^{\star }\left({\mathbf{\mu }}_k,{\mathbf{\Lambda }}_k\right)=\mathcal{N}\left({\mathbf{\mu }}_k|{\mathbf{m}}_k,{\left({\beta }_k{\mathbf{\Lambda }}_k\right)}^{-1}\right)\mathcal{W}\left({\mathbf{\Lambda }}_k|{\mathbf{W}}_k,{\nu }_k\right)$$其中$$\begin{array}{rl}{\beta }_k& ={\beta }_0+N_k\\ {\mathbf{m}}_k& =\frac{1}{{\beta }_k}\left({\beta }_0{\mathbf{m}}_0+N_k{\overline{\mathbf{x}}}_k\right)\\ {\mathbf{W}}_k^{-1}& ={\mathbf{W}}_0^{-1}+N_k{\mathbf{S}}_k+\frac{{\beta }_0{N}_k}{{\beta }_0+N_k}\left({\overline{\mathbf{x}}}_k-{\mathbf{m}}_0\right){\left({\overline{\mathbf{x}}}_k-{\mathbf{m}}_0\right)}^{\mathrm{T}}\\ {\nu }_k& ={\nu }_0+N_k\end{array}$$之前在使用“责任”${\mathbf{r}}_{nk}$，又涉及计算$$\begin{array}{rl}{\mathbb{E}}_{{\mu }_k,{\Lambda }_k}\left[{\left({\mathbf{x}}_n-{\mathbf{\mu }}_k\right)}^{\mathrm{T}}{\mathbf{\Lambda }}_k\left({\mathbf{x}}_n-{\mathbf{\mu }}_k\right)\right]& \\ =& D{\beta }_k^{-1}+{\nu }_k{\left({\mathbf{x}}_n-{\mathbf{m}}_k\right)}^{\mathrm{T}}{\mathbf{W}}_k\left({\mathbf{x}}_n-{\mathbf{m}}_k\right)\\ \ln {\tilde{\Lambda }}_k\equiv \mathbb{E}\left[\ln \left|{\mathbf{\Lambda }}_k\right|\right]=\sum _{i=1}^D\psi \left(\frac{{\nu }_k+1-i}{2}\right)+D\ln 2+\ln \left|{\mathbf{W}}_k\right|& \\ \ln {\tilde{\pi }}_k\equiv \mathbb{E}\left[\ln {\pi }_k\right]=\psi \left({\alpha }_k\right)-\psi (\hat{\alpha })& \end{array}$$这两部分是相互耦合的，当把上面的再代入${\mathbf{r}}_{nk}$$$r_{nk}\propto {\tilde{\pi }}_k{\tilde{\Lambda }}_k^{1/2}\exp \left\{-\frac{D}{2{\beta }_k}-\frac{{\nu }_k}{2}{\left({\mathbf{x}}_n-{\mathbf{m}}_k\right)}^{\mathrm{T}}{\mathbf{W}}_k\left({\mathbf{x}}_n-{\mathbf{m}}_k\right)\right\}$$在最大似然EM中，有$$r_{nk}\propto {\pi }_k{\left|{\mathbf{\Lambda }}_k\right|}^{1/2}\exp \left\{-\frac{1}{2}{\left({\mathbf{x}}_n-{\mathbf{\mu }}_k\right)}^{\mathrm{T}}{\mathbf{\Lambda }}_k\left({\mathbf{x}}_n-{\mathbf{\mu }}_k\right)\right\}$$可以发现变分贝叶斯也是两个步迭代优化。

可以发现优化后，并不会出现GMM的奇点。我们看到，在收敛之后，只有两个分量的混合系数的期望值可以与它们的先验值（先验设置为一个很小的数）区分开。这种效果可以根据贝叶斯模型中数据拟合与模型复杂度之间的折中来定性地理解。这种模型中的复杂度惩罚的来源是参数被推离了它们的先验值。对于解释数据点没有作⽤的分量满⾜$r_{nk}\simeq 0$，从而$N_k\simeq 0$，那么有些退化为先验值${\alpha }_k\simeq {\alpha }_0$。原则上，这些分量会微⼩地适应于数据点，但是对于⼀⼤类先验分布来说，这种微⼩的调整的效果太⼩了，以⾄于⽆法在数值上看出来。对于⾼斯混合模型，后验概率分布中的混合系数的期望值为$$\mathbb{E}\left[{\pi }_k\right]=\frac{{\alpha }_k+N_k}{K{\alpha }_0+N}$$考虑某个退化的分量$N_k\simeq 0$ and ${\alpha }_k\simeq {\alpha }_0$，当${\alpha }_0\to 0,$ then $\mathbb{E}\left[{\pi }_k\right]\to 0$，此时该分量就不起作用，而${\alpha }_0\to \infty ,$ then $\mathbb{E}\left[{\pi }_k\right]\to 1/K$。之前狄利克雷分布的参数表示，当越小，则倾向让某些分量为0，从而本实验设置为${\alpha }_0=10^{-3}$。

10.2.2 Variational lower bound

在实际应⽤中，能够在重新估计期间监视模型的下界是很有⽤的，这可以⽤来检测是否收敛。它也可以为解的数学表达式和它们的软件执⾏提供⼀个有价值的检查，因为在迭代重新估计的每个步骤中，这个下界的值应该不会减⼩。我们可以进⼀步地使⽤变分下界检查更新⽅程的数学推导和它们的软件执⾏的正确性，⽅法是使⽤有限差来检查每次更新确实给出了下界的⼀个（具有限制条件的）极⼤值。$$\begin{array}{rl}\mathcal{L}=& \sum _{\mathbf{Z}}\iiint q(\mathbf{Z},\mathbf{\pi },\mathbf{\mu },\mathbf{\Lambda })\ln \left\{\frac{p(\mathbf{X},\mathbf{Z},\mathbf{\pi },\mathbf{\mu },\mathbf{\Lambda })}{q(\mathbf{Z},\mathbf{\pi },\mathbf{\mu },\mathbf{\Lambda })}\right\}\mathrm{d}\mathbf{\pi }\mathrm{d}\mathbf{\mu }\mathrm{d}\mathbf{\Lambda }\\ =& \mathbb{E}[\ln p(\mathbf{X},\mathbf{Z},\mathbf{\mu },\mathbf{\mu },\mathbf{\Lambda })]-\mathbb{E}[\ln q(\mathbf{Z},\mathbf{\pi },\mathbf{\mu },\mathbf{\Lambda })]\\ =& \mathbb{E}[\ln p(\mathbf{X}|\mathbf{Z},\mathbf{\mu },\mathbf{\Lambda })]+\mathbb{E}[\ln p(\mathbf{Z},\mathbf{\pi },\mathbf{\mu },\mathbf{\Lambda })]\\ & -\mathbb{E}[\ln q(\mathbf{Z})]-\mathbb{E}[\ln q(\mathbf{\pi })]-\mathbb{E}[\ln q(\mathbf{\mu },\mathbf{\Lambda })]\end{array}$$
对于这些期望，由于模型有共轭先验，因此变分后验分布的函数形式是已知的。通过使⽤这些分布的⼀般的参数形式，我们可以推导出下界的形式，将下界作为概率分布的参数的函数。关于这些参数最⼤化下界就会得到所需的重估计⽅程。

10.2.3 Predictive density

为了得到预测分布，即预测$\hat{\mathbf{x}}$，$$p(\hat{\mathbf{x}}|\mathbf{X})=\sum _{\hat{\mathbf{z}}}\iiint p(\hat{\mathbf{x}}|\hat{\mathbf{z}},\mathbf{\mu },\mathbf{\Lambda })p(\hat{\mathbf{z}}|\mathbf{\pi })p(\mathbf{\pi },\mathbf{\mu },\mathbf{\Lambda }|\mathbf{X})\mathrm{d}\mathbf{\pi }\mathrm{d}\mathbf{\mu }\mathrm{d}\mathbf{\Lambda }$$其中$p(\mathbf{\pi },\mathbf{\mu },\mathbf{\Lambda }|\mathbf{X})$是未知的真实后验分布，$$p(\hat{\mathbf{x}}|\mathbf{X})=\sum _{k=1}^K\iiint {\pi }_k\mathcal{N}\left(\hat{\mathbf{x}}|{\mathbf{\mu }}_k,{\mathbf{\Lambda }}_k^{-1}\right)p(\mathbf{\pi },\mathbf{\mu },\mathbf{\Lambda }|\mathbf{X})\mathrm{d}\mathbf{\pi }\mathrm{d}\mathbf{\mu }\mathrm{d}\mathbf{\Lambda }$$可以将真实后验用变分后验代替$$p(\hat{\mathbf{x}}|\mathbf{X})=\sum _{k=1}^K\iiint {\pi }_k\mathcal{N}\left(\hat{\mathbf{x}}|{\mathbf{\mu }}_k,{\mathbf{\Lambda }}_k^{-1}\right)q(\mathbf{\pi })q\left({\mathbf{\mu }}_k,{\mathbf{\Lambda }}_k\right)\mathrm{d}\mathbf{\pi }\mathrm{d}{\mathbf{\mu }}_k\mathrm{d}{\mathbf{\Lambda }}_k$$从而可得一个student分布$$p(\hat{\mathbf{x}}|\mathbf{X})=\frac{1}{\hat{\alpha }}\sum _{k=1}^K{\alpha }_k\mathrm{St}\left(\hat{\mathbf{x}}|{\mathbf{m}}_k,{\mathbf{L}}_k,{\nu }_k+1-D\right)$$其中$${\mathbf{L}}_k=\frac{\left({\nu }_k+1-D\right){\beta }_k}{\left(1+{\beta }_k\right)}{\mathbf{W}}_k$$

10.2.4 Determining the number of components

10.2.5 Induced factorizations

再推导变分更新时，一开始仅仅考虑$$q(\mathbf{Z},\mathbf{\pi },\mathbf{\mu },\mathbf{\Lambda })=q(\mathbf{Z})q(\mathbf{\pi },\mathbf{\mu },\mathbf{\Lambda })$$但是最优的结果又推导出一些新的分解$$q(\mathbf{\pi },\mathbf{\mu },\mathbf{\Lambda })=q(\mathbf{\pi })\prod _{k=1}^{K}q\left({\mathbf{\mu }}_k,{\mathbf{\Lambda }}_k\right)$$这些额外的分解的产⽣原因是假定的分解⽅式与真实分布的条件独⽴性质相互作⽤的结果，正如图模型所告诉的那样

我们会把这些额外的分解⽅式成为诱导分解（induced factorizations），因为它们产⽣于在变分后验分布中假定的分解⽅式与真实联合概率分布的条件独⽴性质之间的相互作⽤。在变分⽅法的数值实现中，考虑这些附加的分解⽅式很重要。例如，对于⼀组变量上的⾼斯分布来说，如果分布的最优形式的精度矩阵总是对⾓矩阵（对应于关于由那个⾼斯分布独⽴描述的变量的分解⽅式），那么在计算过程中始终保留⼀个完整的精度矩阵是⼀种很低效的做法。
使⽤⼀种基于d-划分的简单的图检测⽅法，这种诱导的分解⽅式可以很容易地被检测到。我们将潜在变量划分为三个互斥的组$\mathbf{A},\mathbf{B},\mathbf{C}$，然后让我们假定我们可以在变量$\mathbf{C}$与剩余变量之间进⾏分解，即$$q(\mathbf{A},\mathbf{B},\mathbf{C})=q(\mathbf{A},\mathbf{B})q(\mathbf{C})$$则$$\begin{array}{rl}\ln q^{\star }(\mathbf{A},\mathbf{B})& ={\mathbb{E}}_{\mathbf{C}}[\ln p(\mathbf{X},\mathbf{A},\mathbf{B},\mathbf{C})]+\text{ const }\\ & ={\mathbb{E}}_{\mathbf{C}}[\ln p(\mathbf{A},\mathbf{B}|\mathbf{X},\mathbf{C})]+\text{ const }\end{array}$$那么此处我们可以考虑是否存在$q^{\star }(\mathbf{A},\mathbf{B})=q^{\star }(\mathbf{A})q^{\star }(\mathbf{B})$，这仅仅发生在$\ln p(\mathbf{A},\mathbf{B}|\mathbf{X},\mathbf{C})=\ln p(\mathbf{A}|\mathbf{X},\mathbf{C})+\ln p(\mathbf{B}|\mathbf{X},\mathbf{C})$，也就是说$$\mathbf{A}\perp \mathbf{B}|\mathbf{X},\mathbf{C}$$

10.3 Variational Linear Regression

作为变分推断的第⼆个例⼦，我们回到3.3节的贝叶斯线性回归模型中。在模型证据框架中，我们通过使⽤最⼤化似然函数的⽅法进⾏点估计，从⽽近似了在$\alpha$和$\beta$ 上的积分。⼀个纯粹的贝叶斯⽅法会对所有的超参数和参数进⾏积分。虽然精确的积分是⽆法计算的，但是我们可以使⽤变分⽅法来找到⼀个可以处理的近似。为了简化讨论，我们会假设噪声精度参数$\beta$ 已知，并且固定于它的真实值，虽然这个框架很容易扩展来包含 $\beta$ 上的概率分布。对于线性回归模型来说，可以证明变分⽅法等价于模型证据的框架。$$\begin{array}{rl}p(\mathbf{t}|\mathbf{w})& =\prod _{n=1}^N\mathcal{N}\left(t_n|{\mathbf{w}}^{\mathrm{T}}{\mathbf{\varphi }}_n,{\beta }^{-1}\right)\\ p(\mathbf{w}|\alpha )& =\mathcal{N}\left(\mathbf{w}|0,{\alpha }^{-1}\mathbf{I}\right)\end{array}$$其中${\varphi }_n=\varphi \left({\mathbf{x}}_n\right)$，取$\alpha$的共轭先验$$p(\alpha )=\mathrm{Gam}\left(\alpha |a_0,b_0\right)$$则$$p(\mathbf{t},\mathbf{w},\alpha )=p(\mathbf{t}|\mathbf{w})p(\mathbf{w}|\alpha )p(\alpha )$$

10.3.1 Variational distribution

首要目标是找到一个对后验分布$p(\mathbf{w},\alpha |\mathbf{t})$的近似，采用的变分分解为$$q(\mathbf{w},\alpha )=q(\mathbf{w})q(\alpha )$$对$\alpha$有$$\begin{array}{l}\ln q^{\star }(\alpha )=\ln p(\alpha )+{\mathbb{E}}_{\mathbf{w}}[\ln p(\mathbf{w}|\alpha )]+\text{ const }\\ =\left(a_0-1\right)\ln \alpha -b_0\alpha +\frac{M}{2}\ln \alpha -\frac{\alpha }{2}\mathbb{E}\left[{\mathbf{w}}^{\mathrm{T}}\mathbf{w}\right]+\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}.\end{array}$$刚好为$$q^{\star }(\alpha )=\mathrm{Gam}\left(\alpha |a_N,b_N\right)$$其中$$\begin{array}{rl}{a}_N& =a_0+\frac{M}{2}\\ b_N& =b_0+\frac{1}{2}\mathbb{E}\left[{\mathbf{w}}^{\mathrm{T}}\mathbf{w}\right]\end{array}$$同理有$$\begin{array}{rl}\ln q^{\star }(\mathbf{w})& =\ln p(\mathbf{t}|\mathbf{w})+{\mathbb{E}}_{\alpha }[\ln p(\mathbf{w}|\alpha )]+\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}\\ & =-\frac{\beta }{2}\sum _{n=1}^N{\left\{{\mathbf{w}}^{\mathrm{T}}{\mathbf{\varphi }}_n-t_n\right\}}^2-\frac{1}{2}\mathbb{E}[\alpha ]{\mathbf{w}}^{\mathrm{T}}\mathbf{w}+\text{ const }\\ & =-\frac{1}{2}{\mathbf{w}}^{\mathrm{T}}\left(\mathbb{E}[\alpha ]\mathbf{I}+\beta {\mathbf{\Phi }}^{\mathrm{T}}\mathbf{\Phi }\right)\mathbf{w}+\beta {\mathbf{w}}^{\mathrm{T}}{\mathbf{\Phi }}^{\mathrm{T}}\mathbf{t}+\text{ const. }\end{array}$$刚好为高斯$$q^{\star }(\mathbf{w})=\mathcal{N}\left(\mathbf{w}|{\mathbf{m}}_N,{\mathbf{S}}_N\right)$$其中$$\begin{array}{rl}{\mathbf{m}}_N& =\beta {\mathbf{S}}_N{\mathbf{\Phi }}^{\mathrm{T}}\mathbf{t}\\ {\mathbf{S}}_N& ={\left(\mathbb{E}[\alpha ]\mathbf{I}+\beta {\mathbf{\Phi }}^{\mathrm{T}}\mathbf{\Phi }\right)}^{-1}\end{array}$$注意这个结果与$\alpha$被当成固定参数时得到的后验概率分布${\mathbf{S}}_N^{-1}=\alpha \mathbf{I}+\beta {\mathbf{\Phi }}^{\mathrm{T}}\mathbf{\Phi }$相似。又因$$\begin{array}{rl}\mathbb{E}[\alpha ]& =a_N/b_N\\ \mathbb{E}\left[\mathbf{w}{\mathbf{w}}^{\mathrm{T}}\right]& ={\mathbf{m}}_N{\mathbf{m}}_N^{\mathrm{T}}+{\mathbf{S}}_N\end{array}$$这样就可以迭代地更新变分后验。
当$a_0=b_0=0$时，对应于一个无限宽的$\alpha$先验，那么$$\mathbb{E}[\alpha ]=\frac{a_N}{b_N}=\frac{M/2}{\mathbb{E}\left[{\mathbf{w}}^{\mathrm{T}}\mathbf{w}\right]/2}=\frac{M}{{\mathbf{m}}_N^{\mathrm{T}}{\mathbf{m}}_N+\mathrm{Tr}\left({\mathbf{S}}_N\right)}$$变分⽅法得到的解与使⽤EM算法最⼤化模型证据函数的⽅法得到的解完全相同，唯⼀的区别是$\alpha$的点估计被替换为了它的期望值。由于分布$q(\mathbf{w})$只通过期望$\mathbb{E}[\alpha ]$对$q(\alpha )$产⽣依赖，因此我们看到这两种⽅法对于⽆限宽的先验概率分布会给出相同的结果。

10.3.2 Predictive distribution

$$\begin{array}{rl}p(t|\mathbf{x},\mathbf{t})& =\int p(t|\mathbf{x},\mathbf{w})p(\mathbf{w}|\mathbf{t})\mathrm{d}\mathbf{w}\\ & \simeq \int p(t|\mathbf{x},\mathbf{w})q(\mathbf{w})\mathrm{d}\mathbf{w}\\ & =\int \mathcal{N}\left(t|{\mathbf{w}}^{\mathrm{T}}\mathbf{\varphi }(\mathbf{x}),{\beta }^{-1}\right)\mathcal{N}\left(\mathbf{w}|{\mathbf{m}}_N,{\mathbf{S}}_N\right)\mathrm{d}\mathbf{w}\\ & =\mathcal{N}\left(t|{\mathbf{m}}_N^{\mathrm{T}}\mathbf{\varphi }(\mathbf{x}),{\sigma }^2(\mathbf{x})\right)\end{array}$$其中$${\sigma }^2(\mathbf{x})=\frac{1}{\beta }+\mathbf{\varphi }(\mathbf{x}{)}^{\mathrm{T}}{\mathbf{S}}_N\mathbf{\varphi }(\mathbf{x})$$

10.3.3 Lower bound

另⼀个很重要的量是下界$$\begin{array}{rl}\mathcal{L}(q)=& \mathbb{E}[\ln p(\mathbf{w},\alpha ,\mathbf{t})]-\mathbb{E}[\ln q(\mathbf{w},\alpha )]\\ =& {\mathbb{E}}_{\mathbf{w}}[\ln p(\mathbf{t}|\mathbf{w})]+{\mathbb{E}}_{\mathbf{w},\alpha }[\ln p(\mathbf{w}|\alpha )]+{\mathbb{E}}_{\alpha }[\ln p(\alpha )]\\ & -{\mathbb{E}}_{\alpha }[\ln q(\mathbf{w}){]}_{\mathbf{w}}-\mathbb{E}[\ln q(\alpha )]\end{array}$$使⽤之前章节得到的结果，计算各项的值是很容易的，结果为$$\begin{array}{rl}\mathbb{E}[\ln p(\mathbf{t}|\mathbf{w}){]}_{\mathbf{w}}=& \frac{N}{2}\ln \left(\frac{\beta }{2\pi }\right)-\frac{\beta }{2}{\mathbf{t}}^{\mathrm{T}}\mathbf{t}+\beta {\mathbf{m}}_N^{\mathrm{T}}{\mathbf{\Phi }}^{\mathrm{T}}\mathbf{t}\\ & -\frac{\beta }{2}\mathrm{Tr}\left[{\mathbf{\Phi }}^{\mathrm{T}}\mathbf{\Phi }\left({\mathbf{m}}_N{\mathbf{m}}_N^{\mathrm{T}}+{\mathbf{S}}_N\right)\right]\end{array}$$$$\begin{array}{rl}\mathbb{E}[\ln p(\mathbf{w}|\alpha ){]}_{\mathbf{w},\alpha }=& -\frac{M}{2}\ln (2\pi )+\frac{M}{2}\left(\psi \left(a_N\right)-\ln b_N\right)\\ & -\frac{a_N}{2b_N}\left[{\mathbf{m}}_N^{\mathrm{T}}{\mathbf{m}}_N+\mathrm{Tr}\left({\mathbf{S}}_N\right)\right]\end{array}$$$$\begin{array}{rl}\mathbb{E}[\ln p(\alpha ){]}_{\alpha }=& a_0\ln b_0+\left(a_0-1\right)\left[\psi \left(a_N\right)-\ln b_N\right]\\ & -b_0\frac{a_N}{b_N}-\ln \Gamma \left(a_N\right)\end{array}$$$$\begin{array}{rl}-\mathbb{E}[\ln q(\mathbf{w}){]}_{\mathbf{w}}& =\frac{1}{2}\ln \left|{\mathbf{S}}_N\right|+\frac{M}{2}[1+\ln (2\pi )]\\ -\mathbb{E}[\ln q(\alpha ){]}_{\alpha }& =\ln \Gamma \left(a_N\right)-\left(a_N-1\right)\psi \left(a_N\right)-\ln b_N+a_N\end{array}$$给出了下界$\mathcal{L}(q)$与多项式模型的阶数的关系图像，数据集是从⼀个三阶多项式中⼈⼯⽣成的。这⾥，先验参数被设置为$a_0=b_0=0$，对应于⽆信息先验$p(\alpha )\propto 1/\alpha$。根据2.3.6节的讨论，它是$ln\alpha$ 上的均匀分布。正如我们在10.1节看到的那样，$\mathcal{L}$表⽰模型的对数边缘似然函数$p(M|\mathbf{t})$的下界。因此，变分框架将最⾼的概率赋予了$M=3$的模型。这与最⼤似然的结果相反。最⼤似然⽅法通过增加模型的复杂度尽可能地让误差变⼩，直到误差趋于零，这导致了最⼤似然⽅法倾向于选择具有严重过拟合现象的模型。

10.4 Exponential Family Distributions

考虑独⽴同分布数据的情形。我们将数据的值记作 X = { x n } \mathbf{X}=\left\{\mathbf{x}_{n}\right\} X={xn​}，其相应的隐变量为 Z = { z n } \mathbf{Z}=\left\{\mathbf{z}_{n}\right\} Z={zn​}，它们为指数族分布$$p(\mathbf{X},\mathbf{Z}|\mathbf{\eta })=\prod _{n=1}^{N}h\left({\mathbf{x}}_n,{\mathbf{z}}_n\right)g(\mathbf{\eta })\exp \left\{{\mathbf{\eta }}^{\mathrm{T}}\mathbf{u}\left({\mathbf{x}}_n,{\mathbf{z}}_n\right)\right\}$$引入共轭先验$$p\left(\mathbf{\eta }|{\nu }_0,{\mathbf{\chi }}_0\right)=f\left({\nu }_0,{\mathbf{\chi }}_0\right)g(\mathbf{\eta }{)}^{{\nu }_0}\exp \left\{{\nu }_o{\mathbf{\eta }}^{\mathrm{T}}{\mathbf{\chi }}_0\right\}$$考虑变分分解为$$q(\mathbf{Z},\mathbf{\eta })=q(\mathbf{Z})q(\mathbf{\eta })$$则$$\begin{array}{rl}\ln q^{\star }(\mathbf{Z})& ={\mathbb{E}}_{\eta }[\ln p(\mathbf{X},\mathbf{Z}|\mathbf{\eta })]+\text{ const }\\ & =\sum _{n=1}^N\left\{\ln h\left({\mathbf{x}}_n,{\mathbf{z}}_n\right)+\mathbb{E}\left[{\mathbf{\eta }}^{\mathrm{T}}\right]\mathbf{u}\left({\mathbf{x}}_n,{\mathbf{z}}_n\right)\right\}+\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}.\end{array}$$那么可以诱导分解为$q^{\star }(\mathbf{Z})={\prod }_n{q}^{\star }\left({\mathbf{z}}_n\right)$即$$q^{\star }\left({\mathbf{z}}_n\right)=h\left({\mathbf{x}}_n,{\mathbf{z}}_n\right)g(\mathbb{E}[\mathbf{\eta }])\exp \left\{\mathbb{E}\left[{\mathbf{\eta }}^{\mathrm{T}}\right]\mathbf{u}\left({\mathbf{x}}_n,{\mathbf{z}}_n\right)\right\}$$同样地$$\begin{array}{rl}& \ln q^{\star }(\mathbf{\eta })=\ln p\left(\mathbf{\eta }|{\nu }_0,{\mathbf{\chi }}_0\right)+{\mathbb{E}}_{\mathbf{Z}}[\ln p(\mathbf{X},\mathbf{Z}|\mathbf{\eta })]+\text{ const }\\ =& {\nu }_0\ln g(\mathbf{\eta })+{\mathbf{\eta }}^{\mathrm{T}}{\mathbf{\chi }}_0+\sum _{n=1}^N\left\{\ln g(\mathbf{\eta })+{\mathbf{\eta }}^{\mathrm{T}}{\mathbb{E}}_{{\mathbf{z}}_n}\left[\mathbf{u}\left({\mathbf{x}}_n,{\mathbf{z}}_n\right)\right]\right\}+\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}.\end{array}$$进而$$q^{\star }(\mathbf{\eta })=f\left({\nu }_N,{\mathbf{\chi }}_N\right)g(\mathbf{\eta }{)}^{{\nu }_N}\exp \left\{{\mathbf{\eta }}^{\mathrm{T}}{\mathbf{\chi }}_N\right\}$$其中$$\begin{array}{rl}{\nu }_N& ={\nu }_0+N\\ {\chi }_N& ={\chi }_0+\sum _{n=1}^N{\mathbb{E}}_{{\mathbf{z}}_n}\left[\mathbf{u}\left({\mathbf{x}}_n,{\mathbf{z}}_n\right)\right]\end{array}$$所以可以迭代地更新这两个变分分布。

10.4.1 Variational message passing

我们通过详细讨论⼀个具体的模型来说明变分⽅法的应⽤，这个模型是⾼斯模型的贝叶斯混合。这个模型可以被表⽰为有向图。这⾥我们从更⼀般的⾓度来讨论由有向图描述的模型中对变分⽅法的使⽤，推导出⼀些具有⼴泛适⽤性的结果。
对于有向图有$$p(\mathbf{x})=\prod _{i}p\left({\mathbf{x}}_i|{\mathrm{p}\mathrm{a}}_i\right)$$其中${\mathbf{x}}_i$可能为观测变量也可能为隐变量，现在考虑一个变分近似$$q(\mathbf{x})=\prod _i{q}_i\left({\mathbf{x}}_i\right)$$注意，对于观测结点，在变分分布中没有因⼦$q_i({\mathbf{x}}_i)$。那么$$\ln q_j^{\star }\left({\mathrm{x}}_j\right)={\mathbb{E}}_{i̸=j}\left[\sum _i\ln p\left({\mathrm{x}}_i|{\mathrm{p}\mathrm{a}}_i\right)\right]+\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}.$$$q_j^{\star }\left({\mathrm{x}}_j\right)$所依赖的所有结点组成的集合对应于结点${\mathbf{x}}_j$的马尔科夫毯。也就是更新某一个节点只需要进行局部的计算即可。如果我们现在确定模型的形式，其中所有的条件概率分布都有⼀个共轭-指数族的结构，那么变分推断的过程可以被转化为局部信息传递算法（Winn and Bishop, 2005）。特别地，对于⼀个特定的结点来说，⼀旦它接收到了来⾃所有的⽗结点和所有的⼦结点的信息，那么与这个结点相关联的概率分布就可以被更新。这反过来需要⼦结点从它们的同⽗结点已经接收完毕信息。下界的计算也可以得到简化，因为许多必要的值已经作为信息传递框架的⼀部分计算完毕。分布的信息传递形式有很好的缩放性质，对于⼤的⽹络很合适。

10.5. Local Variational Methods

10.1节和10.2节讨论的变分框架可以被看做“全局”⽅法。之所以这样说，是因为它直接寻找所有随机变量上的完整的后验概率分布的近似。另⼀种“局部”的⽅法涉及到寻找模型中的单独的变量或者变量组上定义的函数的界限。
对于凹函数有 f ( x ) = max ⁡ λ { λ x − g ( λ ) } g ( λ ) = max ⁡ x { λ x − f ( x ) } \begin{aligned} f(x) &=\max _{\lambda}\{\lambda x-g(\lambda)\} \\ g(\lambda) &=\max _{x}\{\lambda x-f(x)\} \end{aligned} f(x)g(λ)​=λmax​{λx−g(λ)}=xmax​{λx−f(x)}​凸函数有 f ( x ) = min ⁡ λ { λ x − g ( λ ) } g ( λ ) = min ⁡ x { λ x − f ( x ) } \begin{aligned} f(x) &=\min _{\lambda}\{\lambda x-g(\lambda)\} \\ g(\lambda) &=\min _{x}\{\lambda x-f(x)\} \end{aligned} f(x)g(λ)​=λmin​{λx−g(λ)}=xmin​{λx−f(x)}​
但对于不凸不凹的函数，只能寻求一种可逆变化，对于$$\sigma (x)=\frac{1}{1+e^{-x}}$$取log之后就称为一个凸函数，从而 g ( λ ) = min ⁡ x { λ x − f ( x ) } = − λ ln ⁡ λ − ( 1 − λ ) ln ⁡ ( 1 − λ ) g(\lambda)=\min _{x}\{\lambda x-f(x)\}=-\lambda \ln \lambda-(1-\lambda) \ln (1-\lambda) g(λ)=xmin​{λx−f(x)}=−λlnλ−(1−λ)ln(1−λ)这样就得到了一个upper bound$$\ln \sigma (x)⩽\lambda x-g(\lambda )$$则$$\sigma (x)⩽\exp (\lambda x-g(\lambda ))$$

我们还可以得到其lower bound$$\begin{array}{rl}\ln \sigma (x)& =-\ln \left(1+e^{-x}\right)=-\ln \left\{e^{-x/2}\left(e^{x/2}+e^{-x/2}\right)\right\}\\ & =x/2-\ln \left(e^{x/2}+e^{-x/2}\right)\end{array}$$其中$f(x)=-\ln \left(e^{x/2}+e^{-x/2}\right)$为$x^2$凹函数，可以通过求二阶导数证明。这样$$g(\lambda )=\underset{x^2}{\max }\left\{\lambda x^2-f\left(\sqrt{x^2}\right)\right\}$$则求导为$$0=\lambda -\frac{dx}{d{x}^2}\frac{d}{dx}f(x)=\lambda +\frac{1}{4x}\tanh \left(\frac{x}{2}\right)$$这是一个关于$x$的值$$\lambda (\xi )=-\frac{1}{4\xi }\tanh \left(\frac{\xi }{2}\right)=-\frac{1}{2\xi }\left[\sigma (\xi )-\frac{1}{2}\right]$$我们将$\xi$看作变分参数，则$$g(\lambda )=\lambda (\xi ){\xi }^2-f(\xi )=\lambda (\xi ){\xi }^2+\ln \left(e^{\xi /2}+e^{-\xi /2}\right)$$从而$$f(x)⩾\lambda x^2-g(\lambda )=\lambda x^2-\lambda {\xi }^2-\ln \left(e^{\xi /2}+e^{-\xi /2}\right)$$则$$\sigma (x)⩾\sigma (\xi )\exp \left\{(x-\xi )/2-\lambda (\xi )\left(x^2-{\xi }^2\right)\right\}$$下面来说说具体怎么使用，当计算$$I=\int \sigma (a)p(a)\mathrm{d}a$$时，$p(a)$往往为高斯分布，这样是没法求得，因此可以使用局部变分$\sigma (a)⩾f(a,\xi )$，则$$I⩾\int f(a,\xi )p(a)\mathrm{d}a=F(\xi )$$

10.6 Variational Logistic Regression

我们现在回到4.5节研究的贝叶斯logistic回归模型，说明局部变分⽅法的应⽤。在4.5节，我们将注意⼒集中于拉普拉斯近似的使⽤，⽽这⾥，我们考虑⼀种贝叶斯的⽅法，本⽅法基于Jaakkola and Jordan（2000）的⽅法。与拉普拉斯⽅法相似，这也会⽣成后验概率分布的⾼斯近似。然⽽，变分⽅法的极⼤的灵活性使得模型的准确率与拉普拉斯相⽐有所提升。此外，与拉普拉斯⽅法不同，变分⽅法最优化⼀个具有良好定义的⽬标函数，这个⽬标函数由模型证据的⼀个严格界限给定。Dybowski and Roberts（2005）也从贝叶斯的⾓度研究了logistic回归问题，使⽤了蒙特卡罗取样的技术。

10.6.1 Variational posterior distribution

首先假设先验超参数固定${\mathbf{m}}_0$ and ${\mathbf{S}}_0$。$$p(\mathbf{t})=\int p(\mathbf{t}|\mathbf{w})p(\mathbf{w})\mathrm{d}\mathbf{w}=\int \left[\prod _{n=1}^{N}p\left(t_n|\mathbf{w}\right)\right]p(\mathbf{w})\mathrm{d}\mathbf{w}$$似然为 p ( t ∣ w ) = σ ( a ) t { 1 − σ ( a ) } 1 − t = ( 1 1 + e − a ) t ( 1 − 1 1 + e − a ) 1 − t = e a t e − a 1 + e − a = e a t σ ( − a ) \begin{aligned} p(t | \mathbf{w}) &=\sigma(a)^{t}\{1-\sigma(a)\}^{1-t} \\ &=\left(\frac{1}{1+e^{-a}}\right)^{t}\left(1-\frac{1}{1+e^{-a}}\right)^{1-t} \\ &=e^{a t} \frac{e^{-a}}{1+e^{-a}}=e^{a t} \sigma(-a) \end{aligned} p(t∣w)​=σ(a)t{1−σ(a)}1−t=(1+e−a1​)t(1−1+e−a1​)1−t=eat1+e−ae−a​=eatσ(−a)​其中$a={\mathbf{w}}^{\mathrm{T}}\mathbf{\varphi }$，为了获得$p(\mathbf{t})$的lower bound，由之前的上一节推导的结果可知$$\sigma (z)⩾\sigma (\xi )\exp \left\{(z-\xi )/2-\lambda (\xi )\left(z^2-{\xi }^2\right)\right\}$$其中$$\lambda (\xi )=\frac{1}{2\xi }\left[\sigma (\xi )-\frac{1}{2}\right]$$则$$p(t|\mathbf{w})=e^{at}\sigma (-a)⩾e^{at}\sigma (\xi )\exp \left\{-(a+\xi )/2-\lambda (\xi )\left(a^2-{\xi }^2\right)\right\}$$注意，由于这个下界分别作⽤于似然函数的每⼀项，因此存在⼀个变分参数${\xi }_n$，对应于训练集的每个观测$\left({\varphi }_n,t_n\right)$。使⽤$a={\mathbf{w}}^{\mathrm{T}}\mathbf{\varphi }$，乘以先验概率分布，我们可以得到下⾯的$\mathbf{t}$和$\mathbf{w}$的联合概率分布。$$p(\mathbf{t},\mathbf{w})=p(\mathbf{t}|\mathbf{w})p(\mathbf{w})⩾h(\mathbf{w},\mathbf{\xi })p(\mathbf{w})$$其中$$\begin{array}{rl}h(\mathbf{w},\mathbf{\xi })=& \prod _{n=1}^N\sigma \left({\xi }_n\right)\exp \{{\mathbf{w}}^{\mathrm{T}}{\mathbf{\varphi }}_n{t}_n-\left({\mathbf{w}}^{\mathrm{T}}{\mathbf{\varphi }}_n+{\xi }_n\right)/2\\ & -\lambda \left({\xi }_n\right)\left({\left[{\mathbf{w}}^{\mathrm{T}}{\mathbf{\varphi }}_n\right]}^2-{\xi }_n^2\right)\}\end{array}$$取ln ln ⁡ { p ( t ∣ w ) p ( w ) } ⩾ ln ⁡ p ( w ) + ∑ n = 1 N { ln ⁡ σ ( ξ n ) + w T ϕ n t n − ( w T ϕ n + ξ n ) / 2 − λ ( ξ n ) ( [ w T ϕ n ] 2 − ξ n 2 ) } \begin{array}{c}{\ln \{p(\mathbf{t} | \mathbf{w}) p(\mathbf{w})\} \geqslant \ln p(\mathbf{w})+\sum_{n=1}^{N}\left\{\ln \sigma\left(\xi_{n}\right)+\mathbf{w}^{\mathrm{T}} \boldsymbol{\phi}_{n} t_{n}\right.} \\ {-\left(\mathbf{w}^{\mathrm{T}} \boldsymbol{\phi}_{n}+\xi_{n}\right) / 2-\lambda\left(\xi_{n}\right)\left(\left[\mathbf{w}^{\mathrm{T}} \boldsymbol{\phi}_{n}\right]^{2}-\xi_{n}^{2}\right) \}}\end{array} ln{p(t∣w)p(w)}⩾lnp(w)+∑n=1N​{lnσ(ξn​)+wTϕn​tn​−(wTϕn​+ξn​)/2−λ(ξn​)([wTϕn​]2−ξn2​)}​将先验代入$$\begin{array}{c}-\frac{1}{2}{\left(\mathbf{w}-{\mathbf{m}}_0\right)}^{\mathrm{T}}{\mathbf{S}}_0^{-1}\left(\mathbf{w}-{\mathbf{m}}_0\right)\\ +{\sum }_{n=1}^N\left\{{\mathbf{w}}^{\mathrm{T}}{\mathbf{\varphi }}_n\left(t_n-1/2\right)-\lambda \left({\xi }_n\right){\mathbf{w}}^{\mathrm{T}}\left({\mathbf{\varphi }}_n{\mathbf{\varphi }}_n^{\mathrm{T}}\right)\mathbf{w}\right\}+\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}.\end{array}$$这刚好为$\mathbf{w}$的二次函数，则说明后验为一个高斯$$q(\mathbf{w})=\mathcal{N}\left(\mathbf{w}|{\mathbf{m}}_N,{\mathbf{S}}_N\right)$$其中$$\begin{array}{rl}{\mathbf{m}}_N& ={\mathbf{S}}_N\left({\mathbf{S}}_0^{-1}{\mathbf{m}}_0+\sum _{n=1}^N\left(t_n-1/2\right){\mathbf{\varphi }}_n\right)\\ {\mathbf{S}}_N^{-1}& ={\mathbf{S}}_0^{-1}+2\sum _{n=1}^N\lambda \left({\xi }_n\right){\mathbf{\varphi }}_n{\mathbf{\varphi }}_n^{\mathrm{T}}\end{array}$$与Laplace framework一样，得到了一个近似后验的高斯分布。然后，变分参数 { ξ n } \left\{\xi_{n}\right\} {ξn​}提供的额外的灵活性使得这个近似的精度更⾼

10.6.2 Optimizing the variational parameters

我们现在得到了后验概率分布的⼀个归⼀化的⾼斯近似。我们稍后会使⽤这个近似得到对于新数据的预测分布。然⽽，⾸先我们需要通过最⼤化边缘似然函数的下界，确定变分参数 { ξ n } \left\{\xi_{n}\right\} {ξn​}。$$\ln p(\mathbf{t})=\ln \int p(\mathbf{t}|\mathbf{w})p(\mathbf{w})\mathrm{d}\mathbf{w}⩾\ln \int h(\mathbf{w},\mathbf{\xi })p(\mathbf{w})\mathrm{d}\mathbf{w}=\mathcal{L}(\mathbf{\xi })$$与3.5节的线性回归模型的超参数$\alpha$的最优化⼀样，有两种⽅法确定 { ξ n } \left\{\xi_{n}\right\} {ξn​}。在第⼀种⽅法中，我们看到函数$\mathcal{L}(\mathbf{\xi })$由$\mathbf{w}$上的积分定义，因此我们可以将$\mathbf{w}$看成⼀个潜在变量，然后使⽤EM算法。在第⼆种⽅法中，我们解析地对$\mathbf{w}$积分，然后直接关于$\xi$进⾏最⼤化。让我们⾸先考虑EM⽅法。
首先初始化 { ξ n } \left\{\xi_{n}\right\} {ξn​}，记作${\mathbf{\xi }}^{\mathrm{o}\mathrm{l}\mathrm{d}}$，在E-step找到关于$\mathbf{w}$的后验$$q(\mathbf{w})=\mathcal{N}\left(\mathbf{w}|{\mathbf{m}}_N,{\mathbf{S}}_N\right)$$在M-step$$Q\left(\mathbf{\xi },{\mathbf{\xi }}^{\mathrm{o}\mathrm{l}\mathrm{d}}\right)=\mathbb{E}[\ln h(\mathbf{w},\mathbf{\xi })p(\mathbf{w})]$$$$Q\left(\mathbf{\xi },{\mathbf{\xi }}^{\text{ old }}\right)=\sum _{n=1}^N\left\{\ln \sigma \left({\xi }_n\right)-{\xi }_n/2-\lambda \left({\xi }_n\right)\left({\mathbf{\varphi }}_n^{\mathrm{T}}\mathbb{E}\left[\mathbf{w}{\mathbf{W}}^{\mathrm{T}}\right]{\mathbf{\varphi }}_n-{\xi }_n^2\right)\right\}+\text{ const }$$求导为0$$0={\lambda }^{\prime }\left({\xi }_n\right)\left({\varphi }_n^{\mathrm{T}}\mathbb{E}\left[\mathbf{w}{\mathbf{w}}^{\mathrm{T}}\right]{\mathbf{\varphi }}_n-{\xi }_n^2\right)$$可以得到$${\left({\xi }_n^{\text{ new }}\right)}^2={\varphi }_n^{\mathrm{T}}\mathbb{E}\left[\mathbf{w}{\mathbf{w}}^{\mathrm{T}}\right]{\mathbf{\varphi }}_n={\mathbf{\varphi }}_n^{\mathrm{T}}\left({\mathbf{S}}_N+{\mathbf{m}}_N{\mathbf{m}}_N^{\mathrm{T}}\right){\mathbf{\varphi }}_n$$这样就可以迭代地更新变分参数。注意到这个lower bound是可以关于$\mathbf{w}$积分的，然后直接对$\xi$求导也能得到和EM一样的结果。$$\begin{array}{rl}\mathcal{L}(\mathbf{\xi })=& \frac{1}{2}\ln \frac{\left|{\mathbf{S}}_N\right|}{\left|{\mathbf{S}}_0\right|}-\frac{1}{2}{\mathbf{m}}_N^{\mathrm{T}}{\mathbf{S}}_N^{-1}{\mathbf{m}}_N+\frac{1}{2}{\mathbf{m}}_0^{\mathrm{T}}{\mathbf{S}}_0^{-1}{\mathbf{m}}_0\\ & +\sum _{n=1}^N\left\{\ln \sigma \left({\xi }_n\right)-\frac{1}{2}{\xi }_n-\lambda \left({\xi }_n\right){\xi }_n^2\right\}\end{array}$$下图为一个具体例子

10.6.3 Inference of hyperparameters

⽬前为⽌，我们将先验概率分布的超参数$\alpha$看成⼀个已知参数。我们现在将贝叶斯logistic回归模型进⾏推⼴，使得这个参数的值可以从数据集中推断出来。这可以通过将全局变分近似和局部变分近似结合到⼀个框架中的⽅式完成，从⽽在每个阶段都保留边缘似然函数的下界
特别地，我们再次考虑⼀个简单的各向同性的⾼斯先验概率分布，形式为$$p(\mathbf{w}|\alpha )=\mathcal{N}\left(\mathbf{w}|0,{\alpha }^{-1}\mathbf{I}\right)$$引入超先验$$p(\alpha )=\mathrm{Gam}\left(\alpha |a_0,b_0\right)$$则边缘似然为$$p(\mathbf{t})=\iint p(\mathbf{w},\alpha ,\mathbf{t})\mathrm{d}\mathbf{w}\mathrm{d}\alpha$$$$p(\mathbf{w},\alpha ,\mathbf{t})=p(\mathbf{t}|\mathbf{w})p(\mathbf{w}|\alpha )p(\alpha )$$然而上式是没法解析处理的，我们会在同⼀个模型中使⽤全局的变分⽅法和局部的变分⽅法来解决这个问题。
⾸先，我们引⼊⼀个变分分布$q(\mathbf{w},\alpha )$，则有$$\ln p(\mathbf{t})=\mathcal{L}(q)+\mathrm{K}\mathrm{L}(q\Vert p)$$其中$$\begin{array}{rl}\mathcal{L}(q)& =\iint q(\mathbf{w},\alpha )\ln \left\{\frac{p(\mathbf{w},\alpha ,\mathbf{t})}{q(\mathbf{w},\alpha )}\right\}\mathrm{d}\mathbf{w}\mathrm{d}\alpha \\ \mathrm{K}\mathrm{L}(q\Vert p)& =-\iint q(\mathbf{w},\alpha )\ln \left\{\frac{p(\mathbf{w},\alpha |\mathbf{t}))}{q(\mathbf{w},\alpha )}\right\}\mathrm{d}\mathbf{w}\mathrm{d}\alpha \end{array}$$然而这个lower bound仍然不可处理，由于$p(\mathbf{t}|\mathbf{w})$。那么可以考虑使用一个局部变分$$\begin{array}{rl}\ln p(\mathbf{t})& ⩾\mathcal{L}(q)⩾\tilde{\mathcal{L}}(q,\mathbf{\xi })\\ & =\iint q(\mathbf{w},\alpha )\ln \left\{\frac{h(\mathbf{w},\mathbf{\xi })p(\mathbf{w}|\alpha )p(\alpha )}{q(\mathbf{w},\alpha )}\right\}\mathrm{d}\mathbf{w}\mathrm{d}\alpha \end{array}$$接下来我们假设变分分布可以在参数和超参数之间进⾏分解，即$$q(\mathbf{w},\alpha )=q(\mathbf{w})q(\alpha )$$按照变分 ln ⁡ q ( w ) = E α [ ln ⁡ { h ( w , ξ ) p ( w ∣ α ) p ( α ) } ] + c o n s t = ln ⁡ h ( w , ξ ) + E α [ ln ⁡ p ( w ∣ α ) ] + c o n s t . \begin{aligned} \ln q(\mathbf{w}) &=\mathbb{E}_{\alpha}[\ln \{h(\mathbf{w}, \boldsymbol{\xi}) p(\mathbf{w} | \alpha) p(\alpha)\}]+\mathrm{const} \\ &=\ln h(\mathbf{w}, \boldsymbol{\xi})+\mathbb{E}_{\alpha}[\ln p(\mathbf{w} | \alpha)]+\mathrm{const.} \end{aligned} lnq(w)​=Eα​[ln{h(w,ξ)p(w∣α)p(α)}]+const=lnh(w,ξ)+Eα​[lnp(w∣α)]+const.​代入后$$\ln q(\mathbf{w})=-\frac{\mathbb{E}[\alpha ]}{2}{\mathbf{w}}^{\mathrm{T}}\mathbf{w}+\sum _{n=1}^N\left\{\left(t_n-1/2\right){\mathbf{w}}^{\mathrm{T}}{\mathbf{\varphi }}_n-\lambda \left({\xi }_n\right){\mathbf{w}}^{\mathrm{T}}{\mathbf{\varphi }}_n{\mathbf{\varphi }}_n^{\mathrm{T}}\mathbf{w}\right\}+\text{ const. }$$可以发现其为高斯$$q(\mathbf{w})=\mathcal{N}\left(\mathbf{w}|{\mathbf{\mu }}_N,{\mathbf{\Sigma }}_N\right)$$其中$$\begin{array}{rl}{\mathbf{\Sigma }}_N^{-1}{\mathbf{\mu }}_N& =\sum _{n=1}^N\left(t_n-1/2\right){\mathbf{\varphi }}_n\\ {\mathbf{\Sigma }}_N^{-1}& =\mathbb{E}[\alpha ]\mathbf{I}+2\sum _{n=1}^N\lambda \left({\xi }_n\right){\mathbf{\varphi }}_n{\mathbf{\varphi }}_n^{\mathrm{T}}\end{array}$$同样的$$\ln q(\alpha )={\mathbb{E}}_{\mathbf{w}}[\ln p(\mathbf{w}|\alpha )]+\ln p(\alpha )+\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}.$$$$\ln q(\alpha )=\frac{M}{2}\ln \alpha -\frac{\alpha }{2}\mathbb{E}\left[{\mathbf{w}}^{\mathrm{T}}\mathbf{w}\right]+\left(a_0-1\right)\ln \alpha -b_0\alpha +\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}.$$$$q(\alpha )=\mathrm{Gam}\left(\alpha |a_N,b_N\right)=\frac{1}{\Gamma \left(a_0\right)}{a}_0^{b_0}{\alpha }^{a_0-1}{e}^{-b_0\alpha }$$其中$$\begin{array}{rl}{a}_N& =a_0+\frac{M}{2}\\ b_N& =b_0+\frac{1}{2}{\mathbb{E}}_{\mathbf{w}}\left[{\mathbf{w}}^{\mathrm{T}}\mathbf{w}\right]\end{array}$$我们还需要最优化变分参数${\xi }_n$，这也可以通过最⼤化下界$\tilde{\mathcal{L}}(q,\mathbf{\xi })$的⽅式得到。略去与$\xi$⽆关的项，对$\alpha$积分，我们有$$\tilde{\mathcal{L}}(q,\mathbf{\xi })=\int q(\mathbf{w})\ln h(\mathbf{w},\mathbf{\xi })\mathrm{d}\mathbf{w}+\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}.$$这个与之前推导的结果一样，则$${\left({\xi }_n^{\text{ new }}\right)}^2={\varphi }_n^{\mathrm{T}}\left({\mathbf{\Sigma }}_N+{\mathbf{\mu }}_N{\mathbf{\mu }}_N^{\mathrm{T}}\right){\mathbf{\varphi }}_n$$这样就能迭代起来了。

10.7 Expectation Propagation

在本章的最后⼀节，我们讨论确定性近似推断的另⼀种形式，被称为期望传播（expectation propagation），或者EP（Minka, 2001a; Minka, 2001b）。与⽬前为⽌讨论的变分贝叶斯⽅法相同，这种⽅法也基于对Kullback-Leibler散度的最⼩化，但是现在形式相反，从⽽得到了性质相当不同的近似结果。
考虑最小化$\mathrm{K}\mathrm{L}(p\Vert q)$，其中$p(\mathbf{z})$固定，且$q(\mathbf{z})$为指数族分布。$$q(\mathbf{z})=h(\mathbf{z})g(\mathbf{\eta })\exp \left\{{\mathbf{\eta }}^{\mathrm{T}}\mathbf{u}(\mathbf{z})\right\}$$则目标函数为$$\mathrm{K}\mathrm{L}(p\Vert q)=-\ln g(\mathbf{\eta })-{\mathbf{\eta }}^{\mathrm{T}}{\mathbb{E}}_{p(\mathbf{z})}[\mathbf{u}(\mathbf{z})]+\mathrm{c}\mathrm{o}\mathrm{n}\mathrm{s}\mathrm{t}$$对自然参数求导$$-\nabla \ln g(\mathbf{\eta })={\mathbb{E}}_{p(\mathbf{z})}[\mathbf{u}(\mathbf{z})]$$之前在第二章有$$-\nabla \ln g(\mathbf{\eta })={\mathbb{E}}_{q(\mathbf{z})}[\mathbf{u}(\mathbf{z})]$$则有$${\mathbb{E}}_{q(\mathbf{z})}[\mathbf{u}(\mathbf{z})]={\mathbb{E}}_{p(\mathbf{z})}[\mathbf{u}(\mathbf{z})]$$我们看到，最优解仅仅对应于将充分统计量的期望进⾏匹配。因此，例如，如果$q(\mathbf{z})$是⼀个⾼斯分布$\mathcal{N}(\mathbf{z}|\mathbf{\mu },\mathbf{\Sigma })$，那么我们通过令$q(\mathbf{z})$的均值$\mathbf{\mu }$等于分布$p(\mathbf{z})$的均值并且令协⽅差$\mathbf{\Sigma }$等于p(z)的协⽅差，即可最⼩化Kullback-Leibler散度。这有时被称为矩匹配（moment matching）。
下面考虑模型为$$p(\mathcal{D},\mathbf{\theta })=\prod _i{f}_i(\mathbf{\theta })$$我们往往需要后验分布$p(\mathbf{\theta }|\mathcal{D})$，则$$p(\mathbf{\theta }|\mathcal{D})=\frac{1}{p(\mathcal{D})}\prod _i{f}_i(\mathbf{\theta })$$其中model evidence为$$p(\mathcal{D})=\int \prod _i{f}_i(\mathbf{\theta })\mathrm{d}\mathbf{\theta }$$但往往这是不可处理的，因此得考虑近似。期望传播基于后验概率分布的近似，这个近似也由⼀组因⼦的乘积给出，即$$q(\mathbf{\theta })=\frac{1}{Z}\prod _i{\tilde{f}}_i(\mathbf{\theta })$$其中每个因子${\tilde{f}}_i(\mathbf{\theta })$近似真实后验$f_i(\mathbf{\theta })$，我们对${\tilde{f}}_i(\mathbf{\theta })$加入一定的限制，我们假设它们为指数族分布，这样一来，$q(\mathbf{\theta })$也是一个指数族分布。
理想情况下，我们想通过最小化KL散度求得变分近似。$$\mathrm{K}\mathrm{L}(p\Vert q)=\mathrm{K}\mathrm{L}\left(\frac{1}{p(\mathcal{D})}\prod _i{f}_i(\mathbf{\theta })\Vert \frac{1}{Z}\prod _i{\tilde{f}}_i(\mathbf{\theta })\right)$$注意与变分推断中使⽤的KL散度相⽐，这个KL的散度恰好相反。通常这个最⼩化是⽆法进⾏的，因为KL散度涉及到关于真实概率分布求平均。作为一个初略的近似，我们最小化$f_i(\mathbf{\theta })$和${\tilde{f}}_i(\mathbf{\theta })$的KL，这个问题容易得多，并且具有算法⽆需迭代的优点。然⽽，由于每个因⼦被各⾃独⽴地进⾏近似，因此因⼦的乘积的近似效果可能很差。
期望传播通过在所有剩余因⼦的环境中对每个因⼦进⾏优化，从⽽取得了⼀个效果好得多的近似。⾸先，这种⽅法初始化因⼦${\tilde{f}}_i(\mathbf{\theta })$，然后在因⼦之间进⾏循环，每次优化⼀个因⼦。这种⽅法的思想类似于之前讨论的变分贝叶斯框架的因⼦更近过程。假设我们希望优化因⼦${\tilde{f}}_j(\mathbf{\theta })$。⾸先，我们将这个因⼦从乘积中移除，得到${\prod }_{i̸=j}{\tilde{f}}_i(\mathbf{\theta })$。从概念上讲，我们要确定因⼦${\tilde{f}}_j(\mathbf{\theta })$的⼀个修正形式，使得乘积$$q^{\text{ new }}(\mathbf{\theta })\propto {\tilde{f}}_j(\mathbf{\theta })\prod _{i̸=j}{\tilde{f}}_i(\mathbf{\theta })$$尽可能接近$$f_j(\mathbf{\theta })\prod _{i̸=j}{\tilde{f}}_i(\mathbf{\theta })$$其中我们保持所有$i̸=j$的因⼦${\tilde{f}}_i(\mathbf{\theta })$固定。这保证了近似在由剩余的因⼦定义的后验概率较⾼的区域最精确。首先定义$$q^{\backslash j}(\mathbf{\theta })=\frac{q(\mathbf{\theta })}{{\tilde{f}}_j(\mathbf{\theta })}$$则$$\frac{1}{Z_j}{f}_j(\mathbf{\theta })q^{\backslash j}(\mathbf{\theta })$$其中$$Z_j=\int f_j(\mathbf{\theta })q^{\backslash j}(\mathbf{\theta })\mathrm{d}\mathbf{\theta }$$然后优化$$\mathrm{K}\mathrm{L}\left(\frac{f_j(\mathbf{\theta })q^{\backslash j}(\mathbf{\theta })}{Z_j}\Vert q^{\mathrm{n}\mathrm{e}\mathrm{w}}(\mathbf{\theta })\right)$$由于$q^{\text{ new }}(\mathbf{\theta })$为一个指数族分布，那么很容易求解上式，利用之前的结果$${\mathbb{E}}_{q(\mathbf{z})}[\mathbf{u}(\mathbf{z})]={\mathbb{E}}_{p(\mathbf{z})}[\mathbf{u}(\mathbf{z})]$$设置为充分统计量的期望即可。我们会假设这是⼀个可以计算的操作。例如，假如$q(\mathbf{\theta })$为高斯分布$\mathcal{N}(\mathbf{\theta }|\mathbf{\mu },\mathbf{\Sigma })$，那么可以设置均值为$\frac{f_j(\mathbf{\theta })q^{\backslash j}(\mathbf{\theta })}{Z_j}$的均值，方差为$\frac{f_j(\mathbf{\theta })q^{\backslash j}(\mathbf{\theta })}{Z_j}$的方差。得到了$q^{\mathrm{n}\mathrm{e}\mathrm{w}}(\mathbf{\theta })$了，则可以得到需要更新的$j$$${\tilde{f}}_j(\mathbf{\theta })=K\frac{q^{\mathrm{n}\mathrm{e}\mathrm{w}}(\mathbf{\theta })}{q^{\backslash j}(\mathbf{\theta })}$$其中$$K=\int {\tilde{f}}_j(\mathbf{\theta })q^{\backslash j}(\mathbf{\theta })\mathrm{d}\mathbf{\theta }$$$K$的值可以通过匹配零阶矩的⽅式得到$$\int {\tilde{f}}_j(\mathbf{\theta })q^{\backslash j}(\mathbf{\theta })\mathrm{d}\mathbf{\theta }=\int f_j(\mathbf{\theta })q^{\backslash j}(\mathbf{\theta })\mathrm{d}\mathbf{\theta }$$则$$K=Z_j$$

具体流程如下

10.7.1 Example: The clutter problem

我们的目标是从观测数据集中推断出多元高斯的均值$\mathbf{\theta }$，其中数据中有背景噪声。

建模为GMM$$p(\mathbf{x}|\mathbf{\theta })=(1-w)\mathcal{N}(\mathbf{x}|\mathbf{\theta },\mathbf{I})+w\mathcal{N}(\mathbf{x}|0,a\mathbf{I})$$其中$w$为一个已经的。设高斯先验$$p(\mathbf{\theta })=\mathcal{N}(\mathbf{\theta }|0,b\mathbf{I})$$其中$a=10,b=100$ and $w=0.5$，则联合分布为$$p(\mathcal{D},\mathbf{\theta })=p(\mathbf{\theta })\prod _{n=1}^{N}p\left({\mathbf{x}}_n|\mathbf{\theta }\right)$$因此后验概率分布由$2^N$个⾼斯分布混合⽽成。从⽽精确解决这个问题的计算代价会随着数据集的规模指数增长，因此对于⼤的$N$值，精确求解是不可⾏的。为了运用EP算法，有$f_0(\mathbf{\theta })=p(\mathbf{\theta })$，$f_n(\mathbf{\theta })=p\left({\mathbf{x}}_n|\mathbf{\theta }\right)$，接下来我们要选择一个指数族分布作为近似分布，为了简单起见，我们选择高斯$$q(\mathbf{\theta })=\mathcal{N}(\mathbf{\theta }|\mathbf{m},v\mathbf{I})$$则可分解为指数二次函数，注意以下并不是一个高斯分布，只是形式上是$${\tilde{f}}_n(\mathbf{\theta })=s_n\mathcal{N}\left(\mathbf{\theta }|{\mathbf{m}}_n,v_n\mathbf{I}\right)$$其中$n=1,\dots ,N$，我们设置${\tilde{f}}_0(\mathbf{\theta })$ equal to the prior $p(\mathbf{\theta })$。${\tilde{f}}_n(\mathbf{\theta }),$ for $n=1,\dots ,N$可以被初始化为1，则对应$s_n={\left(2\pi v_n\right)}^{D/2},v_n\to \infty$和${\mathbf{m}}_n=0$，其中$D$为维度。则此时$q(\mathbf{\theta })$为先验。这样就能开始迭代了，首先从$q(\mathbf{\theta })$除去${\tilde{f}}_n(\mathbf{\theta })$得到$q^{\backslash n}(\mathbf{\theta })$，其中均值和逆方差为$$\begin{array}{rl}{\mathbf{m}}^{\backslash n}& =\mathbf{m}+v^{\backslash n}{v}_n^{-1}\left(\mathbf{m}-{\mathbf{m}}_n\right)\\ {\left(v^{\backslash n}\right)}^{-1}& =v^{-1}-v_n^{-1}\end{array}$$之后计算归一化参数$$Z_n=(1-w)\mathcal{N}\left({\mathbf{x}}_n|{\mathbf{m}}^{\backslash n},\left(v^{\backslash n}+1\right)\mathbf{I}\right)+w\mathcal{N}\left({\mathbf{x}}_n|0,a\mathbf{I}\right)$$然后利用$q^{\backslash n}(\mathbf{\theta })f_n(\mathbf{\theta })$的均值和方差得到$q^{\text{ new }}(\mathbf{\theta })$的均值和方差$$\begin{array}{rl}\mathbf{m}& ={\mathbf{m}}^{\backslash n}+{\rho }_n\frac{v^{\backslash n}}{v^{\backslash n}+1}\left({\mathbf{x}}_n-{\mathbf{m}}^{\backslash n}\right)\\ v& =v^{\backslash n}-{\rho }_n\frac{{\left(v^{\backslash n}\right)}^2}{v^{\backslash n}+1}+{\rho }_n\left(1-{\rho }_n\right)\frac{{\left(v^{\backslash n}\right)}^2{\Vert {\mathbf{x}}_n-{\mathbf{m}}^{\backslash n}\Vert }^2}{D{\left(v^{\backslash n}+1\right)}^2}\end{array}$$其中$${\rho }_n=1-\frac{w}{Z_n}\mathcal{N}\left({\mathbf{x}}_n|0,a\mathbf{I}\right)$$它可以简单地表⽰为点${\mathbf{x}}_n$不在聚类中的概率。接下来就是计算更新后的${\tilde{f}}_n(\theta )$，其参数为$$\begin{array}{rl}{v}_n^{-1}& ={\left(v^{\text{ new }}\right)}^{-1}-{\left(v^{\backslash n}\right)}^{-1}\\ {\mathbf{m}}_n& ={\mathbf{m}}^{\backslash n}+\left(v_n+v^{\backslash n}\right){\left(v^{\backslash n}\right)}^{-1}\left({\mathbf{m}}^{\text{ new }}-{\mathbf{m}}^{\backslash n}\right)\\ s_n& =\frac{Z_n}{{\left(2\pi v_n\right)}^{D/2}\mathcal{N}\left({\mathbf{m}}_n|{\mathbf{m}}^{\backslash n},\left(v_n+v^{\backslash n}\right)\mathbf{I}\right)}\end{array}$$优化过程不断重复，直到满⾜⼀个合适的终⽌准则，例如在对所有因⼦进⾏的⼀次优化迭代中，参数值的最⼤改变量⼩于⼀个阈值。利用此可得model evidence的近似$$p(\mathcal{D})\simeq {\left(2\pi v^{\text{ new }}\right)}^{D/2}\exp (B/2)\prod _{n=1}^N\left\{s_n{\left(2\pi v_n\right)}^{-D/2}\right\}$$其中$$B=\frac{{\left({\mathbf{m}}^{\text{ new }}\right)}^{\mathrm{T}}{\mathbf{m}}^{\text{ new }}}{v}-\sum _{n=1}^N\frac{{\mathbf{m}}_n^{\mathrm{T}}{\mathbf{m}}_n}{v_n}$$注意，因⼦近似可以有⽆穷⼤的或者负数的“⽅差”参数$v_n$。这仅仅对应于曲线向上弯曲⽽不是向下弯曲的情形，并且只要所有的近似后验概率$q(\mathbf{\theta })$有正的⽅差，这种情形就未必有问题。平均场变分

10.7.2 Expectation propagation on graphs

⽬前为⽌在我们对于EP的⼀般的讨论中，我们让概率分布$p(\mathbf{\theta })$中的所有因⼦$f_i(\mathbf{\theta })$是$\mathbf{\theta }$的全部分量的函数，类似地，对于近似分布$q(\mathbf{\theta })$的近似因⼦$\tilde{f}(\mathbf{\theta })$的情形也相同。我们现在考虑下⾯的情形：因⼦只依赖于变量的⼀个⼦集。这种限制可以很⽅便地使⽤第8章讨论的概率图模型的框架来表⽰。这⾥，我们使⽤因⼦图表⽰，因为它同时包含了有向图和⽆向图。
之前对于分解后的$q$，在优化$\mathrm{K}\mathrm{L}(p\Vert q)$时，最优的$q$就是$p$的边缘分布。考虑以下的因子图

$$p(\mathbf{x})=f_a\left(x_1,x_2\right)f_b\left(x_2,x_3\right)f_c\left(x_2,x_4\right)$$目标是近似$$q(\mathbf{x})\propto {\tilde{f}}_a\left(x_1,x_2\right){\tilde{f}}_b\left(x_2,x_3\right){\tilde{f}}_c\left(x_2,x_4\right)$$现在，假设我们将注意⼒集中于近似分布上，其中因⼦本⾝可以关于各个变量进⾏分解，即$$q(\mathbf{x})\propto {\tilde{f}}_{a1}\left(x_1\right){\tilde{f}}_{a2}\left(x_2\right){\tilde{f}}_{b2}\left(x_2\right){\tilde{f}}_{b3}\left(x_3\right){\tilde{f}}_{c2}\left(x_2\right){\tilde{f}}_{c4}\left(x_4\right)$$以上为一个完全分解的分布。现在我们使用EP算法。假设我们已经初始化了所有的因⼦，并且我们选择优化因⼦${\tilde{f}}_b\left(x_2,x_3\right)={\tilde{f}}_{b2}\left(x_2\right){\tilde{f}}_{b3}\left(x_3\right)$。⾸先，我们将这个因⼦从近似分布中移除，
得到$$q^{\backslash b}(\mathbf{x})={\tilde{f}}_{a1}\left(x_1\right){\tilde{f}}_{a2}\left(x_2\right){\tilde{f}}_{c2}\left(x_2\right){\tilde{f}}_{c4}\left(x_4\right)$$之后乘上精确的因子$f_b\left(x_2,x_3\right)$，则$$\hat{p}(\mathbf{x})=q^{\backslash b}(\mathbf{x})f_b\left(x_2,x_3\right)={\tilde{f}}_{a1}\left(x_1\right){\tilde{f}}_{a2}\left(x_2\right){\tilde{f}}_{c2}\left(x_2\right){\tilde{f}}_{c4}\left(x_4\right)f_b\left(x_2,x_3\right)$$然后通过$\mathrm{K}\mathrm{L}\left(\hat{p}\Vert q^{\text{ new }}\right)$优化$q^{\text{ new }}(\mathbf{x})$，之前说过最优的就是边缘分布，这四个边缘分布为$$\begin{array}{ll}\hat{p}\left(x_1\right)& \propto {\tilde{f}}_{a1}\left(x_1\right)\\ \hat{p}\left(x_2\right)& \propto {\tilde{f}}_{a2}\left(x_2\right){\tilde{f}}_{c2}\left(x_2\right){\sum }_{x_3}{f}_b\left(x_2,x_3\right)\\ \hat{p}\left(x_3\right)& \propto {\sum }_{x_2}\left\{f_b\left(x_2,x_3\right){\tilde{f}}_{a2}\left(x_2\right){\tilde{f}}_{c2}\left(x_2\right)\right\}\\ \hat{p}\left(x_4\right)& \propto {\tilde{f}}_{c4}\left(x_4\right)\end{array}$$从而$q^{\text{ new }}(\mathbf{x})$将以上四个相乘得到。为了得到${\tilde{f}}_b\left(x_2,x_3\right)={\tilde{f}}_{b2}\left(x_2\right){\tilde{f}}_{b3}\left(x_3\right)$的更新，首先$\mathrm{divide}{q}^{\mathrm{n}\mathrm{e}\mathrm{w}}(\mathbf{x})$ by $q^{\backslash b}(\mathbf{x})$，则$$\begin{array}{ll}{\tilde{f}}_{b2}\left(x_2\right)& \propto {\sum }_{x_3}{f}_b\left(x_2,x_3\right)\\ {\tilde{f}}_{b3}\left(x_3\right)& \propto {\sum }_{x_2}\left\{f_b\left(x_2,x_3\right){\tilde{f}}_{a2}\left(x_2\right){\tilde{f}}_{c2}\left(x_2\right)\right\}\end{array}$$这些与使⽤置信传播得到的信息完全相同，其中从变量结点到因⼦结点的信息已经被整合到从因⼦结点到变量结点的信息当中。以上分别是${\tilde{f}}_{b2}\left(x_2\right)$对应${\mu }_{f_b\to x_2}\left(x_2\right)$，${\tilde{f}}_{b3}\left(x_3\right)$对应${\mu }_{f_b\to x_3}\left(x_3\right)$这个结果与标准的置信传播稍微有些不同，因为信息同时向两个⽅向传递。我们可以很容易地修改EP步骤，给出加和-乘积算法的标准形式。如果我们每次只优化⼀项，那么我们可以选择我们所希望进⾏的优化的顺序。特别地，对于⼀个树结构的图，我们可以遵循两遍更新的框架，对应于标准的置信传播⽅法，它会产⽣对变量和因⼦的边缘概率分布的精确的推断。这种情况下，近似因⼦的初始化不再重要。
现在，让我们考虑⼀个⼀般的因⼦图，它对应于下⾯的概率分布$$p(\mathbf{\theta })=\prod _i{f}_i\left({\mathbf{\theta }}_i\right)$$其中${\mathbf{\theta }}_i$代表与因子$f_i$相关的变量子集合。考虑使用一个全分解的分布近似$$q(\mathbf{\theta })\propto \prod _i\prod _k{\tilde{f}}_{ik}\left({\theta }_k\right)$$${\theta }_k$代表一个变量，考虑优化${\tilde{f}}_{jl}\left({\theta }_l\right)$，则从$q(\mathbf{\theta })$除去${\tilde{f}}_j\left({\mathbf{\theta }}_j\right)$，则$$q^{\backslash j}(\mathbf{\theta })\propto \prod _{i̸=j}\prod _k{\tilde{f}}_{ik}\left({\theta }_k\right)$$乘上精确因子$f_j\left({\mathbf{\theta }}_j\right)$，则$$q^{\backslash j}(\mathbf{\theta })f_j\left({\mathbf{\theta }}_j\right)$$的在${\tilde{f}}_{jl}\left({\theta }_l\right)$上的边缘分布。则$${\tilde{f}}_{jl}\left({\theta }_l\right)\propto \frac{{\sum }_{{\theta }_{m̸=l}\in {\mathbf{\theta }}_j}{f}_j\left({\mathbf{\theta }}_j\right){\prod }_{k̸=j}{\prod }_{m̸=l}{\tilde{f}}_{km}\left({\theta }_m\right)}{q^{\backslash j}(\mathbf{\theta })}$$