LDA基本介绍以及LDA源码分析(BLEI)

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);
        }

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;

// 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);
}

// 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);
}

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);
}

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);
    }
}