2 Latent Dirichlet Allocation Introduction
LDA是给文本建模的一种方法，它属于生成模型。生成模型是指该模型可以随机生成可观测的数据，LDA可以随机生成一篇由N个主题组成文章。通过对文本的建模，我们可以对文本进行主题分类，判断相似度等。在90年代提出的LSA中，通过对向量空间进行降维，获得文本的潜在语义空间。在LDA中则是通过将文本映射到主题空间，即认为一个文章有若干主题随机组成，从而获得文本间的关系。LDA模型有一个前提：bag of word。意思就是认为文档就是一个词的集合，忽略任何语法或者出现顺序关系。
3 生成模型
LDA的建模过程是逆向通过文本集合建立生成模型，在讨论如何建模时，我们先要理解LDA的生成模型如何生成一篇文档。
假设一个语料库中有三个主题：体育，科技，电影
一篇描述电影制作过程的文档，可能同时包含主题科技和主题电影，而主题科技中有一系列的词，这些词和科技有关，并且他们有一个概率，代表的是在主题为科技的文章中该词出现的概率。同理在主题电影中也有一系列和电影有关的词，并对应一个出现概率。当生成一篇关于电影制作的文档时，首先随机选择某一主题，选择到科技和电影两主题的概率更高；然后选择单词，选择到那些和主题相关的词的概率更高。这样就就完成了一个单词的选择。不断选择N个单词，这样就组成了一篇文档。
具体来说，生成一篇文档按照如下步骤：
1. 选择N，N服从Poisson(ξ)分布，这里N代表文档的长度。
2. 选择θ，θ服从Dirichlet(α)分布，这里θ是列向量，代表的是个主题发生的概率，α是dirichlet分布的参数
3. 对N个单词中的每一个:
a) 选择主题zn，zn服从Multinomial(θ)多项分布。zn代表当前选择的主题
b) 选择wn，根据p(wn | zn; β)：在zn条件下的多项分布。
上式中β是一个K x V的矩阵，βij = P(wj = 1 | zi = 1)，也就是说β记录了某个主题条件下生成某个单词的概率。
观察第二步，这里是LDA和PLSA的区别所在。假设每篇文档由3个主题组成，θ就表明每个主题发生的概率，比如{1/6,2/6,3/6}，这样不同的文档对应的θ也就不同，而θ可以用来判断文档的相似度等。
LDA Graphical model representation：
几乎所有讨论LDA的文章都包括上面的这幅图。它代表的概率模型：
上式计算边缘概率，便可得：
其中D代表一个语料库，M代表语料库中文档的总数。
4 参数估计
通过对LDA生成模型的讨论我们理解到对文本的建模实际上就是要计算α和β两个参数。α和β可以采用极大似然估计，但是这里遇到一个问题，就是似然函数由于α和β的耦合无法直接求出来：
回想前面提到过的variational inference方法，为了估计后验分布，寻找一个似然函数的下界，在这里，这个下界正好可以被用来做为参数估计，因此LDA原始paper[6]选择使用variational inference方法来计算似然函数的下界。这样，分别给定一个α和β的值，就可以计算出一个似然函数的值。极大似然函数的参数估计，就是要找出一对α和β，使得似然函数值最大。这时就用到了EM算法，每次E-STEP输入α和β，计算似然函数，也就是variational inference的过程，M-STEP最大化这个函数，求出α和β。这样不断迭代知道收敛，就求得了最终的α和β值。在variational inference中需要选取Q的分布形式，使得Q容易计算。在LDA原始paper中，作者选取了：
其中γ和Φ为q的参数。这里假设了θ和z相互独立，并丢掉w节点，简化的模型如下：
下面要做的工作就很明显了：
原始paper的作者在其paper的附录做了推导，计算出γ和Φ迭代公式：
其中：
接下来的工作，就是要进行EM迭代，直到α和β收敛。
E-STEP:
对每一篇文档，计算参数γ和Φ
M-STEP：
最大化Variational Inference中的下界，求出此时的 α和β
5 LDA实践
目前LDA的C,java,matlab实现都可以在网络上找到，本节主要讨论LDA的C/C＋＋实现。我在网上找到的C实现一共有两个版本，一个是LDA提出者Blei的C版本，可以在他的主页[7]下载了；另一是Daichi
Mochihashi实现的版本[8]。另外我自己也实现了一个C++版本的LDA，主要参考前面两位作者。三个版本的代码我都上传到http://code.google.com/p/lsa-lda。对比三个版本，我个人认为Blei的C版本最灵活，比如作者使用log形式保存β，这样减少了一定的计算量，再比如保存中间了变量，在计算充分统计量时同时计算了它们的和，方便以后的规范化。Daichi
Mochihashi实现的C版本条例很清晰，很容易读懂，唯一的缺点是其变量名并没有按照LDA原始paper中命名，因此读代码时候经常需要找对应。我实现的版本流程基本和Daichi Mochihashi的版本一样，更容易读懂，但是不灵活，有些地方比如估计α值程序不完整。因此本节讨论Daichi Mochihashi的代码实现。由于其代码的变量和LDA Paper中的变量不对应，因此在开始讨论代码前，我想先把里面关键的变量都罗列出来，这样方便在读代码时候的回头查看。程序涉及到的全局变量有点多，我刚开始读代码时候就因为混淆变量的含义而走了许多弯路。
|
变量名 |
Paper中对应的变量 |
型式 |
|
beta |
β |
二维数组，行代表单词，列代表主题，矩阵单元代表某主题生成某词的概率 |
|
alpha |
α |
数组，对应dirichlet分布的参数 |
|
nclass |
k |
整数，代表主题数，这个是由用户设置的值 |
|
gamma |
γ |
一维数组，variational inference中后验dirichlet分布的参数 |
|
gammas |
充分统计量 |
一维数组，形式同gamma，用于在m-step估计alpha的值 |
|
q |
φ |
二维数组，行代表文档里的单词，列代表主题，矩阵单元代表文档中某主题生成某词的概率 |
|
betas |
充分统计量 |
二维数组，形式同q，该变量用于在e-step统计信息，供m-step估计 β使用。 |
整个程序的步骤如下图（伪代码）：
1. 读入文档
首先需要读入语料库中的文档。Daichi Mochihashi的C实现使用的是SVMLight的文件输入格式，大致来说，就是每个文件行，开头是一个数字，代表有多少单词，接着是id:count的形式，id代表该单词的id号，count代表文档中这个单词出现的次数。读入的文档集合保存在document结构体里。
// feature.h
typedef struct {
int len; //文档长度
int *id; //单词id数组
double *cnt;//各单词出现次数数组
} document;
// lda.c
// function main()
/* open data */
if ((data = feature_matrix(argv[optind], &nlex, &dlenmax)) == NULL) {
fprintf(stderr,
"lda::
cannot open training data.
" 
;
exit(1);
}
2. 调用lda_learn，开始学习
//doc learn.c
//function lda_learn()
void
lda_learn (document *data, double *alpha, double **beta,
int nclass,
int nlex //单词列表长度,
int dlenmax//文档最大长度,
int emmax //最大迭代次数,
int demmax //inference迭代次数,
double epsilon //收敛阈值
{
……
｝
Lda_learn按照以下步骤：
首先随机规范初始化alpha，规范初始化beta
规范化也就是使得行或列之和为1
// learn.c
//function lda_learn()
for (i = 0; i < nclass; i++)
alpha[i] = RANDOM;
for (i = 0, z = 0; i < nclass; i++)
z += alpha[i];
for (i = 0; i < nclass; i++)
alpha[i] = alpha[i] / z;
qsort(alpha, nclass, sizeof(double), // 为alpha排序
(int (*)(const void *, const void *))doublecmp);
//初始化beta
for (i = 0; i < nlex; i++)
for (j = 0; j < nclass; j++)
beta[i][j] = (double) 1 / nlex;
然后初始化充分统计量，variational inference中需要用到的变量等
// learn.c
// function lda_learn()
/*
* initialize posteriors
* gammas和betas相当于gamma和beta的和，gamma和beta代表的是某一个文档，而gammas和betas代表的是所有文档。
*/
if ((gammas = dmatrix(n, nclass)) == NULL) {
fprintf(stderr,
"lda_learn:: cannot allocate gammas.
" ;
return;
}
if ((betas = dmatrix(nlex, nclass)) == NULL) {
fprintf(stderr,
"lda_learn:: cannot allocate betas.
" ;
return;
}
/*
* initialize buffers
*
*/
if ((q = dmatrix(dlenmax, nclass)) == NULL) {
fprintf(stderr,
"lda_learn:: cannot allocate q.
" ;
return;
}
if ((gamma = (double *)calloc(nclass, sizeof(double))) == NULL)
{
fprintf(stderr,
"lda_learn:: cannot allocate gamma.
" ;
return;
}
if ((ap = (double *)calloc(nclass, sizeof(double))) == NULL) {
fprintf(stderr,
"lda_learn:: cannot allocate ap.
" ;
return;
}
if ((nt = (double *)calloc(nclass, sizeof(double))) == NULL) {
fprintf(stderr,
"lda_learn:: cannot allocate nt.
" ;
return;
}
if ((pnt = (double*)calloc(nclass, sizeof(double))) == NULL) {
fprintf(stderr,
"lda_learn:: cannot allocate pnt.
" ;
return;
}
开始EM迭代，如果迭代次数超过设定值则跳出
// learn.c
// function lda_learn
for (t = 0; t < emmax; t++)
{
printf("iteration %d/%d.. ", t + 1, emmax);
fflush(stdout);
/*
* VB-E step
*
*/
/* iterate for data 对每一个文档*/
for (dp = data, i = 0; (dp->len) != -1; dp++, i++)
{
vbem(dp, gamma, q, nt, pnt, ap,
alpha, (const double **)beta, dp->len, nclass, demmax);
accum_gammas(gammas, gamma, i, nclass);
accum_betas(betas, q, nclass, dp);
}
/*
* VB-M step
*
*/
/* Newton-Raphson for alpha */
newton_alpha(alpha, gammas, n, nclass, 0);
/* MLE for beta */
normalize_matrix_col(beta, betas, nlex, nclass);
其中E步如下，具体分析见下面的代码注释：
对每个文档，首先用variatinal inference迭代计算似然函数的下界，并求出此时γ和φ
迭代的具体公式前面已经推导过：
其中
(1)
其中
符号代表log
gamma函数的一阶导数
下面的代码作者用了一个小技巧，首先注意上面的公式(1)
按照公式，在计算时不但要计算![clip_image026[1]](http://olds.blogcn.com/wp-content/uploads/357/35739/2009/11/clip-image02612.jpg "clip_image026[1]"){kind=small}
(gamma)，还要计算![clip_image026[2]](http://olds.blogcn.com/wp-content/uploads/357/35739/2009/11/clip-image02621.jpg "clip_image026[2]"){kind=small}
(sum(gamma))。这样计算量很大。而在代码中我们发现作者没有计算后者，为什么呢？我起初也不明白这个问题，后来写信给作者，作者很快解回信，告诉我其实，通过分析variational
inference的流程我们发现，φ最后是要被规范化的，而公式中的第二项与单个gamma无关，因为它是求和，这样，省略掉后者以后规范化的值等于保留后者规范化的值。因此在作者的程序中省略掉了后者的计算。
//vbem.c
//function vbem()
//初始化nt,由公式gamma＝alpha+sum(phi)
//作者令上式的后半部分sum(phi)为nt
for (k = 0; k < K; k++)
nt[k] = (double) L / K;
//迭代直到超过设定次数
for (j = 0; j < emmax; j++)
{
//计算phi
//首先循环计算公式中phi = eta * exp()的exp部分
//psi为求log gamma函数的一阶导数
//alpha[k] + nt[k]即为gamma，省略了公式中的sum(gamma)具体原因参见上面的分析
for (k = 0; k < K; k++)
ap[k] = exp(psi(alpha[k] + nt[k]));
/* accumulate q */
//计算q
q等于phi = eta * exp()
//其中exp（）就是上面求出的ap［k］
for (l = 0; l < L; l++)
for (k = 0; k < K; k++)
q[l][k] = beta[d->id[l]][k] * ap[k];
/* normalize q */
//规范化q， 也就是使每一行的和为1
for (l = 0; l < L; l++) {
z = 0;
for (k = 0; k < K; k++)
z += q[l][k];
for (k = 0; k < K; k++)
q[l][k] /= z;
}
//上面计算了q，也就是paper中对应的phi
//现在要通过q求gamma
//这里其实是求的nt ，也就是gamma＝alpha+sum(phi)的后半部分
//在函数的最后会加上alpha的
for (k = 0; k < K; k++) {
z = 0;
for (l = 0; l < L; l++)
//注意这里要乘以d->cnt[l]，因为有的单词不只出现一次
z += q[l][k] * d->cnt[l];
nt[k] = z;
}
/* converge? */
//判断是否收敛，这是通过判断nt和上一轮迭代的nt（保存到pnt中了）之差来判断的
if ((j > 0) && converged(nt, pnt, K, 1.0e-2))
break;
//保存这次迭代的nt值
for (k = 0; k < K; k++)
pnt[k] = nt[k];
}
//计算gamma，由公式gamma[k] = alpha[k] + nt[k];
for (k = 0; k < K; k++)
gamma[k] = alpha[k] + nt[k];
这样variational inference的迭代就完成了，我们已经求出了gamma和phi的值，对应paper中后验函数的参数
这是EM迭代的第一步，接下来需要把刚刚得到的gamma和phi（变量q）的值累加到gammas和betas中，因为我们所计算的gamma和beta代表的只是其中一个文档，而gammas和betas代表的是所有文档，用来对alpha和beta进行参数估计。
//learn.c
//function lda_leatrn.c
accum_gammas(gammas, gamma, i, nclass);
accum_betas(betas, q, nclass, dp);
循环每个文档后，E-step就完成了，这时 gammas和betas包含了所有文档对应的gamma和beta值的和。下面就开始M-step
从前面的关于LDA的流程图中我们可以看出来，对beta的估计实际上就是规范化的一个过程，而对alpha的估计是使用牛顿迭代法的一个过程。
// learn.c
// function lda_learn()
/*
* VB-M step
*
*/
/* Newton-Raphson for alpha */
newton_alpha(alpha, gammas, n, nclass, 0);
/* MLE for beta */
normalize_matrix_col(beta, betas, nlex, nclass);
// util.c
// function normalize_matrix_col()
// 规范化矩阵的列，也就是求列的和，然后每个元素除以该和
void
normalize_matrix_col (double **dst, double **src, int rows, int cols)
{
/* column-wise normalize from src -> dst */
double z;
int i, j;
for (j = 0; j < cols; j++) {
for (i = 0, z = 0; i < rows; i++)
z += src[i][j];
for (i = 0; i < rows; i++)
dst[i][j] = src[i][j] / z;
}
}
牛顿迭代法在newton.c中newton_alpha()函数，这里省略
m-step之后，计算likehood似然函数，并判断是否收敛
如果收敛，就完成了LDA的参数估计
最后，想单独提一下blei版本的LDA实现的几个不好理解的地方
在blei版本的C实现中
//lda－estimate。C
//function double doc_e_step()
for (n = 0; n < doc->length; n++)
{
for (k = 0; k < model->num_topics; k++)
{
ss->class_word[k][doc->words[n]] += doc->counts[n]*phi[n][k];
ss->class_total[k] += doc->counts[n]*phi[n][k];
}
}
这里的 ss->class_word就是Daichi Mochihashi版本中的betas，而这里的class_total是什么呢？其实这里class_total是为规范化准备的。我们知道到规范化时需要将行或列求和，作者只不过把求和步骤放到这里了。
而在m-step中：
//lda－model。C
//function void lda_mle
void lda_mle(lda_model* model, lda_suffstats* ss, int estimate_alpha)
{
int k; int w;
for (k = 0; k < model->num_topics; k++)
{
for (w = 0; w < model->num_terms; w++)
{
if (ss->class_word[k][w] > 0)
{
model->log_prob_w[k][w] =
log(ss->class_word[k][w]) -
log(ss->class_total[k]);
作者直接使用 model->log_prob_w[k][w] = log(ss->class_word[k][w]) – log(ss->class_total[k]);而不是规范化所应该的除法，为什么呢？
因为作者使用log保存的beta值，也就是这里的 log_prob_w，而log函数的特性log(a)-log(b)=log(a/b)
6 总结
从LSA到PLSA到LDA，对文本的建模一步步的完善，LDA在document到topic一层引入了dirichlet分布，这是它优于PLSA的地方，使得模型参数的数量不会随着语料库的扩大而增多。LDA建模中最关键的是对参数的估计，原始paper中使用的是variational inference和EM算法，但这不是必须的，实际上有更容易计算的方法：Gibbs Sampling[9]。目前已经有该方法的实现[10]。
7 Reference
[1] Christopher M. Bishop, (2006) Pattern Recognition and Machine Learning, Springer, ISBN 0-387-31073-8.
[2] http://en.wikipedia.org/wiki/Dirichlet_distribution
[3] Sean Borman, The Expectation Maximization Algorithm A short tutorial
[4] ChengXiang Zhai, (2007) A Note on the Expectation-Maximization (EM) Algorithm
[5] J. M. Winn. Variational Message Passing and its Applications. PhD thesis, University of Cambridge, October 2003.
[6] Blei, D.M., Ng, A.Y., Jordan, M.I.: Latent Dirichlet Allocation. Journal of Machine
Learning Research 3 (2003) 993–1022
[7] http://www.cs.princeton.edu/~blei/lda-c/
[8] http://chasen.org/~daiti-m/dist/lda/
[9] G. Heinrich. Parameter estimation for text analysis. Technical report, 2005.
[10] http://gibbslda.sourceforge.net/
8 External Link
[1] http://code.google.com/p/lsa-lda/
本文中程序的代码和相关资料
[2] http://www.cs.princeton.edu/~blei/lda-c/
Lda提出者的LDA实现
[3] http://chasen.org/~daiti-m/dist/lda/
一个比较容易读懂的LDA实现
9 Further Reading
9.1 EM:
[1] Christopher M. Bishop, (2006) Pattern Recognition and Machine Learning, Chapter 9, Springer, ISBN 0-387-31073-8.
[2] Jeff Bilmes, (1998) A Gentle Tutorial of the EM Algorithm and its Application to Parameter Estimation for Gaussian Mixture and Hidden Markov Models
[3] ChengXiang Zhai, (2007) A Note on the Expectation-Maximization (EM) Algorithm
[4] Sean Borman, The Expectation Maximization Algorithm A short tutorial
9.2 Variational Inference:
[1] Christopher M. Bishop, (2006) Pattern Recognition and Machine Learning, Springer, ISBN 0-387-31073-8.第十章
[2] J. M. Winn. Variational Message Passing and its Applications. PhD thesis, University of Cambridge, October 2003.