解决隐马尔科夫模型的学习问题

1. 问题背景与定位

在初识隐马尔可夫模型和隐马尔科夫模型的概率计算算法——解决评估问题中，我们已经了解了 HMM 是什么（模型定义），以及如何计算一个观测序列的概率（评估问题）。现在，我们面临 HMM 三大问题中最核心、也最困难的一个：学习问题 (Learning Problem)。

学习问题描述：在很多现实场景中，我们事先并不知道 HMM 的模型参数 $\lambda =(A,B,\pi )$。我们拥有的，仅仅是一个或多个观测序列 $O$。学习问题的任务就是：如何仅根据给定的观测序列 $O$，来估计出一套最优的模型参数 $\lambda$，使得这套参数能够最好地“解释”我们观测到的数据？

这个问题就好比，我们只知道朋友连续多天的活动记录（观测序列），但对他的城市天气转换规律（状态转移矩阵 $A$）、不同天气下他进行各项活动的偏好（观测概率矩阵 $B$）以及天气初始状态的概率（初始分布 $\pi$）一无所知。我们能从他的活动记录中反推出这些规律吗？

1.1 “有监督”与“无监督”的对比

要理解这个问题的难度，我们可以先从简单的情景即“有监督学习”开始分析。

如果数据是“完整”的（有监督）： 假设我们不仅有观测序列 $O=(o_1,o_2,\dots ,o_T)$，还幸运地拥有与之对应的隐藏状态序列 $Q=(q_1,q_2,\dots ,q_T)$。在这种情况下，估计参数 $\lambda$ 将会变得非常简单，只需要做一些频率统计（即最大似然估计）即可：

• 估计 ${\hat{\pi }}_i$: 初始状态为 $s_i$ 的概率，就是“所有样本中以状态 $s_i$ 开头的次数”除以“总样本数”。
• 估计 ${\hat{a}}_{ij}$: 从状态 $s_i$ 转移到 $s_j$ 的概率，就是“所有样本中从 $s_i$ 转移到 $s_j$ 的次数”除以“所有样本中从 $s_i$ 出发的总次数”。
• 估计 ${\hat{b}}_j(k)$: 在状态 $s_j$ 下观测到 $v_k$ 的概率，就是“所有样本中在状态 $s_j$ 时观测到 $v_k$ 的次数”除以“所有样本中处于状态 $s_j$ 的总次数”。

然而，现实是残酷的（无监督）： 在 HMM 的典型应用中，状态序列 $Q$ 正是那个“隐藏”的变量，我们无法得知。我们只有观测序列 $O$。这就形成了一个“鸡生蛋，蛋生鸡”的困境：

• 如果我们知道模型参数 $\lambda$，我们就能推断隐藏状态（例如，使用维特比算法）。
• 如果我们知道隐藏状态，我们就能通过计数来估计模型参数 $\lambda$。

现在我们两者都不知道，该如何打破这个循环呢？这正是鲍姆-韦尔奇算法要解决的核心矛盾。

2. 鲍姆-韦尔奇算法的核心思想

鲍姆-韦尔奇算法的本质，是此前我们接触的期望最大化 (Expectation-Maximization, EM) 算法（初识EM算法）在 HMM 上的一个具体实现。EM 算法是解决这类含有“隐变量”的最大似然估计问题的通用思想。

我们可以用一个非常直观的“猜测-精炼”策略来理解它：

1. 初始猜测 (Initialization): 我什么都不知道，那就先随机猜测一套模型参数 $\bar{\lambda }=(\bar{A},\bar{B},\bar{\pi })$。这套参数很可能非常糟糕，但它是一个起点。

2. 期望步骤 (Expectation Step / E-Step):

   • 既然有了一套（虽然很差的）模型 $\bar{\lambda }$，我就能“假装”它是对的。
   • 基于这套 $\bar{\lambda }$ 和我的观测序列 $O$，我可以计算出一些期望值。比如，我可以计算在时刻 $t$ 处于状态 $s_i$ 的“期望概率”是多少（这正是我们上次在5. 补充概率概念中学习的 ${\gamma }_t(i)$），以及在时刻 $t$ 从 $s_i$ 转移到 $s_j$ 的“期望概率”是多少（这正是 ${\xi }_t(i,j)$）。
   • 这些“期望”值，可以看作是对缺失的隐藏状态信息的一种“软性”填充或“软计数”。我们不再说“在时刻 $t$ 一定是状态 $s_i$”，而是说“在时刻 $t$ 有 $p$ 的概率是状态 $s_i$”。

3. 最大化步骤 (Maximization Step / M-Step):

   • 现在我有了这些“期望（软）计数”，我就能像在“有监督”情况下一样，来重新估计一套新的、更好的模型参数 $\lambda$。
   • 例如，新的 ${\hat{a}}_{ij}$ 不再是“硬计数”相除，而是“从 $s_i$ 转移到 $s_j$ 的期望总次数”除以“从 $s_i$ 出发的期望总次数”。
   • 通过这一步，我们得到了一套比初始猜测更好的新模型 $\lambda$。

4. 迭代 (Iteration): 将新得到的模型 $\lambda$ 作为下一轮的“旧模型” $\bar{\lambda }$，然后重复第2步和第3步。

这个“猜测-精炼”的过程不断迭代，每一次迭代都会让模型参数变得更好（或者说，至少不会变差），直到模型参数收敛到一个稳定的值。

3. 数学基础与理论推导

现在，我们将上述直观思想转化为严谨的数学推导。

3.1 从最大似然估计出发

我们的目标是找到一组参数 $\lambda$，使得观测数据 $O$ 出现的概率 $P(O|\lambda )$ 最大化。这被称为最大似然估计 (Maximum Likelihood Estimation, MLE)。我们通常处理其对数形式，即最大化对数似然函数:

$$L(\lambda )=\log P(O|\lambda )$$

我们知道，$P(O|\lambda )$ 是通过对所有可能的隐藏状态序列 $Q$ 求和得到的：

$$L(\lambda )=\log \left(\sum _{\text{所有 }Q}P(O,Q|\lambda )\right)$$

这个公式非常棘手，因为“对数的内部是求和”。这导致我们无法直接对其求导并令其为零来解析地求解 $\lambda$。

3.2 期望最大化算法框架与琴生不等式

为了解决这个问题，EM 算法引入了一个巧妙的技巧。假设我们当前的模型参数是 $\bar{\lambda }$，我们希望找到一个更好的新参数 $\lambda$。我们来考察对数似然的增量：

$$L(\lambda )-L(\bar{\lambda })=\log P(O|\lambda )-\log P(O|\bar{\lambda })=\log \left(\frac{P(O|\lambda )}{P(O|\bar{\lambda })}\right)$$

我们可以把 $P(O|\lambda )$ 展开：

$$P(O|\lambda )=\sum _{Q}P(O,Q|\lambda )$$

现在，我们引入一个关于隐藏序列 $Q$ 的任意分布，为了方便，我们就用基于旧参数 $\bar{\lambda }$ 和观测 $O$ 的后验概率 $P(Q|O,\bar{\lambda })$。我们对上式做一个简单的变换（乘以1）：

$$P(O|\lambda )=\sum _{Q}P(Q|O,\bar{\lambda })\frac{P(O,Q|\lambda )}{P(Q|O,\bar{\lambda })}$$

这可以被看作是求 $\frac{P(O,Q|\lambda )}{P(Q|O,\bar{\lambda })}$ 在分布 $P(Q|O,\bar{\lambda })$ 下的期望。

$$P(O|\lambda )=E_{Q\sim P(Q|O,\bar{\lambda })}\left[\frac{P(O,Q|\lambda )}{P(Q|O,\bar{\lambda })}\right]$$

代回到对数似然函数中：

$$L(\lambda )=\log \left(E_{Q\sim P(Q|O,\bar{\lambda })}\left[\frac{P(O,Q|\lambda )}{P(Q|O,\bar{\lambda })}\right]\right)$$

此时，关键的数学工具——琴生不等式 (Jensen’s Inequality)——登场了。对于一个凹函数（如 log 函数），我们有 $\log (E[X])\ge E[\log (X)]$。应用此不等式：

$$\begin{array}{rlrl}L(\lambda )& =\log \left(E_{Q\sim P(Q|O,\bar{\lambda })}\left[\frac{P(O,Q|\lambda )}{P(Q|O,\bar{\lambda })}\right]\right)& & \\ & \ge E_{Q\sim P(Q|O,\bar{\lambda })}\left[\log \left(\frac{P(O,Q|\lambda )}{P(Q|O,\bar{\lambda })}\right)\right]& & \text{（应用琴生不等式）}\\ & =\sum _{Q}P(Q|O,\bar{\lambda })\left(\log P(O,Q|\lambda )-\log P(Q|O,\bar{\lambda })\right)& & \\ & =\underset{Q(\lambda ,\bar{\lambda })}{\underbrace{\sum _{Q}P(Q|O,\bar{\lambda })\log P(O,Q|\lambda )}}-\underset{\text{与 }\lambda \text{ 无关的常数}}{\underbrace{\sum _{Q}P(Q|O,\bar{\lambda })\log P(Q|O,\bar{\lambda })}}& & \end{array}$$

我们发现，对数似然函数 $L(\lambda )$ 有了一个下界 (Lower Bound)。为了让 $L(\lambda )$ 尽可能大，我们只需要最大化这个下界。由于下界的第二项与我们要优化的新参数 $\lambda$ 无关，我们只需要最大化第一项即可。这个第一项，就是大名鼎鼎的 Q 函数。

3.3 E步：构建Q函数

E-Step 的目标就是计算 Q 函数的表达式。

$$Q(\lambda ,\bar{\lambda })=E_{Q|O,\bar{\lambda }}[\log P(O,Q|\lambda )]=\sum _{Q}P(Q|O,\bar{\lambda })\log P(O,Q|\lambda )$$

我们先展开 $\log P(O,Q|\lambda )$。对于一个特定的隐藏路径 $Q=(q_1,\dots ,q_T)$：

$$P(O,Q|\lambda )={\pi }_{q_1}{b}_{q_1}(o_1)\prod _{t=2}^T{a}_{q_{t-1}{q}_t}{b}_{q_t}(o_t)$$

取对数后，乘积变成了求和：

$$\log P(O,Q|\lambda )=\log {\pi }_{q_1}+\sum _{t=1}^{T-1}\log a_{q_t{q}_{t+1}}+\sum _{t=1}^T\log b_{q_t}(o_t)$$

现在，我们将这个表达式代入 Q 函数，并利用期望的线性性质：

$$\begin{array}{rlrl}Q(\lambda ,\bar{\lambda })=& \sum _{Q}P(Q|O,\bar{\lambda })\log {\pi }_{q_1}& & \text{(Part 1: for }\pi \text{)}\\ +& \sum _{Q}P(Q|O,\bar{\lambda })\sum _{t=1}^{T-1}\log a_{q_t{q}_{t+1}}& & \text{(Part 2: for }A\text{)}\\ +& \sum _{Q}P(Q|O,\bar{\lambda })\sum _{t=1}^T\log b_{q_t}(o_t)& & \text{(Part 3: for }B\text{)}\end{array}$$

这三个部分可以分开处理。我们以第二部分为例，交换求和顺序：

$$\text{Part 2}=\sum _{t=1}^{T-1}\sum _{Q}P(Q|O,\bar{\lambda })\log a_{q_t{q}_{t+1}}$$

我们可以把所有状态组合都列出来：

$$\text{Part 2}=\sum _{t=1}^{T-1}\sum _{i=1}^N\sum _{j=1}^N\left(\sum _{Q\text{ s.t. }{q}_t=s_i,q_{t+1}=s_j}P(Q|O,\bar{\lambda })\right)\log a_{ij}$$

括号里的项，正是“给定观测 $O$ 和旧模型 $\bar{\lambda }$，在时刻 $t$ 处于状态 $s_i$ 且在时刻 $t+1$ 处于状态 $s_j$ 的概率”，也就是我们上次定义的 ${\xi }_t(i,j)$！ 同理，Part 1 中的系数是 ${\gamma }_1(i)$，Part 3 中的系数是 ${\gamma }_t(i)$。

因此，Q 函数可以被优美地写成：

$$\begin{array}{rl}Q(\lambda ,\bar{\lambda })=& \sum _{i=1}^{N}P(q_1=s_i|O,\bar{\lambda })\log {\pi }_i\\ +& \sum _{t=1}^{T-1}\sum _{i=1}^N\sum _{j=1}^{N}P(q_t=s_i,q_{t+1}=s_j|O,\bar{\lambda })\log a_{ij}\\ +& \sum _{t=1}^T\sum _{i=1}^{N}P(q_t=s_i|O,\bar{\lambda })\log b_i(o_t)\end{array}$$

使用我们熟悉的 $\gamma$ 和 $\xi$ 符号（注意它们都是基于旧模型 $\bar{\lambda }$ 计算的）：

$$\begin{array}{rl}Q(\lambda ,\bar{\lambda })=& \sum _{i=1}^N{\gamma }_1(i)\log {\pi }_i\\ +& \sum _{t=1}^{T-1}\sum _{i=1}^N\sum _{j=1}^N{\xi }_t(i,j)\log a_{ij}\\ +& \sum _{t=1}^T\sum _{i=1}^N{\gamma }_t(i)\log b_i(o_t)\end{array}$$

E-Step 的任务：就是利用旧参数 $\bar{\lambda }$，通过前向-后向算法计算出所有的 ${\gamma }_t(i)$ 和 ${\xi }_t(i,j)$ 的值。

3.4 M步：最大化求解

M-Step 的任务是最大化 Q 函数，以求解新的参数 $\lambda =(\pi ,A,B)$。由于 Q 函数中的三项分别只与 $\pi ,A,B$ 有关，我们可以将它们分开独立优化。

1. 优化 $\pi$: 我们需要最大化 $Q_{\pi }={\sum }_{i=1}^N{\gamma }_1(i)\log {\pi }_i$，同时满足约束条件 ${\sum }_{i=1}^N{\pi }_i=1$。 这是一个典型的拉格朗日乘子法问题。构建拉格朗日函数：

$$\mathcal{L}(\pi ,\eta )=\sum _{i=1}^N{\gamma }_1(i)\log {\pi }_i+\eta \left(\sum _{i=1}^N{\pi }_i-1\right)$$

对 ${\pi }_i$ 求偏导并令其为 0：

$$\frac{\partial \mathcal{L}}{\partial {\pi }_i}=\frac{{\gamma }_1(i)}{{\pi }_i}+\eta =0⟹{\pi }_i=-\frac{{\gamma }_1(i)}{\eta }$$

将其代入约束条件 ${\sum }_{i=1}^N{\pi }_i=1$：

$$\sum _{i=1}^N\left(-\frac{{\gamma }_1(i)}{\eta }\right)=1⟹-\frac{1}{\eta }\sum _{i=1}^N{\gamma }_1(i)=1$$

我们知道 ${\sum }_{i=1}^N{\gamma }_1(i)={\sum }_{i=1}^{N}P(q_1=s_i|O,\bar{\lambda })=1$，所以 $-\frac{1}{\eta }=1$，即 $\eta =-1$。 代回 ${\pi }_i$ 的表达式，得到更新公式：

$${\hat{\pi }}_i={\gamma }_1(i)$$

2. 优化 $A$: 我们需要最大化 $Q_A={\sum }_{t=1}^{T-1}{\sum }_{i=1}^N{\sum }_{j=1}^N{\xi }_t(i,j)\log a_{ij}$，满足约束 ${\sum }_{j=1}^N{a}_{ij}=1$ 对每一个 $i$ 都成立。 对每一个 $i$，构建拉格朗日函数：

$$\mathcal{L}(a_i,{\eta }_i)=\sum _{j=1}^N\left(\sum _{t=1}^{T-1}{\xi }_t(i,j)\right)\log a_{ij}+{\eta }_i\left(\sum _{j=1}^N{a}_{ij}-1\right)$$

对 $a_{ij}$ 求偏导并令其为 0：

$$\frac{\partial \mathcal{L}}{\partial a_{ij}}=\frac{{\sum }_{t=1}^{T-1}{\xi }_t(i,j)}{a_{ij}}+{\eta }_i=0⟹a_{ij}=-\frac{{\sum }_{t=1}^{T-1}{\xi }_t(i,j)}{{\eta }_i}$$

代入约束 ${\sum }_j{a}_{ij}=1$：

$$\sum _{j=1}^N\left(-\frac{{\sum }_{t=1}^{T-1}{\xi }_t(i,j)}{{\eta }_i}\right)=1⟹-\frac{1}{{\eta }_i}\sum _{j=1}^N\sum _{t=1}^{T-1}{\xi }_t(i,j)=1$$

我们知道 ${\sum }_{j=1}^N{\xi }_t(i,j)={\gamma }_t(i)$，所以 ${\sum }_{j=1}^N{\sum }_{t=1}^{T-1}{\xi }_t(i,j)={\sum }_{t=1}^{T-1}{\gamma }_t(i)$。 因此 ${\eta }_i=-{\sum }_{t=1}^{T-1}{\gamma }_t(i)$。 代回 $a_{ij}$ 的表达式，得到更新公式：

$${\hat{a}}_{ij}=\frac{{\sum }_{t=1}^{T-1}{\xi }_t(i,j)}{{\sum }_{t=1}^{T-1}{\gamma }_t(i)}$$

3. 优化 $B$: 我们需要最大化 $Q_B={\sum }_{t=1}^T{\sum }_{i=1}^N{\gamma }_t(i)\log b_i(o_t)$。我们可以把它按观测值 $v_k$ 重新组织一下：

$$Q_B=\sum _{j=1}^N\sum _{k=1}^M\left(\sum _{t=1,o_t=v_k}^T{\gamma }_t(j)\right)\log b_j(k)$$

约束条件是 ${\sum }_{k=1}^M{b}_j(k)=1$ 对每一个 $j$ 成立。 与优化 $A$ 完全类似，使用拉格朗日乘子法可得更新公式：

$${\hat{b}}_j(k)=\frac{{\sum }_{t=1,o_t=v_k}^T{\gamma }_t(j)}{{\sum }_{t=1}^T{\gamma }_t(j)}$$

这些推导出的更新公式，完美地对应了我们最初的直觉：“期望（软）计数”相除。

4. 鲍姆-韦尔奇算法流程总结

现在，我们可以清晰地总结出完整的鲍姆-韦尔奇算法流程。

输入: 观测序列 $O=(o_1,o_2,\dots ,o_T)$，隐藏状态数 $N$。 输出: HMM 模型参数 $\lambda =(A,B,\pi )$。

算法步骤:

1. 初始化: 随机初始化模型参数 ${\lambda }^{(0)}=(A^{(0)},B^{(0)},{\pi }^{(0)})$。通常需要保证初始概率不为0，且满足行和为1的约束。

2. 迭代循环: 对 $n=0,1,2,\dots$ 进行循环，直到收敛。

   • E步 (Expectation): 根据当前模型 ${\lambda }^{(n)}$ 和观测序列 $O$，利用前向算法和后向算法，计算出对于所有时刻 $t$ 和所有状态 $i,j$ 的：

     • ${\gamma }_t(i)=P(q_t=s_i|O,{\lambda }^{(n)})$
     • ${\xi }_t(i,j)=P(q_t=s_i,q_{t+1}=s_j|O,{\lambda }^{(n)})$

   • M步 (Maximization): 利用 E 步计算出的期望值，更新模型参数，得到新的模型 ${\lambda }^{(n+1)}$：

     • 更新 $\pi$: ${\pi }_i^{(n+1)}={\gamma }_1(i)$
     • 更新 $A$: $a_{ij}^{(n+1)}=\frac{{\sum }_{t=1}^{T-1}{\xi }_t(i,j)}{{\sum }_{t=1}^{T-1}{\gamma }_t(i)}$
     • 更新 $B$: $b_j(k{)}^{(n+1)}=\frac{{\sum }_{t=1,o_t=v_k}^T{\gamma }_t(j)}{{\sum }_{t=1}^T{\gamma }_t(j)}$

3. 终止: 当模型参数 ${\lambda }^{(n+1)}$ 与 ${\lambda }^{(n)}$ 的变化非常小（例如，参数矩阵的范数变化小于一个阈值），或者对数似然 $P(O|{\lambda }^{(n)})$ 的增长非常小时，停止迭代，返回最终的模型参数 ${\lambda }^{(n+1)}$。

5. 总结

本次，我们攻克了 HMM 最核心的学习问题。

• 核心挑战: 在隐藏状态序列未知的情况下，估计模型参数。
• 核心思想: 采用期望最大化 (EM) 的迭代思想，通过“猜测-精炼”的循环来逐步逼近最优解。
• 核心算法: 鲍姆-韦尔奇算法是 EM 思想在 HMM 上的具体实现。
  • E步: 利用现有模型，通过前向-后向算法计算期望（软计数）${\gamma }_t(i)$ 和 ${\xi }_t(i,j)$。
  • M步: 利用这些期望计数，通过简单的公式来更新模型参数，得到一个更好的模型。
• 鲍姆-韦尔奇算法保证了每次迭代后，模型对观测数据的对数似然 $P(O|\lambda )$ 都是非递减的。但需要注意的是，它可能会收敛到局部最优解，而非全局最优解。因此，一个好的初始值或多次不同初始值的尝试，对于获得理想结果非常重要。

至此，我们已经深入探讨了 HMM 的评估问题和学习问题。HMM 的最后一块拼图是解码问题 (Decoding)：给定模型和观测序列，如何找到最可能的那条隐藏状态路径？这个问题将引出另一个优美的动态规划算法——维特比算法 (Viterbi Algorithm)，我们将在后续进行探讨。