第10章 近似推断

10 近似推断

在概率模型的应用中，一个中心任务是在给定观测（可见）数据变量X的条件下，计算潜在变量Z的后验概率分布$p(Z|X)$，以及计算这个概率分布的期望。模型可能是一个纯粹的贝叶斯模型，其中任何未知的参数都有一个先验概率分布，并且被整合到了潜在变量集合中，记作向量Z。

在实际应用中，计算后验概率分布或者计算关于这个后验概率分布的期望是不可行的，这可能是由于潜在空间的维度太高，以至于无法计算，或者由于后验概率分布的形式特别复杂，从而期望无法被解析地计算。

在连续变量的情形中，需要求解的积分可能没有解析解，而空间的维度和被积函数的复杂度可能使得数值积分变得不可行。对于离散变量，求边缘概率的过程涉及对隐含变量的所有可能的配置进行求和。这个过程原则上总是可以计算的，但实际应用中，隐含状态的数量可能有指数多个，从而精确的计算所需代价过高。

在这种情况下，需要借助近似方法。根据近似方法依赖于随机近似还是确定近似，近似方法分为两类：

1. 随机方法：例如马尔可夫链蒙特卡罗方法。这些方法通常具有这样的性质：给定无限多的计算资源，他们可以生成精确的结果，实际应用中，取样方法需要的计算量会相当法，经常将这些方法的应用限制在小规模的问题中。
2. 确定性近似方法：其中有些方法对于大规模数据很实用，这些方法基于对后验概率分布的解析近似。

10.1 变分推断

将所有观测变量的集合记作X，将所有潜在变量和参数组成的集合记作Z。假设有N个独立同分布的数据，其中 X = { x 1 , x 2 , ⋯   , x N }   a n d   Z = { z 1 , ⋯   , z N } X=\{x_1,x_2,\cdots,x_N\} \ and \ Z=\{z_1,\cdots,z_N\} X={x1​,x2​,⋯,xN​} and Z={z1​,⋯,zN​}。概率模型确定了联合概率分布$p(X,Z)$，目标是找到对后验概率分布$p(Z|X)$以及模型证据$p(X)$的近似。与EM算法的讨论相同，将对数边缘概率分解，即
$$lnp(X)=L(q)+KL(q||p)$$
其中相关等式定义如下：
L ( q ) = ∫ Z q ( Z ) l n { p ( Z , X ) q ( Z ) } d Z L(q)=\int_Z{q(Z)ln\{\frac{p(Z,X)}{q(Z)}\}dZ} L(q)=∫Z​q(Z)ln{q(Z)p(Z,X)​}dZ

K L ( q ∣ ∣ p ) = − ∫ Z q ( Z ) l n { p ( Z ∣ X ) q ( Z ) } d Z KL(q||p)=-\int_Zq(Z)ln\{\frac{p(Z|X)}{q(Z)}\}dZ KL(q∣∣p)=−∫Z​q(Z)ln{q(Z)p(Z∣X)​}dZ

***如果隐变量是连续变量，在以上公式中应该使用积分而不是求和。如果某些变量或者全部的变量都是离散变量，那么分析过程不过，只需根据需要把积分替换为求和即可。***

通过关于概率分布$q(Z)$ 的最优化来使得下界$L(q)$达到最大值。这等价于最小化KL散度。如果允许任意选择$q(Z)$ ，那么下界的最大值出现在KL散度等于零的时刻，此时$q(Z)$等于后验概率分布$p(Z|X)$。

然而，假定在需要处理的模型中，对真实的概率分布进行操作是不行的。于是，转而考虑概率分布$q(Z)$的一个受限制的类别，然后寻找这个类别中使得KL散度达到最小值的概率分布。目标是充分限制$q(Z)$可以取得的概率分布的类别范围，使得这个范围中的所有概率分布都是可以处理的概率分布。施加限制条件的唯一目的是为了计算方便，并且在这个限制条件下，我们应该使用尽可能丰富的近似概率分布。特别地，对于高度灵活的概率分布来说，没有“过拟合”现象。使用灵活的近似仅仅使得我们更好地近似真实的后验概率分布。

==限制近似概率分布的范围的一种方法是使用参数概率分布$q(Z|\omega )$，它由参数集合$\omega$控制。==这样，下界$L(q)$变成了$\omega$的函数。可以利用标准的非线形最优化方法确定参数的最优值。

10.1.1 分解概率分布

使用Mean Field方法，在这种方法中，限制概率分布$q(Z)$的范围。假设将$Z$的元素划分成若干个互不相交的组，记作$Z_i$，其中$i=1,\cdots ,M$。假定$q$分布关于这些分组可以进行分解，即
$$q(Z)\prod _{i=1}^M{q}_i(Z_i)$$
注意，这里没有限制各个因子$q_i(Z_i)$的具体函数形式。变分推断的这个分解的形式对应物理学中的一个近似框架，叫做平均场理论（mean field theory）。

在概率分布$q(Z)$中，寻找下界$L(q)$最大的概率分布。于是，希望对$L(q)$关于所有的概率分布$q_i(Z_i)$进行一个自由形式的（变分）最优化，通过关于每个因子进行最优化完成整体的优化过程 。
L ( q ) = ∫ Z q ( Z ) l n { p ( Z , X ) q ( Z ) } d Z = ∫ Z q ( Z ) { ln ⁡ p ( Z , X ) − ln ⁡ q ( Z ) } d Z = ∫ Z ∏ i = 1 M q i { ln ⁡ p ( Z , X ) − ln ⁡ ∏ i = 1 M q i } d Z = ∫ j q j { ∫ − j ∏ i ≠ j q i ⋅ ln ⁡ p ( Z , X ) d Z i } d Z j − ∫ Z ∏ i = 1 M q i ln ⁡ ∏ i = 1 M q i d Z = ∫ j q j { ∫ − j ∏ i ≠ j q i ⋅ ln ⁡ p ( Z , X ) d Z i } d Z j − ∫ j q j { ∫ − j ( ln ⁡ q j + ∑ i ≠ j ln ⁡ q i ) d Z i }   d Z j = ∫ j q j { ∫ − j ∏ i ≠ j q i ⋅ ln ⁡ p ( Z , X ) d Z i } d Z j − ∫ j q j { ln ⁡ q j + c o n s t − j } d Z j = ∫ j q j { ∫ − j ∏ i ≠ j q i ⋅ ln ⁡ p ( Z , X ) d Z i } d Z j − ∫ j q j ln ⁡ q j   d Z j − c o n s t − j = ∫ j q j ln ⁡ p ^ ( X , Z j ) d Z j − ∫ j q j ln ⁡ q j   d Z j − c o n s t − j \begin{aligned} L(q) &= \int_Z{q(Z)ln\{\frac{p(Z,X)}{q(Z)}\}dZ}\\ &=\int_Z{q(Z)\{\ln p(Z,X)-\ln q(Z)\}}dZ\\ &=\int_Z{\prod\limits_{i=1}^M q_i\{\ln p(Z,X)-\ln \prod\limits_{i=1}^M q_i\}}dZ\\ &=\int_jq_j\left\{\int_{-j}\prod\limits_{i\ne j}q_i \cdot\ln p(Z,X)dZ_i\right\}dZ_j-\int_Z\prod\limits_{i=1}^Mq_i\ln\prod\limits_{i=1}^Mq_idZ\\ &=\int_jq_j\left\{\int_{-j}\prod\limits_{i\ne j}q_i \cdot\ln p(Z,X)dZ_i\right\}dZ_j-\int_jq_j\left\{\int_{-j}\left(\ln q_j+\sum_{i\ne j}\ln q_i\right) dZ_i\right \}\ dZ_j\\ &=\int_jq_j\left\{\int_{-j}\prod\limits_{i\ne j}q_i \cdot\ln p(Z,X)dZ_i\right\}dZ_j-\int_jq_j\left\{\ln q_j+const_{-j}\right\}dZ_j\\ &=\int_jq_j\left\{\int_{-j}\prod\limits_{i\ne j}q_i \cdot\ln p(Z,X)dZ_i\right\}dZ_j-\int_jq_j\ln q_j\ dZ_j-const_{-j}\\ &=\int_jq_j\ln \hat{p}(X,Z_j)dZ_j-\int_jq_j\ln q_j\ dZ_j-const_{-j} \end{aligned} L(q)​=∫Z​q(Z)ln{q(Z)p(Z,X)​}dZ=∫Z​q(Z){lnp(Z,X)−lnq(Z)}dZ=∫Z​i=1∏M​qi​{lnp(Z,X)−lni=1∏M​qi​}dZ=∫j​qj​⎩⎨⎧​∫−j​i​=j∏​qi​⋅lnp(Z,X)dZi​⎭⎬⎫​dZj​−∫Z​i=1∏M​qi​lni=1∏M​qi​dZ=∫j​qj​⎩⎨⎧​∫−j​i​=j∏​qi​⋅lnp(Z,X)dZi​⎭⎬⎫​dZj​−∫j​qj​⎩⎨⎧​∫−j​⎝⎛​lnqj​+i​=j∑​lnqi​⎠⎞​dZi​⎭⎬⎫​ dZj​=∫j​qj​⎩⎨⎧​∫−j​i​=j∏​qi​⋅lnp(Z,X)dZi​⎭⎬⎫​dZj​−∫j​qj​{lnqj​+const−j​}dZj​=∫j​qj​⎩⎨⎧​∫−j​i​=j∏​qi​⋅lnp(Z,X)dZi​⎭⎬⎫​dZj​−∫j​qj​lnqj​ dZj​−const−j​=∫j​qj​lnp^​(X,Zj​)dZj​−∫j​qj​lnqj​ dZj​−const−j​​
其中，定义了一个新的概率分布$\hat{p}(X,Z_j)$
$$\ln \hat{p}(X,Z_j)=E_{i\ne j}[\ln p(X,Z)]+const$$
注意，记号$E_{i\ne j}[\cdots ]$表示关于定义在所有$z_i(i\ne j)$上的$q$概率分布的期望，即
$$E_{i\ne j}[\ln p(X,Z)]=\int \ln p(X,Z)\prod _{i\ne j}{q}_{i}d{Z}_i$$
由于
$${\int }_j{q}_j\ln \hat{p}(X,Z_j)d{Z}_j-{\int }_j{q}_j\ln q_{j}d{Z}_j\iff -KL(q_j||\hat{p}(X,Z_j))$$
所以
$$L(q)=-KL(q_j||\hat{p}(X,Z_j))-cons{t}_j$$
最大化$L(q)$的优化任务$\iff$最小化Kullback-Leibler $KL(q_j||\hat{p}(X,Z_j))$，且最小值出现在$q_j(Z_j)=\hat{p}(X,Z_j)$的位置。因此，最优解$q_j^{\ast }(Z_j)$的一般表达式为：
$$\begin{array}{rl}{q}_j(Z_j)& =\hat{p}(X,Z_j)\\ \ln q_j(Z_j)& =\ln \hat{p}(X,Z_j)\\ \because \ln \hat{p}(X,Z_j)& =E_{i\ne j}[\ln p(X,Z)]+常数\\ & \therefore 最优解的一般式如下：\\ \ln q_j^{\ast }(Z_j)& =E_{i\ne j}[\ln p(X,Z)]+常数\end{array}$$
该解表明，为了得到因子$q_j$的最优解的对数，只需考虑所有隐含变量和可见变量上的联合概率分布的对数，然后关于所有其他的因子 { q i } i ≠ j \{q_i\}_{i\ne j} {qi​}i​=j​取期望即可。该公式中的可加性常数通过对概率分布$q_j^{\ast }(Z_j)$进行归一化的方式来设定，取两侧的指数，然后归一化，得
$$q_j^{\ast }(Z_j)=\frac{\exp (E_{i\ne j}[\ln p(X,Z)])}{\int \exp (E_{i\ne j}[\ln p(X,Z)])d{Z}_j}$$
在实际应用中，更方便的做法是直接对归一化之前的一般解进行操作，然后在必要的时候，通过观察的方式恢复出归一化系数。

由最优解的一般式所给定的方程的集合$j=1,2,\dots ,M$，表示在概率能够进行分解这一限制条件下，下界的最大值满足的一组相容的条件。然而，这些方程并没有给出一个显示的解。因为最优化$q_j^{\ast }(Z_j)$一般公式的右侧表达式$E_{i\ne j}[\ln p(X,Z)]$依赖于其它的因子$q_j(Z_j)i\ne j$计算的期望。于是，用以下的方式寻找出一个相容的解：

首先，恰当地初始化所有的因子$q_i(Z_i)$，然后在各个因子上进行循环，每一轮用一个修正后的估计来替换当前因子。这个修正后的估计由$E_{i\ne j}[\ln p(X,Z)]$给出，计算时使用了当前对于所有其他因子的估计。算法保证手链，因为下界关于每个因子$q_i(Z_i)$是一个凸函数。

10.1.2 分解近似的性质

变分推断的方法基于的是真实后验概率分布的分解近似。以下是使用分解概率分布的方式近似一个一般的概率分布的问题，以下以高斯分布问题为例。考虑两个相关的变量$z=(z_1,z_2)$上的高斯分布$p(z)=N(z|u,{\Lambda }^{-1})$，其中均值和精度的元素为
Cannot read properties of undefined (reading 'type')
并且由于精度矩阵的对称性，${\Lambda }_{21}=\Lambda 12$。现在，假设希望使用一个分解的高斯分布$q(z)=q_1(z_1)q_2(z_2)$

首先，直接使用最优解的一般式来寻找最优因子$q_1^{\ast }(z_1)$的表达式，在等式右侧，只需要保留那些与$z_1$有函数依赖关系的项即可，其他的项可以被整合到归一化常数中。
$$\begin{array}{rl}\ln q_1^{\ast }(z_1)& =E_{z_2}[\ln p(z)]+常数\\ & =E_{z_2}\left[-\frac{1}{2}(z_1-u_1{)}^2{\Lambda }_{11}-(z_1-u_1){\Lambda }_{12}(z_2-u_2)\right]+常数\\ & =-\frac{1}{2}{z}_1^2{\Lambda }_{11}+z_1{u}_1{\Lambda }_{11}-z_1{\Lambda }_{12}(E[z_2]-u_2)+常数\end{array}$$

==在PRML中，作者写道“观察这个表达式的右侧是$z_1$的一个二次函数，因此可以将$q^{\ast }(z_1)$看作一个高斯分布。值得强调的是，不事先假设$q^{\ast }(z_1)$是一个高斯分布，而是通过对所有可能的分布$q(z_i)$上的KL散度的变分最优化推导出了这个结果。另外，不需要显示得考虑最优解一般式中的可加性常数，因为它表示归一化常数。如果需要的话，这个常数可以在计算的最后阶段通过观察的方式得到。”==使用配平方可以得到这个高斯分布的均值和方差，有
$$q_1^{\ast }(z_1)=N(z_1|m_1,{\Lambda }_{11}^{-1})$$
其中
$$m_1=u_1-{\Lambda }_{11}^{-1}{\Lambda }_{12}(E[z_2]-u_2)$$
以下是对式[26]进行配平方的过程
$$\begin{array}{rl}\ln q_1^{\ast }(z_1)& =-\frac{1}{2}{z}_1^2{\Lambda }_{11}+z_1{u}_1{\Lambda }_{11}-z_1{\Lambda }_{12}(E[z_2]-u_2)+常数\\ & =-\frac{z_1^2{\Lambda }_{11}+2\cdot z_1(u_1{\Lambda }_{11}-{\Lambda }_{12}(E[z_2]-u_2))}{2}+const\\ & =-\frac{z_1^2+2\cdot z_1(u_1-{\Lambda }_{11}^{-1}{\Lambda }_{12}(E[z_2]-u_2))}{2{\Lambda }_{11}^{-1}}+const\end{array}$$
由此得知，变分后验分布$q_1(z_1)$ 符合高斯分布，其中均值为$u_1-{\Lambda }_{11}^{-1}{\Lambda }_{12}(E[z_2]-u_2)$，标准差为${\Lambda }_{11}^{-1}$。注意，虽然该推导式并不完整，缺少一些必要的项。但作者强调过只需要保留那些与参数$z_1$有函数依赖关系的项即可，其它的项可以被整合到归一化常数中。

根据对称性，$q_2^{\ast }(z_2)$也是一个高斯分布，可以写成
$$q_2^{\ast }(z_2)=N(z_2|m_2,{\Lambda }_{22}^{-1})$$
其中
$$m_2=u_2-{\Lambda }_{22}^{-1}{\Lambda }_{21}(E[z_1]-u_1)$$
注意：这两个解是相互耦合的，即$q^{\ast }(z_1)$依赖于关于$q^{\ast }(z_2)$计算的期望，反之亦然。解决这类问题的方法如下：将变分解看成重估计方程，然后在变量之间循环，更新这些解，知道满足某个收敛准则。

10.1.3 例子：一元高斯分布

目标：在给定$x$的观测值的数据集$D=x1,\cdots ,x_N$的情况下，推断均值$\mu$ 和精度$\tau$的后验概率分布。其中，假设数据是独立地从高斯分布中抽取的。似然函数为
$$p(D|\mu ,\tau )=(\frac{\tau }{2\pi }{)}^{\frac{N}{2}}\exp \left\{-\frac{\tau }{2}\sum _{n=1}^N(x_n-u{)}^2\right\}$$

p ( x i ∣ μ , σ 2 ) = 1 2 π ⋅ σ exp ⁡ { − ( x i − μ ) 2 2 σ 2 } ∵ 数 据 是 独 立 地 从 高 斯 分 布 中 获 取 ∴ p ( D = { x 1 , x 2 , ⋯   , x N } ∣ μ , σ 2 ) = ∏ i N p ( x i ) p ( D ∣ μ , σ 2 ) = ( 1 2 π ⋅ σ ) N exp ⁡ { − 1 2 σ 2 ∑ n = 1 N ( x n − μ ) 2 } = ( 1 2 π ⋅ σ 2 ) N 2 exp ⁡ { − 1 2 σ 2 ∑ n = 1 N ( x n − μ ) 2 } ∵ τ = 1 σ 2 ∴ p ( D ∣ μ , σ 2 ) = ( τ 2 π ) N 2 exp ⁡ { − τ 2 ∑ n = 1 N ( x n − μ ) 2 } \begin{aligned} p(x_i|\mu,\sigma^2) &= \frac{1}{\sqrt{2\pi}\cdot\sigma}\exp\left\{-\frac{(x_i-\mu)^2}{2\sigma^2}\right\}\\ &\because 数据是独立地从高斯分布中获取\\ &\therefore p(D=\{x_1,x_2,\cdots,x_N\}|\mu,\sigma^2)=\prod\limits_i^Np(x_i)\\ p(D|\mu,\sigma^2)&=(\frac{1}{\sqrt{2\pi}\cdot\sigma})^N\exp\left\{-\frac{1}{2\sigma^2}\sum_{n=1}^N(x_n-\mu)^2\right\}\\ &=(\frac{1}{2\pi\cdot\sigma^2})^{\frac{N}{2}}\exp\left\{-\frac{1}{2\sigma^2}\sum_{n=1}^N(x_n-\mu)^2\right\}\\ \because \tau=\frac{1}{\sigma^2}\\ \therefore p(D|\mu,\sigma^2)&=(\frac{\tau}{2\pi})^{\frac{N}{2}}\exp\left\{-\frac{\tau}{2}\sum_{n=1}^N(x_n-\mu)^2\right\} \end{aligned} p(xi​∣μ,σ2)p(D∣μ,σ2)∵τ=σ21​∴p(D∣μ,σ2)​=2π ​⋅σ1​exp{−2σ2(xi​−μ)2​}∵数据是独立地从高斯分布中获取∴p(D={x1​,x2​,⋯,xN​}∣μ,σ2)=i∏N​p(xi​)=(2π ​⋅σ1​)Nexp{−2σ21​n=1∑N​(xn​−μ)2}=(2π⋅σ21​)2N​exp{−2σ21​n=1∑N​(xn​−μ)2}=(2πτ​)2N​exp{−2τ​n=1∑N​(xn​−μ)2}​

引入$\mu$和$\tau$的共轭先验分布，形式如下：
$$\begin{array}{rl}p(u|\tau )& =N(u|u_{0}m({\lambda }_0\tau {)}^{-}1)\\ p(\tau )& =Gam(\tau |a_0,b_0)\end{array}$$
这些分布共同给出了一个高斯- Gamma共轭先验分布。对于这个问题，后验概率可以求出精确解，并且形式是Gaussian- Gamma分布。对后验概率分布分解变分近似：
$$q(\mu ,\tau )=q_{\mu }(\mu )\cdot q_{\tau }(\tau )$$
注意，真实的后验概率分布不可以按照这种形式进行分解，最优的因子$q_{\mu }\mu$ 和$q_{\tau }(\tau )$可以从一般的结果中得到。对于$q_{\mu }(\mu )$
$$\begin{array}{rl}\ln p(D|\mu ,\tau )& =\ln \left((\frac{\tau }{2\pi }{)}^{\frac{N}{2}}\exp \left\{-\frac{\tau }{2}\sum _{n=1}^N(x_n-\mu {)}^2\right\}\right)\\ & =\frac{N}{2}\ln \frac{\tau }{2\pi }-\frac{\tau }{2}\sum _{n=1}^N(x_n-u{)}^2\end{array}$$

$$\begin{array}{rl}\ln p(\mu |\tau )& =\ln \left(\sqrt{\frac{{\lambda }_0\tau }{2\pi }}\cdot \exp \left\{-\frac{{\lambda }_0\tau }{2}(u-u_0{)}^2\right\}\right)\\ & =\frac{1}{2}\ln \frac{{\lambda }_0\tau }{2\pi }-\frac{{\lambda }_0\tau }{2}(u-u_0{)}^2\end{array}$$

$$\begin{array}{rl}\ln q_{\mu }^{\ast }(\mu )& =E_{\tau }[\ln p(D,\mu |\tau )]+const\\ & =E_{\tau }[\ln p(D|\mu ,\tau )\cdot p(\mu |\tau )]+const\\ & =E_{\tau }[\ln p(D|\mu ,\tau )+\ln p(\mu |\tau )]+const\\ & =E_{\tau }[\frac{N}{2}\ln \frac{\tau }{2\pi }-\frac{\tau }{2}\sum _{n=1}^N(x_n-u{)}^2+\frac{1}{2}\ln \frac{{\lambda }_0\tau }{2\pi }-\frac{{\lambda }_0\tau }{2}(u-u_0{)}^2]+const\\ & =\frac{1}{2}{\int }_{\tau }(N+1)\ln \tau d\tau -\frac{1}{2}(\sum _{n=1}^N(x_n-u{)}^2+{\lambda }_0(u-u_0{)}^2)E[\tau ]+const\\ & =-\frac{1}{2}(\sum _{n=1}^N(x_n-u{)}^2+{\lambda }_0(u-u_0{)}^2)E[\tau ]+const\end{array}$$

注意，在PRML中，作者将${\int }_{\tau }\ln \tau d\tau$也视为常数，这是我不太理解的。

从$\ln q_{\mu }^{\ast }\mu$中对参数$\mu$进行配平方，配平方的关键就是保证配出$\frac{{\mu }^2-2\cdot A\cdot \mu +B}{-2\cdot C}$的形式，其中均值为$A$，方差为$C$。
$$\begin{array}{rl}\ln q_{\mu }^{\ast }\mu & =-\frac{E[\tau ]}{2}\left\{{\lambda }_0(\mu -{\mu }_0{)}^2+\sum _{n=1}^N(x_n-\mu {)}^2\right\}+const\\ & =-\frac{E[\tau ]}{2}\left\{{\lambda }_0{\mu }^2+{\lambda }_0{\mu }_0-2{\lambda }_0{\mu }_0\cdot \mu +N(\bar{x}-\mu {)}^2\right\}\\ & =\frac{({\lambda }_0+N)\cdot {\mu }^2-2({\lambda }_0{\mu }_0+N\bar{x})\cdot \mu +({\lambda }_0{\mu }_0+N{\bar{x}}^2)}{-2E[\tau {]}^{-1}}\\ & =\frac{{\mu }^2-2\frac{{\lambda }_0{\mu }_0+N\bar{x}}{{\lambda }_0+N}\cdot \mu +\frac{{\lambda }_0{\mu }_0+N{\bar{x}}^2}{{\lambda }_0+N}}{-2(({\lambda }_0+N)E[\tau ]{)}^{-1}}\end{array}$$
所以对于隐变量$\mu$的变分后验分布是高斯分布$\mu \sim N(\mu |{\mu }_N,{\lambda }_N^{-1})$，其中，均值和方差分别为：
$$\begin{array}{rl}{\mu }_N& =\frac{{\lambda }_0{\mu }_0+N\bar{x}}{{\lambda }_0+N}\\ {\lambda }_N& =({\lambda }_0+N)E[\tau ]\end{array}$$