Latent dirichlet allocation note_0-387-31073-8-CSDN博客

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估计 β使用。

Posted in 信息检索, 机器学习| Tagged Latent dirichlet allocation, lda, 机器学习| 2 Comments

整个程序的步骤如下图（伪代码）：

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

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

其中

(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, (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/