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,主要参考前面两位作者。三个版本的代码我都上传到/p/lsa-lda。对比三个版本,我个人认为Blei的C版本最灵活,比如作者使用log形式保存β,这样减少了一定的计算量,再比如保存中间了变量,在计算充分统计量时同时计算了它们的和,方便以后的规范化。Daichi Mochihashi实现的C版本条例很清晰,很容易读懂,唯一的缺点是其变量名并没有按照LDA原始paper中命名,因此读代码时候经常需要找对应。我实现的版本流程基本和Daichi Mochihashi的版本一样,更容易读懂,但是不灵活,有些地方比如估计α值程序不完整。因此本节讨论Daichi Mochihashi的代码实现。由于其代码的变量和LDA Paper中的变量不对应,因此在开始讨论代码前,我想先把里面关键的变量都罗列出来,这样方便在读代码时候的回头查看。程序涉及到的全局变量有点多,我刚开始读代码时候就因为混淆变量的含义而走了许多弯路。 整个程序的步骤如下图(伪代码): 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() /* * initializeposteriors *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 allocateap.\n" return; } if ((nt = (double *)calloc(nclass, sizeof(double))) == NULL) { fprintf(stderr, "lda_learn:: cannot allocatent.\n" return; } if ((pnt = (double*)calloc(nclass, sizeof(double))) == NULL) { fprintf(stderr, "lda_learn:: cannot allocatepnt.\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-Raphsonfor 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) 按照公式,在计算时不但要计算(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, () Pattern Recognition and Machine Learning, Springer, ISBN 0-387-31073-8. [2]/wiki/Dirichlet_distribution [3] Sean Borman, The Expectation Maximization Algorithm A short tutorial [4] ChengXiang Zhai, () A Note on the Expectation-Maximization (EM) Algorithm [5] J. M. Winn. Variational Message Passing and its Applications. PhD thesis, University of Cambridge, October . [6] Blei, D.M., Ng, A.Y., Jordan, M.I.: Latent Dirichlet Allocation. Journal of Machine Learning Research 3 () 993–1022 [7]http://www.cs.princeton.edu/~blei/lda-c/ [8]/~daiti-m/dist/lda/ [9] G. Heinrich. Parameter estimation for text analysis. Technical report, . [10]/ 8 External Link [1]/p/lsa-lda/ 本文中程序的代码和相关资料 [2]http://www.cs.princeton.edu/~blei/lda-c/ Lda提出者的LDA实现 [3]/~daiti-m/dist/lda/ 一个比较容易读懂的LDA实现 9 Further Reading 9.1 EM: [1] Christopher M. Bishop, () Pattern Recognition and Machine Learning, Chapter 9, Springer, ISBN 0-387-31073-8. [2] Jeff Bilmes, (199A Gentle Tutorial of the EM Algorithm and its Application to Parameter Estimation for Gaussian Mixture and Hidden Markov Models [3] ChengXiang Zhai, () 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, () 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 . Posted in信息检索,机器学习|TaggedLatent dirichlet allocation,lda,机器学习|2 Comments