Slow mixing for Latent Dirichlet allocation

Markov chain Monte Carlo (MCMC) algorithms are ubiquitous in probability theory in general and in machine learning in particular. A Markov chain is devised so that its stationary distribution is some probability distribution of interest. Then one samples from the given distribution by running the Markov chain for a"long time"until it appears to be stationary and then collects the sample. However these chains are often very complex and there are no theoretical guarantees that stationarity is actually reached. In this paper we study the Gibbs sampler of the posterior distribution of a very simple case of Latent Dirichlet Allocation, an attractive Bayesian unsupervised learning model for text generation and text classification. It turns out that in some situations, the mixing time of the Gibbs sampler is exponential in the length of documents and so it is practically impossible to properly sample from the posterior when documents are sufficiently long.


Introduction
Markov chain Monte Carlo (MCMC) is a powerful tool for sampling from a given probability distribution on a very large state space, where direct sampling is difficult, in part because of the size of the state space and in part because of normalizing constants that are difficult to compute.
In machine learning in particular, MCMC algorithms are extremely common for sampling from posterior distributions of Bayesian probabilistic models. The posterior distribution given observed data is then difficult to sample from for the reasons above. One then designs an (irreducible aperiodic) Markov chain whose stationary distribution is precisely the targeted posterior. This is usually fairly easy since the posterior is usually easy to compute up to the normalizing constant (the denominator in Bayes formula). One usually uses Gibbs sampling or the related Metropolis-Hastings algorithm.
Gibbs sampling in these situations can be described as follows. Our state space is a finite set of random variables X = {X a } a∈A , where X a ∈ T for some measurable space T , so that X ∈ T A and the targeted distribution is a given probability measure P on T A . In order to sample from P one starts a Markov chain on T A whose updates are given by first choosing an index a ∈ A at random and then choosing a new value of X a according to the conditional distribution of X a given all X b , b ∈ A \ {a}. Under mild conditions, this chain converges in distribution to P. The Markov chain is then run for a "long time" and then a sample, hopefully approximately from P, is collected. A key question here is for how long the chain actually has to be run, in order for the distribution after that time to be a good approximation of P. Since A is usually large, the number of steps needed should at least be no more than polynomial in the size of A for Gibbs sampling to be feasible. In almost all practical cases, the structure of the sample space and the probability measure P is so complex that is virtually impossible to make a rigorous analysis of the mixing rate. However it may be possible to consider some very simplified special cases. In this paper, we will analyse a special case of Latent Dirichlet allocation, henceforth LDA for short, and demonstrate for such a simple special case that mixing can indeed be a problem.
LDA is model used to classify documents according to their topics. We have a large corpus of documents and we want to determine for each word in each document which topic it comes from. Knowing this we can also classify the documents according to the proportion of words of the different topics it contains. The setup in LDA is that one has a fixed number D of documents of lengths N d a fixed set of topics t 1 , t 2 , . . . , t s and a fixed set of words w 1 , w 2 , . . . , w v . These are specified in advance. The number of topics is usually not large, whereas the number of words is. Next, for each document d = 1, . . . , D, a multinomial distribution θ d = (θ d (1), . . . , θ d (s)) over topics is chosen according to a Dirichlet prior with a known parameter α = (α 1 , . . . , α s ). For each topic t i a multinomial distribution φ i = (φ i (1), . . . , φ i (v)) according to a Dirichlet prior with parameter β = (β 1 , . . . , β v ) independently of each other and of the θ d :s. Given these, the corpus is then generated by for each position p = 1, . . . , N d in each document d picking a topic z d,p according to θ d and then picking the word at that position according to φ z d,p , doing this independently for all positions. Note that the LDA model is a so called "bag of words" model, i.e. it is invariant under permutations within each document.
In this paper, we will consider the very simplified case with D = s = v = 2, N 1 = N 2 = m and α = β = (1, 1) and study the mixing time asymptotics as m → ∞. To simplify the notation, denote the two topics by A and B and the two words by 1 and 2. Define n ij as the number of occurrences of the word j in document i, i, j = 1, 2 and write n i. = n i1 + n i2 (which by assumption equals m) and n .j = n 1j + n 2j and n .. = i,j n ij = 2m. We consider the mixing time for Gibbs sampling of the posterior in a seemingly typical case, namely that the number of 1:s in the first document is 3m/10 and in the second document 6m/10. Our result is the following. Theorem 1.1 Consider the case n 11 = 3m/10 and n 21 = 6m/10. Then there exists a λ > 0 such that for each 0 < κ < 3/4, Remark. The point of this paper is to demonstrate that the mixing time issue can be a real problem for Bayesian inference in machine learning in general and not for LDA in particular. There are also other methods for estimation that seem to work well for LDA, in particular variational inference. Furthermore, experimental results seem to be fairly well in line with what one would expect from a topic classification algorithm.
Before moving on to the proof of Theorem 1.1, we formally introduce the concept of mixing time. Let {X t } ∞ t=0 be a discrete time Markov chain on the finite state space S and for s ∈ S, let P s be the underlying probability measure under X 0 = s. Definition 1.2 Let µ and ν be two probability measures on S. The the total variation distance between µ and ν is given by The essence of Theorem 1.1 is that even after an exponentially long time, the distribution of the state of the Gibbs sampler is concentrated on a set whose probability mass according to the targeted posterior is at most 1/4.

Proof of Theorem 1.1
The following two combinatorial lemmas will be needed. Lemma 2.1 Let X be a standard uniform random variable and 0 ≤ k ≤ n. Then Proof. First recall that for any m = 1, 2, . . ., Next observe that the result holds true for any n with k = n. We want to prove that the claimed result holds for (n, k) = (p, r) for some arbitrary 0 ≤ r ≤ p. We do this by induction. We may then assume that the result is true with (n, k) = (m, l) for any m and l with either m < p or m = p and l > r. Then, taking where the fourth equality follows form the induction hypothesis and the last equality holds because Proof. Both sides of the equality are equal to the multinomial coefficient

✷
Let W = (w 11 , w 12 , . . . , w 1m , w 21 , . . . , w 2m ) be the words in our corpus, let Z = (z 11 , z 12 , . . . , z 2m ) be the latent topics and Z d be the latent topics in document d. Let also θ i be the probability that z d1 = A, d = 1, 2 and let φ t be the conditional probability that w dj = 1 given that z dj = t, t = A, B. In the case under study, these four quantities are all independent standard uniform random variables. We begin by determining the posterior distribution.
. Here we have used the notation k dj for the number words in document d with the topic being A and the word being j and the same dot notation for the k:s as for the n:s. The last equality is Lemma 2.1. For the second factor we have analogously, again using Lemma 2.1, Hence, ignoring factors that do not depend on the k:s and using Lemma 2.2 for the second equality This expression only depends on z via k = k(z) := (k 11 , k 12 , k 21 , k 22 ). Identifying all z having the same k(z), we have (regarding ,with some abuse of notation, a k also as the equivalence class consisting of all z:s having that particular k dj :s) that for any k 1 and k 2 , all z ∈ k 1 have the same probability of transitioning into k 2 . Hence the process where we only record the k:s is a lumped Markov chain and it suffices to regard this chain, whose state space is where of course K = k(Z) and C 0 is the normalizing constant.
For a = 0, the right hand side equals − log(11/7) which is smaller than −9/20. Analogously, the partial derivatives of log g with respect to b, c and d are less than −1/4 (we do not get exactly the same numbers here). It follows that for a sufficiently small ǫ > 0, For such a k, let V k be the number of visits to k of the lumped Markov chain starting at k 0 until it returns to k 0 . Then by an intuitive and well known fact in Markov theory (see e.g. [1], Lemma 6 of Chapter 2) Hence By a Bonferroni bound we get for sufficiently large m (ǫm > 170 is sufficient). Hence, starting at k 0 , the probability of getting more than distance ǫm from k 0 in e Taking the logarithm and differentiating with respect to a gives ∂ ∂a log g(a, 0, 0, 0) = log 1 2 − a 2 )(1 − 10a 3 )( 2 3 + 10a 9 )a ( 1 2 + a 2 ) 10a 3 ( 1 3 − 10a 9 )(1 − a) and taking a = 0, this becomes − log(5/3) < −1/2. Repeating the above argument, it follows for some ν > 0 that starting from k 1 , the probability of getting more than distance νm from k 1 in e ν 9 m steps is at most e − ν 72 m . Combining this with the analogous result for k 0 , this proves Theorem 1.1.