P\'olya Urn Latent Dirichlet Allocation: a doubly sparse massively parallel sampler

TL;DR

This paper introduces a novel, faster, and asymptotically exact parallel sampler for LDA that overcomes memory and convergence issues of previous methods, enabling efficient large-scale topic modeling.

Contribution

A new Pólya-urn-based approximation for LDA sampling that is doubly sparse, parallelizable, and asymptotically exact, improving efficiency over existing methods.

Findings

The proposed sampler is faster than previous Gibbs samplers.

It is asymptotically exact as data size grows.

Partially collapsed samplers can outperform fully collapsed ones.

Abstract

Latent Dirichlet Allocation (LDA) is a topic model widely used in natural language processing and machine learning. Most approaches to training the model rely on iterative algorithms, which makes it difficult to run LDA on big corpora that are best analyzed in parallel and distributed computational environments. Indeed, current approaches to parallel inference either don't converge to the correct posterior or require storage of large dense matrices in memory. We present a novel sampler that overcomes both problems, and we show that this sampler is faster, both empirically and theoretically, than previous Gibbs samplers for LDA. We do so by employing a novel P\'olya-urn-based approximation in the sparse partially collapsed sampler for LDA. We prove that the approximation error vanishes with data size, making our algorithm asymptotically exact, a property of importance for large-scale topic models. In addition, we show, via an explicit example, that - contrary to popular belief in the topic modeling literature - partially collapsed samplers can be more efficient than fully collapsed samplers. We conclude by comparing the performance of our algorithm with that of other approaches on well-known corpora.

Tables

Table 1. TABLE I: Notation for LDA. Sufficient statistics are conditional on algorithm’s current iteration. Bold symbols refer to matrices, bold italic symbols refer to vectors.

Symbol	Description		Symbol	Description
$V$	Vocabulary size		$𝚽:K\times V$	Word-topic probabilities
$D$	Total number of documents		${\varphi }_k:1\times V$	Word probabilities for topic $k$
$N$	Total number of words		$𝜷:1\times V$	Prior concentration vector for ${\varphi }_k$
$K$	Total number of topics		$𝐧:K\times V$	Topic-word sufficient statistic
$v\left(i\right)$	Word type for word $i$		$𝚯:D\times K$	Document-topic probabilities
$d\left(i\right)$	Document for word $i$		${𝜽}_d:1\times K$	Topic probabilities for document $d$
$w_{i,d}$	Word $i$ in document $d$		$𝜶:1\times K$	Prior concentration vector for ${𝜽}_d$
$z_{i,d}$	Topic indicator for word $i$ in $d$		$𝐦:D\times K$	Document-topic sufficient statistic

Table 2. TABLE II: Corpora used in experiments.

Corpus	$V$	$D$	$N$
Enron	27 508	39 860	6 377 365
New York Times	273 405	1 825 116	494 687 810
PubMed	89 987	8 199 999	768 434 972

Equations

$${\beta }_{\text{vbox}\text{scalebox}0.5\bullet }$$

$${\mathbf{\varphi }}_k\sim \mathrm{Dir}\left[{\beta }_{\text{vbox}\text{scalebox}0.5\bullet }+n_{k,\text{vbox}\text{scalebox}0.5\bullet },\frac{{\beta }_{\text{vbox}\text{scalebox}0.5\bullet }}{{\beta }_{\text{vbox}\text{scalebox}0.5\bullet }+n_{k,\text{vbox}\text{scalebox}0.5\bullet }}{\mathbf{F}}_0+\frac{n_{k,\text{vbox}\text{scalebox}0.5\bullet }}{{\beta }_{\text{vbox}\text{scalebox}0.5\bullet }+n_{k,\text{vbox}\text{scalebox}0.5\bullet }}{\hat{\mathbf{F}}}_k\right].$$

$$d_{\mathrm{LP}}(\mathbf{x},{\mathbf{x}}^{\ast })\to 0$$

$${\mathbf{\varphi }}_k=\left[\frac{{\gamma }_1}{{\sum }_{j=1}^V{\gamma }_j},..,\frac{{\gamma }_V}{{\sum }_{j=1}^V{\gamma }_j}\right]$$

$$\mathbf{x}$$

$$\mathbf{x}=\left[\frac{{\tilde{\gamma }}_1}{{\sum }_{j=1}^V{\tilde{\gamma }}_j},..,\frac{{\tilde{\gamma }}_V}{{\sum }_{j=1}^V{\tilde{\gamma }}_j}\right]$$

$$\mathbf{x}$$

$$\frac{{\tilde{\gamma }}_j}{{\sum }_{j=1}^V{\tilde{\gamma }}_j}$$

$$\pi$$

$${\mathbf{\varphi }}_k\sim \mathrm{PPU}\left[{\beta }_{\text{vbox}\text{scalebox}0.5\bullet }+n_{k,\text{vbox}\text{scalebox}0.5\bullet },\frac{{\beta }_{\text{vbox}\text{scalebox}0.5\bullet }}{{\beta }_{\text{vbox}\text{scalebox}0.5\bullet }+n_{k,\text{vbox}\text{scalebox}0.5\bullet }}{\mathbf{F}}_0+\frac{n_{k,\text{vbox}\text{scalebox}0.5\bullet }}{{\beta }_{\text{vbox}\text{scalebox}0.5\bullet }+n_{k,\text{vbox}\text{scalebox}0.5\bullet }}{\hat{\mathbf{F}}}_k\right],$$

$$P(z_{i,d}=k\mid {\mathbf{z}}_{-i,d},\mathbf{\Phi },w_{i,d})$$

$$\propto \underset{a}{\underbrace{{\varphi }_{k,v(i)}{\alpha }_k}}+\underset{b}{\underbrace{{\varphi }_{k,v(i)}{m}_{d,k}^{-i}}}.$$

$${\sigma }_a$$

$$O\left[\sum _{i=1}^N\min \left\{K_{d(i)}^{(\mathbf{m})},K_{v(i)}^{(\mathbf{\Phi })}\right\}\right].$$

$$P(z_{i,d}=k\mid {\mathbf{z}}_{-i,d},\mathbf{\Phi },w_{i,d})$$

$$\propto \underset{a}{\underbrace{{\varphi }_{k,v(i)}{\alpha }_k}}+\underset{b}{\underbrace{{\varphi }_{k,v(i)}{m}_{d,k}^{-i}}}.$$

$$q(z)=\sum _{k=1}^K{\varphi }_{k,v(i)}\left[{\alpha }_k+m_{d,k}^{-i}\right]={\sigma }_a+{\sigma }_b.$$

$${\sigma }_b=\sum _{k=1}^K{\varphi }_{k,v(i)}{m}_{d,k}^{-i}$$

$$O\left[\min \left\{K_{d(i)}^{(\mathbf{m})},K_{v(i)}^{(\mathbf{\Phi })}\right\}\right]$$

$$O\left[\sum _{i=1}^N\min \left\{K_{d(i)}^{(\mathbf{m})},K_{v(i)}^{(\mathbf{\Phi })}\right\}\right].$$

$$\mathbf{x}=\left[\frac{{\gamma }_1}{{\sum }_{j=1}^p{\gamma }_j},..,\frac{{\gamma }_p}{{\sum }_{j=1}^p{\gamma }_j}\right]=\frac{\mathbf{\gamma }}{{\sum }_{j=1}^p{\gamma }_j}$$

$$\mathbf{\gamma }=\sum _{i=1}^{\varpi }{\hat{\mathbf{\gamma }}}_{\mathbf{i}}$$

$$\sqrt{\varpi }\left[\frac{1}{\varpi }\mathbf{\gamma }-\mathbf{F}\right]\stackrel{d}{\to }N(0,\mathbf{\Sigma })$$

$$\sqrt{\varpi}\mathinner{\Big{[}g(\boldsymbol{\gamma}/\varpi)-g(\boldsymbol{F})\Big{]}}\overset{\operatorname{d}}{\rightarrow}\operatorname{N}(\boldsymbol{0},\mathbf{\tilde{\Sigma}})$$

$$|E(x_i)-F_i|$$

$$\left|\mathrm{Cov}(x_i,x_j)-\frac{-F_i{F}_j}{\varpi }\right|\to 0\text{ for }i\ne j.$$

$$|E(x_i)-F_i|$$

$$\left|\mathrm{Cov}(x_i,x_j)-\frac{-F_i{F}_j}{\varpi +1}\right|\to 0\text{ for }i\ne j$$

Full text

††publicationid: pubid: This work has been submitted to the IEEE for possible publication. Copyright may be transferred without notice.

1 Introduction

Latent Dirichlet Allocation (LDA) [1] is a topic model widely used in natural language processing for probabilistically identifying latent semantic themes in large collections of text documents, referred to as corpora. It does this by inferring the latent distribution of topics for each document based only on the word tokens, without any supervised labeling. In LDA, each document $d\in\{1,..,D\}$ is assigned a probability vector $\boldsymbol{\theta}_{d}$, where $\theta_{d,k}=\operatorname{\mathbb{P}}(\text{topic }k\text{ appears in document }d)$ – see Table I for a summary of the notation used in this paper. Then, within document $d$, each word, or token, at position $i$ is assigned a topic indicator $z_{i,d}\in\{1,..,K\}$. If word $w_{i,d}$ is assigned to topic $k$, it is assumed to be randomly drawn from a probability vector $\boldsymbol{\phi}_{k}$ over possible words within that topic.

LDA makes two key exchangeability assumptions.

1. (1)

All documents are exchangeable. 2. (2)

Within each document, word tokens are exchangeable (bag of words).

These assumptions yield a likelihood for the generative process described in Figure 1, which is then combined with conjugate Dirichlet priors on $\boldsymbol{\phi}_{k}$ and $\boldsymbol{\theta}_{d}$ to form a hierarchical Bayesian model. The posterior distribution describes, for every word, which topic the word belongs to.

To train the model, we need to draw samples from the posterior distribution. This can be done in a variety of ways, most of which belong to the following broad classes of methods.

• •

Markov Chain Monte Carlo (MCMC) techniques [2, 3, 4, 5, 6, 7, 8] such as Gibbs sampling can be used, and yield samples from the posterior upon convergence.

• •

Variational Bayesian (VB) methods [1, 9] are also used, which upon convergence yield samples from a variational approximation of the posterior.

• •

Expectation-Maximization (EM) techniques [10] are also applicable, which converge to the max a posteriori approximation of the posterior.

These are by no means exhaustive: other approaches, such as spectral [11] and geometric [12, 13] techniques can also be applied. Some methods [9] focus on rapidly learning the document-topic proportions $\mathbf{\Phi}$ by avoiding batch calculations with the topic indicators $\boldsymbol{z}$, whereas others – including this work – assume $\boldsymbol{z}$ itself is one of the quantities of user interest and attempt to calculate the full posterior as efficiently as possible.

We focus here on MCMC techniques – empirically, their scalability in topic models is comparable with VB and EM, and subject to convergence, they are exact, i.e., they yield the correct posterior.

LDA is increasingly being used with large data sets that are best analyzed in parallel and distributed computational environments. To be efficient in such settings, a sampler must do the following.

1. (1)

Expose sufficient parallelism in a way compatible with the hardware being used. 2. (2)

Take advantage of sparsity found in natural language [14, Ch. 2] to control asymptotic computational complexity and memory requirements.

This makes training models a challenge, because MCMC techniques (along with VB and EM) are iterative algorithms that tend to be inherently sequential. Furthermore, though existing techniques are capable of utilizing sparsity in the sufficient statistics $\mathbf{m}$ and $\mathbf{n}$ for $\mathbf{\Phi}$ and $\mathbf{\Theta}$, current approaches can only use sparsity in one sufficient statistic or the other to reduce computational cost, rather than in both simultaneously. This means that the cost depends on whichever one is less sparse.

To make this precise, we call an algorithm singly sparse if its iterative computational complexity is at least O\mathinner{\big{[}\sum_{i=1}^{N}K_{d(i)}^{(\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptscriptstyle\bullet}}}}})}\big{]}} or O\mathinner{\big{[}\sum_{i=1}^{N}K_{v(i)}^{(\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptscriptstyle\bullet}}}}})}\big{]}}, where $K_{d(i)}^{(\mathbf{m})}$ is the number of nonzero topics in the row of $\mathbf{m}$ corresponding to word $w_{i,d}$, $K_{v(i)}^{(\mathbf{n})}$ is the number of nonzero topics in the column of $\mathbf{n}$ corresponding to word $w_{i,d}$, and $K_{d(i)}^{(\mathbf{\Theta})}$, $K_{v(i)}^{(\mathbf{\Phi})}$ are defined analogously. In contrast, we call an algorithm doubly sparse if its expected iterative complexity is at most O\mathinner{\big{[}\sum_{i=1}^{N}\min\{K_{d(i)}^{(\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptscriptstyle\bullet}}}}})},K_{v(i)}^{(\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptscriptstyle\bullet}}}}})}\}\big{]}} – such an algorithm can use sparsity in both sufficient statistics simultaneously in such a way that its computational cost depends on whichever one is more sparse.

This is particularly the case in the standard MCMC approach – the collapsed Gibbs sampler of [2] modified according to [8] – where $\mathbf{\Phi}$ and $\mathbf{\Theta}$ are integrated out analytically, and each topic indicator $z_{i,d}$ is sampled one by one conditional on all other topic indicators $\boldsymbol{z}_{-i,d}$. In terms of computational cost, the resulting algorithm is singly sparse due to dependence, and exposes very little parallelism, which limits scalability considerably.

We build on the sparse partially collapsed sampler of [3]. This algorithm is exact, parallelizable by documents, and capable of utilizing sparsity in the sufficient statistic $\mathbf{m}$ for $\mathbf{\Theta}$. Further, it possesses favorable iterative complexity, and its convergence rate has been empirically shown to be reasonably fast. Unfortunately, it is singly sparse in the sense that it cannot use sparsity in $\mathbf{n}$, and unlike the standard approach it requires storage of $\mathbf{\Phi}$, a large dense matrix.

Our contribution is a method to bypass this issue by introducing a Pólya urn approximation. This allows us to construct a doubly sparse sampler with lower memory requirements, faster runtime, and more efficient iterative complexity.

2 Previous Work

The standard approach for LDA is the fully collapsed sampler of [2], modified according to [8] to utilize sparsity in the sufficient statistics $\mathbf{n}$ and $\mathbf{m}$. This sampler has iterative sampling complexity O\mathinner{\big{[}\sum_{i=1}^{N}\max\{K_{d(i)}^{(\mathbf{m})},K_{v(i)}^{(\mathbf{n})}\}\big{]}} – singly sparse due to the presence of $\max\{..\}$ – and mixing rate that appears reasonable empirically. A variety of methods for scalable LDA, including in parallel and distributed environments, have been proposed.

Partially Collapsed Samplers. In this class, $\mathbf{\Theta}$ is integrated out analytically but $\mathbf{\Phi}$ is retained. Unlike the fully collapsed sampler, the partially collapsed sampler can be parallelized over documents. This is possible due to the exchangeability assumptions and conjugacy structure present in the LDA model – a general feature of many Bayesian models [15]. In [3] it is shown that this class can utilize sparsity in a fashion similar to [4] to improve sampling speed and yield a singly sparse iterative complexity of O\mathinner{\big{[}\sum_{i=1}^{N}K_{d(i)}^{(\mathbf{m})}\big{]}}.

Partially collapsed samplers were initially studied in the parallel setting by [7] but were quickly dismissed for having slower convergence rates. This claim is incorrect: we exhibit an explicit counterexample in Appendix A in which a partially collapsed Gibbs sampler converges arbitrarily faster than a fully collapsed one. Convergence rate is thus model and algorithm-specific: for LDA, empirical results in [3] show that partially and fully collapsed samplers have comparable convergence behavior.

Asynchronous Distributed LDA. [7] propose modifying the collapsed Gibbs sampler presented by [2] by simply ignoring the sequential requirement. Unfortunately, the resulting algorithm is not exact. In [16], the authors bound the 1-step transition error – their analysis suggests that it is unlikely to accumulate over multiple iterations, and that it is likely to decrease with added parallelism. However, recent theory on Asynchronous Gibbs Sampling in [17] suggests the opposite – that the approximation error will increase with the added parallelism, because the algorithm will miscalculate the posterior’s covariance structure and jump into regions of low probability. AD-LDA has been studied empirically in [3], where it was shown that AD-LDA can stabilize at a lower log-posterior compared with other samplers, suggesting that asynchronicity can cause it to fail to find a good mode of the posterior distribution. We are unsure about how to reconcile these differences and recommend further empirical study of AD-LDA’s performance, which may depend on the posterior’s covariance structure and therefore on the specific corpus being studied.

Metropolis-within-Gibbs Samplers. [4], [5], and [6] propose methods in which Gibbs steps are replaced with Metropolis-Hastings steps. These modifications can reduce the iterative sampling complexity per document to O\mathinner{\big{[}\sum_{i=1}^{N}K_{d(i)}^{(\mathbf{m})}\big{]}}, the number of topics in each document, or even $O(N)$. However, the use of such steps can affect the MCMC algorithm’s mixing rate in ways that are difficult to understand theoretically. As a result, the combined effect on iterative complexity together with convergence rate is unknown and may depend on the specific corpus. These ideas are applicable to both fully collapsed and partially collapsed Gibbs samplers, and can be used in combination with the ideas that we present below.

3 The Algorithm

We begin by introducing the basic partially collapsed Gibbs sampler in [7] and [3]. Let $\mathbf{m}^{-i}$ be the sufficient statistic matrix $\mathbf{m}$ for $\mathbf{\Theta}$ with the portion of the statistic corresponding to word $i$ removed.

Algorithm 3.1** (Partially Collapsed Gibbs Sampler).**

Repeat until convergence:

• •

Sample $\boldsymbol{\phi}_{k}\sim\operatorname{Dir}(\boldsymbol{n}_{k}+\boldsymbol{\beta})$ in parallel for $k=1,..,K$.

• •

Sample $z_{i,d}\propto\phi_{k,v(i)}\,\alpha_{k}+\phi_{k,v(i)}\,m_{d,k}^{-i}$ in parallel for $d=1,..,D$.

3.1 Notation

Define

[TABLE]

The Dirichlet distribution may be parameterized either with a single probability vector $\boldsymbol{\beta}$ or with two inputs, a concentration parameter \beta_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptscriptstyle\bullet}}}}}} and a mean probability vector $\boldsymbol{F}_{0}$ whose product is $\boldsymbol{\beta}$ – this permits us in what follows to write $\operatorname{Dir}(\boldsymbol{\beta})$ and \operatorname{Dir}(\beta_{\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptscriptstyle\bullet}}}}}},\boldsymbol{F}_{0}) interchangeably. We introduce the two-input parametrization here because we consider asymptotics where the concentration parameter approaches $\infty$ while the mean probability vector is held fixed. In Algorithm 1 this permits us to write either $\boldsymbol{\phi}_{k}\sim\operatorname{Dir}(\boldsymbol{n}_{k}+\boldsymbol{\beta})$ or

[TABLE]

Finally, as we consider asymptotic convergence of random variables, we must introduce the necessary notions. We work exclusively with convergence in distribution, which we metrize via the Lévy-Prokhorov metric – though this is not strictly necessary, introducing this metric serves to simplify notation. As the $V$-dimensional probability simplex is a separable space, since it is a subset of $\mathbb{R}^{V}$, convergence in the Lévy-Prokhorov metric is equivalent here to convergence in distribution. Thus we may write

[TABLE]

to denote pointwise convergence of the cumulative distribution functions of random vectors $\boldsymbol{x}$ and $\boldsymbol{x}^{*}$.

3.2 Intuition

We now introduce a small modification to Algorithm 3.1. Instead of sampling from the Dirichlet distribution, we sample from a related distribution – the Pólya urn distribution, which arises in analysis of the Pólya urn model.

Throughout this work, we use ideas originating from the infinite-dimensional Dirichlet process literature to reason about finite-dimensional Dirichlet distributions. These can be viewed as a special case of a Dirichlet process on finite support, and we make use of this connection to develop ideas.

In the Pólya urn model, there is an urn containing a set of balls of different colors. We begin with a finite number of balls, with colors distributed according to some initial distribution. We draw and remove a ball, and then place two balls of the same color back into the urn. We repeat this process iteratively, letting $C_{t}$ be a random variable whose distribution is the distribution of colors inside the urn at time $t$. The sequence $\{C_{t}:t\in\mathbb{N}\}$ can be shown exchangeable, and is called a Pólya sequence. A classical result of [18] states that it admits a Dirichlet process as its de Finetti measure.

This suggests a way to bypass $\mathbf{\Phi}$ being dense: instead of drawing a random probability vector from the Dirichlet distribution, we can instead draw a set of IID samples from a Pólya urn and use them to form a probability vector. Since we are making a finite number of draws from a discrete distribution, some entries will be zero with non-negligible probability, and hence $\boldsymbol{\phi}_{k}$ will be sparse.

To proceed, we need to (a) decide how many draws to make and (b) find a parallel sampling scheme. Recall the representation of the Dirichlet distribution using IID Gamma random variables [19] – letting $\gamma_{j}\sim\operatorname{G}(\varpi_{j},1)$, where $\operatorname{G}$ denotes a Gamma distribution with the shape-scale parametrization, we have that if

[TABLE]

then $\boldsymbol{\phi}_{k}\sim\operatorname{Dir}(\boldsymbol{\varpi})$.

Notice that $\gamma_{j}\in\mathbb{R}^{+}$, $\operatorname{\mathbb{E}}(\gamma_{j})=\varpi_{j}$, and $\operatorname{Var}(\gamma_{j})=\varpi_{j}$. Consider using $\tilde{\gamma}_{j}\sim\operatorname{Pois}(\varpi_{j})$ as a replacement, for which it is also true that $\operatorname{\mathbb{E}}(\tilde{\gamma}_{j})=\varpi_{j}$ and $\operatorname{Var}(\tilde{\gamma}_{j})=\varpi_{j}$. It will be shown below that doing so will precisely yield probability vectors based on sampling from the Pólya urn.

This procedure has been used by [20] to analyze large data sets in the context of $A/B$ testing. It is scalable, as it is just a normalized version of the Poisson sampler described in [21] for use in parallel and distributed environments.

3.3 Convergence

We now prove that the resulting distribution converges to the original Dirichlet with increasing data size. We begin with the necessary definition.

Definition 3.2** (Poisson Pólya Urn).**

For positive $\varpi$ and probability vector $\boldsymbol{F}$, let $\boldsymbol{\varpi}=\varpi\boldsymbol{F}$. We say that

[TABLE]

if we have that

[TABLE]

for $\tilde{\gamma}_{j}\sim\operatorname{Pois}(\varpi_{j})$, provided that the sum is nonzero.

As with the Dirichlet, we consider two parametrizations – one based on a concentration vector, and another based on a concentration parameter together with a mean probability vector. Note that this definition is exactly the same as the Gamma representation of the Dirichlet distribution [19], except for Gammas being replaced with Poissons. We now introduce our main result.

Theorem 3.3** (Poisson Pólya Urn Asymptotic Convergence).**

Let

[TABLE]

Then for all $\boldsymbol{F}$ we have $d_{\operatorname{LP}}(\boldsymbol{x},\boldsymbol{x}^{*})\rightarrow 0$ as $\varpi\rightarrow\infty$.

Proof.

Define a Poisson process $\Pi$ over the unit interval with intensity $\varpi$, and let $\pi\sim\operatorname{Pois}^{+}(\varpi)$ be the number of arrivals over the entire interval – here, $\operatorname{Pois}^{+}$ refers to the Poisson distribution conditioned to be nonzero. Partition the unit interval into a set of subintervals $B_{j}$ with lengths equal to each component of $\boldsymbol{F}$. Let $\tilde{\gamma}_{j}$ be the number of arrivals in each subinterval. Then

[TABLE]

is just the proportion of arrivals of $\Pi$ in the subinterval $B_{j}$. Thus, since $\Pi$ is a Poisson process, $\boldsymbol{x}$ admits the hierarchical representation

[TABLE]

where $\boldsymbol{I}_{\upsilon^{(j)}}$ is a vector that is one at $\upsilon^{(j)}$ and zero everywhere else. A random variable admitting this hierarchical representation can be shown to converge to the desired Dirichlet – proving this is more technical, so we defer the remainder of the proof to Appendix B, noting that the result holds for all $\boldsymbol{F}$, regardless of whether $\boldsymbol{F}$ is fixed or variable. ∎

For LDA we have that under the two-input parametrization

[TABLE]

and, noting that n_{k,\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptscriptstyle\bullet}}}}}} is large, we can now define our algorithm.

Algorithm 3.4** (Pólya Urn Partially Collapsed Gibbs Sampler).**

Repeat until convergence:

• •

Sample $\boldsymbol{\phi}_{k}\sim\operatorname{PPU}(\boldsymbol{n}_{k}+\boldsymbol{\beta})$ in parallel for $k=1,..,K$.

• •

Sample $z_{i,d}\propto\phi_{k,v(i)}\,\alpha_{k}+\phi_{k,v(i)}\,m_{d,k}^{-i}$ in parallel for $d=1,..,D$.

By Theorem 3.3, Algorithm 3.4 converges to Algorithm 3.1 as n_{k,\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptscriptstyle\bullet}}}}}}\rightarrow\infty for all $k$, in the sense that their respective Gibbs steps converge in distribution. Using our method gives a number of improvements: $\mathbf{\Phi}$ is a sparse matrix and can thus be stored inexpensively, $\mathbf{\Phi}$ can be sampled more efficiently, and the resulting algorithm has better iterative computational complexity for sampling $\boldsymbol{z}$. We now describe these improvements in detail.

3.4 Efficient sampling of $\mathbf{\Phi}$

The standard partially collapsed sampler requires the rows of $\mathbf{\Phi}$ to be sampled from a Dirichlet distribution. This is done by drawing Gamma random variables, which can be accomplished with the technique in [22].

On the other hand, the Poisson Pólya urn requires us to draw Poisson random variables instead of Gammas – and, crucially, these Poisson random variables have rate parameters of the form $\beta+l$, $l\in\mathbb{N}_{0}$. For $l=0,..,L$, we can draw these random variables using a precomputed Walker alias table [23]. For $l>L$, we can use the fact that the Poisson distribution is asymptotically Gaussian as its rate parameter increases, and make rounded Gaussian draws. We choose $L=100$ and find that the resulting scheme reduces time spent sampling $\mathbf{\Phi}$ significantly – see Section 4.

3.5 Efficient sampling of $\boldsymbol{z}$

The sparsity of $\mathbf{\Phi}$ can be used to improve sampling speed for $\boldsymbol{z}$. Consider the approach in [3], where each probability is divided into $a$ and $b$, as follows:

[TABLE]

To sample $z_{i}$, we draw a random variable $u\sim\operatorname{U}(0,1)$ and compute $u_{\sigma}=u\mathinner{\left(\sigma_{a}+\sigma_{b}\right)}$, where

[TABLE]

If $u_{\sigma}<\sigma_{a}$, we draw $z_{i}$ from a precomputed alias table, otherwise we iterate over $\phi_{k,v(i)}\,m_{d,k}^{-i}$. Here we can use the fact that both $\mathbf{\Phi}$ and $\mathbf{m}$ are sparse, and iterate over whichever is smaller, to get a doubly sparse algorithm. This yields improved performance – see Section 4 – and also produces improved computational complexity, which we subsequently show. The sparsity in $\mathbf{\Phi}$ also makes it possible to construct sparse alias tables for $a$, further reducing the memory requirements for the algorithm.

3.6 Computational Complexity

We now derive the iterative computational complexity of Pólya Urn LDA. To do so, we need to consider how quickly the number $V$ of unique words grows with the number $N$ of total words. We assume the following.

Assumption 3.5** (Heaps’ Law).**

The number of unique words in a document follows Heaps’ Law $V=\xi N^{\psi}$ with constants $\xi>0$ and $0<\psi<1$.

For most languages, these constants will be in the ranges $5<\xi<50$ and $0.4<\psi<0.6$ – see [24, 25] for more details. Let $K_{d(i)}^{(\mathbf{m})}$ be the number of existing topics in document $d$ associated with word $i$. Let $K_{v(i)}^{(\mathbf{\Phi})}$ be the number of nonzero topics in the row of $\mathbf{\Phi}$ corresponding to word $i$. Note that \operatorname{\mathbb{E}}\mathinner{\big{[}K_{v(i)}^{(\mathbf{\Phi})}\big{]}}=K_{v(i)}^{(\mathbf{n})} at any given iteration. We initially assume that the number of topics is fixed, and subsequently consider the case where it grows with the size of the data.

Theorem 3.6** (Computational Complexity).**

Assuming a vocabulary size following Heaps’ Law, the iterative complexity of Pólya Urn LDA is

[TABLE]

Proof.

The iterative complexity of the sampler is equal to the complexity of sampling all $z_{i,d}$ plus the complexity of sampling $\mathbf{\Phi}$. Provided the number of topics is fixed, we need not consider the latter. Thus, it suffices to show that the complexity of all $z_{i,d}$ is given by (12). Consider a single $z_{i,d}$ with

[TABLE]

To sample $z_{i,d}$ we need to calculate the normalizing constant

[TABLE]

The first term is identical for all $z_{i,d}$, so it can be precomputed, and its iterative complexity is constant. The second term

[TABLE]

is the product of two sparse vectors. Suppose that both are stored in a data structure with $O(1)$ access, such as an array or hash map. By iterating over the smaller vector, the product can be computed with complexity

[TABLE]

for each $z_{i,d}$, and the result follows. ∎

Note the presence of $\min\{..\}$ in the complexity, in contrast with other approaches reviewed in Section 2. This demonstrates that Pólya Urn LDA is a doubly sparse algorithm.

The asymptotic complexity of LDA has also been studied in the setting in which the number $K$ of topics is assumed to grow with the number $N$ of total words [3], under the following assumption.

Assumption 3.7**.**

The number of topics $K$ is assumed to follow the mean of a Dirichlet process mixture, and thus the number of topics is $\gamma\ln\mathinner{\left[1+\frac{N}{\gamma}\right]}$, where $\gamma$ is the concentration parameter of the Dirichlet process.

See [26] for more details. Our result carries forward directly into this setting.

Corollary 3.8**.**

Assuming a vocabulary size following Heaps’ Law, and a number of topics following the mean of a Dirichlet Process mixture, the iterative complexity of Pólya Urn LDA is

[TABLE]

Proof.

As before, the iterative complexity is equal to the complexity of sampling all $z_{i,d}$ plus the complexity of sampling $\mathbf{\Phi}$. As the procedure for sampling $\mathbf{\Phi}$ under the Poisson Pòlya urn is identical to that under the Dirichlet, save for Gamma random variables being replaced with Poisson random variables, it follows immediately from [3] that its contribution to the total complexity is asymptotically negligible under the given assumptions. The result thus follows from Theorem 3.6. ∎

4 Performance Results

To study the empirical performance of Pólya Urn LDA, we implemented our algorithm in Java using the open-source Mallet111See http://mallet.cs.umass.edu and https://github.com/lejon/PartiallyCollapsedLDA [27] framework. We ran the algorithm on the Enron, New York Times, and PubMed corpora222Enron and PubMed can be found at: http://archive.ics.uci.edu/ml/datasets/Bag+of+Words

The New York Times Corpus can be found at: https://catalog.ldc.upenn.edu/ldc2008t19 – see Table II.

Computation for each experiment was performed on resources provided by Linköping University at the National Supercomputer Centre using an HP Cluster Platform 3000 with SL230s Gen8 compute nodes with two 8-core Intel Xeon E5–2660 processors each. We used symmetric priors with hyper-parameters $\alpha=0.1$ and $\beta=0.01$, and ran each experiment five times with different random number seeds. These are typical values and enable comparison with previous results in [3]. Unless specified otherwise, each experiment ran on all 16 cores, and with a rare word limit of 10.

In our first experiment, we compare Pólya Urn LDA with partially collapsed LDA, on which it is based and which was found to be the fastest parallel sampler in [3]. We ran 1,000 iterations with $K=10,100,1000$ with five different random number seeds.

Results can be seen in Figure 2, which gives the marginal unnormalized log-posterior values for $\boldsymbol{z}$ in real time. In all cases, Pólya Urn LDA was faster than standard partially collapsed LDA. Figure 3 shows the marginal log-posterior values per iteration. From this perspective, both algorithms look identical, indicating that Pólya Urn LDA mixes as quickly as partially collapsed LDA. We also ran the algorithm for $K=10$: performance of Pólya Urn LDA and partially collapsed LDA was nearly identical – this is expected, as $\mathbf{\Phi}$ is too small in this case for the Poisson Pólya urn to offer any benefit.

Relative runtime for $\mathbf{\Phi}$ and $\boldsymbol{z}$ separately can be seen in Figure 4 – as the number of topics increases, the relative speedup for $\boldsymbol{z}$ increases. We can also see that sampling $\mathbf{\Phi}$ using the alias tables and the Gaussian approximation in Section 3 is uniformly faster by a factor of four to eight.

In our second experiment, we compared Pólya Urn LDA with partially collapsed LDA [3], sparse AD-LDA [8], and Light AD-LDA [5] on the NYT corpora. Results can be seen in Figure 5, which shows that Pólya Urn LDA is faster for both $K=100$ and $K=1000$.

To ensure that the topics generated by the algorithm are reasonable, we evaluated held-out marginal log-likelihood values for $\boldsymbol{z}$ for Pólya Urn LDA and partially collapsed LDA. For this, we split the NYT corpus into a training and test set at random. Marginal log-likelihood evaluation on the test set was done using the left-to-right algorithm of [28]. We also computed the topic coherence scores of [29] on NYT – all of these diagnostics can be seen in Figure 6. Both algorithms produced similar test set log-likelihoods, and – with the exception of one outlier – similar levels of topic coherence. This indicates that Pólya Urn LDA generates topics of similar quality to partially collapsed LDA, and provides empirical evidence that the approximation error induced by the Poisson Pólya urn is small.

Next, we compared single-core performance of Pólya Urn LDA with partially collapsed LDA and sparse fully collapsed LDA. Pólya Urn LDA was found to outperform fully collapsed sparse LDA and partially collapsed LDA in this setting as well – results, averaged over all five random number seeds, are given in Figure 7. To study parallel scaling, we also looked at performance of Pólya Urn LDA while varying the number of available processor cores. We found that the method scales close to linearly with additional overhead for higher levels of parallelism, as expected for a massively parallel algorithm.

Finally, to better understand where the performance improvement comes from, we compared Pólya Urn LDA with partially collapsed LDA under different data sizes, vocabulary sizes, and rare word thresholds. For each fixed vocabulary size, words were selected according to their TFIDF values. Results can be seen in Figure 8. Here, we see that the relative difference between Pólya Urn LDA and partially collapsed LDA is larger for high vocabulary sizes and lower rare word thresholds. This is because increasing the dimension of the $\mathbf{\Phi}$ matrix, whose number of columns is linear in vocabulary size, doesn’t slow down $\boldsymbol{z}$ as much under Pólya Urn LDA as under partially collapsed LDA. This indicates that at least some of the speedup is obtained from rare words, for which most entries in the corresponding column of $\mathbf{\Phi}$ are zero at any given iteration.

Our algorithm’s performance depends on the particular data set under study. Compared to the standard partially collapsed sampler, the Pólya urn’s speedup comes from the additional sparsity in $\mathbf{\Phi}$. This sparsity will increase with $K$ for a fixed size corpus. In the PubMed example, our sampling time for $\boldsymbol{z}$ is reduced, compared to the partially collapsed sampler, by as much as 20% for $K=1000$ when compared with $K=100$. This matches what is expected from our complexity result in Section 3.6, as there will be a substantial number of words where $K_{v(i)}^{(\mathbf{\Phi})}<K_{d(i)}^{(\mathbf{m})}$, such as rare words or word types that only belong to a few topics.

5 Discussion

In this paper, we introduce the Poisson Pólya urn distribution as an asymptotically exact approximation to the Dirichlet distribution and use this to define Pólya Urn LDA for fast parallel sampling in LDA and similar models. Our algorithm is doubly sparse, in that it can take advantage of sparsity in the sufficient statistics $\mathbf{m}$ and $\mathbf{n}$ for $\mathbf{\Phi}$ and $\mathbf{\Theta}$ simultaneously. It inherits the parallelizability of partially collapsed LDA, which adds the following improvements.

1. (1)

The word-topic probability matrix $\mathbf{\Phi}$ is sparse and can be stored efficiently. 2. (2)

Sampling $\mathbf{\Phi}$ requires Poisson draws, which can be sampled efficiently – see Section 3. 3. (3)

The topic indicators $\boldsymbol{z}$ can be sampled efficiently by utilizing the added sparsity in $\mathbf{\Phi}$ – this improves our iterative complexity to $O\mathinner{\left[\sum_{i=1}^{N}\min\mathinner{\left\{K_{d(i)}^{(\mathbf{m})},K_{v(i)}^{(\mathbf{\Phi})}\right\}}\right]}$.

To our knowledge, this is the lowest computational complexity in any non-Metropolis-Hastings-based LDA sampler. It enables our algorithm to better take advantage of the power-law structure present in natural language. Combined with faster Poisson draws, we obtain significant speedup under standard corpora used for benchmarking – see Section 4.

The Poisson Pólya urn is an asymptotic approximation whose error we prove vanishes with increasing data set size. Our proof is based crucially on the Central Limit Theorem, which suggests that the convergence rate of our approximation is likely to be at least O(\sqrt{n_{k,\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{\displaystyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\textstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptstyle\bullet}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{\scriptscriptstyle\bullet}}}}}}}), though we have not rigorously demonstrated this.

In practice, we find the approximation error to be negligible, as from a convergence point of view, Pólya Urn LDA and partially collapsed LDA behave near-identically. This mirrors the behavior of scalable approximate MCMC methods recently proposed in other areas – see, e.g., [30]. Indeed, if we assume that the standard partially collapsed Gibbs sampler is geometrically ergodic, together with several other regularity conditions, then the theory in [31], as well as [32], provides bounds on the resulting Markov chain’s stationary distribution in terms of the correct posterior. On the other hand, this is also the method’s main drawback: the approximation’s convergence rate is unlikely to be dimension-free and thus may depend on $K$, so we advise caution when considering it for large $K$ with small data sets.

Our analysis – and, indeed, virtually every other complexity analysis in the LDA literature of which we are aware – is limited in that it describes the iterative complexity of the algorithm, not the effective complexity: it does not take into account the rate at which the Markov chain converges to the posterior, which will differ by algorithm. To calculate the effective complexity, we would need to prove that the chain is geometrically ergodic and compute the ergodicity coefficient, or find the spectral gap of the corresponding Markov operator. For LDA, this is currently an open problem. Similar issues occur in comparing MCMC approaches with Variational Bayes, Expectation-Maximization, and other methods, for which posterior approximation error needs to be considered yet is often challenging to assess.

Our approach can be combined with Metropolis-Hastings proposals such as those described in [4] and [5]. We did not consider these, as they significantly alter the convergence rate – this would make the speedup from the Poisson Pólya urn more difficult to understand. To compare these approaches with the standard partially collapsed sampler on a runtime basis see Figure 6 and [3].

Pólya Urn LDA is well-suited to a variety of parallel environments. Here we focus on the multicore setting and find that Pólya Urn LDA outperforms partially collapsed LDA, which in turn was found to outperform other approaches in [3]. Other parallel environments may yield different results: in particular, AD-LDA may perform better in the compute cluster setting, where asynchronous delays may serve to counteract the algorithm’s bias – see [17]. Pólya Urn LDA could also be implemented in such an environment – in fact, it can potentially be implemented using the Exact Asynchronous Gibbs scheme in [17], so that asynchronous delays introduce no approximation error into the posterior. These considerations are left to future work.

Developing massively parallel schemes for topic models that are both practical and theoretically understandable is an area of active current research. Our use of the Poisson Pólya urn approximation is generic: though we focus on LDA in this paper, the technique is equally applicable to any other topic model or language model making use of the Dirichlet distribution. This includes generative document cluster models [33], hidden Markov models used in parts-of-speech tagging [34], mixed membership stochastic blockmodels [35], and others. We hope that our work contributes to greater understanding in these areas.

Acknowledgments

We are thankful to Graham Neubig, Tamara Broderick, Eric P. Xing, Qirong Ho, Wei Dai, Willie Neiswanger, Murat Demirbas, Shawfeng Dong, Thanasis Kottas, Kunal Sarkhel, and Mattias Villani for their thoughts, and to several referees whose comments substantially improved the paper. Membership on this list does not imply agreement with the ideas expressed here, nor responsibility for any errors that may be present. Måns Magnusson was partially financially supported by the Swedish Foundation for Strategic Research (Smart Systems: RIT 15–0097).

Appendix A: Efficiency of Collapsed and Uncollapsed Gibbs Samplers

Here we demonstrate that a collapsed Gibbs sampler can be arbitrarily less efficient than an uncollapsed Gibbs sampler. Our example is a scaled $\operatorname{T}$ distribution on $\mathbb{R}^{2}$ with 5 degrees of freedom, with both components having mean 0 and scale matrix $\mathbf{\Sigma}$ such that $\Sigma_{11}=\Sigma_{22}=1$ and $\Sigma_{21}=\Sigma_{12}=\rho$. We define our variables $\Phi$ and $\boldsymbol{z}$ in a fashion analogous to those in LDA.

Algorithm 5.1** (Collapsed Gibbs Sampler).**

Repeat until convergence:

• •

Sample $z_{1}\mid z_{2}\sim\operatorname{T}[\rho z_{2},(0.8+0.2z_{2}^{2})(1-\rho^{2}),5]$

• •

Sample $z_{2}\mid z_{1}\sim\operatorname{T}[\rho z_{1},(0.8+0.2z_{1}^{2})(1-\rho^{2}),5]$

Algorithm 5.2** (Uncollapsed Gibbs Sampler).**

Repeat until convergence:

• •

Sample $\Phi\sim\operatorname{IG}[3,0.5\boldsymbol{z}^{T}\mathbf{\Sigma}^{-1}\boldsymbol{z}+2]$

• •

Sample $\boldsymbol{z}\sim\operatorname{N}_{2}[\boldsymbol{0},\Phi\mathbf{\Sigma}]$

By taking $\rho$ sufficiently large, Algorithm 5.1 can be made to converge at an arbitrarily slow rate. However, Algorithm 5.2 will converge at a rate independent of $\rho$. This can be seen intuitively by noting that for $\rho=1$, the Markov Chain in Algorithm 1 will be reducible and will thus fail to converge at all, whereas the chain in Algorithm 2 will suffer no such difficulty. This is illustrated in Figure 9, for $\rho=\{0.9,0.99,0.999\}$. This example demonstrates that a collapsed Gibbs sampler may be arbitrarily slower than an uncollapsed or partially collapsed Gibbs sampler.

This serves as an explicit counterexample to the argument made in [7] on the basis of Theorem 1 in [36]. The error in this argument is as follows: [36] assume reversibility and derive regularity conditions under which the correlation of two successive samples $(\boldsymbol{z}^{(i)},\boldsymbol{z}^{(i+1)})$ must strictly increase when some other parameter $\Phi$ is inserted between them. One crucial assumption is that $\boldsymbol{z}$ is sampled as a block whether $\Phi$ is present or not. This assumption does not apply to the example in this appendix, nor to LDA, where $\boldsymbol{z}$ is sampled as a block in one algorithm and component-wise in the other.

More generally, analogous issues have been shown to lead to provably slow mixing in certain large-scale Gibbs samplers – see [37] for details. It may be possible to apply similar arguments to LDA, or to extend the argument in [36] to the case where the dimensionality of $\boldsymbol{z}$ grows – this may show in certain samplers that mixing slows down with data size. This is a concern for any algorithm, as it would result in less favorable effective complexity – see Section 5.

Appendix B: Poisson-Pólya Urn Asymptotic Convergence

In this section we prove the remainder of Theorem 3.3. We first introduce some notation used throughout the rest of Appendix B. Let $\boldsymbol{x},\boldsymbol{x}^{*}\in\mathbb{R}^{p}$ be random vectors, let $\varpi\in\mathbb{R}^{+}$, and let $\boldsymbol{F}=(F_{1},..,F_{p})$ be a probability vector. We allow $\boldsymbol{x}$, $\boldsymbol{x}^{*}$, and $\boldsymbol{F}$ to be implicitly indexed by $\varpi$. The distributions of $\boldsymbol{x}$ and $\boldsymbol{x}^{*}$ will vary by context. As before, we work exclusively with convergence in distribution, which we again metrize via the Lévy-Prokhorov metric to simplify notation. Thus $d_{\operatorname{LP}}(\boldsymbol{x},\boldsymbol{x}^{*})\rightarrow 0$ can be taken to mean that $\boldsymbol{x}$ and $\boldsymbol{x}^{*}$ converge in distribution, i.e., their cumulative distribution functions converge pointwise. With this notation in mind, we now introduce several lemmas.

Lemma 5.3** (Gaussian Asymptotic Distribution of a Dirichlet Random Vector).**

Let $\boldsymbol{x}\sim\operatorname{Dir}(\varpi,\boldsymbol{F})$. Then there exists a sequence of multivariate Gaussian distributions $\boldsymbol{x}^{*}$ indexed by $\varpi$ such that $d_{\operatorname{LP}}(\boldsymbol{x},\boldsymbol{x}^{*})\rightarrow 0$ as $\varpi\rightarrow\infty$.

Proof.

Without loss of generality suppose $\varpi\in\mathbb{N}$. A Dirichlet random vector admits the representation

[TABLE]

with $\gamma_{i}\sim\operatorname{G}(\varpi F_{i},1)$. Since a sum of gamma random variables is gamma, we may write

[TABLE]

where for each $j$ we have $\hat{\gamma}_{ij}\overset{\textrm{iid}}{\sim}\operatorname{G}(F_{j},1)$. Note that $\operatorname{\mathbb{E}}(\boldsymbol{\hat{\gamma}}_{i})=\boldsymbol{F}$. It is now clear that taking $\varpi\in\mathbb{N}$ entails no loss in generality, because we may replace $F_{j}$ with $\varpi F_{j}/\lceil\varpi\rceil$ where $\lceil\cdot\rceil$ is the ceiling function, and take the summation from $1$ to $\lceil\varpi\rceil$. Since $\boldsymbol{\hat{\gamma}}_{i}$ have finite moments, we have by the Central Limit Theorem that

[TABLE]

for some covariance matrix $\mathbf{\Sigma}$ as $\varpi\rightarrow\infty$ where $\overset{\operatorname{d}}{\rightarrow}$ denotes convergence in distribution. Next, we use the delta method to show this carries over to $\boldsymbol{x}$. Define the function $g:\mathbb{R}^{p}\rightarrow\mathbb{R}^{p}$ by $\boldsymbol{x}=g(\boldsymbol{\gamma})$. Since $g$ is continuously differentiable for all strictly positive input, it is continuously differentiable at $\boldsymbol{F}$, as $\boldsymbol{F}$ is a probability vector. Then by Theorem 3.1 of [38], we have

[TABLE]

for a covariance matrix $\mathbf{\tilde{\Sigma}}$ as $\varpi\rightarrow\infty$. But notice that $g(\boldsymbol{\gamma}/\varpi)=g(\boldsymbol{\gamma})=\boldsymbol{x}$. Since the Lévy-Prokhorov metric metrizes convergence in distribution, the result follows. ∎

Lemma 5.4** (Asymptotic Distribution and Moments of the Poisson-Pólya Urn).**

Let $\boldsymbol{x}\sim\operatorname{PPU}(\varpi,\boldsymbol{F})$ be defined via the hierarchical representation (8). Let $\boldsymbol{x}^{*}\sim\operatorname{MN}(\varpi,\boldsymbol{F})$, where $\operatorname{MN}$ denotes the multinomial distribution. Then, as $\varpi\rightarrow\infty$, we have that $d_{\operatorname{LP}}(\varpi\boldsymbol{x},\boldsymbol{x}^{*})\rightarrow 0$ and for all $i,j$ that

[TABLE]

[TABLE]

Proof.

To show convergence in distribution, and convergence of the first two moments, it suffices to show that the probability generating functions of both random variables converge pointwise and are twice differentiable in a neighborhood around $\boldsymbol{1}$. This follows immediately from the argument in [39] by replacing their expression $\frac{1}{N}(x_{1}+..+x_{N})$ with $\sum_{i=1}^{p}F_{i}x_{i}$, and noting that the probability generating functions of $\operatorname{Pois}$ and $\operatorname{Pois}^{+}$ random variables converge pointwise as their means go to $\infty$. ∎

Lemma 5.5** (Gaussian Asymptotic Distribution of a Multinomial Random Vector).**

Let $\boldsymbol{x}\sim\operatorname{MN}(\varpi,\boldsymbol{F})$. Then there exists a sequence of multivariate Gaussian distributions $\boldsymbol{x}^{*}$ indexed by $\varpi$ such that $d_{\operatorname{LP}}(\boldsymbol{x},\boldsymbol{x}^{*})\rightarrow 0$ as $\varpi\rightarrow\infty$.

Proof.

Since a multinomial can be expressed as the sum of IID discrete distributions, the result follows from the Central Limit Theorem. ∎

We can now complete the convergence proof in Section 3.

Theorem 5.6** (Poisson Pólya Urn Asymptotic Convergence).**

Let $\boldsymbol{x}\sim\operatorname{PPU}(\varpi,\boldsymbol{F})$. Let $\boldsymbol{x}^{*}\sim\operatorname{Dir}(\varpi,\boldsymbol{F})$. Then $d_{\operatorname{LP}}(\boldsymbol{x},\boldsymbol{x}^{*})\rightarrow 0$ as $\varpi\rightarrow\infty$ for all $\boldsymbol{F}$.

Proof.

Using Lemma 5.4, Lemma 5.5, and Lemma 5.3, it follows that both random vectors converge in distribution to multivariate Gaussians. Since Gaussians are uniquely characterized by their first two moments, it suffices to show for all $i,j$ that

[TABLE]

[TABLE]

as $\varpi\rightarrow\infty$, since these are the mean and covariance of the Dirichlet. Since $\frac{\varpi}{\varpi+1}\rightarrow 1$, the differences between these moments and those in Lemma 5.4 approach zero, and the result follows. ∎