HMM的应用与前向、后向、Viterbi算法

在上一篇文章机器学习中的隐马尔科夫模型（HMM）详解中，我们已经向大家介绍了HMM的基本概念。本文将通过几个例子来演示HMM可以帮我们做些什么事情。特别地，在这个过程中，我们还将向你介绍HMM中非常重要的前向（Forward）算法以及（基于动态规划思想设计的）维特比（Viterbi）算法。

• 一个引子
• 前向（Forward）算法
• 维特比（Viterbi）算法
• HMM在自然语言处理中的应用举例

欢迎关注白马负金羁的博客 http://blog.csdn.net/baimafujinji，为保证公式、图表得以正确显示，强烈建议你从该地址上查看原版博文。本博客主要关注方向包括：数字图像处理、算法设计与分析、数据结构、机器学习、数据挖掘、统计分析方法、自然语言处理。

---

一个引子

回忆一下，我们通常如何定义一个HMM，上篇文章已经告诉我们定义一个HMM模型，需要三个参数，或者说我们可以把一个HMM记作 $\lambda(A,B,\pi)$。其中，矩阵 $A$ 用来表示转移概率 $P(q_t|q_{t-1})$，$B$ 用来表示测量概率 $P(y_t|q_t)$，最后 $\pi$ 指示了HMM模型的初始状态。

继续上篇文章中最后的那个讨论，假设现在我们有如下这样一个HMM，而且我们已经得到了连续的三个观测值 $y_1=up, y_2=up, y_3=down$，能否据此得到相应的隐状态序列呢？

这其实是HMM可以实现的第一种功能。有时我们称之为Computing Likelihood，其实它所要做的就是 Evaluate $P(Y|\lambda)$。因为我们想知道隐藏在观测值背后的隐状态序列，其实就是要看哪个隐状态序列可以使得上面这个观测值的概率最大。而且在上一篇文章的末尾我也得出了针对上述问题的一个递推公式，即
$$P(y_1,y_2,y_3)=\sum_{q_1}^k\sum_{q_2}^k\sum_{q_3}^kP(y_3|q_3)×P(q_3|q_2)×P(y_2|q_2)×P(q_2|q_1)×P(y_1|q_1)×P(q_1)$$
根据已有的HMM参数以及上面这个公式，当隐状态 $q_1,q_2,q_3$ 取不得不同的值时，显然我们可以算出不同的概率，然后其中取值最大的对应序列就是我们想要找的隐状态。具体来说，我们可以列出下面这张表，其中 $q_1=1,q_2=1,q_3=2$ 就是我们最终要找的隐状态序列。

把这种解法泛化，可得下式（注意其中 $\lambda$ 是已知条件，但是从第二行开始为了书写上的简便，我们将其略去）

$$\begin{split} P(Y|\lambda)&=\sum_Q[P(Y,Q|\lambda)]= \sum_{q_1=1}^k\cdots \sum_{q_T=1}^k[P(y_1,\cdots,y_T, q_1,\cdots,q_T|\lambda)]\\\&=\sum_{q_1=1}^k\cdots \sum_{q_T=1}^kP(q_1)P(y_1|q_1)P(q_2|q_1)\cdots P(q_t|q_{t-1})P(y_t|q_t) \\\&=\sum_{q_1=1}^k\cdots \sum_{q_T=1}^k\pi(q_1)\prod_{t=1}^Ta_{q_{t-1},q_t}b_{q_t}(y_t) \end{split}$$
我们对最后一行的标记稍作解释。因为 $P(q_1)$ 是由初始状态决定的，所以它被记作 $\pi(q_1)$。矩阵 $A$ 用来表示转移概率 $P(q_t|q_{t-1})$，具体来说，其中第 $m$ 行，第 $n$ 列所给出的值就是 $P(q_t=n|q_{t?1}=m)$，所以矩阵中的 $P(q_2|q_1),\cdots,P(q_t|q_{t-1})$ 就可以记作 $\prod_{t=1}^Ta_{q_{t-1},q_t}$ ，其中 $a_{q_{t-1},q_t}$ 就是矩阵中第 $q_{t-1}$ 行，第 $q_t$ 列的项。类似地，（注意我们现在讨论的观测值服从离散分布），$B$ 用来表示测量概率 $P(y_t|q_t)$，所以矩阵$B$ 中的项 $b_{q_t}(y_t)$ 就对应概率 $P(y_t|q_t)$。

上这种做法确实可以得出我们想要的结果，但是这种做法的计算量过大。你很容易发现，它的复杂度是 $O(k^T)$ ，实际应用中需要一种更高效的算法，这就是本文后续要讨论的内容。

---

前向（Forward）算法

首先针对下图所示的情况，我们定义一个新的记号：

$$\alpha_i(t)=P(y_1,y_2,\cdots, y_t, q_{t}=i\ |\ \lambda)$$
其注意这里的 $\alpha$ 和 $a$ 是不同的。我们的新记号 $\alpha_i(t)$ 给出了当状态 $q_{t}=i$ 时，它与到时刻 $t$ 为止所有的观测值的联合分布的概率。

借助上面这个定义，我们便可以简化原来相当繁杂的运算。不妨来看看 $\alpha_i(1)$ 所表示的内容：
$$\alpha_i(1)=P(y_1, q_{1}=i)=P(y_1|q_1=i)\cdot P(q_1)=b_i(y_1)\cdot \pi(q_1)$$
同理，我们还可以写出$\alpha_j(2)$ 所表示的内容：
$$\begin{split} \alpha_j(2)=P(y_1,y_2, q_{2}=j)&=\sum_{i=1}^kP(y_1,y_2,q_1=i,q_2=j)\&=\sum_{i=1}^kP(y_2|q_2=j)\cdot P(q_2=j|q_1=i)\cdot P(y_1,q_1=i)\&=\sum_{i=1}^kP(y_2|q_2=j)\cdot P(q_2=j|q_1=i)\cdot \alpha_i(1)\&=P(y_2|q_2=j)\sum_{i=1}^kP(q_2=j|q_1=i)\cdot \alpha_i(1)\&=b_j(y_2) \sum_{i=1}^k a_{i,j} \cdot \alpha_i(1) \end{split}$$
可见 $\alpha_i(2)$ 和$\alpha_i(1)$ 之间产生了递推关系，所以我们可以写出$\alpha_j(T)$ 的表达式：
$$\alpha_j(T)=b_j(y_T) \sum_{i=1}^k a_{i,j} \cdot \alpha_i(T-1)$$
到此为止，其实我们已经得到了一个相当高效的计算 $P(y_1,y_2,\cdots,y_T)$ 的方法。因为根据定义，我们可以知道 $\alpha_j(T)=（y_1,y_2,\cdots,y_T,q_t=j)$ 。而我们最终要求的问题就可以通过下式得到：
$$P(y_1,y_2,\cdots,y_T)=\sum_{j=1}^k\alpha_j(T)$$
而且现在问题的复杂度已经降为了 $O(K^2T)$。

---

维特比（Viterbi）算法

我们已经知道，前向算法的目的在于求 $P(y_1,y_2,\cdots,y_t)$，之所以在求解公式前面会出现很多的求和，这是因为，当给定一个观测序列（例如UP，UP，DOWN）时，最终的概率值其实是所有可能的情况的概率加总。或者说，针对“UP，UP，DOWN”这里例子而言，前向算法求出的其实是本文开头给出的那个表中所有行的概率之和。而这个表中的概率一共有 $K^T$ 行（这个计算量太大），如果采用前向算法，我们就可以大大简化这个计算量。

但现在我们要回到最开始抛出的那个问题，如果给定一个观测序列，例如“UP，UP，DOWN”，那最可能的隐状态序列应该是什么？注意这个问题的答案不再整张概率表，而其实是要求其中（概率值最大的）一行。当然，我们的困难可能在于如果不把整张表列出来，那又如何知道哪一行的概率值最大呢？其实根据前向算法已经给了我们启示。

对于任何一个包含有隐变量的模型来说，例如HMM，当我们有一组观测值时，确定相应的最有可能之隐变量的方法通常又称为 decoding。

HMM中用来decoding的算法就是所谓的维特比算法。它是一种基于动态规划思想设计的算法，最初用于数字通信链路中的卷积码decoding。它的提出者 Andrew Viterbi 被称为是现代CDMA技术之父，他同时还是高通公司的联合创始人和IEEE院士。

经由最可能的状态序列 $q_0, q_1, \cdots, q_{t-1}$ 后进入状态 $q_t=j$ 时，前面这些状态序列与到 $t$ 时刻为止所有观测值的联合分布之概率，我们记为

$$v_j(t)=\max_{q_0,q_1,\cdots,q_{t-1}}P(q_1,q_2,\cdots,q_{t-1},y_1,y_2,\cdots,y_t,q_t=j \ | \ \lambda)$$

注意，我们表示“最可能的状态序列”的方法就是在之前所有的状态序列里找到使上述概率却得最大的那个状态序列，这也同当前我们所面对的问题相一致。

$$\begin{split} v_j(t)&=\max_{q_1,\cdots,q_{t-1}}P(q_1,q_2,\cdots,q_{t-1},y_1,y_2,\cdots,y_t,q_t=j)\&=\max_{q_1,\cdots,q_{t-1}}P(q_1,q_2,\cdots,q_{t-1},y_1,y_2,\cdots,y_{t-1})\cdot P(y_t\ |\ q_t=j)\cdot P(q_t=j\ |\ q_{t-1})\&=\max_{q_1,\cdots,q_{t-1}}P(q_1,q_2,\cdots,q_{t-1},y_1,y_2,\cdots,y_{t-1})\cdot b_j(y_t)\cdot a_{q_{t-1},\ j}\&=b_j(y_t) \max_i^k v_i(t-1)\cdot a_{i,\ j} \end{split}$$

于是我们同样得到了一个递归的公式（注意这个递归公式的base case由初始状态给定）。这也就是维特比算法的核心原理。此外，由于我们最终要确定最大可能的隐状态序列，所以记 $\psi_j(t)$ 是在 $t-1$ 时刻的一个状态，且从该状态转换到状态 $j$ 的概率最大，即

$$\psi_j(t)= b_j(y_t) \arg\max_i^k v_i(t-1)\cdot a_{i,\ j}$$
接下来，我们回到开篇给出的那个例子来演示一下维特比算法的执行过程。首先，我们根据已知条件 $\pi$ 来初始化在 $t=1$ 时刻的概率，则有
$v_{j=1}(t=1)=0.5\times 0.7=0.35$，$v_{j=2}(t=1)=0.2\times 0.1=0.35$， $v_{j=3}(t=1)=0.3\times 0.3=0.09$，如下图中最左边一列所示。接下来，当 $t=2$ 时，要求确定 $v_{j=1}(t=2)$ 的值时，对于 $i=1,2,3$，我们必须考虑三个值，即 $v_i(1)a_{i,1}$ , 其中第一个最大，所以根据算法执行流程，我们得到 $v_{j=1}(t=2)=0.147$ 。也就是下图中红色实线所标示的转移路径，其他两条在后续的步骤中不会再参与计算我们将其用虚线标出。而且，我们还可以得到 $\psi_1(2)=1$，这表明当 $t=2$ 时，转入状态 $j=1$ 前的的状态是 $j=1$ 。

同理，当要确定 $v_j=2(t=2)$ 的值时，对于 $i=1,2,3$，我们必须考虑三个值，即 $v_i(1)a_{i,2}$ , 于是得到$v_{j=2}(t=2)=0.007$ 。要确定 $v_{j=3}(t=2)$ 的值时，对于 $i=1,2,3$，我们必须考虑三个值，即 $v_i(1)a_{i,3}$ , 于是得到$v_{j=3}(t=2)=0.021$ , 我们同样用红色实线在下图中标出。而且，还有$\psi_2(2)=1$，$\psi_3(2)=1$ 。



在接下来的步骤里，我们要确定 $v_{j=1}(t=3)$，$v_{j=2}(t=3)$，$v_{j=3}(t=3)$，然后在所有的值中，我们会发现$v_{j=2}(t=3)=0.01764$最大，如下图所示，用红色实线标出（当 $t=2$时，由$j=2,3$ 发出的路径未画出）。



到此为止，因为我们仅有三个观测值，因此可以结束我们的算法执行，从而知道当观测值序列为 $y_1=up, \ y_2=up, \ y_3=down$ 时，最可能的隐状态序列为 $q_1=1, \ q_2=1,\ q_3=2$（这个序列由$\psi$给 出），且概率为 $0.01764$。

---

HMM在自然语言处理中的应用举例

前面给出的例子都比较简单，不免让人产生一种错觉，以为HMM在实践中难当大任。然而，HMM在语音识别和自然语言处理中的用处非常广泛，而且准确率也相当之高（POS Tagging的准确率可以达到96%以上）。

例如我们现在有一个句子 “time flies like an arrow”. 然后我们可以为每个标注词性，则有“time/Noun flies/Verb like/Prep an/Art arrow/Noun/”。显然，当我们单独来看单词 flies，它可能是一个Noun，也可能是一个Verb，但是如果它前面那个词是形容词，那么它显然是名词的概率更高，而如果它前面的那个词是名词，则它是动词的可能性更高。所以我们就可以建立一个HMM模型，来根据观测到的状态序列（及历史）来预测具体某个词到底是什么词性。其中，单词的词性就是HMM中的隐状态，而我们所看到的具体的每个词就是观测状态。对此有兴趣的读者可以参阅文献[4]以了解HMM在NLP中的更多应用细节。

---

参考文献

[1] Stuart Russell and Peter Norvig. Artificial Intelligence: A Modern Approach. 3rd Edition.

[2] The University of Melbourne “Web Search and Text Analysis” 课程之幻灯片素材

[3] 同时推荐悉尼科大徐亦达博士的机器学习公开课中关于HMM的部分

[4] Daniel Jurafsky & James H. Martin. Speech and Language Processing, 3rd Edition. Chapter 6