EM Algorithm in speech progressing

前言

最近准备东大IME项目面试，正在读Automatic Speech Recognition，在HMM的部分中遇到了一些比较复杂的数学过程，这里简单写一些笔记。

背景知识

前向概率:

$\alpha_t(i)=P(q_t=i,o^t_1), t=1,...,T$

后向概率:

$\beta_t(i)=P(o^T_{t+1}|q_t=i), t=1,...,T-1$

其中$o_a^b$表示观测序列$(o_a,o_{a+1},...,o_b)$

利用前向和后向概率可以递推的计算马尔科夫链中的$P(o_1^T)$，否则几乎其时间复杂度几乎不可接受。

其递推式为分别为

$\alpha_t(j)=\sum^N_{i=1}\alpha_{t-1}(i)a_{ij}b_j(o_t), t=2,3,...,T; j=1,2,...,N$

$\beta_t(i)=\sum^N_{j=1}\beta_{t+1}(j)a_{ij}b_j(o_{t+1}), t=T-1, T-2,...; i=1,2,...,N$

另外，由定义我们不难得到初值。对于$\alpha$，有

$\alpha_1(i) = P(q_1=i,o_1)=P(q_1=i)P(o_1|q_1)=\pi_ib_i(o_1)$

对于$\beta$,因为在定义域之外，我们选择$\beta_T(i)=1$

其中，$b_i(o_1)$表示当前状态为$i$时观测结果为$o_1$的概率，而$\pi_i$表示马尔科夫链的稳定分布(stationary distribution)，即当$t\to \infty$时$p_i(t)$的值

而定义了前向和后向概率后，我们可以将$P(o_1^T)$分解，有

$\begin{aligned} P(q_t&=i,o_1^T)=P(q_t=i,o_1^t,o_{t+1}^T)\\ &=P(q_t=i,o_1^t)P(q_t=i,o^T_{t+1}|o_1^t,q_t=i)\\ &=P(q_t=i,o_1^t)P(q_t=i,o^T_{t+1},q_t=i)\\ &=\alpha_t(i)\beta_t(i) \end{aligned}$

注意到第二行到第三行成立是因为在特定状态下观测是独立同分布(IID)的。

由此有

$P(o_1^T)=\sum^N_{i=1}P(q_t=i,o^T_1)=\sum^N_{i=1}\alpha_t(i)\beta_t(i)$

注意到这里的t，即$\alpha$和$\beta$的分断点是任意选择的，那么我们只要令$t=T$就有

$P(o_1^T)=\sum^N_{i=1}\alpha_T(i)$

EM算法

背景

考虑全数据$y={o,h}$，其中$o$是观测到的数据，而$h$是隐藏的变量。那么要找到参数$\theta$，我们需要最大化其对数似然函数$log p(o;\theta)$，但是这个概率分布函数是相对难以计算的，因此EM算法试图迭代地处理全数据$y$，再通过建立$y$到$o$的映射$o=g(y)$来解决问题。而全数据的选择对于不同的问题通常具有其独特性。

而通过合理的选择隐数据$h$,可以配合观测数据$o$来构造全数据$y$，这样可以简化似然估计的过程。

考虑以下条件期望

$Q(\theta|\theta_0)=E_{h|o}[\log p(y;\theta)o;\theta_0] = E[\log p(o,h;\theta)|o;\theta_0]$

对$h$在其所有的可能性上求和(如果是连续的则为积分)就可以得到似然函数的值。注意到这个似然函数的值是与$\theta_0$相关的，相当于在初值$\theta_0$的基础上计算得到了一个新的$\theta$,而可以保证$\theta$是收敛的，即$$Q(\theta|\theta_{k+1})\ge Q(\theta|\theta_k)$$，仅在$\theta_k$已是一个最大似然估计时取到等号。

Baum-Welch 算法

如在上一节中所述，对于EM算法我们需要一个观测数据和一个隐藏数据来构造全数据，对于HMM来说，分别为观测序列$o^T_1$和隐马尔科夫链状态序列$q_1^T$，代入公式则有

$Q(\theta|\theta_0)= E[\log p(o^T_1,q^T_1;\theta)|o^T_1;\theta_0]$

E步骤

在E步骤中，需要将$Q(\theta|\theta_0)$化简到便于在M步骤中最大化的形式。将期望展开则有

$\begin{aligned} Q(\theta|\theta_0)&= E[\log p(o^T_1,q^T_1;\theta)|o^T_1;\theta_0]\\ &=\sum_{q_1^T}P(q_1^T|o_1^T,\theta_0)\log P(o_1^T,q_1^T|\theta) \end{aligned}$

，其中$\theta$和$\theta_0$表示当次和上一次迭代中的参数

由于在高斯混合HMM中，$P(o_1^T|q_1^T)$服从高斯分布，但其形式过于复杂，为了方便简写，我们令

$N_t(i) = -\frac{D}{2}\log(2\pi)-\frac12\log|\Sigma_i|-\frac12(o_t-\mu_i)^\mathrm{T}\Sigma^-1_i(o_t-\mu_i)$

则

$\log P(o_1^T|q_1^T) = \sum_{t=1}^T N_t(q_t)$

同时，由转移概率$a_{ij}$的定义可以得到$P(q_1^T)=\prod_{t=1}^{T-1}a_{q_tq_{t+1}}$，因此有

$\log P(o_1^T,q_1^T|\theta) = \log P(o_1^T|q_1^T)P(q_1^T) = \sum_{t=1}^T N_t(q_t) + \sum_{t=1}^{T-1}\log a_{q_tq_{t+1}}$

回代到$Q$的定义式中有

$Q(\theta|\theta_0)=\sum_{q_1^T}P(q_1^T|o_1^T,\theta_0)\sum_{t=1}^T N_t(q_t) + \sum_{q_1^T}P(q_1^T|o_1^T,\theta_0)\sum_{t=1}^{T-1}\log a_{q_tq_{t+1}}$

首先引入克罗内克$\delta$函数(Kronecker delta)

$\delta_{ij} = \begin{cases} 1 & (i=j),\\ 0 & (i\neq j) \end{cases}$

由于$\delta$仅在一个点上取到非$0$值，因此乘上它并不会改变函数的值。所以可以利用构造一个合适的$\delta$函数并调换乘项位置来化简$Q$.

在本例中，对$Q$的第一项有

$\begin{aligned} Q_1(\theta|\theta_0)&=\sum^N_{i=1}\left\{\sum_{q^T_1}P(q_1^T|o_1^T,\theta_0)\sum^T_{t=1}N_t(q_t) \right\}\delta_{q_t,i}\\ &=\sum^N_{i=1}\left\{\sum_{q^T_1}P(q_1^T|o_1^T,\theta_0)\delta_{q_t,i}\sum^T_{t=1}N_t(q_t) \right\}\\ &=\sum^N_{i=1}\sum^T_{t=1}P(q_t=i|o_1^T,\theta_0)N_t(q_t)\\ \end{aligned}$

同理，对$Q$的第二项乘上$\delta_{q_t,i}\delta_{q_{t+1},j}$可得

$Q_2(\theta|\theta_0)=\sum^N_{i=1}\sum^N_{j=1}\sum^{T-1}_{t=1}P(q_t=i,q_{t+1}=j|o_1^T,\theta_0)\log a_{ij}$

由此，不难注意到Q的两项分别只与高斯分布($N_t(I)$)和马尔科夫链($a_{ij}$)有关，因此可以独立的最大化这两项。

令$Q_1$和$Q_2$中连加时使用的权分别为

$\begin{aligned} \gamma_t(i) &= P(q_t=i|o_1^T,\theta_0)\\ \xi_t(i,j) &= P(q_t=i, q_{t+1}=j|o_1^T,\theta_0) \end{aligned}$

则由其定义，不难得出

$\xi_t(i,j) = \frac{\alpha_t(i)\beta_{t+1}(j)a_{ij}\exp(N_{t+1}(j))}{P(o_1^T|\theta_0)}, t=1,2,...,T-1.$

由定义式，显然的，对$\xi_t$在j上求和就可以得到$\gamma$

$\gamma_t(i)=\sum^N_{j=1}\xi_t(i,j), t=1,2,...,T-1.$

又由定义可以直接得到

$\gamma_T(i)=P(q_T=i|o_1^T,\theta_0) = \frac{P(q_t=i,o_1^T|\theta_0)}{P(o_1^T|\theta_0)}=\frac{\alpha_T(i)}{P(o_1^T|\theta_0)}$

由此将$Q$化简为了只要对状态$i$或状态对$i,j$遍历的便于计算的形式，以在M步骤中处理。

M步骤

在M步骤中需要将E步骤中得到的$Q_1$和$Q_2$的值在$\theta_0$下最大化。对于$Q_2$，令$\frac{\partial Q_2}{\partial a_{ij}}=0$。由于有边界条件$\sum^N_{j=1}a_{ij}=1$，使用标准拉格朗日乘数法可以得到

$\hat a=\frac{\sum^{T-1}_{t=1}\xi_t(i,j)}{\sum^{T-1}_{t=1}\gamma_t(i)}$

具体计算步骤略去。

而对$Q_1$，我们省掉与优化无关的常数项和常数系数后可以得到等价的目标函数

$Q_1(\mu_i,\Sigma_i)=\sum^N_{i=1}\sum^{Tr}_{t=1}\gamma_t(i)(o_t-\mu_i)^\mathrm{T}\Sigma^-1_i(o_t-\mu_i)-\frac12\log|\Sigma_i|$

省略的部分包括高斯分布中与$\Sigma$无关的项，以及最后一项的系数(显然的有$\sum^N_{i=1}\gamma_t(i)=1$,那么再过一次求和之后仍然是常数，求导之后会被消去)

则相等于解

$\frac{\partial Q_1}{\partial\Sigma_i}=0, i=1,2,...,N.$

令$K=\Sigma^-1$，则对$\log|\Sigma_i|$项，$k_{lm}$的导数就是$\Sigma$中对应的项，即$\sigma_{lm}$，于是单项的导数$\frac{\partial Q_1}{\partial k_{lm}}=0$,即

$\sum^T_{t=1}\gamma_t(i)\left\{\frac12\sigma_{lm}-\frac12(o_t-\mu_i)_l(o_t-\mu_i)_m\right\}=0$

将$\sigma_{lm}$解出并写成矩阵形式有

$\hat\Sigma_i = \frac{\sum^T_{t=1}\gamma_t(i)(o_t-\hat\mu_i)(o_t-\hat\mu_i)^\mathrm{T}}{\sum^T_{t=1}\gamma_t(i)}$

从而有

$\hat\mu_i=\frac{\sum^T_{t=1}\gamma_t(i)o_t}{\sum^T_{t=1}\gamma(i)}$

Viterbi算法

Viterbi算法用于解码HMM状态，即，已知一个观测$o_1^T$，求出最可能的HMM状态序列$q_1^T$

其核心思想是动态规划，由于这个算法在OI中相当常见就不再赘述其背景了。

在Viterbi算法中，我们要最大化的量是

$\delta_i(t) = \mathop{max}\limits_{q1,q2,...,q_{t-1}} P(o_1^t,q_1^{t-1},q_t=i)$

也就是对于观测序列$o^t_1$,假设t时刻之前的HMM状态序列为$q_1^{t-1}$的情况下，$t$时刻状态为$i$的概率。

而有了$t$时刻的状态i的目标函数则可以简单的得到下一时刻状态j下的目标函数表达式，即

$\delta_j(t+1) = \mathop{max}\limits_{i}\delta_i(t)a_{ij}b_j(o_{t+1})$

这个式子的三项分别为：转移前的状态i的概率，从前状态i转移到该状态j的概率，以及对于相应的观测$o_{t+1}$，该时刻处于状态j的概率。

由此迭代到T时刻就可以得到最大概率$P^*$，以及所有时刻的状态路径$q^*(t)$

如果曾经有做过DP相关的题目这部分应该不难理解。

至此，利用EM算法和Viterbi算法，我们就完成了对HMM参数的估计以及隐参数HMM状态序列的解码。