Memory and speed for an LDA model

[Julian Brooke]

Hello folks,

I'm trying to test some adaptions of the latent Dirichlet allocation (LDA) model, adaptions where collapsed Gibbs sampling isn't going to work right out of the box; I thought I'd try using a general purpose sampler, and STAN looked like it might be the best option in terms of speed. To begin, I coded up a version of standard LDA in stan (see below) and it seems to be working in very small data sets. However, if I push the number of text/size of the vocabulary to anything approaching what I'd need for a real world application, STAN quickly uses up all the memory on my server (92 Gigs!) and crashes, before the sampling processes even begins. I can scale it down a bit to the point where it has enough memory to start sampling, but I started a 50 text corpus a few days ago and its still running. This is much worse than the other general-purpose sampler I've tried (pymc). Anyway, I thought I'd check to see if maybe the problem was in my code, otherwise I have to reluctantly conclude that STAN (and general-purpose sampling in general) is not what I need, and I might as well go back to trying to hack existing LDA implementations for my purposes, which is not ideal:(. If anyone has any advice, specific or general, it would be much appreciated.

(Note: I don't sample the hyperparameters, instead I pass them in. Beta is a matrix of vectors so that I can use it to seed topics with specific words; Also, the entire corpus is a single vector, with switches between documents stored in a separate matrix (doc_ind); oh, and I have to add 1 to topic because categorical gives a range starting with zero whereas arrays start with 1, this mismatch seems like a bad design choice BTW)

data {

 int<lower=1> K; //topics

 int<lower=1> V; //vocabulary

 int<lower=1> D; //number of documents

 int<lower=1> T; // total number of words

 vector[K] alpha;

 vector[V] beta[K];

 

 int words[T];

 int doc_ind[D];

}

parameters {

simplex[V] psi[K];

simplex[K] theta[D];

}

model {

int start;

int topic;

for (n in 1:K) {

psi[n] ~ dirichlet(beta[n]);

}

for (n in 1:D) {

theta[n] ~ dirichlet(alpha);

}

start <- 1;

for (i in 1:D) {

for (j in start:doc_ind[i] - 1) {

topic ~ categorical(theta[i]);

words[j] ~ categorical(psi[topic + 1]);

}

start <- doc_ind[i];

}

 

}

[Jiqiang Guo]

Please see this related post 

https://groups.google.com/forum/?fromgroups=#!topic/stan-users/MJf_SwNmGFo

The answer is that you Stan for now does not have a sampler for discrete parameters so that your model is not supported in Stan now. 

--

Jiqiang

--
 
 

[Bob Carpenter]

We don't really recommend Stan for models like LDA.
The problem is that the multi-modality of the posterior
makes Bayesian inference intractable, not just for Stan,
but for any known technique.

That being the case, you might as well go for an optimizer.
Like the very scalable online one Matt (the inventor of
our sampler, NUTS) wrote for Vowpal Wabbit:

http://hunch.net/~vw/

and described in this paper:

http://www.cs.princeton.edu/~blei/papers/HoffmanBleiBach2010b.pdf

The reason Stan's consuming so much memory is that we're
building up the log probability expression tree in memory to
use to compute auto-diff.

I believe PyMC uses Metropolis, which may run and provide
samples, but isn't going to give you very efficient sampling,
even around a single mode.

Finally, the Stan model isn't well formed because
we don't do integer sampling. So the "topic ~ categorical(...)"
isn't going to do what you think -- what it will do is
add log categorical(topic|theta[i]) to the log probabilty
where topic is uninitialized.

What you need to do is marginalize out the discrete
variables because Stan doesn't sample them. Here's
a note on how to do that for basic LDA:

http://lingpipe-blog.com/2011/04/14/summing-out-topic-assignments-ld/

This can also be a useful move if you're applying some
kind of optimization.

> model {
> int start;
> int topic;
> for (n in 1:K) {
> psi[n] ~ dirichlet(beta[n]);
> }
> for (n in 1:D) {
> theta[n] ~ dirichlet(alpha);
> }
> start <- 1;
> for (i in 1:D) {
> for (j in start:doc_ind[i] - 1) {
> topic ~ categorical(theta[i]);
> words[j] ~ categorical(psi[topic + 1]);
> }
> start <- doc_ind[i];
> }
> }

Also, I don't follow the use of start here. Why not
write this loop as in the Blei et al. LDA paper? That is,
as:

for (d in 1:D) {
theta[d] ~ dirichlet(alpha);
for (j in 1:N[d]) {
topic[j] ~ categorical(theta[d]);
words[j] ~ categorical(psi[topic[j]]);
}
}

The indexing in Stan is from 1, by the way, so I doubt
you wanted the +1 inside of psi.

You'd also need an index, as in topic[j], if we
supported discrete sampling. At least the above is how
I'd write it in BUGS/JAGS (modulo their matrix slicing
syntax which I can never remember).

- Bob

> --
>
>

[Matt Hoffman]

I agree with everything Bob said. And would add that if you want to do
sampling for LDA, it's likely that (for now) you'll be better off with
the collapsed Gibbs sampler of Griffiths and Steyver (2004). There's
an implementation in Jonathan Chang's lda R package, which for some
reason was removed from CRAN but can still be downloaded here:
http://cran.r-project.org/src/contrib/Archive/lda/

Matt

> --
>
>

[Bob Carpenter]

There's also a Java implementation of the collapsed
sampler in my company's open-source (but not free)
software along with a tutorial on usage:

http://alias-i.com/lingpipe/demos/tutorial/cluster/read-me.html

The downside is the non-standard, restrictive, open-source
license:

http://alias-i.com/lingpipe/web/licensing.html

It's very very easy to implement the collapsed sampler. I
have a full derivation of integrating out the continuous
parameters linked from here:

http://lingpipe-blog.com/2010/07/13/collapsed-gibbs-sampling-for-lda-bayesian-naive-bayes/

- Bob

[Julian Brooke]

Ah, so a number of issues. Oddly, the compiler was giving me errors that made me think that categorical was doing what I wanted (e.g. saying topic was in the wrong range, e.g. 0-9 instead of 1-10, when I added 1 it stopped complaining) but I guess not. But even if I get rid of that topic selection categorical (which I actually did in my pymc implementation, now that I think about it), is the other one (the one that picks the words) kosher? Can't do LDA without something categorical, I wouldn't think...

Am I wrong in thinking that STAN can't deal with 2-dimensional objects where the second dimension is of varying length? That's why I don't just iterate over documents, instead I iterate over document end indices, and then over a range in the single corpus vector they point to. Start is really just there so I don't get an out-of-range error on the first iteration, though now that I think about it I could get rid of it by fixing my data.

Anyway, the bigger issue is performance, I guess STAN isn't the right thing from that perspective. Which is too bad, because it would be perfect otherwise; the problem with an off-the-shelf LDA solution is that I don't care about LDA specifically, I just want to test out some models which are similar to see if they work for my particular application, and I'd love to be able to fiddle at the exactly the level of abstraction provided in STAN (or other samplers, though I was particularly taken with STAN's detailed, well-written docs, thumbs up on that). I'm not much of a low-level programmer nor a statistics guru, so making collapsed gibbs sampling or any other optimization work for the actual models I have in mind (if even possible) is a slightly daunting prospect. But I think I'm finally convinced that's the only way to proceed.

BTW, Bob, we're actually acquainted, I met you briefly when you visited the UofT CL group a while back, unfortunately Gibbs sampling wasn't on my radar yet then, or you might have warned me away then. Thanks, everyone, for the info.

Julian

[Bob Carpenter]

On 10/8/12 6:44 PM, Julian Brooke wrote:
> Ah, so a number of issues. Oddly, the compiler was giving me errors that made me think that categorical was doing what I
> wanted (e.g. saying topic was in the wrong range, e.g. 0-9 instead of 1-10, when I added 1 it stopped complaining) but I
> guess not. But even if I get rid of that topic selection categorical (which I actually did in my pymc implementation,
> now that I think about it), is the other one (the one that picks the words) kosher? Can't do LDA without something
> categorical, I wouldn't think...

Yes you can. Just marginalize out the categoricals.
Here's the derivation:

http://lingpipe-blog.com/2011/04/14/summing-out-topic-assignments-ld/

> Am I wrong in thinking that STAN can't deal with 2-dimensional objects where the second dimension is of varying length?

Actually, we can handle ragged arrays inside Stan, but
it was unintentional and all the external interfaces assume
rectangular arrays.

> That's why I don't just iterate over documents, instead I iterate over document end indices, and then over a range in
> the single corpus vector they point to. Start is really just there so I don't get an out-of-range error on the first
> iteration, though now that I think about it I could get rid of it by fixing my data.

I see. Other people's index fiddles are always hard to read.

> Anyway, the bigger issue is performance, I guess STAN isn't the right thing from that perspective.

Very unlikely.

> Which is too bad,
> because it would be perfect otherwise; the problem with an off-the-shelf LDA solution is that I don't care about LDA
> specifically, I just want to test out some models which are similar to see if they work for my particular application,
> and I'd love to be able to fiddle at the exactly the level of abstraction provided in STAN (or other samplers, though I
> was particularly taken with STAN's detailed, well-written docs, thumbs up on that).

You might be happier with Hal Daume's HMC, which in part
is tailored to not only LDA but some process generalizations
thereof:

http://www.umiacs.umd.edu/~hal/HBC/

> I'm not much of a low-level
> programmer nor a statistics guru, so making collapsed gibbs sampling or any other optimization work for the actual
> models I have in mind (if even possible) is a slightly daunting prospect. But I think I'm finally convinced that's the
> only way to proceed.

Summing out discrete vars is easier. But the continuous
ones tend to do the obvious thing, once you have enough
experience that is, to know what the obvious thing is!

> BTW, Bob, we're actually acquainted, I met you briefly when you visited the UofT CL group a while back, unfortunately
> Gibbs sampling wasn't on my radar yet then, or you might have warned me away then. Thanks, everyone, for the info.

If by T you mean Toronto, then you've got lots of
in-house expertise with sampling and approximate
sampling. But not being a statistics expert may make it
hard to communicate. I know it was hard for me when
making the transition from CS/machine learning.

- Bob

[Julian Brooke]

On Monday, October 8, 2012 7:44:59 PM UTC-4, Bob Carpenter wrote: 

You might be happier with Hal Daume's HMC, which in part
is tailored to not only LDA but some process generalizations
thereof:

http://www.umiacs.umd.edu/~hal/HBC/

Yeah, actually tried HBC already, got it complied with some effort but ran into a completely inexplicable error trying to run it, and decided I didn't want to bother anyone about 4 year old code. Maybe I should reconsider that, but after doing a couple of reads of Matt's online Variational Bayes paper, I'm now thinking that adapting his python code for online LDA might be the sweet spot, i.e. code that I can hope to make sense of (yay python!), modify to my needs (no collapsing required, as far as I can tell), and actually finish running on a big web corpus before the heat death of the universe. Of course, a few days ago Stan looked like the sweet spot, so what do I know? I'll give it a try, anyway. Thanks again.

 

If by T you mean Toronto, then you've got lots of
in-house expertise with sampling and approximate
sampling.  But not being a statistics expert may make it
hard to communicate.  I know it was hard for me when
making the transition from CS/machine learning.

 

Indeed. I talked to Radford Neal once shortly after I arrived, thoroughly embarrassed myself, and I haven't seen him since...

Julian 

[Bentley Coffey]

Bob and Matt,

Do you still recommend against using Stan to estimate a LDA model. The reference manual (version 2.2.0) contains on page 139 the following Stan code for performing LDA:

data {

int<lower=2> K; // num topics

int<lower=2> V; // num words

int<lower=1> M; // num docs

int<lower=1> N; // total word instances

int<lower=1,upper=V> w[N]; // word n

int<lower=1,upper=M> doc[N]; // doc ID for word n

vector<lower=0>[K] alpha; // topic prior

vector<lower=0>[V] beta; // word prior

}

parameters {

simplex[K] theta[M]; // topic dist for doc m

simplex[V] phi[K]; // word dist for topic k

}

model {

for (m in 1:M)

theta[m] ~ dirichlet(alpha); // prior

for (k in 1:K)

phi[k] ~ dirichlet(beta); // prior

for (n in 1:N) {

real gamma[K];

for (k in 1:K)

gamma[k] <- log(theta[doc[n],k]) + log(phi[k,w[n]]);

increment_log_prob(log_sum_exp(gamma)); // likelihood

}

}

It seems to me that this code is missing a crucial piece of the likelihood, which appears as the last term in your linked source write-up of the collapsed likelihood: 

And also on the bottom of page 138 of the Stan reference manual... Also noteworthy is that this could be re-written to be more efficient by using sufficient statistics (count of each vocabulary word in each document, obviating the next for stacking a ragged array). 

I'm aware of the multimodality issues but was planning to overcome that labeling-like identification concern with prior knowledge (I happen to know the 50 different topic names, so I know that some of the words in the vocabulary should be more probable for particular topics,eg the words "oil" and "gas" should be more likely in the topic "oil and gas"). Should I really abandon Stan as my platform of choice for this kind of application?

Bentley

[Bob Carpenter]

There wasn't anything after the ":" ? What's missing? It
looks like it matches what we have in the doc.

> And also on the bottom of page 138 of the Stan reference manual... Also noteworthy is that this could be re-written to be more efficient by using sufficient statistics (count of each vocabulary word in each document, obviating the next for stacking a ragged array).

For cases where the word usage is dense and the counts not mostly ones,
this would indeed be much more efficient. I'm used to dealing with
natural language data with LDA which doesn't have this property.

Nevertheless, calculating sufficient statistics on a per-doc level
in a sparse way could still be worthwhile as an efficiency boost.

> I'm aware of the multimodality issues but was planning to overcome that labeling-like identification concern with prior knowledge (I happen to know the 50 different topic names, so I know that some of the words in the vocabulary should be more probable for particular topics,eg the words "oil" and "gas" should be more likely in the topic "oil and gas"). Should I really abandon Stan as my platform of choice for this kind of application?

Stan's good for exploring variant model structures and running one-off
model fits. If you need something running at scale in a production environment,
you probably want to use some kind of optimization tailored for your problem,
such as a variational approach.

- Bob

> --
> You received this message because you are subscribed to the Google Groups "Stan users mailing list" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to stan-users+...@googlegroups.com.
> For more options, visit https://groups.google.com/d/optout.

[Bentley Coffey]

Sorry, my attempt to paste an image of the math somehow failed so I'll try to write it in plain text...

The likelihood on the bottom of page 138 is:

ln p(theta,phi|w,alpha,beta = sum_{m=1}^M ln (Dirichlet(theta_m|alpha)) + sum_{k=K}^M ln (Dirichlet(phi_k|beta)) + sum_{m=1}^M sum_{n=1}^{N[m]} ln (sum_{z=1}^K Categorical(z|theta_m) * Categorical(w_{m,n}|phi_z))

which reassuringly matches the math in the link that you gave in your 2012 post (http://lingpipe-blog.com/2011/04/14/summing-out-topic-assignments-ld/)

Yet in the Stan code provided in the reference manual, unless I'm missing something, the likelihood constructs the first two pieces and not the third (suppressing summations):

ln (Dirichlet(theta_m|alpha)) <=> log(theta[doc[n],k])

ln (Dirichlet(phi_k|beta)) <=> log(phi[k,w[n]]) 

ln (sum_{z=1}^K Categorical(z|theta_m) * Categorical(w_{m,n}|phi_z)) <=> ???

Bentley

You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/3aCWlCguL_s/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Bob Carpenter]

The increment_log_prob(

On Aug 28, 2014, at 2:46 PM, Bentley Coffey <bentley...@gmail.com> wrote:

> Sorry, my attempt to paste an image of the math somehow failed so I'll try to write it in plain text...
>
> The likelihood on the bottom of page 138 is:
>
> ln p(theta,phi|w,alpha,beta = sum_{m=1}^M ln (Dirichlet(theta_m|alpha)) + sum_{k=K}^M ln (Dirichlet(phi_k|beta)) + sum_{m=1}^M sum_{n=1}^{N[m]} ln (sum_{z=1}^K Categorical(z|theta_m) * Categorical(w_{m,n}|phi_z))
>
> which reassuringly matches the math in the link that you gave in your 2012 post (http://lingpipe-blog.com/2011/04/14/summing-out-topic-assignments-ld/)
>
> Yet in the Stan code provided in the reference manual, unless I'm missing something, the likelihood constructs the first two pieces and not the third (suppressing summations):
>
> ln (Dirichlet(theta_m|alpha)) <=> log(theta[doc[n],k])
> ln (Dirichlet(phi_k|beta)) <=> log(phi[k,w[n]])
> ln (sum_{z=1}^K Categorical(z|theta_m) * Categorical(w_{m,n}|phi_z)) <=> ???

That last line is translated as follows:

> for (k in 1:K)
> gamma[k] <- log(theta[doc[n],k]) + log(phi[k,w[n]]);
> increment_log_prob(log_sum_exp(gamma));

The increment_log_prob() adds the term to the total log probability,
and log_sum_exp(v) =def= log(SUM_n exp(v[n])).

- Bob

[Matt Hoffman]

Stan can be used to fit LDA. But you may be disappointed if you compare its ability to find good parameters quickly with a standard collapsed Gibbs sampler.

CGS is a really surprisingly good stochastic search algorithm for LDA—it seems to be much easier to optimize the log-likelihood in "z-space" (the word-topic assignments) than in "theta-phi-space" (the topics and topic weights). One reason is sparsity—if you use Dirichlet priors with small hyperparameters and sum out the z variables (which you need to in Stan), then you wind up with a highly nonconvex sparsity-inducing prior. Gradient-based optimization on these kinds of nonconvex problems is very, very hard, and since NUTS relies very heavily on the use of gradient information to tell it where to explore, it's challenging for it to figure out what topics should be active in what documents. (Although in principle it'll figure it out eventually.)

What NUTS is good at is mixing quickly and delivering large numbers of effectively independent samples. Nobody really knows what the mixing speed of CGS is like, and most people don't care—for most topic modeling applications people aren't really interested in quantifying the uncertainty in their parameter estimates.

Matt

[Janne Sinkkonen]

We have been using hca with good results:

https://github.com/wbuntine/topic-models

http://mloss.org/software/view/527/

http://delivery.acm.org/10.1145/2630000/2623691/p881-buntine.pdf

I would vote for it if one wants to estimate LDA or its nonparametric variants by collapsed Gibbs. 

[Bentley Coffey]

Bob: yes, I see what I was missing now... thanks!

Matt: what you've written makes complete sense to me... but...

Unfortunately, the implementations of CGS available in r packages (lda and topicmodels) seem to fall rather short of my needs. I'm a researcher, so uncertainty matters to me because producing just a point estimate wouldn't be acceptable for publication (for good reason, IMHO). It looks like I could "hack" either of those packages to call them repeatedly but the mixing across such iterations may just be so poor that NUTS would be faster. More importantly, neither lda or topicmodels is designed to accommodate prior information on how the word distribution differs across topics -- both have just a single parameter for beta (as well as alpha) but my prior knowledge is in how beta differs across topics. Recall that this was my strategy for navigating the multi-modality problems that Bob mentioned (I should be able to keep increasing the strength of my priors until the posterior becomes unimodal). To give you more specific detail on my particular problem: my purpose is to estimate the share of each topic discussed in each of 50 documents. I have 50 to 100 pre-determined topics (they are industry names taken from NAICS codes, such as "Oil and Gas Extraction") to which I add an additional topic I call "misc" for none of the above. From these topic descriptions, I have constructed my 1000 to 4000 word vocabulary (e.g. it contains "oil", "gas", "extraction", "oil extraction", and "gas extractrion"). In setting the prior on beta, I can make it so that the word "oil" is more probable than other words when the topic is "oil and gas extraction" (I might have to experiment with how strong it needs to be to ensure that the term "oil" is most associated with "oil and gas extraction", which is what I meant by making the posterior unimodal). So, unless there is some way to rehabilitate the lda and/or topicmodels packages, calling stan from r appears to be my only hope short of a radical departure from the comfort of my r platform... well, unless there is some model other than LDA that would be more appropriate here (as an economist, I'm somewhat new to topic models and LDA seemed to be sufficiently appropriate and accepted as a standard approach but I would certainly welcome a better suggestion). Unfortunately, after letting my machine crunch overnight (approximately 12 hours), stan had made 0 progress on my problem (still on the first iteration). Here's my stan code (which is slightly modified, from the version in manual that Bob and I discussed, to accommodate a different prior on beta for each topic):

data{

  int <lower = 2> T;                // Num Topics

  int <lower = 2> V;                // Num Words in Vocabulary

  int <lower = 1> M;                // Num Docs

  int <lower = 0> K[V,M];           // search term counts in each of the M docs

  vector <lower = 0> [T] alpha[M];  // topic prior

  vector <lower = 0> [V] beta[T];   // word prior

}

parameters{

  simplex [T] theta[M]; // topic distribution for each of the M docs

  simplex [V] phi[T]; // word distribution for each of the T topics

}

model{

  for (m in 1:M)

    theta[m] ~ dirichlet(alpha[m]);  //prior

  for (t in 1:T)

    phi[t] ~ dirichlet(beta[t]);     //prior

  for (v in 1:V) {

    for (m in 1:M) {

      real gamma [T];

      for (t in 1:T)

        gamma[t] <- log(theta[m,t]) + log(phi[t,v]); //likelihood

      increment_log_prob(K[v,m]*log_sum_exp(gamma));

    }

  }

}

Here's my r code that calls stan:

library(rstan)

library(reshape2)

##### Import data files

Search.Term.Counts <- read.csv("string_title_year_count.csv")

Industry.Terms <- read.csv("search_terms_and_industries.csv")

# data cleaning: restrict to year 2012 and remove a few negative counts (probably due to a coding mistake upstream)

Search.Term.Counts.t <- subset(Search.Term.Counts,year==2012)

Search.Term.Counts.t$count <- ifelse(Search.Term.Counts.t$count < 0, 0, Search.Term.Counts.t$count)

# extract Num.Words.In.Vocabulary and Num.Docs from Search.Term.Counts 

Num.Words.In.Vocabulary <- length(levels(Search.Term.Counts.t$string))

Num.Docs <- length(levels(as.factor(Search.Term.Counts.t$title)))

# extract Num.Topics [= Num.Industries + 1] from Industry.Terms

Num.Topics <- length(levels(as.factor(Industry.Terms$code)))+1

# extract Search.Term.Counts

melted.Search.Term.Counts <- melt(Search.Term.Counts.t[,c("title","string","count")], id=c("title", "string"))

Search.Term.Counts.Array <- acast(melted.Search.Term.Counts, string ~ title)

# Create the simply flat topic prior 

# alpha here represents the virtual count of times we've seen this topic for this doc (akin to the beta distribution's parameters so that 1 is a flat uniform)

# Steyvers n Griffiths (2007) suggest 50/T, which is actually around 1 here anyway

Topic.Prior <- matrix(1,nrow=Num.Topics, ncol=Num.Docs)

rownames(Topic.Prior) <- c(levels(as.factor(Industry.Terms$code)),"000")

colnames(Topic.Prior) <- levels(as.factor(Search.Term.Counts.t$title))

# Create the topic prior which is flat but to which I'll add additional weight allocated for search terms known to be associated with particular industries

# beta here represents the virtual count of times we've seen this word used in this topic 

# Steyvers n Griffiths (2007) suggest 0.01, which is improper, but they acknowledge that it should depend on vocabulary, num docs, etc... so I'll stick with 1 that corresponds to a flat uniform

Word.Prior <- matrix(1,nrow=Num.Words.In.Vocabulary, ncol=Num.Topics)

rownames(Word.Prior) <- levels(Search.Term.Counts.t$string)

colnames(Word.Prior) <- c(levels(as.factor(Industry.Terms$code)),"000")

# add information to the word prior

# loop thru each industry (not including the misc topic category)

# find the maximum number of words in that industry's search terms

# loop thru that industry's search terms to appropriately increment the prior for cases of individual terms in a search string of multiple words:

# the individual word should be more associated with misc and the whole search string of multiple words should be more associated with the industry 

for (Current.Industry.Index in 1:(ncol(Word.Prior)-1)) {

  Industry.Code <- colnames(Word.Prior)[Current.Industry.Index]

  Current.Industry.Terms <- subset(Industry.Terms,code == Industry.Code)

  Current.Industry.Terms$Num.Terms <- sapply(gregexpr("\\W+", Current.Industry.Terms$string), function(x) sum(x>0) ) + 1

  Max.Num.Terms <- max(Current.Industry.Terms$Num.Terms)

  for (Current.String.Index in 1:nrow(Current.Industry.Terms)) {

    if (as.character(Current.Industry.Terms$string[Current.String.Index]) %in% rownames(Word.Prior)) {

      Word.Prior[as.character(Current.Industry.Terms$string[Current.String.Index]),Industry.Code] <- Word.Prior[as.character(Current.Industry.Terms$string[Current.String.Index]),Industry.Code] + Current.Industry.Terms$Num.Terms[Current.String.Index]

      Word.Prior[as.character(Current.Industry.Terms$string[Current.String.Index]),"000"] <- Word.Prior[as.character(Current.Industry.Terms$string[Current.String.Index]),"000"] + (Max.Num.Terms - Current.Industry.Terms$Num.Terms[Current.String.Index])

    }

  }

}

##########################################################################################

# Analyze with Stan

##########################################################################################

# reshaping the data as required by Stan -- packed into a list: counts as a VxM array, topic prior as a TxM array, and word prior as a VxM array

StanInputData <- list(T = Num.Topics, V = Num.Words.In.Vocabulary, M = Num.Docs, K = Search.Term.Counts.Array, alpha = t(Topic.Prior), beta = t(Word.Prior))

# run stan on the file, compile to C [stanc()], and iteratively draw samples from the posterior [sampling()]

stan.fit <- stan(file="LDA.stan", data=StanInputData, seed = 12345, iter=1000, chains=4)

--
You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/3aCWlCguL_s/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Bob Carpenter]

On Aug 31, 2014, at 11:09 PM, Bentley Coffey <bentley...@gmail.com> wrote:

> Bob: yes, I see what I was missing now... thanks!
> Matt: what you've written makes complete sense to me... but...
>
> Unfortunately, the implementations of CGS available in r packages (lda and topicmodels) seem to fall rather short of my needs. I'm a researcher, so uncertainty matters to me because producing just a point estimate wouldn't be acceptable for publication (for good reason, IMHO).

Using LDA for anything other than exploratory data analysis
seems dangerous. What do you think the uncertainty is going
to be characterizing that isn't highly dependent on your prior
strength in this case?

> It looks like I could "hack" either of those packages to call them repeatedly but the mixing across such iterations may just be so poor that NUTS would be faster.

It's pretty much hopeless trying to mix the multiple modes of LDA.

> More importantly, neither lda or topicmodels is designed to accommodate prior information on how the word distribution differs across topics -- both have just a single parameter for beta (as well as alpha) but my prior knowledge is in how beta differs across topics. Recall that this was my strategy for navigating the multi-modality problems that Bob mentioned (I should be able to keep increasing the strength of my priors until the posterior becomes unimodal).

It still won't be convex, but with enough prior strength
you can push down the modes you don't want. You still probably
won't get proper exploration of the posterior.

> To give you more specific detail on my particular problem: my purpose is to estimate the share of each topic discussed in each of 50 documents. I have 50 to 100 pre-determined topics (they are industry names taken from NAICS codes, such as "Oil and Gas Extraction") to which I add an additional topic I call "misc" for none of the above. From these topic descriptions, I have constructed my 1000 to 4000 word vocabulary (e.g. it contains "oil", "gas", "extraction", "oil extraction", and "gas extractrion"). In setting the prior on beta, I can make it so that the word "oil" is more probable than other words when the topic is "oil and gas extraction" (I might have to experiment with how strong it needs to be to ensure that the term "oil" is most associated with "oil and gas extraction", which is what I meant by making the posterior unimodal).

> So, unless there is some way to rehabilitate the lda and/or topicmodels packages, calling stan from r appears to be my only hope short of a radical departure from the comfort of my r platform... well, unless there is some model other than LDA that would be more appropriate here (as an economist, I'm somewhat new to topic models and LDA seemed to be sufficiently appropriate and accepted as a standard approach but I would certainly welcome a better suggestion).

It depends what you're trying to do. If you are assuming each document
is made up out of several topics, then LDA's a reasonable model.

But you have to be careful with the prior on the document topic proportions
(theta in your model and the original LDA papers). The strength of that prior
will determine how many topics get non-negligible probability in a document,
and this in turn plays off against the topic distributions themselves.

> Unfortunately, after letting my machine crunch overnight (approximately 12 hours), stan had made 0 progress on my problem (still on the first iteration).

Do you know if you ran out of memory? I'd think this would be slow,
but not 12 hours per iteration slow.

> Here's my stan code (which is slightly modified, from the version in manual that Bob and I discussed, to accommodate a different prior on beta for each topic):
>
> data{
> int <lower = 2> T; // Num Topics
> int <lower = 2> V; // Num Words in Vocabulary
> int <lower = 1> M; // Num Docs
> int <lower = 0> K[V,M]; // search term counts in each of the M docs
> vector <lower = 0> [T] alpha[M]; // topic prior

Should alpha really vary by document? Do you have prior information
on how many topics area in each document?

Given your strong priors, you could try initializing theta and phi
to the prior means. alpha / sum(alpha) and beta / sum(beta).

> vector <lower = 0> [V] beta[T]; // word prior
> }
> parameters{
> simplex [T] theta[M]; // topic distribution for each of the M docs
> simplex [V] phi[T]; // word distribution for each of the T topics
> }
> model{
> for (m in 1:M)
> theta[m] ~ dirichlet(alpha[m]); //prior
> for (t in 1:T)
> phi[t] ~ dirichlet(beta[t]); //prior
> for (v in 1:V) {
> for (m in 1:M) {
> real gamma [T];
> for (t in 1:T)
> gamma[t] <- log(theta[m,t]) + log(phi[t,v]); //likelihood
> increment_log_prob(K[v,m]*log_sum_exp(gamma));
> }
> }
> }

Presumably what you have is very sparse. So I doubt this will work well.

You might be able to define a collapsed Gibbs sampler for your
model in JAGS. Or just write the Gibbs sampler in R directly --- it's
a pretty easy one.

> Here's my r code that calls stan:

I don't speak R well enough to get through the rest.

- Bob

[Bentley Coffey]

Mechanically, other than the prior, the uncertainty must come from the data's likelihood. What exactly in the data describes each model parameter is always tough to decipher, particularly with mixtures of distributions... but here goes: the density of each word across documents resembles a mixture of exponential densities and there is a clear covariation in word density across documents (i.e. documents where "gas" appears quite frequently also have "oil" appearing quite frequently); this is what I expect to identify the theta parameters and hence their uncertainty should come from how much the data varies around that signal...

I'm rather certain that each document includes more than 1 topic, although perhaps not every topic appears in each document... So, it still seems to me that LDA is the most appropriate model among those of which I'm aware...

I didn't run out of memory. This takes up about 1GB of memory, of which I still have 10GB free. Interestingly, when I look at the CPU usage of the core handling the process, it reads as only 80% busy...

Theta is the parameter of interest for me; the beta parameters are basically just "nuisance" parameters for my purposes. I want the prior on theta to be completely flat, hence me setting alpha = 1 for all topics. In response to your question, I have revised my model construction in Stan so that the value of alpha fed to Stan no longer varies across documents. I'm not sure why that construction would save computation time since either way alpha = 1 for all topics and documents...

The other change that I made was that initialization that you suggested. Perhaps this produced a serious speed-up? I'm not sure... but, somehow, it ran the first 20 iterations in about 20 minutes. The next 20 iterations took about an hour... and it's still chugging along now...

Yeah, I may have to resort to coding the [fast] CGS myself... or grab someone else's C code and embed it in-line within my r script...

> CGS is a really surprisingly good stochastic search algorithm for LDA--it seems to be much easier to optimize the log-likelihood in "z-space" (the word-topic assignments) than in "theta-phi-space" (the topics and topic weights). One reason is sparsity--if you use Dirichlet priors with small hyperparameters and sum out the z variables (which you need to in Stan), then you wind up with a highly nonconvex sparsity-inducing prior. Gradient-based optimization on these kinds of nonconvex problems is very, very hard, and since NUTS relies very heavily on the use of gradient information to tell it where to explore, it's challenging for it to figure out what topics should be active in what documents. (Although in principle it'll figure it out eventually.)
>
> What NUTS is good at is mixing quickly and delivering large numbers of effectively independent samples. Nobody really knows what the mixing speed of CGS is like, and most people don't care--for most topic modeling applications people aren't really interested in quantifying the uncertainty in their parameter estimates.

[Matt Hoffman]

> I'm a researcher, so uncertainty matters to me because producing just a point estimate wouldn't be acceptable for publication (for good reason, IMHO).

Given that mixing across the combinatorially explosive number of modes for a model like LDA is impossible, any estimate of posterior uncertainty you wind up with is likely to be at best a very loose lower bound. As Bob said, the source of this multimodality is not just label-switching, but also the sparsity induced by the Dirichlet distribution. Reporting posterior-based confidence intervals may therefore be more misleading than not including any confidence intervals at all.

Perhaps you can describe the specific question you're trying to answer and get uncertainty estimates about?

Theta is the parameter of interest for me; the beta parameters are basically just "nuisance" parameters for my purposes. I want the prior on theta to be completely flat, hence me setting alpha = 1 for all topics.

FYI, this isn't an "uninformative" prior. Like all priors on the simplex, it corresponds to a very specific set of assumptions about how many topics you expect to see per document.

A less rigid approach would be to learn the prior on theta, as suggested by Wallach et al. in "Rethinking LDA: Why Priors Matter":

https://people.cs.umass.edu/~wallach/publications/wallach09rethinking.pdf

Matt

[Bentley Coffey]

You're right, I'd want to be clear about that. I thought that my prior scheme would also rule out the sparsity-induced multi-modality. If every Dirichlet parameter is greater than or equal to 1 then it should be unimodal in the simplex. Again, I could be missing something...

Matt, I sent you a private email about my particular application, which is unlikely to be of interest to the broader user group...

Yeah, if the Dirichlet parameters are all equal to 1, then the prior is uniform in the simplex. Such an assumption increases the appeal of Bayesian work to an audience dominated by frequentists because it makes the posterior proportional to the data's likelihood (with which they are comfortable). However, you're right that any prior corresponds to a very specification set of assumptions about how many topics we'd expect to see (a priori) in a document. Hence, flat isn't exactly uninformative, which can be seen via a reparameterization of the model. 

I found the approach proposed by Wallach et al. (2009), ie optimization with respect to the prior parameters, to be rather appealing; however, I didn't immediately understand the intuition behind their results. I would have thought that the (abstract) symmetry in the problem would result in the symmetric-asymmetric prior being just as good as the asymmetric-symmetric prior. That said, if I tried to wholely adopt their recommended approach, then I would lose the ability to incorporate my prior information (which, intuitively, is the whole purpose of a prior) and I may not be able to cleanly assign the pre-determined topic names to the latent topics generated... By the way, I noticed that this Wallach et al. (2009) approach had been implemented in the c code that Janne Sinkkonen linked to (hca: http://mloss.org/software/view/527/)... that site also references a new paper (Li et al. 2014) which seems very promising for speeding up LDA (see: http://www.sravi.org/pubs/fastlda-kdd2014.pdf)...

A last note, for those of you keeping score, Stan is generating samples at a rate of about 10 every 2 hours...

--

[Bentley Coffey]

You're right, I'd want to be clear about that. I thought that my prior scheme would also rule out the sparsity-induced multi-modality. If every Dirichlet parameter is greater than or equal to 1 then it should be unimodal in the simplex. Again, I could be missing something...

Matt, I sent you a private email about my particular application, which is unlikely to be of interest to the broader user group...

Yeah, if the Dirichlet parameters are all equal to 1, then the prior is uniform in the simplex. Such an assumption increases the appeal of Bayesian work to an audience dominated by frequentists because it makes the posterior proportional to the data's likelihood (with which they are comfortable). However, you're right that any prior corresponds to a very specification set of assumptions about how many topics we'd expect to see (a priori) in a document. Hence, flat isn't exactly uninformative, which can be seen via a reparameterization of the model. 

I found the approach proposed by Wallach et al. (2009), ie optimization with respect to the prior parameters, to be rather appealing; however, I didn't immediately understand the intuition behind their results. I would have thought that the (abstract) symmetry in the problem would result in the symmetric-asymmetric prior being just as good as the asymmetric-symmetric prior. That said, if I tried to wholely adopt their recommended approach, then I would lose the ability to incorporate my prior information (which, intuitively, is the whole purpose of a prior) and I may not be able to cleanly assign the pre-determined topic names to the latent topics generated... By the way, I noticed that this Wallach et al. (2009) approach had been implemented in the c code that Janne Sinkkonen linked to (hca: http://mloss.org/software/view/527/)... that site also references a new paper (Li et al. 2014) which seems very promising for speeding up LDA (see: http://www.sravi.org/pubs/fastlda-kdd2014.pdf)...

A last note, for those of you keeping score, Stan is generating samples at a rate of about 10 every 2 hours...

--

[Bob Carpenter]

On Sep 1, 2014, at 12:11 PM, Bentley Coffey <bentley...@gmail.com> wrote:

> Mechanically, other than the prior, the uncertainty must come from the data's likelihood. What exactly in the data describes each model parameter is always tough to decipher, particularly with mixtures of distributions... but here goes: the density of each word across documents resembles a mixture of exponential densities and there is a clear covariation in word density across documents (i.e. documents where "gas" appears quite frequently also have "oil" appearing quite frequently); this is what I expect to identify the theta parameters and hence their uncertainty should come from how much the data varies around that signal...

I'd suggest doing a sensitivity analysis with different priors.
They make a huge difference in terms of the implicit correlation
among words being modeled. High count priors in alpha spread
out each document, which in effect shrinks all the topic models
and lessens the correlation effect in the posterior.

LDA's pretty good at modeling the burstiness of words across
documents, but you'll need to tune the priors heavily if you
expect it to be well calibrated to new data in the language
modeling sense of making good predictions (high probability) for
new documents.

> I'm rather certain that each document includes more than 1 topic, although perhaps not every topic appears in each document... So, it still seems to me that LDA is the most appropriate model among those of which I'm aware...
>
> I didn't run out of memory. This takes up about 1GB of memory, of which I still have 10GB free. Interestingly, when I look at the CPU usage of the core handling the process, it reads as only 80% busy...
>
> Theta is the parameter of interest for me; the beta parameters are basically just "nuisance" parameters for my purposes. I want the prior on theta to be completely flat, hence me setting alpha = 1 for all topics.

I often see people use alpha << 1 to avoid the following issue.

Suppose you have 100 words in a document and 100 topics.
And you see 10 words in one topic, 30 in a second, and
60 in a third. Then you get 11/200, 31/200, and 61/200
posterior mean probabilities and about half the posterior
mass assigned to the 97 "out-of-topic" topics.

> In response to your question, I have revised my model construction in Stan so that the value of alpha fed to Stan no longer varies across documents. I'm not sure why that construction would save computation time since either way alpha = 1 for all topics and documents...

It won't save time. I just wondered why it was varying.

> The other change that I made was that initialization that you suggested. Perhaps this produced a serious speed-up?

It can if a random init is way out in the tails.

> I'm not sure... but, somehow, it ran the first 20 iterations in about 20 minutes. The next 20 iterations took about an hour... and it's still chugging along now...

It goes over a hump in the early iterations when it gets serious about
setting both the mass matrix and the step size. Then it should speed up.

> Yeah, I may have to resort to coding the [fast] CGS myself... or grab someone else's C code and embed it in-line within my r script...

You probably don't even need to do it in C. But
you'd need to customize it to deal with varying priors.

- Bob

[Stefano]

Hello,

Would this be a way to implement such sufficient statistics?

  data {

    int<lower=2> K;               // num topics

    int<lower=2> W;               // num words

    int<lower=1> M;               // num docs

    int<lower=1> N;               // total word instances //UNIQUE FOR EVERY DOCUMENT/WORD

    int<lower=0> tot_w_count;     //upper value for word counts

    vector<lower=0, upper=tot_w_count>[N] w_count;   // word count per document (e.g., cat in document 1 is present 5 times)

    int<lower=1,upper=W> w[N];    // word n //UNIQUE FOR EVERY DOCUMENT/WORD

    int<lower=1,upper=M> doc[N];  // doc ID for word n

    vector<lower=0>[K] alpha;     // topic prior

    vector<lower=0>[K] beta;     // word prior

    

  }

  parameters {

    simplex[K] theta[M];   // topic dist for doc m     

    simplex[W] phi[K];            // word dist for topic k

  } 

  model {

    for (m in 1:M) theta[m] ~ dirichlet(alpha);  // prior

    for (k in 1:K) phi[k] ~ dirichlet(beta);  // prior

    for (n in 1:N) {

      real gamma[K];

      for (k in 1:K) gamma[k] <- (log(theta[doc[n],k]) + log(phi[k,w[n]]));

      increment_log_prob(log_sum_exp(gamma) * w_count[n]);  // likelihood

    }

  }

Thanks.

[Bob Carpenter]

I think that should work. Did you try it?

It only gives you as much speedup as there is
duplication of words.

And you should only compute log(theta) and log(phi)
once --- you can define them as local variables once and
then reuse them.

And I think the likelihood can be simplified by using
vector addition on a slightly rearranged log_theta and
log_phi.

> for (n in 1:N) {
> real gamma[K];
> for (k in 1:K) gamma[k] <- (log(theta[doc[n],k]) + log(phi[k,w[n]]));
> increment_log_prob(log_sum_exp(gamma) * w_count[n]); // likelihood
> }

model {
vector[K] log_phi[W];
vector[K] log_theta[M];
... define these locals once ...

for (n in 1:N)

increment_log_prob(w_count[n] * log_sum_exp(log_theta[m] + log_phi[w[n]]));

I don't see how to do much better than that.

- Bob

> To post to this group, send email to stan-...@googlegroups.com.

[Stefano]

On Monday, January 11, 2016 at 7:09:18 PM UTC+11, Bob Carpenter wrote:

I think that should work.  Did you try it?

Yes it seems to work
 

It only gives you as much speedup as there is
duplication of words.

In some applications such as biostatistics, the number of words (RNA sequences) count can be in the order of 10^4

And you should only compute log(theta) and log(phi)
once --- you can define them as local variables once and
then reuse them.

And I think the likelihood can be simplified by using
vector addition on a slightly rearranged log_theta and
log_phi.

>     for (n in 1:N) {
>       real gamma[K];
>       for (k in 1:K) gamma[k] <- (log(theta[doc[n],k]) + log(phi[k,w[n]]));
>       increment_log_prob(log_sum_exp(gamma) * w_count[n]);  // likelihood
>     }


model {
  vector[K] log_phi[W];
  vector[K] log_theta[M];
  ... define these locals once ...

  for (n in 1:N)
    increment_log_prob(w_count[n] * log_sum_exp(log_theta[m] + log_phi[w[n]]));

Thanks a lot, I will apply these modifications

[Bob Carpenter]

> On Jan 11, 2016, at 4:56 AM, Stefano <mangiol...@gmail.com> wrote:
>
>
>
> On Monday, January 11, 2016 at 7:09:18 PM UTC+11, Bob Carpenter wrote:
> I think that should work. Did you try it?
>
> Yes it seems to work
>
>
> It only gives you as much speedup as there is
> duplication of words.
>
> In some applications such as biostatistics, the number of words (RNA sequences) count can be in the order of 10^4

Other than for playing around with models, I'd think a custom
LDA fitter would make more sense. It's not like you can do
real Bayesian inference with LDA anyway --- the posterior is
combinatorially intractable.

But the speedup from aggregation will depend on how many times
"words" are duplicated in "documents". If each word occurs once,
it may even be slower.

- Bob

[stefano]

Stefano Mangiola

Data analyst @ 

Melbourne University 

Royal Melbourne Hospital - WEHI

+61 (0)450135954

2016-01-13 7:14 GMT+11:00 Bob Carpenter <ca...@alias-i.com>:

> On Jan 11, 2016, at 4:56 AM, Stefano <mangiol...@gmail.com> wrote:
>
>
>
> On Monday, January 11, 2016 at 7:09:18 PM UTC+11, Bob Carpenter wrote:
> I think that should work.  Did you try it?
>
> Yes it seems to work
>
>
> It only gives you as much speedup as there is
> duplication of words.
>
> In some applications such as biostatistics, the number of words (RNA sequences) count can be in the order of 10^4

Other than for playing around with models, I'd think a custom
LDA fitter would make more sense.  It's not like you can do
real Bayesian inference with LDA anyway --- the posterior is
combinatorially intractable.

So, are you saying that Stan, BUGS, JAGS are not suitable for any LDA implementation? What could be an example of a custom

LDA fitter?

 

You received this message because you are subscribed to a topic in the Google Groups "Stan users mailing list" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/stan-users/3aCWlCguL_s/unsubscribe.
To unsubscribe from this group and all its topics, send an email to stan-users+...@googlegroups.com.

[Bob Carpenter]

> On Jan 12, 2016, at 8:58 PM, stefano <mangiol...@gmail.com> wrote:

> Other than for playing around with models, I'd think a custom
> LDA fitter would make more sense. It's not like you can do
> real Bayesian inference with LDA anyway --- the posterior is
> combinatorially intractable.

...

> So, are you saying that Stan, BUGS, JAGS are not suitable for any LDA implementation? What could be an example of a custom
> LDA fitter?

What I'm saying is that you can't do Bayesian inference
for LDA. Nobody can do Bayesian inference for LDA --- there
are too many modes in the posterior.

The best an MCMC algorithm can do for LDA is explore around a single
mode. So if you start multiple chains, they'll never converge to
the same answer with realistic data (at least any kind of natural
language data I've ever seen). Now there's some stability in
the major clusters found, as analyzed by Griffiths and Steyvers
for their Gibbs sampler (which marginalizes out the continuous
parameters).

So what I'm saying is that you might as well use something like
variational inference and be content with a single mode. Or
just a point estimate. You'll find a super-scalable stochastic
variational inference implemented by Matt Hoffman (who also invented
the NUTS algorithm we use in Stan's MCMC) in Vowpal Wabbit.
Here's the writeup:

http://jmlr.csail.mit.edu/papers/v14/hoffman13a.html

- Bob