引言
LDA相关的资料已经非常丰富，比较好的如[1], [2], [3], [4] 和 [5]。
[1]是一篇中文博文，核心的点都提到了，基本参考了[2], [3], [4]的内容。
[2]是一篇科普佳作，结合数学史，把LDA涉及的很多函数及背后的数据原理通俗地讲了一遍，适合慢慢品读。
由于原始的LDA训练很慢，[3]主要讲了LDA的并行化实现，但是文献对LDA使用Gibbs采样涉及的公式进行了严密的推导，适合了解算法推导的每个环节。
[4]是较早一篇比较详细讲解LDA的论文，基本上[1], [2], [3]都参考了该文献, 并且GibbsLDA++ 的实现，也是主要参考该论文。
[5]讲述了LDA涉及的概率图模型的更多相关知识，但是还是需要有一定概率图理论基础。
总的来说，LDA算法的推导还是需要一定数理基础，但是Gibbs采样的实现却非常简单，伪代码一看就懂。
模型理解
简单地再次总结下要点。
概率基础：
- 多项分布
\begin{align*} \displaystyle Mult(\overrightarrow{n} |\overrightarrow{p},N) = \binom{N}{\overrightarrow{n}}\prod_{k=1}^K p_k^{n_k} \end{align*}
- 狄利克雷分布：容易发现先验是狄利克雷分布，似然是多项分布，两者相乘得到的后验也是狄利克雷分布
\begin{align*} Dir(\overrightarrow{p}|\overrightarrow{\alpha}) &= \frac{1}{\Delta(\overrightarrow{\alpha})} \prod_{k=1}^K p_k^{\alpha_k-1} \end{align*}\begin{align*} \int_{\overrightarrow{p}}^{} \prod_{k=1}^K p_k^{\alpha_k-1} = \Delta(\overrightarrow{\alpha}) \end{align*}
- Delta函数/Gama函数： 狄利克雷分布函数中有一个Delta函数，有的地方也说Beta函数，可以表示成Gama函数的形式，而Gama(x+1) = x * Game(x)
\begin{align*} \Delta(\overrightarrow{\alpha}) = \frac{\prod_{k=1}^K\Gamma(\alpha_k)} {\Gamma(\sum_{k=1}^K\alpha_k)} \end{align*}
Gibbs Sampling 推导过程：
1) 写出主题和词的联合概率公式，利用狄利克雷分布连乘部分积分等于Delta函数的性质，可以将联合分布只用Delta函数表示
\begin{align*}
p(\overrightarrow{\mathbf{w}},\overrightarrow{\mathbf{z}} |\overrightarrow{\alpha}, \overrightarrow{\beta})
&= p(\overrightarrow{\mathbf{w}} |\overrightarrow{\mathbf{z}}, \overrightarrow{\beta}) p(\overrightarrow{\mathbf{z}} |\overrightarrow{\alpha}) \\
&= \prod_{k=1}^K p(\overrightarrow{w}_{(k)} | \overrightarrow{z}_{(k)}, \overrightarrow{\beta}) \prod_{m=1}^M p(\overrightarrow{z}_m |\overrightarrow{\alpha}) \\
&= \prod_{k=1}^K \frac{\Delta(\overrightarrow{n}_k+\overrightarrow{\beta})}{\Delta(\overrightarrow{\beta})} \prod_{m=1}^M \frac{\Delta(\overrightarrow{n}_m+\overrightarrow{\alpha})}{\Delta(\overrightarrow{\alpha})}
\end{align*}
2) 写出排除当前词主题的条件概率公式，其实就是1中得到的两个联合概率的比值，全部都是Delta函数，而Delta函数又都可以表示成Gama函数，利用Gama(x+1) = x * Game(x)的性质，得到最终结果，这个就是Gibbs采样公式
\begin{align*}
p(z_i = k|\overrightarrow{\mathbf{z}}_{\neg i}, \overrightarrow{\mathbf{w}}) \propto
\frac{n_{m,\neg i}^{(k)} + \alpha_k}{\sum_{k=1}^K (n_{m,\neg i}^{(k)} + \alpha_k)}
\cdot \frac{n_{k,\neg i}^{(t)} + \beta_t}{\sum_{t=1}^V (n_{k,\neg i}^{(t)} + \beta_t)}
\end{align*}
3) 采样结束后，得到每篇文档每个词的主题，根据Dirichlet分布的期望，得到doc-topic这个多项分布的参数theta，以及topic-word这个多项分布的参数phi
GibbsLDA++代码剖析
相关编译和运行可以直接参考官方文档，这里只剖析一下代码。
整个代码的实现偏C风格，结构很清晰，注释也很丰富。
- lda.cpp 主函数
- model.h 定义了类model，核心实现
- dataset.h 定义了类document和类dataset，存储输入文档数据
- constants.h 定义了BUFF值和模型状态值
- strtokenizer.h 定义类strtokenizer，分割文本并保存token
- utils.h 定义类 utils，包含解析参数，生成模型保存名称，排序等函数
主要是理解类model的实现，LDA训练相关参数如下，其中alpha和beta对K维都是一样的：
// --- model parameters and variables ---
int M; // dataset size (i.e., number of docs)
int V; // vocabulary size
int K; // number of topics
double alpha, beta; // LDA hyperparameters
int niters; // number of Gibbs sampling iterations
int liter; // the iteration at which the model was saved
int savestep; // saving period
int twords; // print out top words per each topic
int withrawstrs;
double * p; // temp variable for sampling
int ** z; // topic assignments for words, size M x doc.size()
int ** nw; // cwt[i][j]: number of instances of word/term i assigned to topic j, size V x K
int ** nd; // na[i][j]: number of words in document i assigned to topic j, size M x K
int * nwsum; // nwsum[j]: total number of words assigned to topic j, size K
int * ndsum; // nasum[i]: total number of words in document i, size M
double ** theta; // theta: document-topic distributions, size M x K
double ** phi; // phi: topic-word distributions, size K x V
关键函数如下：
// init for estimation
int init_est();
int init_estc();
// estimate LDA model using Gibbs sampling
void estimate();
int sampling(int m, int n);
void compute_theta();
void compute_phi();
init_est()从头开始训练，init_estc()继续训练。文档为M，词表为V，主题数为K，初始状态下，每个词对主题k的数目nw[w][k]都为0，每个文档对主题k的数目nd[m][k]也都为0，每个主题的数目nwsum[k]均为0，每篇文档的词数ndsum[m]均为0。然后对于每篇文档的每个词，随机采样一个主题，更新数组nw, nd, nwsum, ndsum。theta数组为M*K， phi数组为K*V。
z = new int*[M];
for (m = 0; m < ptrndata->M; m++) {
int N = ptrndata->docs[m]->length;
z[m] = new int[N];
// initialize for z
for (n = 0; n < N; n++) {
int topic = (int)(((double)random() / RAND_MAX) * K);
z[m][n] = topic;
// number of instances of word i assigned to topic j
nw[ptrndata->docs[m]->words[n]][topic] += 1;
// number of words in document i assigned to topic j
nd[m][topic] += 1;
// total number of words assigned to topic j
nwsum[topic] += 1;
}
// total number of words in document i
ndsum[m] = N;
}
estimate()进行训练，对于每篇文档的每个词，采样一个主题，然后保存模型，并且计算theta和phi：
// for all z_i
for (int m = 0; m < M; m++) {
for (int n = 0; n < ptrndata->docs[m]->length; n++) {
// (z_i = z[m][n])
// sample from p(z_i|z_-i, w)
int topic = sampling(m, n);
z[m][n] = topic;
}
}
sampling(int m, int n)函数使用Gibbs采样公式计算当前词的主题分布，然后随机选择一个主题，用到了累计概率值高于随机值的方式，类似轮盘赌：
int model::sampling(int m, int n) {
// remove z_i from the count variables
int topic = z[m][n];
int w = ptrndata->docs[m]->words[n];
nw[w][topic] -= 1;
nd[m][topic] -= 1;
nwsum[topic] -= 1;
ndsum[m] -= 1;
double Vbeta = V * beta;
double Kalpha = K * alpha;
// do multinomial sampling via cumulative method
for (int k = 0; k < K; k++) {
p[k] = (nw[w][k] + beta) / (nwsum[k] + Vbeta) *
(nd[m][k] + alpha) / (ndsum[m] + Kalpha);
}
// cumulate multinomial parameters
for (int k = 1; k < K; k++) {
p[k] += p[k - 1];
}
// scaled sample because of unnormalized p[]
double u = ((double)random() / RAND_MAX) * p[K - 1];
for (topic = 0; topic < K; topic++) {
if (p[topic] > u) {
break;
}
}
// add newly estimated z_i to count variables
nw[w][topic] += 1;
nd[m][topic] += 1;
nwsum[topic] += 1;
ndsum[m] += 1;
return topic;
}
compute_theta()函数调用了推导公式：
\begin{align*}
\hat{\theta}_{mk} &= \frac{n_{m}^{(k)} + \alpha_k}{\sum_{k=1}^K (n_{m}^{(k)} + \alpha_k)}
\end{align*}
for (int m = 0; m < M; m++) {
for (int k = 0; k < K; k++) {
theta[m][k] = (nd[m][k] + alpha) / (ndsum[m] + K * alpha);
}
}
compute_phi()函数调用了推导公式：
\begin{align*}
\hat{\varphi}_{kw} &= \frac{n_{k}^{(w)} + \beta_t}{\sum_{w=1}^V (n_{k}^{(w)} + \beta_w)}
\end{align*}
for (int k = 0; k < K; k++) {
for (int w = 0; w < V; w++) {
phi[k][w] = (nw[w][k] + beta) / (nwsum[k] + V * beta);
}
}
对新的数据推导主题分布时，重新走了一遍Gibbs采样，相应的计算theta和phi的公式不变，但有略微调整，这里不再赘述。
LDA开源工具包
相关参考
- [1] 概率语言模型及其变形系列(2)-LDA及Gibbs Sampling
- [2] LDA数学八卦 (pdf)
- [3] Distributed Gibbs Sampling of Latent Topic Models: The Gritty Details (pdf)
- [4] Parameter estimation for text analysis (pdf)
- [5] Graphical Representation, Generative Model and Gibbs Sampling (pdf)