贝叶斯网专题9：变分推理-朴素平均场方法

第一部分：贝叶斯网基础

1.1 信息论基础

1.2 贝叶斯网基本概念

1.3 变量独立性的图论分析

第二部分：贝叶斯网推理

2.1 概率推理中的变量消元方法

2.2 团树传播算法

2.3 近似推理

2.3.1 蒙特卡洛方法

2.3.1.1 重要性抽样法

2.3.1.2 马尔可夫蒙特卡洛抽样法(MCMC)

2.3.2 变分推理

上一节我们介绍了利用MCMC抽样进行贝叶斯网后验概率推理的方法，其是一种随机近似算法，存在的缺点是抽样收敛速度慢，马尔可夫链往往需要花费较多步迭代才能趋于平稳分布，且当贝叶斯网中存在极端概率0或1时，无法保证马尔可夫链存在平稳分布。
相对于基于抽样的随机近似方法，变分法是另一类具有确定性的近似方法。该方法借助一个更简单的分布去逼近待推理的复杂后验分布，比如朴素平均场方法利用一个变量间都相互独立的联合分布来逼近后验分布，朴素则意味着变量间的独立性。
如何求解这样的变分分布呢？我们知道，概率分布本质是一类满足归一化条件的函数，我们可以通过施加独立性条件使得变分分布具有易于推理的简单形式，并从这些具有简单形式的函数中寻找一个最逼近后验分布的函数，这其实是泛函问题。其中有两个问题需要解决，一是如何度量变分分布与后验分布之间的距离，二是如何求解该距离的极值点和极值。在介绍朴素平均场这一具体方法之前，我们先介绍泛函中具有更普遍意义的变分原理。

2.3.2.1 变分原理

2.3.2.1.1 问题来源

变分与微分相对应的，变分针对的是泛函，而微分针对的是普通函数。泛函和普通函数的区别在于，普通函数是从数域到数域的映射，而泛函是从函数域到数域的映射。
泛函最早起源于最速降曲线问题，比该问题更简单的是求两点之间最短曲线的问题。我们以后者为例，来介绍泛函求极值的变分原理。
两点间的最短曲线显然是直线，这是我们的常识，但如何严格推导这一结论呢？
首先我们对这一问题进行建模。设$P_1(x_1,y_1),P_2(x_2,y_2)$是平面上给定的两点，$y(x)$为连接两点的任意曲线。于是，该曲线的长度为：
$$L(y)={\int }_{x_1}^{x_2}\sqrt{(dx{)}^2+(dy{)}^2}={\int }_{x_1}^{x_2}\sqrt{1+y^{\prime 2}}dx$$
连接$P_1,P_2$两点的曲线有无数条，每条曲线都有一个$L$值与之对应。在上式中，$y(x)$称为自变函数，$L(y)$是依赖于自变函数的函数，称为泛函。满足边界条件$y(x_1)=y_1,y(x_2)=y_2$的函数$y$称为容许函数，我们需要从这些$y$中找出使得$L$最小的那条。
对于这类问题，我们可以概括出更一般的框架，即求解泛函：
$$J(y)={\int }_{x_1}^{x_2}F(x,y,y^{\prime })dx$$
在边界条件$y(x_1)=y_1,y(x_2)=y_2$的极小值问题。其中，$J(y)$称为最简单的积分形式泛函，简称最简泛函，被积函数$F$称为拉格朗日函数。

2.3.2.1.2 变分分析

微分是解决函数求极值问题，变分则是解决泛函求极值问题。下面给出变分的定义。
假设函数$y(x)$是泛函$L(y)$取极小值时对应的极小值点。邻近于$y(x)$的任意容许函数$\hat{y}(x)$可表示为如下形式：
$$\hat{y}(x)=y(x)+ϵ\eta (x)$$
其中$ϵ$是任意小参数，而$\eta (x)$为边界值为0的任意函数。当$x$固定时，容许函数$\hat{y}(x)$与$y(x)$之差$\delta y$称为泛函自变函数的变分，即：
$$\delta y=\hat{y}(x)-y(x)=ϵ\eta (x)$$
类似地，容许函数$\hat{y}(x)$的斜率与$y(x)$的斜率之差$\delta y^{\prime }$称为泛函自变函数斜率的变分，即：
$$\delta y^{\prime }={\hat{y}}^{\prime }(x)-y^{\prime }(x)=ϵ{\eta }^{\prime }(x)$$
应该注意$\delta y$与函数$y(x)$的微分$dy$之间的差别，$dy$是自变量$x$的改变量$dx$引起的函数$y(x)$的无穷小增量，对应于泛函变分$\delta J(y)$。而变分$\delta y$是$y(x)$的任意微小改变量，其应该与$dx$对应。
设泛函增量为：
$$\Delta J={\int }_{x_1}^{x_2}F(x,\hat{y},{\hat{y}}^{\prime })dx-{\int }_{x_1}^{x_2}F(x,y,y^{\prime })dx$$
对上式第一项按泰勒级数展开：
$$\begin{array}{rl}{\int }_{x_1}^{x_2}F(x,\hat{y},{\hat{y}}^{\prime })dx=& {\int }_{x_1}^{x_2}F(x,y+\delta y,y^{\prime }+\delta y^{\prime })dx\\ =& {\int }_{x_1}^{x_2}F(x,y,y^{\prime })dx\\ & +{\int }_{x_1}^{x_2}(\frac{\partial F}{\partial y}\delta y+\frac{\partial F}{\partial y^{\prime }}\delta y^{\prime })dx\\ & +\frac{1}{2!}{\int }_{x_1}^{x_2}(\frac{{\partial }^{2}F}{\partial y^2}(\delta y{)}^2+\frac{{\partial }^{2}F}{\partial y\partial y^{\prime }}\delta y\delta y^{\prime }+\frac{{\partial }^{2}F}{\partial y^{\prime 2}}(\delta y^{\prime }{)}^2)dx\\ & +\cdots \end{array}$$
记泛函的各阶变分为：
$$\begin{array}{rl}\delta J& ={\int }_{x_1}^{x_2}(\frac{\partial F}{\partial y}\delta y+\frac{\partial F}{\partial y^{\prime }}\delta y^{\prime })dx\\ {\delta }^{2}J& ={\int }_{x_1}^{x_2}(\frac{{\partial }^{2}F}{\partial y^2}(\delta y{)}^2+\frac{{\partial }^{2}F}{\partial y\partial y^{\prime }}\delta y\delta y^{\prime }+\frac{{\partial }^{2}F}{\partial y^{\prime 2}}(\delta y^{\prime }{)}^2)dx\\ & \cdots \\ {\delta }^{k}J& ={\int }_{x_1}^{x_2}(\delta y\frac{\partial }{\partial y}+\delta y^{\prime }\frac{\partial }{\partial y^{\prime }}{)}^{(k)}Fdx\end{array}$$
从而，$\Delta J=\delta J+\frac{1}{2!}{\delta }^{2}J+\cdots$
由于一阶变分$\delta y$与小参数$ϵ$成正比，二阶变分${\delta }^{2}y$与${ϵ}^2$成正比，以此类推。因此，当$ϵ$趋于无穷小时，二阶以上各项变分都是一阶变分的高阶无穷小，可以忽略不记。所以泛函增量$\Delta J$趋于0的条件是一阶变分$\delta J=0$，此时达到$J$的极值点。若要进一步分析泛函的极值是极大值还是极小值，则还需要考虑泛函的二阶变分${\delta }^{2}J$，因为处于极值点时：
$$\Delta J=\frac{1}{2!}{\delta }^{2}J+o({\delta }^{2}J)$$
当${\delta }^{2}J>0$时，$\Delta J>0$，泛函$J(y)$为极小值，反之则为极大值。

2.3.2.1.3 泛函极值的必要条件-欧拉方程

泛函取极值时，其一阶变分为0，即：
$$\delta J={\int }_{x_1}^{x_2}(\frac{\partial F}{\partial y}\delta y+\frac{\partial F}{\partial y^{\prime }}\delta y^{\prime })dx=0$$
由于$\delta y$与$\delta y^{\prime }$不是独立无关的，因此上式第二项可写成（推导过程用到了分部积分法）：
$$\begin{array}{rl}{\int }_{x_1}^{x_2}\frac{\partial F}{\partial y^{\prime }}\delta y^{\prime }dx& ={\int }_{x_1}^{x_2}\frac{\partial F}{\partial y^{\prime }}\frac{d(\delta y)}{dx}dx\\ & ={\int }_{x_1}^{x_2}\frac{\partial F}{\partial y^{\prime }}d(\delta y)\\ & ={\frac{\partial F}{\partial y^{\prime }}\delta y|}_{x_1}^{x_2}-{\int }_{x_1}^{x_2}(\delta y)d\frac{\partial F}{\partial y^{\prime }}\\ & ={\frac{\partial F}{\partial y^{\prime }}\delta y|}_{x_1}^{x_2}-{\int }_{x_1}^{x_2}\frac{d}{dx}(\frac{\partial F}{\partial y^{\prime }})\delta ydx\end{array}$$
由于函数$y(x)$在边界处的值已知，$\delta y$在边界处为0，从而：
$${\frac{\partial F}{\partial y^{\prime }}\delta y|}_{x_1}^{x_2}=0$$
回代则得：
$$\delta J={\int }_{x_1}^{x_2}(\frac{\partial F}{\partial y}-\frac{d}{dx}(\frac{\partial F}{\partial y^{\prime }}))\delta ydx=0$$
由于$\delta y$为在区间$[x_1,x_2]$内关于$x$的任意函数，因此积分为0则需要满足：
$$\frac{\partial F}{\partial y}-\frac{d}{dx}(\frac{\partial F}{\partial y^{\prime }})=0$$
该式称为欧拉方程。
求解该方程，并利用边界条件，则可以确定泛函$J(y)$的极值点$y(x)$。
回到一开始提到的求解两点之间最短曲线的问题。待求曲线为$y(x)$，总弧长为泛函$J(y)={\int }_{x_1}^{x_2}\sqrt{1+y^{\prime 2}}dx$。
将拉格朗日函数$F=\sqrt{1+y^{\prime 2}}$，代入欧拉方程可得：
$$0-\frac{d}{dx}(\frac{y^{\prime }}{\sqrt{1+y^{\prime 2}}})=0$$
从而可得：
$$\frac{y^{\prime }}{\sqrt{1+y^{\prime 2}}}=C\Rightarrow y^{\prime }=\frac{C}{\sqrt{1-C^2}}$$
由此可见斜率$y^{\prime }$为一个常数，该曲线是一条直线。从而证明了两点之间直线最短的常识。进一步地，由边界的$P_1,P_2$条件，可解出直线的具体方程。

2.3.2.1.4 泛函变分的基本运算法则

泛函变分法则和微分运算法则基本相同。设$F_1,F_2$均为$x,y,y^{\prime },\cdots$的函数，则：
$$\begin{array}{rl}& \delta (F_1+F_2)=\delta F_1+\delta F_2\\ & \delta (F_1\cdot F_2)=F_2\delta F_1+F_1\delta F_2\\ & \delta (\frac{F_1}{F_2})=\frac{1}{F_2^2}(F_2\delta F_1-F_1\delta F_2)\\ & \delta (F^n)=n{F}^{n-1}\delta F\end{array}$$

2.3.2.2 变分贝叶斯网推理

上一节我们给出了利用变分法分析泛函极值的框架，并推导出针对最简泛函『拉格朗日函数』的极值必要条件：欧拉方程。
在贝叶斯网后验概率推理问题中，由于实际问题中贝叶斯网往往比较复杂。我们知道，一个分布本质就是一个满足归一化条件的函数。我们借鉴变分分析的思路，通过设定一个简化的分布形式，并尝试从符合这种形式的容许分布中寻找最接近目标后验概率分布的那一个。
为了便于推理，我们所设定的简化分布具有如下朴素形式（变量互相独立）：
$$Q(Z)=\prod _{i=1}^n{Q}_i(Z_i)\text{\hspace{1em}(1)}$$
利用这类具有变量独立性的分布去逼近复杂的后验概率分布$P(Z|E)$，就称为朴素平均场方法。
在本专题第1篇『信息论基础』中介绍了相对熵（KL散度）的概念。我们知道，可以通过KL散度来度量两个概率分布的距离。于是，我们就将问题转化为寻找满足下式条件的具有朴素形式的$Q(Z)$：
$$\begin{array}{rl}Q(Z)& =\underset{Q(Z)}{\mathrm{argmin}}KL(Q(Z)||P(Z|E))\\ & =\underset{Q(Z)}{\mathrm{argmin}}\sum _{Z}Q(Z)\log \frac{Q(Z)}{P(Z|E)}\end{array}\text{\hspace{1em}(2)}$$
我们对KL散度进行变形：
$$\begin{array}{rl}KL(Q(Z)||P(Z|E))& =-\sum _{Z}Q(Z)\log \frac{P(Z|E)}{Q(Z)}\\ & =-\sum _{Z}Q(Z)\log \frac{P(Z,E)}{Q(Z)P(E)}\\ & =-\sum _{Z}Q(Z)\log P(Z,E)+\sum _{Z}Q(Z)(\log Q(Z)+\log P(E))\\ & =-\sum _{Z}Q(Z)\log P(Z,E)+\sum _{Z}Q(Z)\log Q(Z)+\log P(E)\end{array}$$
令$J(Q)=-{\sum }_{Z}Q(Z)\log P(Z,E)+{\sum }_{Z}Q(Z)\log Q(Z)$，则上式可简写为：
$$KL(Q||P)=J(Q)+\log P(E)$$
由于$\log P(E)$为确定值，因此式(2)的优化问题就转化为求解使得$J(Q)$最小的最优变分分布$Q(Z)$。
我们将式(1)关于$Q(Z)$的变量独立条件代入$J(Q)$，并将$J(Q)$分为前后两部分来分析。首先来看后半部分：
$$\begin{array}{rl}{H}_Q(Z)& =-\sum _{Z}Q(Z)\log Q(Z)\\ & =-\sum _{Z_1}\cdots \sum _{Z_n}\left(\prod _{i=1}^n{Q}_i(Z_i)\sum _{j=1}^n\log Q_j(Z_j)\right)\\ & =-\sum _{j=1}^n\sum _{Z_1}\cdots \sum _{Z_n}\left(\prod _{i=1}^n{Q}_i(Z_i)\log Q_j(Z_j)\right)\\ & =-\sum _{j=1}^n\sum _{Z_1}\cdots \sum _{Z_n}\left(\prod _{i\ne j}{Q}_i(Z_i)Q_j(Z_j)\log Q_j(Z_j)\right)\\ & =-\sum _{j=1}^n\underset{except{Z}_j}{\underbrace{\sum _{Z_1}\cdots \sum _{Z_n}}}\prod _{i\ne j}{Q}_i(Z_i)\sum _{Z_j}\left(Q_j(Z_j)\log Q_j(Z_j)\right)\\ & =-\sum _{j=1}^n\sum _{Z_j}\left(Q_j(Z_j)\log Q_j(Z_j)\right)\underset{except{Z}_j}{\underbrace{\sum _{Z_1}{Q}_1(Z_1)\sum _{Z_2}{Q}_2(Z_2)\cdots \sum _{Z_n}{Q}_n(Z_n)}}\\ & =-\sum _{j=1}^n\sum _{Z_j}\left(Q_j(Z_j)\log Q_j(Z_j)\right)\\ & =\sum _{j=1}^n{H}_{Q_j}(Z_j)\end{array}$$
可以看到，后半部分实际就是分布$Q(Z)$的信息熵，由于$Z$的各边缘分布$Q_i(Z_i)$互相独立，因此$Q(Z)$的信息熵等于各边缘分布$Q_i(Z_i)$的信息熵之和。上面给出的是严密的数学推导过程，从信息熵的物理意义更容易直观理解：$Z$的随机程度等于这$n$个独立的$Z_i$的随机程度之和。正是因为有了朴素平均场理论，我们将原本高维的分布$Q(Z)$拆分成单维度形式。
再来看前半部分：
$$\begin{array}{rl}\sum _{Z}Q(Z)\log P(Z,E)& =\sum _{Z_1}\cdots \sum _{Z_n}\prod _{i=1}^n{Q}_i(Z_i)\log P(Z,E)\\ & =\sum _{Z_j}{Q}_j(Z_j)\underset{except{Z}_j}{\underbrace{\sum _{Z_1}\cdots \sum _{Z_n}}}\prod _{i\ne j}{Q}_i(Z_i)\log P(Z,E)\end{array}$$
我们可以将$\underset{except{Z}_j}{\underbrace{\sum _{Z_1}\cdots \sum _{Z_n}}}{\prod }_{i\ne j}{Q}_i(Z_i)\log P(Z,E)$看作函数$\log P(Z,E)$在分布${\prod }_{i\ne j}{Q}_i(Z_i)$上的期望，记作$E_Q\left[\log P(Z,E)|Z_j\right]$，代入上式有：
$$\sum _{Z}Q(Z)\log P(Z,E)=\sum _{Z_j}{Q}_j(Z_j)E_Q\left[\log P(Z,E)|Z_j\right]$$
将前后两部分都代入$J(Q)$可得：
$$\begin{array}{rl}J(Q)=& -\sum _{Z_j}{Q}_j(Z_j)E_Q\left[\log P(Z,E)|Z_j\right]\\ & -\sum _{j=1}^n{H}_{Q_j}(Z_j)\end{array}$$
我们对每一个$Q_j$轮流优化，在优化$Q_j$时，其它 { Q i ∣ i ≠ j } \{Q_i|i\neq j\} {Qi​∣i​=j}都当做常数处理。
除了以上等式，我们还有概率归一化条件：${\sum }_{Z_j}{Q}_j(Z_j)=1$。为了优化$J(Q_j)$，引入拉格朗日乘数$\lambda$，并求解如下方程：
$$\frac{\partial }{\partial Q_j(Z_j=z)}\left[J(Q_j)+\lambda (\sum _{Z_j}{Q}_j(Z_j)-1)\right]=0,\forall z\in {\Omega }_{Z_j}$$
即：
$$-E_Q\left[\log P(Z,E)|Z_j=z\right]+\log Q_j(Z_j=z)+1+\lambda =0$$
从而可解得用于更新$Q_j(Z_j=z)$的平均场方程：
$$Q_j(Z_j=z)=\frac{1}{a}{e}^{E_Q\left[\log P(Z,E)|Z_j=z\right]}$$
其中$a$是归一化常数，
$$a=\sum _{Z_j}{e}^{E_Q\left[\log P(Z,E)|Z_j\right]}$$
依次轮流对$Q_j$进行优化，每一次迭代都使$J(Q)$减小，因此最终算法会收敛。但值得注意的是，最终不一定收敛到全局最优，初始点的选取和迭代更新的顺序都会影响到解的质量。

朴素平均场法使用朴素分布直接对目标后验概率分布进行逼近，概念清晰，算法收敛速度比随机抽样法块。当网络中变量数目多，连接稠密，但变量之间的依赖关系较弱时，该方法可以给出较好的结果。但是在实际中，由于变分平均场的假设过于简单，在变量间依赖关系较强的情况下，不能很好地逼近真实后验分布。