HMM基本算法

Posted 2022-06-27 liuerdou

隐马尔科夫模型（Hidden Markov Model，以下简称HMM）作为语音信号的一种统计模型，在语音处理的各个领域中获得了广泛的应用。当然，随着目前深度学习的崛起，尤其是RNN，LSTM等神经网络序列模型的火热，HMM的地位有所下降。但是作为一个经典的模型，学习HMM的模型和对应算法，对我们解决问题建模的能力提高以及算法思路的拓展还是很好的。本篇将介绍HMM应用到语音处理中经常会面临的3大基本问题的解决方案，也就是给出三种基本算法：前向-后向算法、Baum_Welch算法以及Viterbi解码算法。

前向-后向算法

? 该算法用来计算给定一个观察值序列 \(O=o_1,o_2,\ldots,o_T\) 以及一个模型 \(\lambda=(\pi,A,B)\) 时，由模型 \(\lambda\) 计算出O的概率 $ P(O| \lambda ) $ 。
HMM组成如下图所示：

\(P(O|\lambda)\) 最直接的求取方式如下：

对一个固定的状态序列 \(Q=q_1,q_2,\ldots,q_T\) ,有
\[ P(O|Q,\lambda)=\mathop\Pi_t=1^TP(o_t|q_t,\lambda) = b_q_1(o_1)b_q_2(o_2) \ldots b_q_T(o_T)\tag1-1 \]
其中，
\[ b_q_t(o_t)=\left. b_jk \right|_q_t = \theta_j,o_t=V_k， 1 \leq t \leq T \tag1-2 \]
而对给定\(\lambda\)，产生Q的概率为
\[ P（Q|\lambda）= \pi_q_1a_q_1q_2 \cdots a_q_T-1q_T\tag1-3 \]
因此，所求概率为
\[ \begineqnarray* P(O|\lambda) & = & \sum_all QP(O|Q,\lambda)P(Q|\lambda) \& = & \sum_q_1,q_2,\ldots,q_T \pi_q_1b_q_1(o_1)a_q_T-1q_Tb_q_T(o_T) \endeqnarray*\tag1-4 \]
? 显而易见，上式的计算量是惊人的，大约为\(2TN^T\)数量级，当\(N=5,T=100\)时，计算量达\(10^72\)，这是完全不能接受的。在此情况下，要想求出\(P(O|\lambda)\)，就必须寻求更为有效的算法，前向后向算法就是解决这一问题的一种有效算法。
? 后面算法中为方便表示，对状态\(\theta_i\)的形式简记为\(i\)。

前向算法

定义前向变量为
\[ \alpha_t(i) = P(o_1o_2 \cdots o_t,q_t=i|\lambda)\tag1-5 \]
那么有，

• 初始化：对 \(1 \leq i \leq N\) ,有

\[ \alpha_1(j) = \pi_ib_i(o_1)\tag1-6 \]

• 递推：对 \(1 \leq t \leq T-1,1 \leq j \leq N\),有

\[ \alpha_t+1(j) = \left[\sum_i=1^N \alpha_t(i)a_ij\right]b_j(o_t+1)\tag1-7 \]

• 终止：

\[ P(O|\lambda)=\sum_i=1^N \alpha_T(i)\tag1-8 \]
式（1-7）中
\[ b_j(o_t+1)=\left. b_jk \right|_o_t+1=V_k\tag1-9 \]
这种算法计算量大为减少，只需要\(N^2T\)次运算。对\(N=5,T=100\)时，只需要大约3000次乘法计算。它是一种格形结构，图给出前向算法示意图：

后向算法

与前向算法类似，定义后向变量为
\[ \beta_t(i)=P(o_t+1o_t+2 \cdots o_T|q_t=i,\lambda)\tag1-10 \]
类似前向算法，有

• 初始化：对\(1 \le i \le N\)，有

\[ \beta_T(i)=1\tag1-11 \]

• 递推：对\(t=T-1,T-2,\cdots,1,1 \leq i \leq N\),有
  \[ \beta_t(i) = \sum_j=1^N a_ijb_j(o_t+1) \beta_t+1(j)\tag1-12 \]
• 终止：

\[ P(O|\lambda)=\pi_ib_i(o_1)\beta_1(i)\tag1-13 \]
后向算法的计算量大约也在\(N^2T\)数量级，它也是一种格形结构，HMM后向算法中计算后向变量\(\beta_t(i)\)的递推过程如下图所示。

Viterbi算法

这个算法解决对给定一个观察值序列\(O=o_1o_2 \cdots o_T\)和一个模型\(\lambda=(\pi,A,B)\)，如何确定一个最佳状态序列\(Q^*=q_1^*,q_2^*,\cdots,q_T^*\)的问题。
“最佳”的意义有很多种，由不同的定义可得到不同的结论，这里讨论的最佳意义上的状态序列\(Q^*\)是指使\(P(Q,O|\lambda)\)最大时确定的状态序列\(Q^*\)。这可用Viterbi算法来实现，其描述如下：
定义\(\delta_t(i)\)为时刻t时沿一条路径\(q_1,q_2,\cdots,q_t\)，且q_t=i,产生出\(o_1o_2\cdots o_t\)的最大概率，即有
\[ \delta_t(i)=\mathopmax_q_1,q_2,\cdots,q_t-1P(q_1q_2 \cdots q_t,q_t=i,o_1o_2\cdots o_t|\lambda)\tag2-1 \]
那么，求取最佳状态序列$Q^* $的过程为

• 初始化：对$1 \le i \le N $,有
  \[ \delta_t(i)=\pi_i b_i(o_t)\tag2-2 \]
  \[ \varphi_1(i)=0\tag2-3 \]
• 递推：对$2 \le t \le T，1 \le j \le N $，有
  \[ \delta_t(j)=\mathopmax_1 \le i \le N[\delta_t-1(i)a_ij]b_j(o_t)\tag2-4 \]
  \[ \varphi_t(j)=\mathopargmax_1 \le i \le N[\delta_t-1(i)a_ij]\tag2-5 \]
• 终止：
  \[ P^*=\mathopmax_1\le i \le N[\delta_T(i)]\tag2-6 \]
  \[ q_T^*=\mathopargmax_1 \le i \le N[\delta_T(i)]\tag2-7 \]
• 路径回溯，确定最佳状态序列：
  \[ q_t^*=\varphi_t+1(q_t+1^*),t=T-1,T-2,\cdots,1\tag2-8 \]
  对语音处理应用而言，\(P(Q,O|\lambda)\)动态范围很大，或者说不同的Q使\(P(Q,O|\lambda)\)的值差别很大，而\(\mathop\max_q P(Q,O|\lambda)\)事实上是\(\sum_qP(Q,O|\lambda)\)中举足轻重的唯一成分，因此常常等价的使用\(max_qP(Q,O|\lambda)\)近似\(\sum_qP(Q,O|\lambda)\),那么，Viterbi 算法也就能用来计算\(P(O|\lambda)\)。
  此外，上述的Viterbi算法也是一种格形结构，而且类似于前向算法。同样，由后向算法的思想出发，亦可推导出Viterbi算法的另一种实现方式。

Baum-Welch算法

这个算法用于解决HMM训练问题，即HMM参数估计问题。可描述为：给定一个观察虚列\(O=o_1,o_2,\cdots,o_T\)，确定一个\(\lambda=(\pi,A,B)\)，使\(P(O|\lambda)\)最大。
显然，由前向算法和后向算法定义的前向和后向变量，即式(1-5)、(1-10)，有
\[ P(O|\lambda)=\sum_i=1^N \sum_j=1^N \alpha_t(i) a_ij b_j(o_t+1) \beta_t+1(j),1 \le t \le T-1 \tag3-1 \]
求取\(\lambda\)使\(P(O|\lambda)\)最大是一个泛函极值问题。但由于给定的训练序列有限，因而不存在一个估计\(\lambda\)的最佳方法。在这种情况下，Baum-Welch算法利用递归的思想，使\(P(O|\lambda)\)局部极大，最后得到模型参数\(\lambda=(\pi,A,B)\)。此外，用梯度方法也可以达到类似的目的。
定义\(\xi_t(i,j)\)为给定训练序列O和模型\(\lambda\)时，HMM模型在t时刻处于i状态，t+1时刻处于j状态的概率，即
\[ \xi_t(i,j)=P(q_t=i,q_t+1=j|O,\lambda)\tag3-2 \]
可以推导出
\[ \xi_t(i,j)=[\alpha_t(i)a_ijb_j(o_t+1)\beta_t+1(j)]\tag3-3 \]
那么，HMM模型在时刻 t 处于 i 状态的概率为
\[ \gamma_t(i)=P(q_t=i|O,\lambda)=\sum_j=1^N \xi_t(i,j)=\alpha_t(i)\beta_t(i)/P(O|\lambda)\tag3-4 \]
因此，\(\mathop\sum_t=1^T-1 \gamma_t(i)\)表示从 i 状态转移出去次数的期望值 ，而\(\sum_t=1^T-1\xi_t(i,j)\)表示从 i 状态转移到状态 j 次数的期望值。由此，导出Baum-Welch算法中著名的重估公式，即
\[ \mathop\pi_i^-=\gamma_1(i)\tag3-5 \]

\[ \mathopa_ij^-=\mathop\sum_t=1^T-1\xi_t(i,j) / \mathop\sum_t=1^T-1\gamma_t(i)\tag3-6 \]
\[ \mathopb_jk^-=\mathop\sum_t=1,o_k=V_k^T\gamma_t(i) / \mathop\sum_t=1^T\gamma_t(i)\tag3-7 \]
那么，HMM参数\(\lambda=(\pi,A,B)\)的求取过程为，根据观察值序列O和选取的初始模型\(\lambda=(\pi,A,B)\)，由重估公式(3-5)、(3-6)和(3-7)求得一组新参数\(\mathop\pi_i^-\)，\(\mathopa_ij^-\)，\(\mathopb_jk^-\)，亦即得到了一个新模型\(\mathop\lambda^-=(\mathop\pi^-,\mathopA^-,\mathopB^-)\)。由下面的证明可以看出，\(P(O|\mathop\lambda^-)>P(O|\lambda)\)，即由重估公式得到的\(\mathop\lambda^-\)比\(\lambda\)在表示观察值序列O方面要好，那么重复这个过程，逐步改进模型参数，直到\(P(O|\mathop\lambda^-)\)收敛，即不再明显增大，此时的\(\mathop\lambda^-\)即为所求的模型。
应当指出，HMM训练，或者称为参数估计问题，是HMM在语音处理中应用的关键问题，与前面讨论的两个问题相比，这是最困难的问题。Baum-Welch算法只是得到广泛应用的解决这一问题的经典方法，但并不是唯一的，也远不是最完善的方法。

重估算法证明

Baum-Welch算法一般的证明方式是构造Q函数作为辅助函数，这和EM算法中的Q函数构造是一样的。这里简要介绍EM算法，这样在介绍Q函数构造原则的同时，也将证明Baum-Welch算法的有效性。

EM算法是一种从“不完全数据”中求解模型分布参数的最大似然估计方法。所谓“不完全数据”一般是指两种情况：一种情况是由于观测过程本身的限制或者错误造成观测数据成为有错漏的“不完全”数据；另一种情况是对参数的似然函数直接优化十分困难，而引入额外的参数（隐含的或丢失的）后优化就会比较容易。于是定义原始观测数据加上额外参数组成“完全数据”，原始观测数据自然就成为“不完全数据”。在模式识别及相关领域中，后一种情况更常见。由于所要优化的似然函数是一个随机变量，直接对其求最大化不好计算，但它的期望却是一个确定性的函数，优化起来相对容易。这就是EM算法的基本思路。该算法中包含两个主要方面：一是求期望，用E来表示，一是最大化，用M来表示；这也是这种方法名字的由来。EM算法在矢量优化和HMM模型的参数估计中有着很重要的作用。下面以HMM模型训练为例介绍EM算法。

在HMM模型中，数据是由两部分组成，一部分是可以观测到的数据，如观测特征序列O，称为可观测数据；另一部分数据无法观测到，如状态序列Q，称为隐含序列。由这两部分数据可以共同构成一个完全数据集（O,Q）。EM算法的目的是通过迭代的将完全数据集的对数似然度期望最大化，来实现对可观测数据的对数似然度的最大化。

根据Bayes公式，完全数据集的似然度和可观测数据集的似然度之间存在以下关系：
\[ P(O,Q|\lambda)=P(Q|O,\lambda)P(O|\lambda)\tag4-1 \]
那么，观测数据的对数似然度可以表示为
\[ logP(O|\lambda)=logP(O,Q|\lambda)-logP(Q|O,\lambda)\tag4-2 \]
对于两个参数集\(\lambda\)和\(\mathop\lambda^-\)，在已知O和\(\lambda\)的情况下，对公式（4-2）在完全数据集上求期望，则
\[ E[logP(o|\mathop\lambda^-)|O,\lambda]=E[logP(O,Q|\mathop\lambda^-)|O,\lambda]-E[logP(Q|O,\mathop\lambda^-)|O,\lambda]\tag4-3 \]
其中式（4-3）的左半部分为
\[ \begineqnarray* E[logP(o|\mathop\lambda^-)|O,\lambda] & = & \int logP(O|\mathop\lambda^-)P(Q|O,\lambda)dQ \& = & logP(O|\mathop\lambda^-) \& = & L(O,\mathop\lambda^-)\tag4-4 \endeqnarray* \]
令
\[ Q(\lambda,\mathop\lambda^-)=E[logP(O,Q|\mathop\lambda^-)|O,\lambda]=\int logP(O,Q|\mathop\lambda^-)P(Q|O,\lambda)dQ\tag4-5 \]
而
\[ H(\lambda,\mathop\lambda^-)=E[logP(Q|O,\mathop\lambda^-)|O,\lambda]=\int logP(Q|\mathop\lambda^-)P(Q|O,\lambda)dQ\tag4-6 \]
则式（4-4）变为
\[ L(O,\mathop\lambda^-)=Q(\lambda,\mathop\lambda^-)-H(\lambda,\mathop\lambda^-)\tag4-7 \]
由Jensen不等式可以知道，\(H(\lambda,\mathop\lambda^-) \le H(\lambda,\lambda)\)。EM算法的基本原理在于，如果可以保证\(Q(\lambda,\mathop\lambda^-) \ge Q(\lambda,\lambda)\)不等式成立，那么不等式\(L(O,\mathop\lambda^-) \ge L(O,\lambda)\)一定成立。从上面分析可以看出，Q()函数实质上就是完全数据集的对数似然度的期望，这样，通过将Q()函数最大化就可以实现观测数据的对数似然度的最大化。
当隐藏数据为离散时，Q()函数一般表示为
\[ Q(\lambda,\mathop\lambda^-)=\sum_所有Q \frac P(O,Q|\lambda)P(O|\lambda) logP(O,Q|\mathop\lambda^-)\tag4-8 \]
EM算法一般描述——给定一个当前的参数集\(\lambda\)，可以通过如下方式获得新的参数集\(\mathop\lambda^-\):

• 选择初始参数$ \lambda $；
• 求期望，即在给定的参数集$ \lambda $上求Q()函数；
• 最大化，选择$ \lambda=\mathopargmax_\mathop\lambda^-Q(\lambda,\mathop\lambda^-)$。