【统计学习方法死磕系列Ⅱ】潜在狄利克雷分配

一、引言：从文档到主题

想象你面前摆着 10000 篇新闻文章。你想知道它们大约在讨论哪些主题——政治、体育、科技、娱乐？每个主题由哪些关键词定义？每篇文章又涉及哪些主题、各占多少比例？

潜在狄利克雷分配（Latent Dirichlet Allocation, LDA）是解决这个问题的经典模型。它由 Blei、Ng 和 Jordan 于 2003 年提出，是主题模型（Topic Model）中最具代表性的工作，也是贝叶斯层级模型在自然语言处理中最成功的应用之一。

LDA 的核心假设优雅而深刻：每篇文档是关于若干主题的分布，每个主题是关于若干词语的分布。这一简单的三层贝叶斯模型，能够从无标签的文档集合中自动发现”主题”及其对应的词语分布。

二、从 LSA 到 pLSA 到 LDA 的演变

在进入 LDA 的细节之前，有必要理解主题模型的发展脉络。

2.1 三个模型的对比

模型	年份	数学基础	概率	文档-主题	可生成
LSA	1988	SVD	非概率	确定值	否
pLSA	1999	EM 算法	概率	确定值	否（对新文档不行）
LDA	2003	贝叶斯推断	概率	随机变量	是

LDA 相对于 pLSA 的关键改进是：在文档-主题分布和主题-词分布上加了 Dirichlet 先验，使模型成为完全的贝叶斯生成模型，能够对新文档进行推断，也避免了 pLSA 的参数爆炸问题。

三、LDA 的核心假设与生成过程

3.1 模型假设

LDA 假设一个包含 $K$ 个主题和 $V$ 个词汇的语料库 $D$（含 $M$ 篇文档）按以下方式生成：

每个主题 $k$ 有一个词分布 $\varphi_k \sim \text{Dirichlet}(\beta)$（$V$ 维向量）
每篇文档 $d$ 有一个主题分布 $\theta_d \sim \text{Dirichlet}(\alpha)$（$K$ 维向量）
文档 $d$ 中的每个词通过两步生成：
- 先按 $\theta_d$ 选一个主题 $z_{d,n} \sim \text{Multinomial}(\theta_d)$
- 再按该主题的词分布 $\varphi_{z_{d,n}}$ 选一个词 $w_{d,n} \sim \text{Multinomial}(\varphi_{z_{d,n}})$

其中 $\alpha$ 和 $\beta$ 是 Dirichlet 先验的浓度参数。

3.2 生成过程的伪代码

Algorithm: LDA Generative Process
----------------------------------------------
Input: M, K, V, α, β

1. For each topic k = 1..K:
       φ_k ~ Dirichlet(β)       // topic-word distribution

2. For each document d = 1..M:
       θ_d ~ Dirichlet(α)       // document-topic distribution
       
       For each word position n = 1..N_d:
           z_{d,n} ~ Categorical(θ_d)     // topic assignment
           w_{d,n} ~ Categorical(φ_{z_{d,n}})  // word generation

3.3 盘子表示法（Plate Notation）

LDA 的图模型可以用以下盘子表示：

外层盘子（M）：重复 M 次（每篇文档）
内层盘子（N）：重复 N_d 次（文档中的每个词）
外层 $K$：重复 K 次（每个主题）

变量的依赖关系：

$\alpha \to \theta_d \to z_{d,n}$（文档的主题分布影响每个词的主题）
$\beta \to \varphi_k \to w_{d,n}$（每个主题的词分布影响词的生成）
$z_{d,n} \to w_{d,n}$（词的主题决定词的生成）

3.4 Dirichlet 分布的直观理解

Dirichlet 分布是 Beta 分布的多元推广。一个 $K$ 维 Dirichlet 分布的概率密度为：

$$ p(\theta | \alpha) = \frac{\Gamma(\sum_k \alpha_k)}{\prod_k \Gamma(\alpha_k)} \prod_{k=1}^{K} \theta_k^{\alpha_k - 1} $$

其中 $\theta_k \geq 0, \sum_k \theta_k = 1$（因此在 $K-1$ 维单纯形上）。

对称 Dirichlet 先验（$\alpha_1 = \alpha_2 = \dots = \alpha_K = \alpha$）：

$\alpha < 1$：倾向于稀疏的分布——文档倾向于聚焦在少数几个主题上
$\alpha = 1$：均匀先验——等价于在单纯形上均匀采样
$\alpha > 1$：倾向于平滑的分布——文档倾向于混合多个主题

对于主题模型，通常取 $\alpha < 1$（如 0.1）来鼓励每个文档只涉及少量主题，同理 $\beta < 1$（如 0.01）来鼓励每个主题只涉及少量词汇。

四、LDA 的数学与推断

4.1 联合概率分布

LDA 的完全联合概率为：

$$ p(W, Z, \Theta, \Phi | \alpha, \beta) = \prod_{k=1}^{K} p(\varphi_k | \beta) \prod_{d=1}^{M} p(\theta_d | \alpha) \prod_{n=1}^{N_d} p(z_{d,n} | \theta_d) p(w_{d,n} | z_{d,n}, \varphi) $$

我们观测到的只有词 $W$，需要推断的是隐变量 $Z$（主题分配）、$\Theta$（文档-主题分布）和 $\Phi$（主题-词分布）。

4.2 推断目标

我们想要计算后验分布：

$$ p(Z, \Theta, \Phi | W, \alpha, \beta) = \frac{p(W, Z, \Theta, \Phi | \alpha, \beta)}{p(W | \alpha, \beta)} $$

分母 $p(W|\alpha,\beta)$ 涉及对 $\Theta$、$\Phi$ 和 $Z$ 的高维积分，无法解析计算。因此需要近似推断方法。

4.3 推断方法概述

方法	核心思想	优缺点
变分推断（VB）	用简单的变分分布近似后验	Blei 原文使用，快速但可能偏差
坍缩 Gibbs 采样	积分掉 $\Theta, \Phi$，仅对 $Z$ 采样	实现简单，结果较好，但慢
坍缩变分推断（CVB）	融合 Gibbs 采样和 VB	速度和精度的折中
在线变分推断	每次用部分数据更新	适合流数据或大数据集

本文将重点介绍坍缩 Gibbs 采样——这是最直观、也最容易手写实现的方法。

五、坍缩 Gibbs 采样推导

5.1 Dirichlet-Multinomial 共轭

Gibbs 采样的关键是利用 Dirichlet 与 Multinomial 的共轭关系：

如果 $\theta \sim \text{Dirichlet}(\alpha)$，且 $z_i \sim \text{Multinomial}(\theta)$，则 $\theta$ 的后验也是 Dirichlet：

$$ p(\theta | z_1, \dots, z_N) = \text{Dirichlet}(\alpha + c_1, \alpha + c_2, \dots, \alpha + c_K) $$

其中 $c_k$ 是观察到 $z_i = k$ 的次数。

这一共轭性质允许我们把 $\Theta$ 和 $\Phi$ 积分掉（”坍缩”），大幅简化采样问题。

5.2 坍缩 Gibbs 采样公式

将 $\Theta$ 和 $\Phi$ 积分掉后，只剩下主题分配 $Z$ 需要采样。对于文档 $d$ 中的第 $n$ 个词（设位置为 $i$），其主题分配 $z_{d,n}$ 的条件分布为：

$$ p(z_{d,n} = k | Z_{\setminus(d,n)}, W, \alpha, \beta) \propto \frac{n_{d,k}^{\setminus(d,n)} + \alpha_k}{\sum_{j}(n_{d,j}^{\setminus(d,n)} + \alpha_j)} \cdot \frac{m_{k, w_{d,n}}^{\setminus(d,n)} + \beta}{\sum_{v}(m_{k,v}^{\setminus(d,n)} + \beta)} $$

其中：

$n_{d,k}$：文档 $d$ 中分配给主题 $k$ 的词数
$m_{k,v}$：词汇 $v$ 被分配给主题 $k$ 的次数（在所有文档中）
$\setminus(d,n)$：除去当前位置的计数

直观理解：

第一项（文档-主题）：文档 $d$ 中主题 $k$ 的比例。话题 $k$ 在文档中越”火”，当前位置越倾向于选它
第二项（主题-词）：词 $w_{d,n}$ 在主题 $k$ 中出现的比例。该词在主题 $k$ 中越常见，当前位置越倾向选它

两者相乘，体现了 LDA 的”富者愈富”（rich get richer）效应。

5.3 算法伪代码

Algorithm: Collapsed Gibbs Sampling for LDA
----------------------------------------------
Input: corpus (documents of word indices), K, α, β, iterations
Output: topic assignments Z, φ_k, θ_d

// Initialization
For each document d:
    For each word position n in document d:
        Randomly assign topic z_{d,n} ∈ {1..K}
        Increment counters:
            n_{d,k} += 1
            m_{k,w_{d,n}} += 1

// Gibbs sampling iterations
For iter = 1 to iterations:
    For each document d:
        For each word position n in document d:
            // Remove current assignment
            Decrement counters for (d, n)
            
            // Compute conditional probability for each topic k
            For k = 1 to K:
                p(k) ∝ (n_{d,k} + α) * (m_{k,word} + β) 
                        / (Σ_j n_{d,j} + Kα) / (Σ_v m_{k,v} + Vβ)
            
            // Sample new topic from p(k)
            z_{d,n} ~ Categorical(softmax(log_p))
            
            // Update counters

// Estimate φ and θ after burn-in
For k = 1..K:
    φ_kv = (m_{k,v} + β) / (Σ_v' m_{k,v'} + Vβ)
For d = 1..M:
    θ_dk = (n_{d,k} + α) / (Σ_k' n_{d,k'} + Kα)

5.4 Python 实现

import numpy as np
from collections import defaultdict

class LDA:
    def __init__(self, n_topics=10, alpha=0.1, beta=0.01, n_iters=1000):
        self.K = n_topics
        self.alpha = alpha
        self.beta = beta
        self.n_iters = n_iters
    
    def fit(self, docs, vocab_size):
        """
        docs: list of list of word indices
        vocab_size: int
        """
        M = len(docs)
        V = vocab_size
        K = self.K
        
        # Counters
        ndk = np.zeros((M, K))       # doc-topic counts
        nkv = np.zeros((K, V))       # topic-word counts
        nk = np.zeros(K)             # total words per topic
        
        # Topic assignments
        z = []
        for d, doc in enumerate(docs):
            zd = []
            for w in doc:
                k = np.random.randint(K)
                zd.append(k)
                ndk[d, k] += 1
                nkv[k, w] += 1
                nk[k] += 1
            z.append(np.array(zd))
        
        # Gibbs sampling
        for it in range(self.n_iters):
            for d, doc in enumerate(docs):
                for n, w in enumerate(doc):
                    old_k = z[d][n]
                    # Remove current assignment
                    ndk[d, old_k] -= 1
                    nkv[old_k, w] -= 1
                    nk[old_k] -= 1
                    
                    # Compute probabilities
                    p = np.zeros(K)
                    for k in range(K):
                        # Doc-topic prob
                        theta = (ndk[d, k] + self.alpha)
                        # Topic-word prob
                        phi = (nkv[k, w] + self.beta) / (nk[k] + V * self.beta)
                        p[k] = theta * phi
                    
                    # Normalize and sample
                    p /= p.sum()
                    new_k = np.random.choice(K, p=p)
                    
                    # Update
                    z[d][n] = new_k
                    ndk[d, new_k] += 1
                    nkv[new_k, w] += 1
                    nk[new_k] += 1
        
        # Estimate distributions
        self.phi_ = np.zeros((K, V))
        for k in range(K):
            self.phi_[k] = (nkv[k] + self.beta) / (nk[k] + V * self.beta)
        
        self.theta_ = np.zeros((M, K))
        for d in range(M):
            self.theta_[d] = (ndk[d] + self.alpha) / (ndk[d].sum() + K * self.alpha)
        
        return self
    
    def get_top_words(self, vocab, n_words=10):
        """Return top n_words for each topic"""
        top_words = []
        for k in range(self.K):
            top_idx = np.argsort(self.phi_[k])[::-1][:n_words]
            top_words.append([vocab[i] for i in top_idx])
        return top_words

五-B、变分推断视角（Blei et al. 原始方法）

虽然坍缩 Gibbs 采样是最容易手写实现的方法，但 Blei、Ng 和 Jordan 在 2003 年的原始论文中使用的推断方法是变分推断（Variational Inference, VB）。它值得了解，因为它是后续大量扩展工作（如在线 LDA、动态主题模型）的基础。

变分推断的核心思想

LDA 的后验 $p(\Theta, \Phi, Z | W, \alpha, \beta)$ 由于分母（evidence）不可计算，无法直接使用。变分推断的策略是：用一族容易处理的变分分布 $q(\Theta, \Phi, Z | \gamma, \phi, \lambda)$ 来近似真实后验，然后通过最小化两者之间的 KL 散度来找到最佳近似：

$$ q^* = \arg\min_q \text{KL}(q \| p) $$

Blei et al. 使用平均场假设（mean-field assumption），将变分分布分解为：

$$ q(\Theta, \Phi, Z | \gamma, \phi, \lambda) = \prod_{d=1}^{M} q(\theta_d | \gamma_d) \prod_{k=1}^{K} q(\varphi_k | \lambda_k) \prod_{d=1}^{M} \prod_{n=1}^{N_d} q(z_{d,n} | \phi_{d,n}) $$

其中：

$\gamma_d$：文档 $d$ 的变分 Dirichlet 参数（近似 $\theta_d$ 的后验）
$\lambda_k$：主题 $k$ 的变分 Dirichlet 参数（近似 $\varphi_k$ 的后验）
$\phi_{d,n}$：词 $w_{d,n}$ 的变分 Multinomial 参数（近似 $z_{d,n}$ 的后验）

变分 EM 算法

变分推断通过坐标上升（coordinate ascent）迭代更新三个参数：

变分 E 步（对每篇文档）：

初始化 $\gamma_d = \alpha + N_d/K$，$\phi_{d,n,k} = 1/K$
迭代更新直到收敛：
- 更新 $\phi_{d,n,k} \propto \exp(\Psi(\gamma_{dk}) + \Psi(\lambda_{k,w_{d,n}}) - \Psi(\sum_v \lambda_{kv}))$
- 更新 $\gamma_{dk} = \alpha_k + \sum_{n=1}^{N_d} \phi_{d,n,k}$

其中 $\Psi$ 是 digamma 函数（$\Gamma$ 函数的对数导数）。

变分 M 步（对所有文档）：

$$ \lambda_{kv} = \beta_v + \sum_{d=1}^{M} \sum_{n: w_{d,n}=v} \phi_{d,n,k} $$

变分推断 vs Gibbs 采样

方面	变分推断	坍缩 Gibbs 采样
性质	优化问题（坐标上升）	采样问题（MCMC）
确定性	确定性	随机
收敛检测	ELBO 单调递增	迹图和 $\hat{R}$ 诊断
速度	通常更快（特别是在线变分）	通常较慢（但容易并行化）
理论基础	KL 散度最小化	平稳分布收敛
扩展性	容易导出在线版本	容易扩展到非共轭模型

两者在实践中都给出可接受的结果。选择取决于具体场景：需要处理流式数据选在线变分，模型非共轭选 Gibbs 采样。

六、LDA 的应用与实践

6.1 应用场景

领域	应用
文本分析	新闻主题发现，科技文章分类
推荐系统	基于主题相似度推荐内容
社交网络	用户兴趣画像，社区发现
生物信息学	基因表达模式发现
计算机视觉	视觉词袋模型，场景分类
数字人文	文学作品主题分析，历史文献挖掘

6.2 sklearn 实现

from sklearn.decomposition import LatentDirichletAllocation
from sklearn.feature_extraction.text import CountVectorizer

# Prepare data
vectorizer = CountVectorizer(max_features=10000, stop_words='english')
X = vectorizer.fit_transform(documents)  # n_samples x n_features

# Fit LDA
lda = LatentDirichletAllocation(
    n_components=10,     # number of topics
    doc_topic_prior=0.1,  # alpha
    topic_word_prior=0.01, # beta
    learning_method='batch',
    max_iter=100,
    random_state=42
)
lda.fit(X)

# Get top words for each topic
feature_names = vectorizer.get_feature_names_out()
for topic_idx, topic in enumerate(lda.components_):
    top_words = [feature_names[i] for i in topic.argsort()[:-10-1:-1]]
    print(f"Topic {topic_idx}: {', '.join(top_words)}")

# Get document-topic distributions
doc_topic_dist = lda.transform(X)  # shape: (n_samples, n_topics)

6.3 困惑度评估

困惑度（Perplexity）是评估 LDA 质量的标准指标，越低越好：

$$ \text{Perplexity}(D_{\text{test}}) = \exp\left(-\frac{\sum_{d \in D_{\text{test}}} \log p(w_d)}{\sum_{d \in D_{\text{test}}} N_d}\right) $$

# Compute perplexity on held-out data
perplexity = lda.perplexity(X_test)
print(f"Perplexity: {perplexity:.2f}")

七、面试高频问题

Q1：LDA 中的 $\alpha$ 和 $\beta$ 如何影响模型？如何选择？

$\alpha$ 控制文档-主题分布的稀疏性：

$\alpha < 1$：每篇文档主要涉及少数主题（更符合实际）
$\alpha > 1$：文档倾向于混合所有主题
典型值：$\alpha = 0.1$ 或 $50/K$

$\beta$ 控制主题-词分布的稀疏性：

$\beta < 1$：每个主题由少数特征词定义
$\beta > 1$：主题的词分布更平滑
典型值：$\beta = 0.01$（让每个主题用少量词就能刻画）

两者都可以通过超参数搜索（grid search 在 perplexity 上）来选择。或使用非对称的 $\alpha$（允许某些主题在先验上更常见）。

Q2：LDA 和 pLSA 的核心区别是什么？

pLSA 的文档-主题分布 $\theta_d$ 和主题-词分布 $\varphi_k$ 是模型参数（确定性值），直接用 MLE 估计。这导致两个问题：

过拟合：参数个数随文档数线性增长（$M \times K + K \times V$）
无法处理新文档：没有先验，无法为新文档推断主题分布

LDA 将 $\theta_d$ 和 $\varphi_k$ 视为随机变量，加 Dirichlet 先验。好处：

正则化：先验防止过拟合，参数个数固定（$K \times V$）
生成性：新文档可以直接推断 $\theta$ 而不需要重新训练

Q3：坍缩 Gibbs 采样中的”坍缩”是什么意思？

“坍缩”（collapsed）指的是通过共轭关系将连续变量 $\Theta$ 和 $\Phi$ 解析积分掉，只留下离散变量 $Z$ 进行采样。

积分后的采样空间大幅缩小（从连续参数空间 + 离散主题空间变为仅离散主题空间），收敛更快。采样后 $\Theta$ 和 $\Phi$ 可以通过计数平均值恢复：

$$ \hat{\theta}_{dk} = \frac{n_{dk} + \alpha}{\sum_j n_{dj} + K\alpha} $$

Q4：如何确定主题数 K？

困惑度：在不同 K 下计算 perplexity（hold-out），选困惑度最低的 K
一致性（Coherence）：衡量主题中 top words 的语义一致性（如使用 NPMI），越高越好
人类评估：人工检查主题的可解释性（topic interpretability）
先验知识：根据业务场景预设主题数

实践中通常尝试若干 K 值，综合考虑 perplexity、coherence 和人类可解释性。

Q5：LDA 假设词袋模型，忽略了词序信息。这会导致什么问题？如何改进？

词袋假设意味着 LDA 认为”我爱中国”和”中国爱我”是相同的文档。这在以下场景可能有问题：

情感分析（”not good” vs “good”）
需要词序的语义（短语、命名实体等）

改进方法：

N-gram 特征：使用 bigram 或 trigram 作为”词”
短语检测：先用算法检测短语，再将短语作为 token
HMM-LDA：结合 HMM 建模词序（如 Griffiths et al. 2005）
Dynamic Topic Model：建模主题随时间的变化
Contextualized Topic Model：使用 BERT 嵌入代替词频

Q6：LDA 的结果如何评估？Perplexity 和 Topic Coherence 哪个更好？

LDA 的评估分为自动评估和人工评估两个层面。

困惑度（Perplexity）：衡量模型对 hold-out 文档的预测能力。数学定义为：

$$ \text{Perplexity} = \exp\left(-\frac{\sum_{d \in D_{test}} \log p(w_d|\Phi, \alpha)}{\sum_{d \in D_{test}} N_d}\right) $$

困惑度越低，模型对未见文档的预测能力越强。但它有一个致命缺陷：困惑度与人类对主题质量的判断常常不一致。低困惑度的模型未必生成语义连贯的主题——模型可能学到了高频词搭配的统计规律，但没学到有意义的语义结构。

主题一致性（Topic Coherence）：衡量主题的 top-N 词在外部语料（如 Wikipedia）或训练语料中的共现程度。常用指标：

UMass Coherence（Mimno et al., 2011）：用训练语料中的词共现计算

$$ C_{\text{UMass}}(t; V^{(t)}) = \frac{2}{N(N-1)} \sum_{i=1}^{N} \sum_{j=1}^{i-1} \log \frac{D(w_i, w_j) + \epsilon}{D(w_j)} $$

其中 $D(w_i, w_j)$ 是词 $w_i$ 和 $w_j$ 的共文档频率，$D(w_j)$ 是 $w_j$ 的文档频率。
NPMI（Normalized Pointwise Mutual Information）：使用外部大语料（如 Wikipedia）计算

$$ \text{NPMI}(w_i, w_j) = \frac{\log \frac{P(w_i, w_j)}{P(w_i)P(w_j)}}{-\log P(w_i, w_j)} $$

主题一致性与人类判断的相关性远高于困惑度。实践建议：用 Topic Coherence 选择 K 和评估主题质量，困惑度仅作为参考。

人工评估（黄金标准）：

Topic Intrusion：在每个主题的 top words 中注入一个不相关的词，让标注者找出——标注者越容易找出，主题越连贯
Topic Interpretability：直接让标注者给主题的”可解释性”打分

最终，LDA 是否成功取决于主题能否为下游任务（搜索、推荐、分析）提供价值。数字指标是指南，人类判断才是终点。

八、总结

潜在狄利克雷分配（LDA）是主题模型的标杆。它通过 Dirichlet-Multinomial 共轭这一优雅的贝叶斯构造，实现了从无标签文本中自动发现语义主题。Blei、Ng 和 Jordan 在 2003 年提出的这个模型，至今仍然是文本挖掘、信息检索和推荐系统中的重要工具。

LDA 的美丽在于它的简单性：仅用三个变量（$\alpha, \beta, K$）就能从数以百万计的文档中发现有意义的主题结构。而这种”无监督发现隐藏结构”的能力，正是机器学习最令人着迷的特质之一。

在 LDA 之外，主题模型近年来也在进化——从词袋到词嵌入（如 lda2vec），从 Gibbs 采样到变分自编码器（如 NVDM），从静态话题到动态话题追踪。但无论怎样演变，LDA 的贝叶斯层级思想都将是理解这一切的出发点。

参考文献：

Blei, D. M., Ng, A. Y., & Jordan, M. I. (2003). Latent Dirichlet Allocation. Journal of Machine Learning Research, 3, 993-1022.
Griffiths, T. L., & Steyvers, M. (2004). Finding scientific topics. Proceedings of the National Academy of Sciences, 101(suppl 1), 5228-5235.
Heinrich, G. (2005). Parameter estimation for text analysis. Technical report.
李航. (2019). 《统计学习方法》（第2版）. 第17章.
Hoffman, M. D., Blei, D. M., & Bach, F. (2010). Online learning for latent Dirichlet allocation. NIPS.