关闭

Latent dirichlet allocation note

标签: allocation文档algorithmmatrixfunctionnull
11962人阅读 评论(5) 收藏 举报
分类:

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:

clip_image002

几乎所有讨论LDA的文章都包括上面的这幅图。它代表的概率模型:

clip_image004

上式计算边缘概率,便可得:

clip_image006

其中D代表一个语料库,M代表语料库中文档的总数。

4 参数估计

通过对LDA生成模型的讨论我们理解到对文本的建模实际上就是要计算α和β两个参数。α和β可以采用极大似然估计,但是这里遇到一个问题,就是似然函数由于α和β的耦合无法直接求出来:

clip_image008

回想前面提到过的variational inference方法,为了估计后验分布,寻找一个似然函数的下界,在这里,这个下界正好可以被用来做为参数估计,因此LDA原始paper[6]选择使用variational inference方法来计算似然函数的下界。这样,分别给定一个α和β的值,就可以计算出一个似然函数的值。极大似然函数的参数估计,就是要找出一对α和β,使得似然函数值最大。这时就用到了EM算法,每次E-STEP输入α和β,计算似然函数,也就是variational inference的过程,M-STEP最大化这个函数,求出α和β。这样不断迭代知道收敛,就求得了最终的α和β值。在variational inference中需要选取Q的分布形式,使得Q容易计算。在LDA原始paper中,作者选取了:

clip_image010

其中γ和Φ为q的参数。这里假设了θ和z相互独立,并丢掉w节点,简化的模型如下:

clip_image012

下面要做的工作就很明显了:

clip_image014

原始paper的作者在其paper的附录做了推导,计算出γ和Φ迭代公式:

clip_image016

其中:

clip_image018

接下来的工作,就是要进行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估计 β使用。

整个程序的步骤如下图(伪代码):

clip_image020

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&quot ;) ;

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&quot ;) ;

return;

}

if ((betas = dmatrix(nlex, nclass)) == NULL) {

fprintf(stderr, "lda_learn:: cannot allocate betas.\n&quot ;) ;

return;

}

/*

* initialize buffers

*

*/

if ((q = dmatrix(dlenmax, nclass)) == NULL) {

fprintf(stderr, "lda_learn:: cannot allocate q.\n&quot ;) ;

return;

}

if ((gamma = (double *)calloc(nclass, sizeof(double))) == NULL)

{

fprintf(stderr, "lda_learn:: cannot allocate gamma.\n&quot ;) ;

return;

}

if ((ap = (double *)calloc(nclass, sizeof(double))) == NULL) {

fprintf(stderr, "lda_learn:: cannot allocate ap.\n&quot ;) ;

return;

}

if ((nt = (double *)calloc(nclass, sizeof(double))) == NULL) {

fprintf(stderr, "lda_learn:: cannot allocate nt.\n&quot ;) ;

return;

}

if ((pnt = (double*)calloc(nclass, sizeof(double))) == NULL) {

fprintf(stderr, "lda_learn:: cannot allocate pnt.\n&quot ;) ;

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迭代计算似然函数的下界,并求出此时γ和φ

迭代的具体公式前面已经推导过:

clip_image022

其中

clip_image024(1)

其中clip_image026符号代表log gamma函数的一阶导数

下面的代码作者用了一个小技巧,首先注意上面的公式(1)

按照公式,在计算时不但要计算clip_image026[1](gamma),还要计算clip_image026[2](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 8) 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.

Posted in 信息检索机器学习| Tagged | 2 Comments

0
0

查看评论
* 以上用户言论只代表其个人观点,不代表CSDN网站的观点或立场
    个人资料
    • 访问:1116620次
    • 积分:8658
    • 等级:
    • 排名:第2291名
    • 原创:32篇
    • 转载:304篇
    • 译文:0篇
    • 评论:126条
    最新评论