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.\n" ;
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.\n" ;
return;
}
if ((betas = dmatrix(nlex, nclass)) == NULL) {
fprintf(stderr, "lda_learn:: cannot allocate betas.\n" ;
return;
}
/*
* initialize buffers
*
*/
if ((q = dmatrix(dlenmax, nclass)) == NULL) {
fprintf(stderr, "lda_learn:: cannot allocate q.\n" ;
return;
}
if ((gamma = (double *)calloc(nclass, sizeof(double))) == NULL)
{
fprintf(stderr, "lda_learn:: cannot allocate gamma.\n" ;
return;
}
if ((ap = (double *)calloc(nclass, sizeof(double))) == NULL) {
fprintf(stderr, "lda_learn:: cannot allocate ap.\n" ;
return;
}
if ((nt = (double *)calloc(nclass, sizeof(double))) == NULL) {
fprintf(stderr, "lda_learn:: cannot allocate nt.\n" ;
return;
}
if ((pnt = (double*)calloc(nclass, sizeof(double))) == NULL) {
fprintf(stderr, "lda_learn:: cannot allocate pnt.\n" ;
return;
}
开始EM迭代,如果迭代次数超过设定值则跳出
// learn.c
// function lda_learn
for (t = 0; t < emmax; t++)
{
printf("iteration %d/%d..\t", 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)
按照公式,在计算时不但要计算(gamma),还要计算
(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 = \beta * 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 = \beta * 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, (199 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.