9 概率图模型(三)：推断

推断(Inference) 这个词，对于有一定机器学习基础的同学来说，一定是听说过，这也是贝叶斯方法中一个非常重要的理论性研究。那么什么是推断呢？推断说白了，就是求概率。比如，对于一个联合概率密度函数$p(x)=p(x_1,x_2,\cdots ,x_p)$。我们需要求的有哪些呢？

1. 边缘概率：$p\left(x_i\right)={\sum }_{x_1}\cdots {\sum }_{x_{i-1}}\cdots {\sum }_{x_{i+1}}\cdots {\sum }_{x_p}p(x)$
2. 条件概率：$p\left(x_A|x_B\right),$ 今 $x=x_A\cup x_B$
3. MAP Inference：也就是$\hat{z}=\arg {\max }_{z}p(z|x)\propto \mathrm{argmax}p(x,z)$ 因为, $p(z|x)=\frac{p(x,z)}{p(x)}\propto p(x,z)$

因为我们的目标是求一个最优的参数z，并不需要知道具体的数值是多少，只要知道谁大谁小就行，所以$p(x)$ 可以直接不看了。
现在我们知道了，Inference 在求什么？下一步，我们要总结Inference 有哪些方法。

1 背景

1.1 Inference 求解方法

1.1.1 精准推断(Deterministic Inference)

• Variable Elimination (VE)，变量消除法；
• Belief Propagation (BP) 信念传播，这个可不是我们之前学习的反向传播算法，这里需要注意。同时这个算法衍生出的Sum Product Algorithm，这就是推断的核心，这是一种树结构的；
• Junction Tree Algorithm，这是一种普通图结构。

1.1.2 近似推断(Approximate Inference)

• 典型的有有向环(Loop Belief Propagation)；
• 采样方法，包括Mente Carlo Inference：Importance Sampling，MCMC；
• 变分推断(Variational Inference)。

1.2 隐马尔可夫模型(Hidden Markov Model)

Hidden Markov Model (HMM) 算法将在后面的章节中做详细的描述，在这一小节中，我们主要做一点概述性的描述。HMM 的模型可视化为上图所示，其中O 是隐变量，也就是我们的观测变量。
我们主要考虑三个问题，在三个问题为Inference 的问题。1. Evaluation，也就是求一个边缘密度$p(O)={\sum }_{I}P(I,O)$。2. Learning，也就是寻找${\hat{\lambda }}_{\circ }$ 3. Decoding: $\hat{I}={\mathrm{argmax}}_{I}P(I|O)$，包括Vitebi Algorithm，这是一个动态规划算法。而隐马尔可夫模型实际上一种动态规划模型(Dynamic Bayesian Network)。

2 变量消除：Variable Elimination

在上一小节中，我们简单的介绍了推断的背景和分类，我们知道了大致有哪些推断的方法。推断的任务可以被我们介绍为：给定已知的$p(x)=p(x_1,x_2,\cdots ,x_p)$，我们需要求的有三个：

1. 边缘概率：$p\left(x_i\right)={\sum }_{x_1,\cdots ,x_{i-1},x_{i+1},\cdots ,x_p}p\left(x_1,x_2,\cdots ,x_p\right)$
2. 条件概率：$p(x_A|x_B)$，也就是在已知$x_B$ 集合的情况下，如何求得xA 集合的概率。
3. 最大后验概率(MAP)：${\hat{x}}_A={\mathrm{argmax}}_{x_A}p\left(x_A|x_B\right)=\mathrm{argmax}p\left(x_A,x_B\right)$。
   下面我们要介绍最简单的一个精确推断中中的东西，名为变量消除法(Variable Elimination)。这是一种最简单的推断方法，也是我们学习推断法的核心概念之一。下面我们做详细的解释。

2.1变量消除法(Variable Elimination Algorithm)

假如我们有一个马氏链：

那么我们怎么来求$p(d)$ 呢？根据公式我们可以得到：
$$p(d)=\sum _{a,b,c}p(a,b,c,d)$$
然后使用因子分解，我们可以得到：$$p(d)=\sum _{a,b,c}p(a)p(b|a)p(c|b)p(d|c)$$
假定，$a,b,c,d$ 都为均匀离散的二值random variable，所以$a,b,c,d\in$ {0,1}。
所以，
$$\begin{array}{rl}p(d)& =p(a=0)\cdot p(b=0|a=0)\cdot p(c=0|b=0)\cdot p(d|c=0)\\ & +p(a=1)\cdot p(b=0|a=1)\cdot p(c=0|b=0)\cdot p(d|c=0)\\ & +\cdots \\ & +p(a=1)\cdot p(b=1|a=1)\cdot p(c=1|b=1)\cdot p(d|c=1)\end{array}$$
实际上，这里有8 个因子的积。那么我们来做进一步的分析，我们可以令
$$\begin{array}{rl}p(d)& =\sum _{a,b,c}p(a)p(b|a)p(c|b)p(d|c)\\ & =\sum _{b,c}p(c|b)p(d|c)\sum _{a}p(a)p(b|a)\end{array}$$
而$p(a)p(b|a)=p(a,b)$，而${\sum }_{a}p(a)p(b|a)=p(b)$。我们可以将$a$ 看成$\varphi (a)$ 这是一个和$a$ 相关的函数，同理$p(b|a)$看成$\varphi (a,b)$。所以，我们可以将${\sum }_{a}p(a)p(b|a)$ 看成${\varphi }_a(b)$，这样就相当于一个关于$b$ 的函数，并且是从$a$ 中导出的。所以，我们做如下替换可得：
$$\sum _{b,c}p(c|b)p(d|c)\sum _{a}p(a)p(b|a)=\sum _{b,c}p(c|b)p(d|c){\varphi }_{a}b=\sum _{c}p(d|c)\sum _{b}p(c|b){\varphi }_{a}b$$
同理，我们将${\sum }_{b}p(c|b){\varphi }_{a}b$ 看成${\varphi }_b(c)$。所以，原始将被改写为：$$\sum _{c}p(d|c){\varphi }_b(c)={\varphi }_b(c)$$
这个算法的核心就是乘法对加法的分配律。那我们怎么类比到乘法的分配律呢？首先先来简单的回顾一下乘法的分配律，也就是$ac+ab=a(b+c)$。那么我们仔细的来看看这个计算p(d) 的过程。这是不是就是一个不断的提取公因子，进行计算的过程？有没有觉得和分配律很像？先提取$a$ 的部分，计算$a$ 的部分，然后再依次的提取$b$ 的部分，$c$ 的部分，最后剩下的就是$d$ 的部分。那么，我们就可以把这么一长串的公式进行逐步化简了，这就是变量消元的思想。同样，在无向图中，我们也可以使用到马尔可夫网络中。
$$p(a,b,c,d)=\frac{1}{z}\prod _{i=1}^k{\varphi }_{c_i}(x_{c_i})$$
写成因子分解的形式就是$$p(x)=\prod _{x_i}{\varphi }_i(x_i)$$
这实际上就是分配律，一个变量一个变量的提取，然后进行分解计算。同时这种算法的缺点也非常的明显。
首先，就是重复计算的问题，无论计算那个变量的概率都要重复的计算一遍所有的概率。这个原因就会导致算法的计算难度非常的大。第二个就是计算次序的问题，我们举的例子还比较的简单，所以我们可以一眼就看出来，按$a-b-c-d$ 的次序开始算。但是，实际上，并没有这么容易就得到计算的次序，而且计算次序不一样会导致计算的难度有很大的区别。而有数学家已经证明了，确定最优的计算顺序，本身就是一个NP hard 的问题，非常难求解。

3 Belief Propagation:信念传播

在上一小节中，我们已经介绍了变量消除(Variable Elimination)，Variable Elimination 的思想是 Probability Graph 中的核心思想之一。上一节中我们就已经介绍了，这实际上就是乘法对加法的分配律。但是，Variable Elimination 中有很多的问题，比如重复计算和最优计算次序不好确定的问题。所以，我们这一节来介绍Belief Propagation 来解决重复计算的问题。

3.1 Forward and Backward Algorithm

假设，我们现在有一个马氏链模型：

联合概率可以被我们表示为：$p(a,b,c,d,e)=p(a)\cdot p(b|a)\cdot p(c|b)\cdot p(d|c)\cdot p(e|d)$。
如果，我们想要求的是$p(e)$，那么：
$$\begin{array}{l}p(e)={\sum }_{a,b,c,d}p(a,b,c,d,e)\\ \\ ={\sum }_{a}p(e|d)\cdot {\sum }_{c}p(d|c)\cdot {\sum }_{b}p(c|b)\cdot {\sum }_{a}p(b|a)\cdot p(a)\end{array}$$
为了简化表达，这里我们需要定义一个很重要的符号。因为${\sum }_{a}p(b|a)\cdot p(a)$，是一个关于b 的表达式，也就是相当于把a 给约掉了。所以我们可以把${\sum }_{a}p(b|a)\cdot p(a)$记为$m_{a\to b}(b)$。同理，我们也可以将${\sum }_{c}p(d|c)\cdot {\sum }_{b}p(c|b)\cdot m_{a\to b}(b)$记为$m_{b\to c}(c)$。那么为了求得p(e)，我们依次的求解顺序为$m_{a\to b}(b)$，$m_{b\to c}(c)$，$m_{c\to d}(d)$和$m_{d\to e}(e)$。也就相当于沿着这个链这个马氏链一直往前走，也就是前向算法(Forward Algorithm)。我们用公式表达即为：

如果是要求$p(c)$，那么我们的传递过程为$a\to b\to c\leftarrow d\leftarrow e$。这里，我们就不能只用前向算法来解决了，需要用到Forward-Backward 算法来解决了。也就是同时使用Forward 和Backward的方法，那么我们来看看$p(c)$怎么求？
$$\begin{array}{rl}p(c)& =\sum _{a,b,d,e}p(a,b,c,d,e)\\ & \\ & =(\sum _{b}p(c|b)\sum _{a}p(b|a)p(a))\cdot (\sum _{d}p(d|c)\sum _{e}p(e|d))\end{array}$$
对比上面的计算$p(e)$ 的过程，我们就可以发现，${\sum }_{b}p(c|b){\sum }_{a}p(b|a)p(a)$部分的计算也就是$m_{b\to c}(c)$的计算是一模一样的。所以说，Variable Elimination 里面有大量的重复计算。Belief 的想法很简单，也就是将$m_{i\to j}(j)$，全部事先计算好，就像一个个积木一样，然后再用这个积木来搭建运算。那么也就是，我们事先将方向全部定义好，正向和反向的全部都求了再说。为了进一步探究Belief Background，我们需要讨论一些更加Generalize 的情况。也就是从$Chain\to Tree$，有向$\to$无项的情况。

3.2 Belief Propagation 的扩展

我们的Generalize 的后，分析了一个树形的无向图结构。图的网络结构如下所示：

下面第一步，我们把上面那个模型的设置写出来。所以，我们需要进行因式分解，我们用$\phi (i)$来表示和$i$ 有关的部分。所以，我们可以将联合概率密度写为：$$p(a,b,c,d)=\frac{1}{Z}{\phi }_a(a){\phi }_b(b){\phi }_c(c){\phi }_d(d)\cdot {\phi }_{a,b}(a,b){\phi }_{b,c}(b,c)\phi (b,d)(b,d)$$
我们要求$p(a)={\sum }_{b,c,d}p(a,b,c,d)$和$p(b)={\sum }_{a,c,d}p(a,b,c,d)$，其间一定会出现大量的重复计算。这个模型中有四个点，三条边，每条边都有两个方向，所以我们要求的是6 个“积木”。我们来一步步
的看看，如何可以得到想要的$p(a)$。

1. 首先，我们需要求的是$c\to b$和$d\to b$两个过程。其中，$c\to b$的过程也就是$m_{c\to b}(b)$可以被我们表达为${\sum }_c{\phi }_c\cdot {\phi }_{b,c}$。同理，$m_{d\to b}(b)$可以被我们表达为${\sum }_d{\phi }_b\cdot {\phi }_{b,d}$
2. 第二步，我们需要求$b\to a$的过程，也就是$m_{b\to a}(a)$。它等于$m_{c\to b}(b)$，$m_{d\to b}(b)$乘上$b$己的部分求和得到，我们可以写为：$$m_{b\to a}(a)=\sum _b{m}_{c\to b}(b)\cdot {\phi }_b\cdot {\phi }_{a,b}\cdot m_{d\to b}(b)$$
3. 最后，$m_{b\to a}(a)$乘上a 自己的部分就得到了$p(a)$，也就是：$p(a)={\phi }_a\cdot m_{b\to a}(a)$。
   所以，我们总结一下：
   $$m_{b\to a}(x_a)=\sum _{x_b}{\phi }_{a,b}{\phi }_b{m}_{c\to b}(x_b)m_{d\to b}(x_b)$$
   而$$p(a)={\phi }_a\cdot m_{b\to a}(a)$$。
   我相信到这里，大家应该是可以理解这个意思的，会有点抽象，但并不是很难。下一步我想做个Generalize 为了便于大家进行理解，我这里尽量不跳步：
   m b → a ( x a ) = ∑ x b φ a , b φ b ∏ { k ∈ N B ( b ) } − a m k → b ( x b ) m _ { b \rightarrow a } ( x _ { a } ) = \sum _ { x _ { b } } \varphi _ { a , b } \varphi _ { b } \prod _ { \{ k \in N B ( b ) \} - a } m _ { k \rightarrow b } ( x _ { b } ) mb→a​(xa​)=xb​∑​φa,b​φb​{k∈NB(b)}−a∏​mk→b​(xb​)
   这里的NB(b) 代表的是所有节点b 的邻接节点。我们可以进一步表示为：

m j → i ( x i ) = ∑ x j φ i , j φ j ∏ { k ∈ N B ( j ) } − i m k → j ( x j ) m _ { j \rightarrow i } ( x _ { i } ) = \sum _ { x _ { j } } \varphi _ { i , j } \varphi _ { j } \prod _ { \{ k \in N B ( j ) \} - i } m _ { k \rightarrow j} ( x _ { j } ) mj→i​(xi​)=xj​∑​φi,j​φj​{k∈NB(j)}−i∏​mk→j​(xj​)
而
$$p(x_i)={\phi }_i\prod _{k\in NB(i)}{m}_{k\to i}(x_i)$$
通过对上面表达式的观察，我们是不是发现了一个很有意思的现象。也就是这些概率都是由$m_{i\to j}$这样的小积木拼接起来的。所以我们可以get a conclusion：我们不要一上来就直接去求边缘概率密度，比如$p(a),p(b),p(c),p(d)$这些的。我们可以先建立一个Cache，把$m_{i\to j}$全部算出来。然后，要求什么的话，直接进行搭建和拼接就可以了。从这里，我们就引出了Belief Propagation。

3.2.1 Sequential Implementation

顺序计算的思路，我们需要借助一个队列来实现：

1. Get Root: 首先我们需要假设一个节点为根节点。
2. Collect Message: 对于每一个在根节点的邻接点中的节点$x_i$，Collect Message ($x_i$)。对应的就是图2 中蓝色的线条。
3. Distribute Message: 对于每一个在根节点的邻接点中的节点$x_i$Distribute Message ($x_i$)。对应的就是图2 中红色的线条。

经过这三个步骤以后，我们可以得到所有的$i,j\in v$，从而计算出$p(x_k),k\in v$。

3.2.2 Parallel Implementation

这种思想在图结构的网络中经常使用，大致也就是随意选一个节点，然后向四周其他的节点辐射信息。其他节点收到信息之后，更新自己的状态，然后反馈自己的信息给源节点，源节点再更新自己的信息。不断地重复这个工作，直到最后收敛为止。

3.3 小结

实际上，Belief Propagation 就是一种Variable Elimination。但是，我们发现了Variable Elimination中有很多的重复计算。所以，我们想到了提取出来先算好，要用的时候直接放进去就行了。所以，Belief Propagation 中分解了传递的过程，先计算消息传递的机制，再来组装出计算边缘概率。其实本质还是Variable Elimination 算法，不过就是使表达更加的规范了，通过拆解的方法来消除重复计算。
BP算法的优点：

• 缓存$m_{i\to j}$，减少重复计算，时间复杂度低
• 计算$m_{i\to j}$的过程中不想干的分支可以进行并行计算

4 Max Product Algorithm

我们在这里再总结一下概率图模型有什么用。对于一个图，Graph ={$X,E$}，其中$X$ 代表的是普通变量，$E$ 代表的是Evidence，也就是观测变量。

1. 首先要解决的是边缘变量的问题，也就是已知：$E=${$e_1,e_2,\cdots ,e_k$}，如何求$p(E)$的问题，其中$E$ 为一个变量或者为一个子集。实际上就是一个likelihood 的问题。
2. 条件概率，也就是一个求后验概率的问题，目标概率为$X=(Y;Z)$。而$p(Y|E)={\sum }_{z}p(X|E)$。
3. 最大后验估计(MAP)，也被我们称为Decoding 的问题。也就是我们希望找到一个隐序列，使得：$\hat{x}=argma{x}_{x}p(X|E)，\hat{y}=argma{x}_{y}p(Y|E)$。

这里的Max-Product 算法和隐马尔可夫模型(HMM) 中的Viterbi 算法非常的类似。其实，从算法上讲它就是Belief Propagation 算法的一种改进，从模型上讲是Viterbi 算法的推广。在这里我们求的不是概率值了，而是一个最优的概率序列$(\hat{a},\hat{b},\hat{c},\hat{d})={\mathrm{argmax}}_{a,b,c,d}p(x_a,x_b,x_c,x_d|E)$

4.1 Max Product Algorithm

下面展示一个树的拓扑结构图。

在这个树中，我们将$m_{b\to a}$看成是能使$p(x_a,x_b,x_c,x_d|E)$ 联合概率达到最大的值。每一个节点代表的是到这个节点为止的路径联合概率达到最大的值。我们表达为：$$m_{j\to i}=\underset{x_j}{\max }{\varphi }_j{\varphi }_{ij}\prod _{k\in Neighbour(j)-i}{m}_{k\to j}$$
那么，在图一所示的概率图模型中，$m_{c\to b}$ 可以表示为：$$m_{c\to b}=\underset{x_c}{\max }{\phi }_c\cdot {\phi }_{bc}$$
其中，${\phi }_c\cdot {\phi }_{bc}$ 可以表示为和$c$ 相关的函数。
$$m_{d\to b}=\underset{x_d}{\max }{\phi }_d\cdot {\phi }_{bd}$$
其中，${\phi }_d\cdot {\phi }_{bd}$ 可以表示为和$d$ 相关的函数。
$$m_{b\to a}=\underset{x_b}{\max }{\phi }_b\cdot {\phi }_{ab}\cdot m_{c\to b}\cdot m_{d\to b}$$
最终，我们将得到的是：$$\max p(a,b,c,d)=\underset{x_a}{\max }{\phi }_a\cdot m_{b\to a}$$
而${\phi }_a\cdot m_{b\to a}$ 就可以看成是一个关于a 的函数。这里我们再提一下Belief Propagation，这里的Max-Product 实际上就是Belief Propagation 的一个变形。

4.2 Belief Propagation

实际上这个算法的提出时因为，多次求边缘概率密度会发现中间有很多的步骤是重复的。我们用$m_{i\to j}$ 记录每一个边缘概率，最后进行组合就行。所以，
m j → i ( x i ) = ∑ x j φ i , j ( x i , x j ) φ j ( x j ) ∏ { k ∈ N B ( j ) } − i m k → j ( x j ) \left. \begin{array} { c } { m _ { j \rightarrow i } ( x _ { i } ) = \sum _ { x _ { j } } \varphi _ { i , j } ( x _ { i } , x _ { j } ) \varphi _ { j } ( x _ { j } ) \prod _ { \{ k \in N B ( j ) \} -i } m _ { k \rightarrow j } ( x _ { j } ) } \end{array} \right. mj→i​(xi​)=∑xj​​φi,j​(xi​,xj​)φj​(xj​)∏{k∈NB(j)}−i​mk→j​(xj​)​
$$\begin{array}{c}p(x_i)=\phi (x_i){\prod }_{k\in NB(x_i)}{m}_{k\to i}|(x_i)\end{array}$$

4.3 Compare

其实，我们一比较就可以很简单的看出，Max-product 和Belief Propagation 只有一个地方不一样。那就是前者是求最大，后者是求和。也就是Max-product 到Sum-product。在求得了$maxp(a,b,c,d)=ma{x}_{x_a}{\phi }_a{m}_{b\to a}$ 之后，我们利用回溯法我们比较就可以比较简单的得到$x_a^{\ast },x_b^{\ast },x_c^{\ast },x_d^{\ast }$了。在这个算法中，我们就不需要事先计算$m_{i\to j}$ 了，直接在迭代中进行计算就可以了，也不会存在什么重复计算的问题。