机器学习数学推导（十五）：隐马尔可夫模型

雾里传来脚步声，有人在你身后。你看不见人，只能听到短促、轻快、急促的声音。从节奏和音调判断，对方是在走路、跑步，还是跛行？如果听到一整段声音序列，哪种步态最有可能产生它？再进一步，根据我对“走路”建立的模型，这段声音本身出现的概率有多大？

这正是 HMM 的三个核心问题。有意思的是，它们最终都归结为一个技巧：把联合概率 $P(\mathbf{O}, \mathbf{I})$ 按时间分解成局部因子的乘积，然后用动态规划在时间维度上复用子计算结果。暴力枚举需要 $$O(N^T)$$ 的复杂度，而前向后向算法、Viterbi 算法和 Baum-Welch 算法则全都降到 $$O(N^2 T)$$ 。指数级复杂度之所以能塌缩成多项式，全靠马尔可夫假设——给定当前状态，未来与过去条件独立。

你将学到什么#

HMM 的联合分布和两个条件独立假设
前向 / 后向算法：计算 $P(\mathbf{O}\mid\lambda)$ 和后验平滑 $\gamma_t(i), \xi_t(i,j)$
Viterbi：在动态规划中把 sum 替换为 max，直接求解 MAP
Baum-Welch：EM 算法在 HMM 中的特化版本；M 步公式为什么对应“期望计数”
数值问题（下溢、归一化）以及现代改进方法（CRF、RNN、CTC）

前置知识#

马尔可夫链和随机矩阵
动态规划 / 记忆化
EM 算法与 ELBO（详见第十三篇）

模型#

HMM 是序列建模中最简单但又非平凡的隐变量模型。它包含两条随时间演化的平行链：

隐状态链 $\mathbf{I}=(i_1,\dots,i_T)$ ，其中 $i_t\in\{1,\dots,N\}$ ；
观测链 $\mathbf{O}=(o_1,\dots,o_T)$ ，其中 $o_t\in\{v_1,\dots,v_M\}$ 。

整个模型由三组参数 $\lambda=(\boldsymbol{\pi}, \mathbf{A}, \mathbf{B})$ 完全定义：

参数模块	符号	含义
初始分布	$\pi_i = P(i_1 = i)$	隐状态链如何开始
转移矩阵	$a_{ij} = P(i_{t+1}=j \mid i_t=i)$	状态之间如何转移
发射矩阵	$b_j(k) = P(o_t = v_k \mid i_t = j)$	状态如何生成观测

两个条件独立性假设让计算变得可行：

状态的一阶马尔可夫性： $P(i_{t+1}\mid i_{1:t}) = P(i_{t+1}\mid i_t)$ 。
观测独立性： $$o_t$$ 只依赖当前状态 $$i_t$$ ，与其他状态或观测无关。

P(\mathbf{O}, \mathbf{I} \mid \lambda) = \pi_{i_1}\,b_{i_1}(o_1) \prod_{t=2}^{T} a_{i_{t-1},i_t}\,b_{i_t}(o_t).

本文中所有算法的核心目标，就是在不枚举 $$N^T$$ 条隐路径的情况下，对这个乘积进行求和或求最大值。

三状态天气模型#

全文用一个三状态天气模型——晴天、雨天、多云——作为示例，转移概率如图所示。“晴天”具有较强的黏性（ $a_{\text{晴晴}}=0.70$ ），而“多云”是连接其他状态的枢纽。这个模型规模小到可以手动推导，结构却足够复杂，能展示算法的所有细节。

三个问题，一个联合分布#

给定 $\lambda$ ，我们可以提出以下三个核心问题：

编号	问题	输入	输出	算法
1	概率计算	$\lambda, \mathbf{O}$	$P(\mathbf{O}\mid\lambda)$	前向 / 后向
2	解码	$\lambda, \mathbf{O}$	$\arg\max_{\mathbf{I}} P(\mathbf{I}\mid\mathbf{O},\lambda)$	Viterbi
3	学习	$\mathbf{O}$ （与 $$N$$ ）	$\hat\lambda = \arg\max_\lambda P(\mathbf{O}\mid\lambda)$	Baum-Welch（EM）

如果直接枚举所有隐路径来解决第一个问题，当 $$N{=}50, T{=}100$$ 时，计算量会达到约 $10^{170}$ 次。接下来三节将把这三个问题的复杂度降到多项式时间。

前向算法：从左到右计算#

前向算法的网格图：α 值从左到右流动，每个节点将入射路径概率求和后乘以本地发射概率

\alpha_t(i) = P(o_1, o_2, \dots, o_t,\; i_t = i \mid \lambda),

它表示“到时刻 $$t$$ 为止生成了所有观测，并且当前处于状态 $$i$$ ”的联合概率。

\alpha_1(i) = \pi_i\, b_i(o_1).

\begin{aligned} \alpha_t(j) &= P(o_1,\dots,o_t,\, i_t=j) \\ &= \sum_{i=1}^N P(o_1,\dots,o_{t-1},\, i_{t-1}=i)\,P(i_t=j\mid i_{t-1}=i)\,P(o_t\mid i_t=j)\\ &= \left[\sum_{i=1}^N \alpha_{t-1}(i)\, a_{ij}\right] b_j(o_t). \end{aligned}

方括号中的求和是唯一跨越时间边界的计算，其余部分都是局部操作。这是动态规划最简洁的形式。

P(\mathbf{O}\mid\lambda) = \sum_{i=1}^N \alpha_T(i).

复杂度
总共有 $$T$$ 个时间步，每步涉及 $$N$$ 个目标状态，每个目标状态需要对 $$N$$ 个前驱状态求和，因此复杂度为 $$O(N^2 T)$$ 。当 $$N=50, T=100$$ 时，总共需要约 $2.5\times 10^{5}$ 次运算，比暴力枚举快了大约 $10^{165}$ 倍。

下溢问题
由于概率是几何级数相乘， $\alpha_t$ 最终会因为数值过小而发生下溢。有两种常用解决方法：

对数空间 + logsumexp：

\log\alpha_t(j) = \mathrm{lse}_i\!\big(\log\alpha_{t-1}(i) + \log a_{ij}\big) + \log b_j(o_t).

缩放前向算法：
在每一步用 $c_t = \sum_j \alpha_t(j)$ 对 $\alpha_t$ 进行归一化，并累加 $\log P(\mathbf{O}) = \sum_t \log c_t$ 。

后向算法与后验平滑#

前向计算只能求出 $P(\mathbf{O})$ 。如果想得到每个时刻隐状态的后验分布，还得从右往左再来一轮。

\beta_t(i) = P(o_{t+1}, o_{t+2}, \dots, o_T \mid i_t = i, \lambda).

注意条件的区别： $\beta$ 是在已知当前状态下，未来观测的条件概率；而 $\alpha$ 是过去观测和状态的联合概率。两者互为镜像。

边界条件：时间超过 $$T$$ 后没有观测了，所以对所有 $$i$$ 都有 $\beta_T(i) = 1$ 。

\beta_t(i) = \sum_{j=1}^N a_{ij}\, b_j(o_{t+1})\, \beta_{t+1}(j).

P(\mathbf{O}\mid\lambda) = \sum_i \pi_i\,b_i(o_1)\,\beta_1(i) = \sum_i \alpha_T(i).

驱动学习的两个后验#

把 $\alpha$ 和 $\beta$ 都算好之后，下面两个量可以用 $$O(1)$$ 的额外开销直接得到：

\gamma_t(i) = P(i_t=i \mid \mathbf{O},\lambda) = \frac{\alpha_t(i)\,\beta_t(i)}{P(\mathbf{O}\mid\lambda)}.

\xi_t(i,j) = P(i_t=i, i_{t+1}=j \mid \mathbf{O},\lambda) = \frac{\alpha_t(i)\,a_{ij}\,b_j(o_{t+1})\,\beta_{t+1}(j)}{P(\mathbf{O}\mid\lambda)}.

这两个统计量正是 Baum-Welch 算法 E 步需要的核心结果。

Viterbi：从求和到求最大#

Viterbi 网格：max-product 动态规划高亮出唯一一条最可能的隐路径

\mathbf{I}^* = \arg\max_{\mathbf{I}}\, P(\mathbf{I}\mid\mathbf{O},\lambda) = \arg\max_{\mathbf{I}}\, P(\mathbf{O},\mathbf{I}\mid\lambda).

分母 $P(\mathbf{O})$ 和 $\mathbf{I}$ 无关，可以直接忽略。

\delta_t(j) = \max_{i_1,\dots,i_{t-1}} P(i_1,\dots,i_{t-1}, i_t=j,\, o_1,\dots,o_t\mid \lambda),

\boxed{\;\delta_t(j) = \max_{i}\big[\delta_{t-1}(i)\, a_{ij}\big]\, b_j(o_t).\;}

\psi_t(j) = \arg\max_i\big[\delta_{t-1}(i)\, a_{ij}\big].

前向计算完成后，取 $i_T^* = \arg\max_i \delta_T(i)$ ，然后通过回溯恢复路径： $i_t^* = \psi_{t+1}(i_{t+1}^*)$ 。

为什么把 sum 换成 max 是可行的？ 因为两种操作符在时间分解的联合分布上都满足分配律——“max-times” 和 “sum-times” 都构成可交换半环。动态规划的框架完全相同，只是局部归约方式不同。

\log\delta_t(j) = \max_i\big[\log\delta_{t-1}(i) + \log a_{ij}\big] + \log b_j(o_t).

所有运算都变成加法，彻底避免了下溢问题。

Baum-Welch：HMM 中的 EM 算法#

如果 $\lambda$ 是未知的，我的目标是通过最大化 $\log P(\mathbf{O}\mid\lambda)$ 来学习它。隐状态 $\mathbf{I}$ 是不可观测的，所以用 EM 算法来解决。这就是 Baum-Welch 算法——它比通用的 EM 框架早了好几年（Baum & Petrie, 1966），也是 EM 在实际应用中最漂亮的例子之一。

E 步：后验统计量#

假设当前参数是 $\lambda^{(k)}$ ，运行一次前向后向算法，得到 $\gamma_t(i)$ 和 $\xi_t(i,j)$ 。它们分别是指示变量 $\mathbb{1}[i_t=i]$ 和 $\mathbb{1}[i_t=i, i_{t+1}=j]$ 在当前模型下的期望值。

\log P(\mathbf{O},\mathbf{I}\mid\lambda) = \log\pi_{i_1} + \sum_{t=1}^{T-1}\log a_{i_t i_{t+1}} + \sum_{t=1}^{T}\log b_{i_t}(o_t).

Q(\lambda;\lambda^{(k)}) = \sum_i \gamma_1(i)\log\pi_i + \sum_{t=1}^{T-1}\sum_{i,j}\xi_t(i,j)\log a_{ij} + \sum_{t=1}^{T}\sum_{j,k}\gamma_t(j)\,\mathbb{1}[o_t=v_k]\log b_j(k).

这三项分别对应 $\boldsymbol{\pi}, \mathbf{A}, \mathbf{B}$ 的更新，彼此独立——M 步就是解三个独立的带约束优化问题。

M 步：三组闭式更新公式#

\text{参数} = \frac{\text{事件的期望计数}}{\text{条件上下文的期望计数}}.

\boxed{\; \hat\pi_i = \gamma_1(i),\qquad \hat a_{ij} = \frac{\sum_{t=1}^{T-1}\xi_t(i,j)}{\sum_{t=1}^{T-1}\gamma_t(i)},\qquad \hat b_j(k) = \frac{\sum_{t:\,o_t=v_k}\gamma_t(j)}{\sum_{t=1}^{T}\gamma_t(j)}. \;}

每个比值都可以理解为：从 $$i$$ 转移到 $$j$$ 的期望次数除以离开 $$i$$ 的总期望次数，发射概率同理。形式上和基于可观测数据的最大似然估计一致，只是把“实际观测”换成了“后验期望”。

收敛性#

EM 算法保证 $P(\mathbf{O}\mid\lambda^{(k+1)}) \geq P(\mathbf{O}\mid\lambda^{(k)})$ ——图中的曲线理论上一定是单调上升的。但问题是，似然函数是多峰的，Baum-Welch 只能找到局部最优解。常用的解决方法包括随机重启、用 k-means 或 segmental k-means 初始化，或者加入信息性先验。

应用：词性标注#

词性标注：Viterbi 选出最可能的标签序列；下方热力图展示了每个标签下各词的发射概率

词性标注（POS tagging）是 HMM 的经典应用。隐状态对应词性标签，比如 PRON、VERB、ADJ、NOUN 等，观测值则是具体的单词。转移概率建模语法规则，比如限定词后面通常跟着形容词或名词；发射概率建模词汇表信息，比如 love 大多时候是动词，偶尔是名词。

对于句子 “I love natural language processing”，Viterbi 算法给出的最佳路径是 PRON / VERB / ADJ / NOUN / NOUN。尽管 processing 在词汇上有歧义，但结果仍然正确，因为 ADJ -> NOUN 的转移概率非常高。

同样的引擎还用在语音识别中，状态是音素，观测是 MFCC 帧；也用在生物信息学的 profile HMM 中，状态是匹配、插入或删除列；还用在手势识别中。只要场景符合“离散隐状态生成带噪声的观测流”，HMM 就能发挥作用。

数值缩放：为什么简单的前向-后向算法会失效#

前向变量 $\alpha_t(i) = P(\mathbf{o}_{1:t}, q_t = s_i \mid \lambda)$ 会以几何级数的速度衰减。每一步都要乘以转移概率和发射概率，而这两者的值都不超过 1。假设 $$T = 200$$ 个 token，平均乘子为 $$0.1$$ ，那么 $\alpha_T \sim 10^{-200}$ ，远远低于 float32 的下溢阈值 $\approx 10^{-38}$ ，甚至比 float64 的 $\approx 10^{-308}$ 还要小。

c_t = \sum_i \tilde\alpha_t(i), \qquad \hat\alpha_t(i) = \tilde\alpha_t(i) / c_t,

\log P(\mathbf{o}_{1:T} \mid \lambda) = \sum_{t=1}^{T} \log c_t.

后向计算时，直接使用相同的 $$c_t$$ 进行归一化，这样平滑后验 $\gamma_t(i) = \hat\alpha_t(i)\hat\beta_t(i)$ 就能自动保持正确，不需要额外记录任何信息。

\log\alpha_{t+1}(j) = \log b_j(\mathbf{o}_{t+1}) + \mathrm{logsumexp}_i\big(\log\alpha_t(i) + \log a_{ij}\big).

这种方法每一步需要多做一次 exp 和 log 运算，但能干净地处理发射概率已经是 log 密度的情况（比如连续 Gaussian 发射）。像 hmmlearn 和 pomegranate 这些主流 HMM 库，默认都采用 log-space 实现。Viterbi 算法天然就在 log-space 中运行——因为 $\max$ 和单调变换可以交换顺序，所以不存在缩放问题。

成本、复杂度与 HMM 的适用场景#

三个核心算法的渐进时间复杂度都是 $$O(N^2 T)$$ ，其中 $$N$$ 是状态数， $$T$$ 是序列长度。 $$N^2$$ 来自状态转移的求和操作， $$T$$ 则来自序列扫描。EM（Baum-Welch）每轮迭代需要一次完整的前向-后向计算，因此单条序列的训练复杂度是 $O(N^2 T \cdot I)$ ，而整个语料库则是 $O(N^2 \sum_s T_s \cdot I)$ ，其中 $$I$$ 是迭代次数。

什么时候这个复杂度会成为问题？以词性标注为例，当 $$N = 45$$ 、 $$T = 200$$ 时， $N^2 T = 4 \times 10^5$ ，计算几乎是瞬时完成的。但在语音识别中， $N \approx 5000$ （每个上下文相关的三音子对应一个状态），此时 $N^2 = 2.5 \times 10^7$ ，每个 token 都要乘一次。这正是为什么传统语音系统会在 Viterbi 算法上使用 beam search（每步只保留 top- $$K$$ 状态），也是为什么这条赛道最终被 RNN 和 Transformer 接管——它们完全避开了 $$N^2$$ 的开销。

到 2026 年，什么情况下还应该选择 HMM 而不是 RNN 或 Transformer？

标注数据少，结构先验强： 如果只有 50 条句子和 10 个隐状态，“每步一个离散状态、Markov 转移”这种假设能比神经网络更好地防止过拟合。
需要精确的后验概率： Viterbi 能给出全局 MAP 路径，而 Transformer 上的 beam search 做不到。生物信息学中的基因预测和蛋白质二级结构分析至今仍在用 HMM，原因就在于此。
流式处理且内存有限： 前向滤波的内存需求是 $$O(N)$$ ，与序列长度 $$T$$ 无关；而 Transformer 的 attention 是 $$O(T)$$ 。
可解释性强： $A_{ij}$ 可以直接解读为“从状态 $$i$$ 转移到状态 $$j$$ 的概率”，并用领域知识验证。相比之下，Transformer 的某个 attention head 就很难解释了。

其他场景——比如长序列、大规模标注语料、复杂的发射分布——现代序列模型已经全面占优。HMM 仍然是很好的教学案例，因为它清晰地区分了“算什么”（滤波、平滑、解码、学习）和“怎么算”。

深入问答#

Q1：前向算法和 Viterbi 算法有什么区别？为什么算子的替换很重要？
前向算法计算的是边缘概率 $P(\mathbf{O}\mid\lambda) = \sum_{\mathbf{I}} P(\mathbf{O},\mathbf{I})$ ，而 Viterbi 算法计算的是 $\max_{\mathbf{I}} P(\mathbf{O},\mathbf{I})$ 。两者的动态规划框架相同，但使用的半环不同（前向用 sum-product，Viterbi 用 max-product）。前向算法回答“这些观测数据有多可信”，Viterbi 算法则回答“哪种状态序列最能解释这些数据”。

Q2：为什么 Viterbi 算法最大化联合概率而不是后验概率？
因为 $P(\mathbf{I}\mid\mathbf{O}) = P(\mathbf{O},\mathbf{I})/P(\mathbf{O})$ ，分母 $P(\mathbf{O})$ 对所有状态序列 $\mathbf{I}$ 都是常数。所以最大化联合概率等价于最大化后验概率，还能省去一次归一化计算。

Q3：Baum-Welch 算法在什么情况下会失效？
主要有三种典型问题：(a) 初始化不好——容易陷入平坦或退化的局部最优；(b) 标签置换——状态只能通过排列来区分，无法固定顺序；(c) 观测坍缩——某个状态的发射分布只集中在训练数据中出现过的符号上，未见符号的概率为零。对于 (c)，可以通过加一个小的 Dirichlet 先验来平滑发射分布。

Q4：为什么序列标注任务更常用 CRF 而不是 HMM？
CRF 是判别式模型，直接建模 $P(\mathbf{I}\mid\mathbf{O})$ ，可以利用观测序列 $\mathbf{O}$ 的全局特征（如大小写、后缀模式、上下文信息），不受条件独立假设的限制。而 HMM 在需要生成序列或训练数据稀缺时仍然有优势。

Q5：深度学习时代，HMM 还有用吗？
作为端到端的独立建模方法，HMM 基本被 RNN 和 Transformer 取代了。但它的推断算法依然活跃：现代语音系统中的 CTC 解码本质上是对齐网格上的前向算法；序列级蒸馏用的是 Viterbi 算法；结构化输出 Transformer 借鉴了 CRF，而 CRF 又借鉴了 HMM 的思想。

Q6：如何选择状态数 $$N$$ ？
先从小值开始尝试，然后用信息准则（ $\text{AIC} = -2\log L + 2|\lambda|$ ， $\text{BIC} = -2\log L + |\lambda|\log T$ ）或留出数据的似然值来评估。如果使用贝叶斯非参数方法（如带层次 Dirichlet 过程先验的 iHMM），可以直接把 $$N$$ 当作随机变量，让数据决定最佳值。

Q7：如何处理连续观测数据？
把发射矩阵换成概率密度函数即可：可以是单个高斯分布、混合高斯分布（语音领域经典的 GMM-HMM），或者神经网络密度（如混合密度网络 -> “neural HMM”）。前向后向递推公式保持不变，只是 $$b_j(o_t)$$ 变成了对观测值的似然评估。

练习题#

E1. 设 $$N=2$$ ， $\boldsymbol{\pi}=(0.6, 0.4)$ ， $\mathbf{B} = \begin{pmatrix}0.5 & 0.5 \\ 0.4 & 0.6\end{pmatrix}$ ，且 $$o_1 = v_1$$ ，计算 $\alpha_1$ 。

解答

$\alpha_1(1) = 0.6 \cdot 0.5 = 0.30$ ， $\alpha_1(2) = 0.4 \cdot 0.4 = 0.16$ 。

E2. 当 $$N=50$$ 、 $$T=100$$ 时，比较前向算法和暴力枚举的复杂度。

解答

前向算法需要 $N^2 T = 2.5\times 10^{5}$ 次乘法。暴力枚举则有 $N^T \approx 10^{170}$ 条路径——完全不可行。

E3. 解释公式 $\hat a_{ij} = \frac{\sum_t \xi_t(i,j)}{\sum_t \gamma_t(i)}$ 的含义。

解答

从状态 $$i$$ 转移到状态 $$j$$ 的期望次数，除以从状态 $$i$$ 离开的总期望次数。这是对硬性 MLE 计数的一种软性推广。

E4. 证明：先运行前向算法再运行后向算法，得到的 $P(\mathbf{O}\mid\lambda)$ 和单独运行前向算法的结果相同。

解答

对于任意 $$t$$ ，有 $\sum_i \pi_i b_i(o_1)\beta_1(i) = \sum_i \alpha_1(i)\beta_1(i)/1 = \sum_i \alpha_t(i)\beta_t(i)$ （链式法则的结果）。当 $$t = T$$ 时， $\beta_T \equiv 1$ ，最终结果为 $\sum_i \alpha_T(i)$ 。

E5. 为什么 Viterbi 算法的时间复杂度是 $$O(N^2 T)$$ ，而回溯指针的空间复杂度是 $$O(NT)$$ ？

解答

Viterbi 算法的正向扫描和前向算法类似，复杂度为 $$O(N^2 T)$$ 。每个 (状态，时刻) 单元额外存储一个整数回溯指针，总共需要 $$NT$$ 个槽位。回溯阶段只需线性时间读取这些指针。

下一步#

HMM 是生成式的——它建模 $$P(x, z)$$ ，再通过贝叶斯反推 $P(z\mid x)$ 。这条路在监督序列标注（命名实体识别、词性标注）里有一个明显的浪费：我其实只关心 $P(z\mid x)$ ，却被迫把 $$P(x)$$ 也建模了一遍，而 $$P(x)$$ 的形式（特征独立性、马尔可夫性）经常和真实数据冲突。下一章是这个问题的判别式答案：条件随机场（CRF）。

CRF 直接建模 $P(z\mid x)$ ，对 $$x$$ 不做任何分布假设，因此可以自由地把任意特征（前后窗口、词形、外部词典）塞进特征函数里。它的训练目标是条件对数似然，凸的；推断仍然用前向-后向和 Viterbi——和 HMM 几乎一模一样的动态规划。理解 HMM 到 CRF 的转变，是理解"生成式 → 判别式"在序列模型里的完整复刻，也是后面 BiLSTM-CRF、Transformer-CRF 这些深度结构里那个 CRF 头到底在做什么的关键。

参考文献#

Rabiner, L. R. (1989). A Tutorial on Hidden Markov Models and Selected Applications in Speech Recognition. Proc. IEEE 77(2), 257-286.
Baum, L. E., & Petrie, T. (1966). Statistical Inference for Probabilistic Functions of Finite State Markov Chains. Annals of Math. Statist. 37(6).
Viterbi, A. J. (1967). Error Bounds for Convolutional Codes and an Asymptotically Optimum Decoding Algorithm. IEEE Trans. IT 13(2).
Murphy, K. P. (2012). Machine Learning: A Probabilistic Perspective, Ch. 17. MIT Press.
Bishop, C. M. (2006). Pattern Recognition and Machine Learning, Ch. 13. Springer.
Jelinek, F. (1997). Statistical Methods for Speech Recognition. MIT Press.
Lafferty, J., McCallum, A., & Pereira, F. (2001). Conditional Random Fields. ICML.
Eddy, S. R. (1998). Profile Hidden Markov Models. Bioinformatics 14(9).