机器学习笔记之隐马尔可夫模型(五)学习问题——EM算法

引言

前面两节分别介绍了使用前向算法(Forward Algorithm)和后向算法(Backward Alogrithm)对求值问题(Evaluation)进行求解。本节将介绍学习问题(Learning)并使用EM算法进行求解。

关于学习问题的介绍

学习问题(Learning)我们并不陌生，其本质上是 给定数据集合$\mathcal{X}$，通过$\mathcal{X}$求解模型的参数信息$\theta$：
之前介绍的模型中：

• 极大似然估计(Maximum Likelihood Estimate,MLE)。如：线性回归(Linear Regression)、一些线性分类(Linear Classification)算法如线性判别分析(Linear Discriminant Analysis,LDA)、逻辑回归(Logistic Regression)等等。
  $${\hat{\theta }}_{MLE}=\underset{\theta }{\mathrm{argmax}}P(\mathcal{X}\mid \theta )$$
• 狭义EM算法(Expectation Maximization,EM)。暂时仅介绍了高斯混合模型(Gaussain Mixture Model,GMM)一种。
  $${\theta }^{(t+1)}=\underset{\theta }{\mathrm{argmax}}{\int }_{\mathcal{Z}}P(\mathcal{X},\mathcal{Z}\mid \theta )\cdot P(\mathcal{Z}\mid \mathcal{X},{\theta }^{(t)})d\mathcal{Z}$$

本节介绍的隐马尔可夫模型同样可以采用狭义EM算法进行求解。

场景介绍

概念的数学符号表示

针对隐马尔可夫模型，对场景与相关数学符号进行介绍：

• 隐马尔可夫模型中的变量由状态序列$\mathcal{I}$与观测序列$\mathcal{O}$构成。假设观测序列长度为$T$，状态序列与观测序列分别表示如下：
  I = { i 1 , i 2 , ⋯   , i T } O = { o 1 , o 2 , ⋯   , o T } \mathcal I = \{i_1,i_2,\cdots,i_T\} \\ \mathcal O = \{o_1,o_2,\cdots,o_T\} I={i1​,i2​,⋯,iT​}O={o1​,o2​,⋯,oT​}
• 隐马尔可夫模型中状态序列中的 任意状态变量$i_t\in \mathcal{I}$的取值结果均是离散的。即：
  ∀ i t ∈ I , i t ∈ Q = { q 1 , q 2 , ⋯   , q K } \forall i_t \in \mathcal I, \quad i_t \in \mathcal Q = \{q_1,q_2,\cdots,q_{\mathcal K}\} ∀it​∈I,it​∈Q={q1​,q2​,⋯,qK​}
  我们同样定义观测序列中任意观测变量$o_t\in \mathcal{O}$的取值结果是离散的。即：
  ∀ o t ∈ O , o t ∈ V = { v 1 , v 2 , ⋯   , v M } \forall o_t \in \mathcal O,\quad o_t \in \mathcal V = \{v_1,v_2,\cdots,v_{\mathcal M}\} ∀ot​∈O,ot​∈V={v1​,v2​,⋯,vM​}
  我们称$\mathcal{Q}$为状态值集合，称$\mathcal{V}$为观测值集合。

模型参数的数学符号表示

隐马尔可夫模型的参数记作$\lambda$，模型参数共包含三个部分：

• 初始概率分布$\pi$：
  即 状态变量$i_1$ 选择任意状态值$q_i(i\in 1,2,\cdots ,\mathcal{K})$的概率结果组成的向量形式。即：
  $$\pi ={\left[\begin{array}{c}P(i_1=q_1)\\ P(i_1=q_2)\\ ⋮\\ P(i_1=q_{\mathcal{K}})\end{array}\right]}_{\mathcal{K}\times 1}\sum _{k=1}^{\mathcal{K}}P(i_1=q_k)=1$$
• 状态转移矩阵$\mathcal{A}$：矩阵$\mathcal{A}$中的任意一个元素$a_{ij}$表示当前时刻的状态变量$i_t=q_i$的条件下，下一时刻的状态变量$i_{t+1}=q_j$的后验概率结果。即：
  $$\mathcal{A}=[a_{ij}{]}_{\mathcal{K}\times \mathcal{K}},a_{ij}=P(i_{t+1}=q_j\mid i_t=q_i)$$
• 发射矩阵$\mathcal{B}$：矩阵中的任意一个元素$b_j(k)$表示当前时刻的状态变量$i_t=q_i$的条件下，对应时刻的观测变量$o_t=v_k$的后验概率结果。即：
  $$\mathcal{B}=[b_j(k){]}_{\mathcal{K}\times \mathcal{M}},b_j(k)=P(o_t=v_k\mid i_t=q_j)$$

HMM模型的特殊性质

• 齐次马尔可夫假设：任一时刻状态变量$i_t$的后验概率，只和其前一时刻的状态变量$i_{t-1}$相关，与其他变量无关。即：
  $$P(i_t\mid i_{t-1},\cdots ,i_1,o_{t-1},\cdots ,o_1)=P(i_t\mid i_{t-1})$$
• 观测独立性假设：任一时刻观测变量$o_t$的后验概率，只和对应时刻的状态变量$i_t$相关，与其他变量无关。
  $$P(o_t\mid i_t,i_{t-1},\cdots ,i_1,o_{t-1},\cdots ,o_1)=P(o_t\mid i_t)$$

模型参数$\lambda$求解

我们在前向算法(Forward Algorithm)中介绍过，如果直接求解$P(\mathcal{O}\mid \lambda )$，它的表达式结果如下：
$$\begin{array}{rl}P(\mathcal{O}\mid \lambda )& =\sum _{\mathcal{I}}P(\mathcal{O},\mathcal{I}\mid \lambda )\\ & =\sum _{i_1}\cdots \sum _{i_T}\left[\pi \cdot \prod _{t=2}^T{a}_{i_t,i_{t+1}}\cdot \prod _{t=1}^T{b}_{i_t}(o_t)\right]\end{array}$$
如果对$P(\mathcal{O}\mid \lambda )$直接使用极大似然估计，即：
$${\hat{\lambda }}_{MLE}=\underset{\lambda }{\mathrm{argmax}}P(\mathcal{O}\mid \lambda )$$
直接求解的方式存在许多问题：

• 首先，观测序列$\mathcal{O}$中的个观测变量$o_1,\cdots ,o_T$之间不是独立同分布，不能使用$\underset{\lambda }{\mathrm{argmax}}{\prod }_{i=1}^{T}P(o_i\mid \lambda )$进行求解；
• 并且直接求解$P(\mathcal{O}\mid \lambda )$的时间复杂度是$O({\mathcal{K}}^2\times T)$，因此求解过程十分复杂。本节将介绍使用EM算法通过迭代的方式求解模型参数$\lambda$。

EM算法求解模型参数

E步操作

EM算法公式表示如下：
$${\theta }^{(t+1)}=\underset{\theta }{\mathrm{argmax}}{\int }_{\mathcal{Z}}P(\mathcal{X},\mathcal{Z}\mid \theta )\cdot P(\mathcal{Z}\mid \mathcal{X},{\theta }^{(t)})d\mathcal{Z}$$
将EM算法中出现的变量与隐马尔可夫模型中出现的概念进行映射：

• $\mathcal{X}$：观测数据(Observed Data) $\to$ 观测序列$\mathcal{O}$；
• $\mathcal{Z}$：隐变量(Latent Data) $\to$ 状态序列$\mathcal{I}$；
• $\theta$：参数(Parameter) $\to$ 模型参数$\lambda$；

经过映射后，基于隐马尔可夫模型的EM算法表示如下：
'状态序列'$\mathcal{I}$是离散的，因此积分过程中改成$\sum$;
$${\lambda }^{(t+1)}=\underset{\lambda }{\mathrm{argmax}}\sum _{\mathcal{I}}\left[\log P(\mathcal{O},\mid \lambda )\cdot P(\mathcal{I}\mid \mathcal{O},{\lambda }^{(t)})\right]$$
为了简化运算，对上述公式进行变形：
条件概率公式~
$${\lambda }^{(t+1)}=\underset{\lambda }{\mathrm{argmax}}\sum _{\mathcal{I}}\left[\log P(\mathcal{O},\mid \lambda )\cdot \frac{P(\mathcal{I},\mathcal{O}\mid {\lambda }^{(t)})}{P(\mathcal{O},{\lambda }^{(t)})}\right]$$

观察后项的分母部分$P(\mathcal{O},{\lambda }^{(t)})$：

• ${\lambda }^{(t)}$是上一次迭代步骤的参数结果，是已知项；
• 观测序列$\mathcal{O}$(就是样本数据)，它和$\lambda$的取值无关；

至此，EM算法可化简为如下形式：
需要注意的点：化简后的结果已经不是期望形式了~
$$\begin{array}{rl}{\lambda }^{(t+1)}& =\underset{\lambda }{\mathrm{argmax}}\sum _{\mathcal{I}}\left[\log P(\mathcal{I},\mathcal{O}\mid \lambda )\cdot P(\mathcal{I},\mathcal{O}\mid {\lambda }^{(t)})\right]\end{array}$$
并且${\lambda }^{(t)}$以及迭代求解后的${\lambda }^{(t+1)}$本身不是一个参数，而是由三个参数组成：
$${\lambda }^{(t)}=({\pi }^{(t)},{\mathcal{A}}^{(t)},{\mathcal{B}}^{(t)});{\lambda }^{(t+1)}=({\pi }^{(t+1)},{\mathcal{A}}^{(t+1)},{\mathcal{B}}^{(t+1)})$$
将EM算法的公式部分表示为关于$\lambda ,{\lambda }^{(t)}$的函数：
$$\mathcal{Q}(\lambda ,{\lambda }^{(t)})=\sum _{\mathcal{I}}\left[\log P(\mathcal{I},\mathcal{O}\mid \lambda )\cdot P(\mathcal{I},\mathcal{O}\mid {\lambda }^{(t)})\right]$$
将直接求解得到的$P(\mathcal{O},\mathcal{I}\mid \lambda )=\pi \cdot {\prod }_{t=2}^T{a}_{i_t,i_{t+1}}\cdot {\prod }_{t=1}^T{b}_{i_t}(o_t)$带入$\mathcal{Q}(\lambda ,{\lambda }^{(t)})$中：
$$\begin{array}{rl}\mathcal{Q}(\lambda ,{\lambda }^{(t)})& =\sum _{\mathcal{I}}\left[\log \left(\pi \cdot \prod _{t=2}^T{a}_{i_t,i_{t+1}}\cdot \prod _{t=1}^T{b}_{i_t}(o_t)\right)\cdot P(\mathcal{I},\mathcal{O}\mid {\lambda }^{(t)})\right]\\ & =\sum _{\mathcal{I}}\left[\left(\log \pi +\log \prod _{t=2}^T{a}_{i_t,i_{t+1}}+\log \prod _{t=1}^T{b}_{i_t}(o_t)\right)\cdot P(\mathcal{I},\mathcal{O}\mid {\lambda }^{(t)})\right]\\ & =\sum _{\mathcal{I}}\left[\left(\log \pi +\sum _{t=2}^T\log a_{i_t,i_{t+1}}+\sum _{t=1}^T\log b_{i_t}(o_t)\right)\cdot P(\mathcal{I},\mathcal{O}\mid {\lambda }^{(t)})\right]\end{array}$$
这里以求解${\pi }^{(t+1)}$为例，进行求解。
观察上式中的小括号部分，只有$\log \pi$和$\pi$有关，而剩余两项均和$\pi$无关联，看做常数。因此，针对求解${\pi }^{(t+1)}$的$\mathcal{Q}(\lambda ,{\lambda }^{(t)})$表达如下：
$$\begin{array}{rl}{\pi }^{(t+1)}& =\underset{\pi }{\mathrm{argmax}}\mathcal{Q}(\lambda ,{\lambda }^{(t)})\\ & =\underset{\pi }{\mathrm{argmax}}\sum _{\mathcal{I}}[\log \pi \cdot P(\mathcal{O},\mathcal{I}\mid {\lambda }^{(t)})]\\ & =\underset{\pi }{\mathrm{argmax}}\sum _{i_1}\cdots \sum _{i_T}[\log \pi \cdot P(\mathcal{O},i_1,\cdots ,i_T\mid {\lambda }^{(t)})]\end{array}$$
观察上述展开结果，$\pi$根据定义，状态序列$\mathcal{I}$ 的初始概率分布$\pi$只和 第一个状态变量$i_1$相关，和其他变量无关。因此，将上式继续化简成如下形式：
$i_1$的取值范围是离散的，是状态值集合$\mathcal{Q}$中的一个元素。因此将${\sum }_{i_1}$改为${\sum }_{k=1}^{\mathcal{K}}$。
$$\begin{array}{rl}{\pi }^{(t+1)}& =\underset{\pi }{\mathrm{argmax}}\sum _{i_1}[\log \pi \cdot P(\mathcal{O},i_1\mid {\lambda }^{(t)})]\\ & =\underset{\pi }{\mathrm{argmax}}\sum _{k=1}^{\mathcal{K}}[\log P(i_1=q_k)\cdot P(\mathcal{O},i_1=q_k\mid {\lambda }^{(t)})]\end{array}$$
需要注意的是，将$\pi$写成概率组成向量的形式，是存在约束条件的。即：
$$\sum _{k=1}^{\mathcal{K}}P(i_1=q_k)=1$$
至此，将${\pi }^{(t+1)}$的求问题转化为 带一个约束的优化问题：
$$\{\begin{array}{l}\underset{\pi }{\mathrm{argmax}}{\sum }_{k=1}^{\mathcal{K}}[\log P(i_1=q_k)\cdot P(\mathcal{O},i_1=q_k\mid {\lambda }^{(t)})]\\ s.t.{\sum }_{k=1}^{\mathcal{K}}P(i_1=q_k)=1\end{array}$$

M步操作

使用拉格朗日乘数法求解该问题：
定义拉格朗日函数$\mathcal{L}(\pi ,\eta )$表示如下：
${\sum }_{k=1}^{\mathcal{K}}P(i_1=q_k)$和$1$之间怎么去减，需要视情况而定。
$$\mathcal{L}(\pi ,\eta )=\sum _{k=1}^{\mathcal{K}}[\log P(i_1=q_k)\cdot P(\mathcal{O},i_1=q_k\mid {\lambda }^{(t)})]+\eta \left(\sum _{k=1}^{\mathcal{K}}P(i_1=q_k)-1\right)$$
令${\pi }_k=P(i_1=q_k)$，即$\pi$向量中的其中一个分量。令$\mathcal{L}(\pi ,\eta )$对${\pi }_k$求偏导：
$$\frac{\partial \mathcal{L}}{\partial {\pi }_k}=\frac{1}{{\pi }_k}P(\mathcal{O},i_1=q_k\mid {\lambda }^{(t)})+\eta (1-0)$$
令$\frac{\partial \mathcal{L}}{\partial {\pi }_k}\triangleq 0$，有：
$$P(\mathcal{O},i_1=q_k\mid {\lambda }^{(t)})+{\pi }_k\cdot \eta =0$$
基于上式，则有：
$$\sum _{k=1}^{\mathcal{K}}\left[P(\mathcal{O},i_1=q_k\mid {\lambda }^{(t)})+{\pi }_k\cdot \eta \right]=0$$
该式可从两种角度解释：

• 从积分角度解释：上式左右两端均对$i_1$求积分：常数0的积分依然是常数。
• 从偏导角度解释：对$i_1$可选择的每一种可能对应的概率求偏导，并令其均为0。则有：
  $$\sum _{k=1}^{\mathcal{K}}\left[P(\mathcal{O},i_1=q_k\mid {\lambda }^{(t)})+{\pi }_k\cdot \eta \right]=0+0+\cdots +0=0$$

将${\sum }_{k=1}^{\mathcal{K}}\left[P(\mathcal{O},i_1=q_k\mid {\lambda }^{(t)})+{\pi }_k\cdot \eta \right]=0$展开：
边缘概率分布~
约束条件~

• ${\sum }_{k=1}^{\mathcal{K}}P(\mathcal{O},i_1=q_k\mid {\lambda }^{(t)})={\sum }_{i_1}P(\mathcal{O},i_1=q_k\mid {\lambda }^{(t)})=P(\mathcal{O}\mid {\lambda }^{(t)})$
• ${\sum }_{k=1}^{\mathcal{K}}{\pi }_k\cdot \eta =\eta \cdot \left({\sum }_{k=1}^{\mathcal{K}}{\pi }_k\right)=\eta \cdot 1=\eta$

最终结果有：
$$\begin{gathered} P(\mathcal{O}\mid {\lambda }^{(t)})+\eta =0 \\ \to \eta =-P(\mathcal{O}\mid {\lambda }^{(t)}) \end{gathered}$$

将$\eta =-P(\mathcal{O}\mid {\lambda }^{(t)})$带回$P(\mathcal{O},i_1=q_k\mid {\lambda }^{(t)})+{\pi }_k\cdot \eta =0$中，求得$p{i}_k$的最终结果为：
$${\pi }_k^{(t+1)}=\hat{P}(i_1=q_k)=\frac{P(\mathcal{O},i_1=q_k\mid {\lambda }^{(t)})}{P(\mathcal{O}\mid {\lambda }^{(t)})}$$

同理，可以求解出其他分量的迭代结果，从而得到迭代后的初始概率分布${\pi }^{(t+1)}$：
$${\pi }^{(t+1)}=({\pi }_1^{(t+1)},{\pi }_2^{(t+1)},\cdots ,{\pi }_{\mathcal{K}}^{(t+1)}{)}^T$$
此时，就已经将${\pi }^{(t+1)}$和${\lambda }^{(t)}$的迭代方式找出，同理${\mathcal{A}}^{(t+1)},{\mathcal{B}}^{(t+1)}$也使用相似方式。

下一节将介绍解码问题(Decoding)。
相关参考：
机器学习-隐马尔可夫模型5-Learning问题-Baum Welch算法(EM)