基本介紹：
topic model，主題模型介紹：
http://www.cnblogs.com/lixiaolun/p/4455764.html 以及 (
http://blog.csdn.net/hxxiaopei/article/details/7617838)
topic model本質上就一個套路，在doc-word user-url user-doc等關系中增加topic層，擴充為2層結構，一方面可以降維，另一方面挖掘深層次的關系，用戶doc word user url的聚類。
LDA的理論知識不介紹太多，基本就講了原理以及推導兩個內容，原理比較簡單，推導過程貌似很簡單，就一個變分加上一些參數估計的方法就搞定了，但是具體的細節還沒明白，以后慢慢研究。簡單介紹下基本原理以及有意義的幾個公式：
plsa作為topic-model ，每篇文檔對應一系列topics,每個topic對應一批terms，有如下問題：
1.每篇文檔及其在topic上的分布都是模型參數，也就是模型參數隨着文檔的數目增加而增加，這樣容易導致overfitting
2.對於new doc，如何確定其topic 分布
LDA解決這個問題，沒必要把每個doc-topic分布作為模型參數，為doc-topic分布增加一個先驗概率，限制整體上文檔的topic分布，具體先驗分布的作用，之前已經介紹過。
doc-topic分布服從多項分布，狄利克雷分布是其共軛先驗。這樣參數的個數就變成K+N*K, N為詞個數，K為topic個數，與文檔個數無關。如果我想知道一個文檔的topic分布怎么辦？下面介紹下train以及predic的方法。作者采用了varitional inference進行推導，過程就免了，列出來幾個重要的公式：
論文中幾個重要公式：
概率模型：
D 表示文檔集合，最后就是保證P(D|α，β)最大。
phi的迭代公式，表示文檔中單詞n在topic i上的分布:
gamma的迭代公式，文檔在topic上的分布
Beta的迭代公式，model中topic-word分布：
alpha的迭代公式，model中文檔-topic分布的先驗參數，利用梯度下降法即可求解：
LDA最核心的迭代公式，針對每一篇文檔，計算每個詞的topic分布，從而計算文檔的topic分布：
基本邏輯：
1.初始模型參數，開始迭代，執行2,3,4,直至收斂
2.針對每一篇文檔，初始gamma以及phi參數，迭代該參數，直至收斂
2.1.計算文檔中每個單詞在topic上的分布，利用model中beta以及文檔-topic分布（2.2）
2.2.計算文檔-topic分布，利用模型參數alpha以及word-topic分布（2.1結果)
3.update模型參數beta，利用word-topic分布。
4.update模型參數alpha，利用2.2的結果gamma
源碼介紹
模型參數：
var_gamma:文檔-topic分布，每篇文檔都會計算其topic分布
phi: word-topic分布，針對每篇文檔，計算文檔中每個word的topic分布
lda_mode: lad的模型參數，里面包括beta以及alpha
lda_suffstats: 記錄統計信息，比如每個topic上每個word出現的次數，這是為了計算lda_model而存在
corpus： 全部文檔信息
document： 文檔的具體信息，包括word信息
corpus_initialize_ss：初始化 ss，隨機選擇一些文檔，利用文檔詞信息update ss 里面topic上term的頻度
| for (k = 0; k < num_topics; k++) { for (i = 0; i < NUM_INIT; i++) { d = floor(myrand() * c->num_docs); printf("initialized with document %d\n", d); doc = &(c->docs[d]); for (n = 0; n < doc->length; n++) { ss->class_word[k][doc->words[n]] += doc->counts[n]; } } for (n = 0; n < model->num_terms; n++) { ss->class_word[k][n] += 1.0; ss->class_total[k] = ss->class_total[k] + ss->class_word[k][n]; } } |
random_initialize_ss：隨機選取一些值初始化ss
run_em：執行EM算法，針對每篇文檔，利用其單詞以及初始化的β和α信息 更新模型，直至收斂。
該函數是一個框架性質的，具體的見下。
| void run_em(char* start, char* directory, corpus* corpus) { int d, n; lda_model *model = NULL; double **var_gamma, **phi; // allocate variational parameters var_gamma = malloc(sizeof(double*)*(corpus->num_docs)); for (d = 0; d < corpus->num_docs; d++) var_gamma[d] = malloc(sizeof(double) * NTOPICS); int max_length = max_corpus_length(corpus); phi = malloc(sizeof(double*)*max_length); for (n = 0; n < max_length; n++) phi[n] = malloc(sizeof(double) * NTOPICS); // initialize model char filename[100]; lda_suffstats* ss = NULL; if (strcmp(start, "seeded")==0) { model = new_lda_model(corpus->num_terms, NTOPICS); ss = new_lda_suffstats(model); corpus_initialize_ss(ss, model, corpus); lda_mle(model, ss, 0); model->alpha = INITIAL_ALPHA; } else if (strcmp(start, "random")==0) { model = new_lda_model(corpus->num_terms, NTOPICS); ss = new_lda_suffstats(model); random_initialize_ss(ss, model); lda_mle(model, ss, 0); model->alpha = INITIAL_ALPHA; } else { model = load_lda_model(start); ss = new_lda_suffstats(model); } sprintf(filename,"%s/000",directory); save_lda_model(model, filename); // run expectation maximization int i = 0; double likelihood, likelihood_old = 0, converged = 1; sprintf(filename, "%s/likelihood.dat", directory); FILE* likelihood_file = fopen(filename, "w"); while (((converged < 0) || (converged > EM_CONVERGED) || (i <= 2)) && (i <= EM_MAX_ITER)) { i++; likelihood = 0; zero_initialize_ss(ss, model); // e-step //這里是核心，針對每篇文檔計算相關模型參數 for (d = 0; d < corpus->num_docs; d++) { likelihood += doc_e_step(&(corpus->docs[d]), var_gamma[d], phi, model, ss); } // m-step lda_mle(model, ss, ESTIMATE_ALPHA); // check for convergence converged = (likelihood_old - likelihood) / (likelihood_old); if (converged < 0) VAR_MAX_ITER = VAR_MAX_ITER * 2; likelihood_old = likelihood; |
doc_e_step：執行EM中的E-step，干了兩件事情
1.基於文檔的單詞，update φ以及γ，得到doc-topic 分布以及 word-topic分布(lda_inference）
2.利用γ信息，更新α，這里的α只有一個參數，所以計算方式不同。
本來應該是αi = αi +
，但是實際α只有一個，所以作者通過在所有topic上的分布計算出α。
，但是實際α只有一個，所以作者通過在所有topic上的分布計算出α。
每篇文檔迭代一次，遍歷完所有文檔后，就計算出最后的α。
| double doc_e_step(document* doc, double* gamma, double** phi, lda_model* model, lda_suffstats* ss) { double likelihood; int n, k; // posterior inference likelihood = lda_inference(doc, model, gamma, phi); // update sufficient statistics double gamma_sum = 0; for (k = 0; k < model->num_topics; k++) { gamma_sum += gamma[k]; ss->alpha_suffstats += digamma(gamma[k]); } ss->alpha_suffstats -= model->num_topics * digamma(gamma_sum); 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->num_docs = ss->num_docs + 1; return(likelihood); } |
lad_inference：這是最核心的code，基於每個doc，計算對應γ和φ，也就是doc-topic以及word-topic分布。也就是如下代碼：
利用topic個數以及word個數，初始化γ以及φ。
嚴格實現這個過程，工程做了優化，對φ取了對數logφ，這樣降低計算復雜度，同時利用
log_sum接口
，計算log(φ1) log(φ2)...log(φk)計算出log(φ1+φ2.....+φk)，這樣利用(logφ1)-log(φ1+φ2.....+φk)即可完成歸一化。
利用：
解釋下代碼：
針對文檔中的每一個詞 n:
針對每個topic i：
單詞n在topic i上的分布為：φni = topic i在word n上的分布 × exp(Digamma(該文檔在toic i上的分布))
歸一化 φni
文檔在topic上的分布γi = 其先驗分布αi + 文檔內所有詞在topic i上的分布值和
lda-infernce，不只是在train的時候使用，對於一片新的文檔，同樣是通過該函數/方法計算文檔在tpoic上的分布，只依賴於模型參數α以及β
利用：
| double lda_inference(document* doc, lda_model* model, double* var_gamma, double** phi) { double converged = 1; double phisum = 0, likelihood = 0; double likelihood_old = 0, oldphi[model->num_topics]; int k, n, var_iter; double digamma_gam[model->num_topics]; // compute posterior dirichlet //init gama and php for (k = 0; k < model->num_topics; k++) { //初始化 γ以及φ var_gamma[k] = model->alpha + (doc->total/((double) model->num_topics)); digamma_gam[k] = digamma(var_gamma[k]); for (n = 0; n < doc->length; n++) phi[n][k] = 1.0/model->num_topics; } var_iter = 0; while ((converged > VAR_CONVERGED) && ((var_iter < VAR_MAX_ITER) || (VAR_MAX_ITER == -1))) { var_iter++; for (n = 0; n < doc->length; n++) { phisum = 0; for (k = 0; k < model->num_topics; k++) { oldphi[k] = phi[n][k]; //對於每個word，更新對應的topic，也就是公式中  的對數結果phi[n][k] = digamma_gam[k] + model->log_prob_w[k][doc->words[n]]; if (k > 0) //為歸一化做准備，通過log(a) +log(b)計算log(a+b) phisum = log_sum(phisum, phi[n][k]); else phisum = phi[n][k]; // note, phi is in log space } for (k = 0; k < model->num_topics; k++) { phi[n][k] = exp(phi[n][k] - phisum);//歸一化 //update γ，這里面沒有用到α，原始公式不同 var_gamma[k] =var_gamma[k] + doc->counts[n]*(phi[n][k] - oldphi[k]); // !!! a lot of extra digamma's here because of how we're computing it // !!! but its more automatically updated too. digamma_gam[k] = digamma(var_gamma[k]); printf("%d:%d: gmama: %f php: %f\n", n, k, var_gmama[k], php[n][k]); } } //計算似然結果，觀察是否收斂，計算采用公式 likelihood = compute_likelihood(doc, model, phi, var_gamma); assert(!isnan(likelihood)); converged = (likelihood_old - likelihood) / likelihood_old; likelihood_old = likelihood; // printf("[LDA INF] %8.5f %1.3e\n", likelihood, converged); } return(likelihood); } |
compute
likehood：這個函數實現的是blei LDA 公式15，也就是定義的似然函數，比較復雜，但是嚴格按照這個實現。
需要注意的是，blei的代碼，k個α值相同，計算時 包含α的計算進行了簡化。
利用：
| double compute_likelihood(document* doc, lda_model* model, double** phi, double* var_gamma) { double likelihood = 0, digsum = 0, var_gamma_sum = 0, dig[model->num_topics]; int k, n; for (k = 0; k < model->num_topics; k++) { dig[k] = digamma(var_gamma[k]); var_gamma_sum += var_gamma[k]; } digsum = digamma(var_gamma_sum); lgamma(α*k) - k*lgamma(alpha) likelihood = lgamma(model->alpha * model -> num_topics) - model -> num_topics * lgamma(model->alpha) - (lgamma(var_gamma_sum)); for (k = 0; k < model->num_topics; k++) { likelihood += (model->alpha - 1)*(dig[k] - digsum) + lgamma(var_gamma[k]) - (var_gamma[k] - 1)*(dig[k] - digsum); for (n = 0; n < doc->length; n++) { if (phi[n][k] > 0) { likelihood += doc->counts[n]* (phi[n][k]*((dig[k] - digsum) - log(phi[n][k]) + model->log_prob_w[k][doc->words[n]])); } } } return(likelihood); } |
lda_mle：針對每個文檔執行do_e_step后，更新了ss,也就是計算模型所需要數據，topic-word對的次數
| 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]); } else model->log_prob_w[k][w] = -100; } } if (estimate_alpha == 1) { model->alpha = opt_alpha(ss->alpha_suffstats, ss->num_docs, model->num_topics); printf("new alpha = %5.5f\n", model->alpha); } } |
轉自：http://blog.csdn.net/hxxiaopei/article/details/8034308