Scalable Inference for Nested Chinese Restaurant Process Topic Models

TL;DR

This paper introduces a scalable inference algorithm for nested Chinese restaurant process topic models, enabling efficient processing of massive text corpora while maintaining high model quality.

Contribution

A novel partially collapsed Gibbs sampling algorithm combined with an efficient distributed implementation for large-scale nCRP topic models.

Findings

111 times more efficient than previous methods

Can extract 1,722 topics from 131 million documents

Operates on a 28-billion-token corpus using 50 machines in 7 hours

Abstract

Nested Chinese Restaurant Process (nCRP) topic models are powerful nonparametric Bayesian methods to extract a topic hierarchy from a given text corpus, where the hierarchical structure is automatically determined by the data. Hierarchical Latent Dirichlet Allocation (hLDA) is a popular instance of nCRP topic models. However, hLDA has only been evaluated at small scale, because the existing collapsed Gibbs sampling and instantiated weight variational inference algorithms either are not scalable or sacrifice inference quality with mean-field assumptions. Moreover, an efficient distributed implementation of the data structures, such as dynamically growing count matrices and trees, is challenging. In this paper, we propose a novel partially collapsed Gibbs sampling (PCGS) algorithm, which combines the advantages of collapsed and instantiated weight algorithms to achieve good scalability as well as high model quality. An initialization strategy is presented to further improve the model quality. Finally, we propose an efficient distributed implementation of PCGS through vectorization, pre-processing, and a careful design of the concurrent data structures and communication strategy. Empirical studies show that our algorithm is 111 times more efficient than the previous open-source implementation for hLDA, with comparable or even better model quality. Our distributed implementation can extract 1,722 topics from a 131-million-document corpus with 28 billion tokens, which is 4-5 orders of magnitude larger than the previous largest corpus, with 50 machines in 7 hours.

Tables

Table 1. Table 1. Statistics of the datasets.

Dataset	$D$	# tokens	$V$
NYTimes (subset)	$3\times {10}^3$	$7.23\times {10}^5$	101635
NIPS	$1.5\times {10}^3$	$1.93\times {10}^6$	12375
NYTimes	$2.93\times {10}^5$	$9.7\times {10}^7$	101635
PubMed	$8.2\times {10}^6$	$7.38\times {10}^8$	141043
ClueWeb12 (small)	$1.5\times {10}^7$	$5.6\times {10}^9$	100000
ClueWeb12 (large)	$1.31\times {10}^8$	$2.8\times {10}^{10}$	100000

Equations

$$\{\begin{array}{ll}p(c_{dl}=t_{ik})=\frac{m_{t_{ik}}}{{\gamma }_l+m_i},& k=1,\dots ,K_i\\ p(c_{dl}=T+1)=\frac{{\gamma }_l}{{\gamma }_l+m_i},& \end{array}$$

$$p(\mathbf{w},\mathbf{z},\mathbf{c})=$$

$$f_C(d,t)=\frac{B({\mathbf{C}}_t^{\neg d}+{\mathbf{C}}_t^d+{\beta }_{l_t})}{B({\mathbf{C}}_t^{\neg d}+{\beta }_{l_t})}.$$

$$f_I(d,t)=\prod _{v=1}^V{\varphi }_{tv}^{C_{tv}^d}.$$

$$\propto \text{mbox}nCRP(\mathbf{c};\gamma ,{\mathbf{c}}_{\neg d})\prod _{l=1}^L\{\begin{array}{ll}{f}_I(d,c_l)& c_l\in \mathcal{I},\\ f_C(d,c_l)& c_l\in \mathcal{C}.\end{array}$$

$$\log f_I(d,t)=\log \prod _{v=1}^V({\varphi }_{tv}{)}^{C_{tv}^d}=\sum _{v\in {\mathcal{W}}_{d,l_t}}\log {\varphi }_{tv},$$

$$\log f_C(d,t)=$$

$$p({\mathbf{w}}^{\text{mbox}h}|{\mathbf{w}}^{\text{mbox}o},\varphi )\approx \prod _{d=1}^D\frac{1}{S}\sum _{s=1}^S\prod _{n=1}^{L_d}\sum _{l=1}^{L}p({\mathbf{w}}_{dn}^{\text{mbox}h},{\mathbf{z}}_{dn}^{\text{mbox}h}=l|{\mathbf{c}}_d^{(s)},{\theta }_d^{(s)},\varphi ),$$

$$\text{mbox}perplexity=\exp (-\text{mbox}loglikelihood/\text{mbox}numberoftokens),$$

Full text

In this paper, we propose a novel partially collapsed Gibbs sampling (PCGS) algorithm, which combines the advantages of collapsed and instantiated weight algorithms to achieve good scalability as well as high model quality. An initialization strategy is presented to further improve the model quality. Finally, we propose an efficient distributed implementation of PCGS through vectorization, pre-processing, and a careful design of the concurrent data structures and communication strategy.

Empirical studies show that our algorithm is 111 times more efficient than the previous open-source implementation for hLDA, with comparable or even better model quality. Our distributed implementation can extract 1,722 topics from a 131-million-document corpus with 28 billion tokens, which is 4-5 orders of magnitude larger than the previous largest corpus, with 50 machines in 7 hours.

distributed computing, Bayesian inference, nested Chinese Restaurant Process, topic models

††copyright: rightsretained††doi: 10.475/123_4††isbn: 123-4567-24-567/08/06††conference: ACM SIGKDD conference; 2017; Halifax, Nova Scotia Canada††journalyear: 2017††price: 15.00††ccs: Mathematics of computing Bayesian nonparametric models††ccs: Mathematics of computing Gibbs sampling††ccs: Computing methodologies Topic modeling††ccs: Computing methodologies Distributed algorithms

1. Introduction

Topic models are popular tools in the machine learning toolbox. They extract a set of latent topics from an input text corpus. Each topic is a unigram distribution over words, and the high-probability words often present strong semantic correlation. Topic models have been widely used in information retrieval (Wei and Croft, 2006), text analysis (Boyd-Graber et al., 2007; Zhu et al., 2012), information visualization (Wang et al., 2016), and many other application areas for feature extraction and dimensionality reduction.

However, the traditional topic models, such as Latent Dirichlet Allocation (LDA) (Blei et al., 2003), are flat. They do not learn any relationships between topics by assuming that the probabilities of observing all the topics are independent. On the other hand, topics are naturally organized in a hierarchy (Paisley et al., 2015). For example, when a topic on “unsupervised learning” is observed in a document, it is likely to also observe the more general topics containing the topic, such as “machine learning” and “computer science” in the same document. By capturing such relationships, hierarchical topic models can achieve deeper understanding and better generalization (Ahmed et al., 2013; Paisley et al., 2015) of the corpus than the flat models.

There are many different approaches to learning the topic hierarchy. For example, Google’s Rephil (Murphy, 2012) puts a hierarchical noisy-or network on the documents; the super-topic approach learns topics of topics (Li and McCallum, 2006; Pujara and Skomoroch, 2012); and the nested Chinese Restaurant Process (nCRP) (Blei et al., 2010; Ahmed et al., 2013; Paisley et al., 2015) approach utilizes the nCRP as a prior on topic hierarchies. Hierarchical topic models have been successfully applied to document modeling (Paisley et al., 2015), online advertising (Murphy, 2012) and microblog location prediction (Ahmed et al., 2013), outperforming flat models. Amongst these approaches, the nCRP method has a non-parametric prior on the topic hierarchy structure, which leads to a natural structure learning algorithm with Gibbs sampling, avoiding the slow-mixing Metropolis-Hastings proposals or neural rules (Murphy, 2012). The hierarchical Latent Dirichlet Allocation (hLDA) model is a popular instance of nCRP topic models (Blei et al., 2010). In hLDA, topics form a tree with an nCRP prior, while each document is assigned with a path from the root topic to a leaf topic, and words in the document are modeled with an admixture of topics on the path.

However, due to the lack of scalable algorithms and implementations, hLDA has only been evaluated at a small scale, e.g., with thousands of documents and tens of topics (Blei et al., 2010; Wang and Blei, 2009), which limits its wider adoption in real-life applications. The training of topic models can be accelerated via distributed computing, which has been successfully applied to LDA to handle hundreds of billions of tokens and millions of topics (Ahmed et al., 2012; Chen et al., 2016). Unfortunately, the previous algorithms for hLDA are unsuitable for distributed computing. Specifically, the collapsed Gibbs sampler (Blei et al., 2010) is difficult to parallelize because collapsing the topic distributions breaks the conditional independence between document-wise latent variables; on the other side, the instantiated weight variational inference algorithm (Wang and Blei, 2009) has inferior model quality because of label switching and local optimum, as we will analyze in Sec. 3.2. Moreover, the data structures used by hLDA, such as the dynamically growing count matrices and trees, are much more sophisticated than the data structures for LDA, and their efficient distributed implementations are challenging.

In this paper, we propose a novel partially collapsed Gibbs sampling (PCGS) algorithm, which combines the advantages of the collapsed Gibbs sampler (Blei et al., 2010) and the instantiated weight variational inference method (Wang and Blei, 2009) to achieve a good trade-off between the scalability and the quality of inference. We keep most topic distributions as instantiated to maximize the degree of parallelism; while we integrate out some rapid changing topic distributions to preserve the quality of inference. To further improve the model quality, we propose an initialization strategy. Finally, we present an efficient distributed implementation of PCGS through vectorization, pre-processing, and a careful design of the concurrent data structures and the communication strategy.

We design a set of experiments to extensively examine the model quality of our PCGS algorithm as well as its efficiency and scalability. The experimental results show that our single-thread PCGS is 111 times faster than the previous state-of-the-art implementation, hlda-c (Blei, 2009); and our distributed PCGS can extract 1,722 topics from a 131-million-document corpus with 28-billion tokens, which is 4-5 orders of magnitude larger than the previous largest corpus, with 50 machines in 7 hours. To the best of our knowledge, this is the first time to scale up hLDA for large-scale datasets.

2. Hierarchical LDA

We first review nCRP and hLDA for learning a topic hierarchy.

2.1. Nested Chinese Restaurant Process

Nested Chinese Restaurant process (nCRP) (Blei et al., 2010) represents a powerful nonparametric Bayesian method to learn a tree structure, whose width and depth are unbounded. Suppose there is a truncated tree with $L$ levels, where each node except the leaves has an infinite number of children. An unique ID is assigned to each node, where the root node has the ID 1. nCRP defines a probability distribution on a series of paths $(\bm{\mathbf{c}}_{1},\bm{\mathbf{c}}_{2},\dots)$ on the tree, where each path $\bm{\mathbf{c}}_{d}\in\mathbb{N}_{+}^{L}$ consists of $L$ node IDs from the root to a certain leaf. Given $\bm{\mathbf{c}}_{1},\dots,\bm{\mathbf{c}}_{d-1}$, we mark a node $i$ as visited if any of the paths passes through it, and the next path $\bm{\mathbf{c}}_{d}$ is generated as follows: (1) let $c_{d1}=1$; (2) for each level $l=2,\dots,L$, denote $i$ as a shortcut for $c_{d,l-1}$. Assume that there are already $T$ visited nodes, where the children of $i$ are denoted as $t_{i1},\dots,t_{iK_{i}}$. The next node of the path $c_{dl}$ can be generated as

[TABLE]

where $m_{i}:=\#\{(d,l)|c_{dl}=i\}$ is the number of times that node $i$ is visited, $\gamma_{l}$ is a hyper-parameter, and $\#\{\cdot\}$ denotes the cardinality of a set. If $c_{dl}=T+1$, the path goes through a child node of $c_{d,l-1}$, which is not visited before, we assign it the ID $T+1$. We refer this operation as the generation of a new child, although in fact it is just visiting a node that is never visited before. The above procedure is denoted as $\bm{\mathbf{c}}_{d}\sim\mbox{nCRP}(\bm{\mathbf{c}}_{d};\bm{\mathbf{\gamma}},\bm{\mathbf{c}}_{<d})$, where the subscript $<d$ stands for all the possible indices that are smaller than $d$, i.e., $\bm{\mathbf{c}}_{<d}=\{\bm{\mathbf{c}}_{1},\dots,\bm{\mathbf{c}}_{d-1}\}$.

Intuitively, nCRP puts a CRP (Teh, 2011) on each parent node, where the probability of visiting each node is proportional to its previous times of visit. Due to this fact, we can easily extend the stick-breaking formulation for CRP (Sethuraman, 1994) to nCRP (Wang and Blei, 2009). The generative procedure is:

• •

For each node $t\in\mathbb{N}_{+}$ on each level $l$, draw $\bm{\mathbf{\pi}}^{t}\sim\mbox{GEM}(\gamma_{l})$, where $\bm{\mathbf{\pi}}^{t}$ is a distribution over the children of $t$, and $\mbox{GEM}(\cdot)$ is the stick-breaking distribution (Sethuraman, 1994). A sample $\bm{\mathbf{\pi}}\sim\mbox{GEM}(\gamma)$ can be obtained as follows:

• –

For $i=1,\dots,\infty$, draw $V_{i}\sim\mbox{Beta}(1,\gamma)$, and let $\pi_{i}=V_{i}\prod_{j=1}^{i-1}(1-V_{j})$.

• •

For each path $\bm{\mathbf{c}}_{d}$, let $c_{d1}=1$, and for $l=2,\dots,L$, select a child $k_{dl}\sim\mbox{Mult}(\bm{\mathbf{\pi}}^{c_{d,l-1}})$, and let the corresponding id $c_{dl}=t_{c_{d,l-1},k_{dl}}$.

The original nCRP is recovered by integrating out $\{\bm{\mathbf{\pi}}^{t}\}$. In the stick-breaking formulation, the probability of visiting each child is explicitly instantiated, and the paths $\{\bm{\mathbf{c}}_{d}\}$ are conditionally independent given the probabilities $\{\bm{\mathbf{\pi}}^{t}\}$.

2.2. Hierarchical Latent Dirichlet Allocation

Given a corpus of $D$ bag-of-words documents $\bm{\mathbf{W}}=\{\bm{\mathbf{w}}_{d}\}_{d=1}^{D}$, where each document $\bm{\mathbf{w}}_{d}=\{w_{dn}\}_{n=1}^{N_{d}}$ has $N_{d}$ tokens, and each token is represented by its word id $w_{dn}\in\{1,\dots,V\}$ in the vocabulary of $V$ unique words. hLDA is an nCRP-based topic model to learn a topic hierarchy (Blei et al., 2010). In hLDA, topics form a $L$-level tree, i.e., each tree node $t$ is a topic, and is associated with a distribution over words $\bm{\mathbf{\phi}}_{t}\in\Delta^{V-1}$, where $\Delta^{V-1}$ is the $(V-1)$-simplex. Since nodes and topics have one-to-one correspondence, we do not distinguish them in the sequel.

In hLDA, each document is assigned with a path $\bm{\mathbf{c}}_{d}$, and its words are modeled with a mixture of the topics in $\bm{\mathbf{c}}_{d}$, with the document-specific mixing proportion $\bm{\mathbf{\theta}}_{d}$. The generative process for the corpus is:

• •

For each node $t$, draw $\bm{\mathbf{\phi}}_{t}\sim\mbox{Dir}(\beta_{l_{t}}\mathbf{1})$, where $l_{t}$ is the level of node $t$ and $\mathbf{1}=(1,\dots,1)$ is an all-one vector;

• •

For each document $d$:

• –

Draw $\bm{\mathbf{c}}_{d}\sim\mbox{nCRP}(\bm{\mathbf{c}}_{d};\bm{\mathbf{\gamma}},\bm{\mathbf{c}}_{<d})$;

• –

Draw $\bm{\mathbf{\theta}}_{d}\sim\mbox{Dir}(\bm{\mathbf{\alpha}})$;

• –

For each position $n$, draw $z_{dn}\sim\mbox{Mult}(\bm{\mathbf{\theta}}_{d})$, and draw $w_{dn}\sim\mbox{Mult}(\bm{\mathbf{\phi}}_{c_{d,z_{dn}}})$,

where $\mbox{Dir}(\cdot)$ is the Dirichlet distribution, and $\bm{\mathbf{\alpha}}$ and $\bm{\mathbf{\beta}}$ are Dirichlet hyper-parameters. There are two special cases of hLDA. When the tree degenerates to a chain, hLDA recovers the vanilla LDA with $L$ topics, and when the tree has two levels and the probability of assigning to the first level $\theta_{d1}$ is close to zero, hLDA recovers the Dirichlet Process Mixture Model (DPMM) (Neal, 2000).

3. Inference for HLDA

There are two classes of algorithms for the posterior inference in hLDA—the collapsed weight algorithm that integrates out the mixing weights $\bm{\mathbf{\pi}}$ and parameters $\bm{\mathbf{\phi}}$, and the instantiated weight algorithm that explicitly infers these parameters. These algorithms present a trade-off between scalability and the quality of inference.

In this section, we first introduce the collapsed Gibbs sampler (Blei et al., 2010), which is not scalable. To address this problem, we present an instantiated weight block Gibbs sampler, which is based on the same stick-breaking formulation as the variational inference algorithm (Wang and Blei, 2009) and has an excellent scalability. However, it suffers from local optima. To tackle this issue, we propose a partially collapsed Gibbs sampler that has good scalability as well as high-quality inference. We also present an initialization strategy to find better local optima.

3.1. Collapsed Gibbs Sampling (CGS)

CGS is a collapsed weight algorithm that is based on the nCRP formulation. The generative process of hLDA (Sec. 2.2) defines a joint distribution $p(\bm{\mathbf{w}},\bm{\mathbf{z}},\bm{\mathbf{\theta}},\bm{\mathbf{c}},\bm{\mathbf{\phi}})$= $\prod_{t=1}^{\infty}$ $p(\bm{\mathbf{\phi}}_{t})$ $\prod_{d=1}^{D}$ $p(\bm{\mathbf{c}}_{d}|\bm{\mathbf{c}}_{<d})$ $p(\bm{\mathbf{\theta}}_{d})$ $\prod_{d=1}^{D}$ $\prod_{n=1}^{N_{d}}$ $p(z_{dn}|\bm{\mathbf{\theta}}_{d})$ $p(w_{dn}|\bm{\mathbf{\phi}}_{c_{d,z_{dn}}}).$ Based on the conjugacy between Dirichlet and multinomial, $\bm{\mathbf{\theta}}$ and $\bm{\mathbf{\phi}}$ are integrated out to get the collapsed distribution:

[TABLE]

where $\bm{\mathbf{C}}_{d}=(C_{d1},\dots,C_{dL})$ are the document-level counts, i.e., $C_{dl}=\#\{n|z_{dn}=l\}$, $\bm{\mathbf{C}}_{t}$ are the topic-word counts, i.e., $C_{tv}=\#\{(d,n)|w_{dn}=v\wedge c_{d,z_{dn}}=t\}$, and $B(\cdot)$ is the multivariate beta function, $B(\bm{\mathbf{\alpha}})=\prod_{k}\Gamma(\alpha_{k})/\Gamma(\sum_{k}\alpha_{k})$.

CGS alternatively samples $\bm{\mathbf{z}}$ and $\bm{\mathbf{c}}$ from their conditional distributions: $p(z_{dn}=l|w_{dn}=v,\bm{\mathbf{w}},\bm{\mathbf{z}}_{\neg dn},\bm{\mathbf{c}})\propto(C_{dl}^{\neg dn}+\alpha_{l})\frac{C_{c_{dl},v}^{\neg dn}+\beta_{l}}{C_{c_{dl}}^{\neg dn}+V\beta_{l}}$, and $p(\bm{\mathbf{c}}_{d}=\bm{\mathbf{c}}|\bm{\mathbf{w}},\bm{\mathbf{z}},\bm{\mathbf{c}}_{\neg d})\propto\mbox{nCRP}(\bm{\mathbf{c}};\bm{\mathbf{\gamma}},\bm{\mathbf{c}}_{\neg d})\prod_{l=1}^{L}f_{C}(d,c_{l})$, where $\neg d$ represents excluding document $d$, and $\neg dn$ means excluding the token $(d,n)$, e.g., $C_{tv}^{\neg d}=\#\{(d^{\prime},n)|w_{d^{\prime}n}=v\wedge c_{d^{\prime},z_{d^{\prime}n}}=t\wedge d^{\prime}\neq d\}$. Furthermore, the superscript d denotes only considering document $d$, e.g., $C_{tv}^{d}=\#\{n|w_{dn}=v\wedge c_{d,z_{dn}}=t\}$. Finally, $C_{t}=\sum_{v=1}^{V}C_{tv}$ and

[TABLE]

A straightforward approach for parallelizing CGS is to let each thread work on a disjoint set of documents and synchronize the counts $\{\bm{\mathbf{C}}_{t}\}$ between threads and machines. One possible solution for the synchronization is the parameter server (Ahmed et al., 2012), which maintains a local copy of $\{\bm{\mathbf{C}}_{t}\}$ on each worker machine as well as on a parameter server. Each machine periodically synchronizes its local copy with the parameter server by pushing its change to the server and fetching the latest parameters. Multiple worker threads read and write the local copy concurrently, as illustrated in Fig. 1(a).

However, this approach has several disadvantages, limiting its scalability and efficiency. Firstly, due to the limited network bandwidth, the period of synchronization can be relatively long, e.g., minutes. While this is acceptable for LDA, the stale state can potentially harm the quality of inference for hLDA, which is much more sensitive to local optima, as we will analyze in Sec. 3.2. Secondly, the local copy of $\{\bm{\mathbf{C}}_{t}\}$ needs to support concurrent reads, updates, and resizes from the worker threads, which is notoriously difficult and much more expensive than a serial version (Dechev et al., 2006). Finally, even in a serial setting, the computational cost of CGS is high because $f_{C}(d,t)$ involves the computation of the multivariate beta function, which is computed with gamma functions (see the appendix), that are much more expensive to compute than simple arithmetic.

3.2. Block Gibbs Sampling

To address the scalability and efficiency issues of CGS, we begin with a block Gibbs sampler (BGS), which is an instantiated weight algorithm that is based on the same model formulation of the variational inference algorithm (Wang and Blei, 2009), but the per-iteration time complexity is made lower by replacing expectation with sampling.

The BGS is based on the stick-breaking formulation of nCRP (defined in Sec. 2.1), which defines a joint distribution $p(\bm{\mathbf{w}},\bm{\mathbf{z}},\bm{\mathbf{\theta}},\bm{\mathbf{c}},\bm{\mathbf{\phi}},\bm{\mathbf{\pi}})$. Integrating out $\bm{\mathbf{\theta}}$, BGS samples in the posterior distribution of $(\bm{\mathbf{z}},\bm{\mathbf{c}},\bm{\mathbf{\phi}},\bm{\mathbf{\pi}})$ by alternatively sampling $\bm{\mathbf{z}}$, $\bm{\mathbf{c}}$, $\bm{\mathbf{\phi}}$ and $\bm{\mathbf{\pi}}$ given the others. The resultant updates are as follows:

Sample $\bm{\mathbf{z}}$: for each token, draw a level assignment from the distribution $p(z_{dn}=l|w_{dn}=v,\bm{\mathbf{w}},\bm{\mathbf{c}},\bm{\mathbf{\phi}},\bm{\mathbf{\pi}})\propto(C_{dl}^{\neg dn}+\alpha_{l})\phi_{c_{dl},v}.$

Sample $\bm{\mathbf{c}}$: for each document, sample a path from the conditional distribution $p(\bm{\mathbf{c}}_{d}=\bm{\mathbf{c}}|\bm{\mathbf{w}},\bm{\mathbf{z}},\bm{\mathbf{c}}_{\neg d},\bm{\mathbf{\phi}},\bm{\mathbf{\pi}})\propto\prod_{l=2}^{L}\pi_{c_{l-1}\rightarrow c_{l}}\prod_{l=1}^{L}$ $f_{I}(d,c_{l}),$ where $\pi_{i\rightarrow j}$ is the probability of going from node $i$ to its child $j$, i.e., $\pi^{i}_{k}=\pi_{i\rightarrow t_{ik}}$, and

[TABLE]

Sample $\bm{\mathbf{\pi}}$: Draw the stick-breaking weights $\pi^{t}_{k}=V^{t}_{k}\sum_{j=1}^{k-1}(1-V^{t}_{j})$, where $V^{t}_{k}\sim\mbox{Beta}(1+m^{t}_{k},\gamma+m^{t}_{>k})$, $m^{t}_{k}$ is the number of times that a path go through $t$ and its $k$-th child, and $m^{t}_{>k}=\sum_{j=k+1}^{\infty}m^{t}_{j}$.

Sample $\bm{\mathbf{\phi}}$: Draw the topic distribution $p(\bm{\mathbf{\phi}}_{t})=\mbox{Dir}(\beta_{l_{t}}+\bm{\mathbf{C}}_{t})$.

A subtlety here is on sampling $\bm{\mathbf{\pi}}$ and $\bm{\mathbf{\phi}}$, where $\bm{\mathbf{\pi}}^{t}$ is infinite-dimensional, and there are infinite $\bm{\mathbf{\phi}}_{t}$’s. We approximate the sampling by truncating $\bm{\mathbf{\pi}}^{t}$ to be finite-dimensional, i.e., there are finite children for each node, so that the whole tree has a finite number of nodes. This approximation can be avoided with slice sampling (Ge et al., 2015), but truncation is not the main reason affecting the model quality, as there are some more severe issues as we shall see soon.

Due to conditional independence, the document-specific random variables $\bm{\mathbf{z}}$ and $\bm{\mathbf{c}}$ can be sampled in parallel for each document. This fits in a bulk-synchronous parallel (BSP) pattern. In each iteration, $\bm{\mathbf{\pi}}$ and $\bm{\mathbf{\phi}}$ are sampled and broadcast-ed to each worker, and then the workers sample $\bm{\mathbf{z}}$ and $\bm{\mathbf{c}}$ without any communication, as illustrated in Fig. 1(b). BSP has been successfully adopted for LDA (Zhai et al., 2012; Zaheer et al., 2016; Chen et al., 2016) and achieved superior throughput than the asynchronous version (e.g., parameter server) for being lock free (Zaheer et al., 2016). Moreover, the BGS update for $\bm{\mathbf{c}}$ is also cheaper than that of CGS because the computation of $f_{I}(d,t)$ only involves the logarithm of $\bm{\mathbf{\phi}}$, which remains invariant during the sampling of $\bm{\mathbf{z}}$ and $\bm{\mathbf{c}}$, and therefore, can be pre-processed.

Unfortunately, while BSP works well for LDA, its quality of inference is unsatisfactory for hLDA in practice. We provide a number of explanations for this phenomenon:

$\bm{\mathbf{C}}_{t}$** is not slow-changing.** BSP works well for LDA because the topic-word count $\bm{\mathbf{C}}_{t}$ is slow-changing. Therefore, a stale $\bm{\mathbf{C}}_{t}$ is close to the fresh version, and the result is not much affected. However, this assumption is not true for hLDA because a topic can have very few documents assigned to it. For instance, if the sampler assigns a document to a topic that do not have any assigned documents, the topic-word count $\bm{\mathbf{C}}_{t}$ for that topic will suddenly change from a zero vector to non-zero, which differs significantly with its stale version.

Label switching. If two different workers generate two new topics, it is not clear whether they are the same one. For example, in a certain iteration, documents $d_{1}$ and $d_{2}$ should be assigned to two different new topics $t_{1}$ and $t_{2}$. But in BGS, two different workers may decide to assign $d_{1}$ and $d_{2}$ to a same topic $t$, because both workers do not know the changes made by the other worker, and just regard $t$ as a topic that no document is assigned to it. As the result, instead of learning two different topics $t_{1}$ and $t_{2}$, BGS learns one topic $t$ that is a combination of $t_{1}$ and $t_{2}$. For flat models such as DPMM or hierarchial Dirichlet process (HDP) (Teh et al., 2004), label switching is sometimes (approximately) resolved by running an algorithm that matches the new topics from different workers (Campbell et al., 2015). However, it is not clear how to match topics on trees.

Local optima. In flat models, even when label switching happens, e.g., two topics $t_{1}$ and $t_{2}$ are mixed as one topic, the algorithm may gradually separate them by generating a new topic and assigning the documents that should belong to $t_{2}$ to the new topic (Yang et al., 2016). However, these moves are more difficult for hLDA because it is more sensitive to local optima. For instance, if two topics $t_{1}$ and $t_{2}$ are incorrectly mixed as one topic $t$, and $t$ has a sub-tree. To correctly separate $t_{1}$ and $t_{2}$, the sampler needs to create a new brother of $t$, and move some decedents of $t$ to its brother. These operations can hardly be achieved with local moves. Wang and Blei (Wang and Blei, 2009) attempted to make this kind of moves by split-and-merge operations, whose time complexity is typically quadratic with the number of topics, and does not scale to a large number of topics.

3.3. Partially Collapsed Gibbs Sampling

It can be seen from the aforementioned discussion that there is a trade-off between scalability and the quality of inference. CGS learns good models but is not scalable, while BGS is very scalable but sacrifices the quality of inference. To combine their advantages, we propose a partially collapsed Gibbs sampler (PCGS).

Intuitively, if a topic $t$ has lots of assigned documents, its count $\bm{\mathbf{C}}_{t}$ changes slowly. Based on this insight, we categorize the topics as slow changing topics (SCTs) $\mathcal{I}$ and rapid changing topics (RCTs) $\mathcal{C}$, such that the number of SCTs dominates. Then, we perform CGS for the RCTs $\mathcal{C}$ and BGS for the SCTs $\mathcal{I}$. The quality of inference is not greatly affected because we perform CGS for the RCTs, and the scalability and efficiency is good because for most topics we perform the scalable and efficient BGS. We define a topic $t$ to be rapid changing if it is assigned to less than $M$ (a user-defined constant) documents, and slow changing otherwise.

Formally, let $\bm{\mathbf{\phi}}_{\mathcal{C}}=\{\bm{\mathbf{\phi}}_{t}|t\in\mathcal{C}\}$, $\bm{\mathbf{\phi}}_{\mathcal{I}}=\{\bm{\mathbf{\phi}}_{t}|t\in\mathcal{I}\}$, $\bm{\mathbf{C}}_{\mathcal{C}}=\{\bm{\mathbf{C}}_{t}|t\in\mathcal{C}\}$ and $\bm{\mathbf{C}}_{\mathcal{I}}=\{\bm{\mathbf{C}}_{t}|t\in\mathcal{I}\}$, we derive the following Gibbs sampling updates, where the details can be found in the appendix:

Sample $\bm{\mathbf{z}}$: draw the level assignment for each token from $p(z_{dn}=l|w_{dn}=v,\bm{\mathbf{w}},\bm{\mathbf{z}}_{\neg dn},\bm{\mathbf{c}},\bm{\mathbf{\phi}}_{\mathcal{I}})$ $\propto(C_{dl}^{\neg dn}+\alpha_{l})\begin{cases}\phi_{c_{dl},v}&c_{dl}\in\mathcal{I},\\ \frac{C_{c_{dl},v}^{\neg dn}+\beta_{l_{t}}}{C_{c_{dl}}^{\neg dn}+V\beta_{l_{t}}}&c_{dl}\in\mathcal{C}.\end{cases}$

Sample $\bm{\mathbf{c}}$: sample the path from $p(\bm{\mathbf{c}}_{d}=\bm{\mathbf{c}}|\bm{\mathbf{w}},\bm{\mathbf{z}},\bm{\mathbf{c}}_{\neg d},\bm{\mathbf{\phi}}_{\mathcal{I}})$

[TABLE]

Sample $\bm{\mathbf{\phi}}^{\mathcal{I}}$: For $t\in\mathcal{I}$, draw $\bm{\mathbf{\phi}}_{t}\sim\mbox{Dir}(\beta_{l_{t}}+\bm{\mathbf{C}}_{t})$.

These PCGS updates just combine the update rules of CGS and BGS, which utilizes CGS rule for $t\in\mathcal{C}$ and BGS rule for $t\in\mathcal{I}$.

Since the document visit counts $\bm{\mathbf{m}}$ (defined in Sec. 2.1) only requires $O(T)$ space, where $T$ is the number of topics, it is cheap to synchronize. We keep the entire tree weight $\bm{\mathbf{\pi}}$ collapsed out, and periodically synchronize the counts across machines. PCGS creates new topics in the same way as CGS. Therefore, PCGS does not require truncation and has the correct stationary distribution.

For distributed computing, PCGS performs asynchronous updates for the rapid-changing $(\bm{\mathbf{C}}_{\mathcal{C}},\bm{\mathbf{m}})$ and perform BSP updates for the slow-changing counts $\bm{\mathbf{C}}_{\mathcal{I}}$, as illustrated in Fig. 1(c). Since there are few rapid-changing topics, the amount of asynchronous updates of PCGS is much smaller than that of CGS, which needs to update all the counts $(\{\bm{\mathbf{C}}_{t}\},\bm{\mathbf{m}})$ asynchronously. Thanks to the small size of asynchronous updates, network bandwidth is not a bottleneck for PCGS, and PCGS can update the counts more frequently than CGS. In the sequel, the PCGS counts are more fresh than CGS in distributed setting. Because the number of slow-changing topics dominates, PCGS enjoys similar scalability and efficiency as BGS.

3.4. Initialization Strategy

hLDA is sensitive to local optima, so a proper initialization is crucial for obtaining good results. We adopt the progressive online initialization strategy (Blei et al., 2010; Wang and Blei, 2009), which begins with an empty collection of documents, and gradually adds documents by inferring the posterior of document-specific variables $(\bm{\mathbf{c}}_{d},\bm{\mathbf{z}}_{d})$ given all the previously observed documents. The documents are organized into mini-batches, and $\bm{\mathbf{\phi}}_{\mathcal{I}}$ is sampled per mini-batch.

To further improve the model quality, we noticed that all the aforementioned algorithms update $\bm{\mathbf{c}}$ and $\bm{\mathbf{z}}$ while keeping the other fixed, which can severely trap the sampler in local optima. For example, after a document $d$ is assigned to a certain path $\bm{\mathbf{c}}_{d}$, its words are assigned to the levels $\bm{\mathbf{z}}_{d}$ of the current path. In the next iteration, even if there is another path $\bm{\mathbf{c}}^{\prime}_{d}$ such that $p(\bm{\mathbf{w}}_{d}|\bm{\mathbf{c}}^{\prime}_{d})$ is larger than $p(\bm{\mathbf{w}}_{d}|\bm{\mathbf{c}}_{d})$, $p(\bm{\mathbf{w}}_{d}|\bm{\mathbf{c}}^{\prime}_{d},\bm{\mathbf{z}}_{d})$ is not likely to be larger than $p(\bm{\mathbf{w}}_{d}|\bm{\mathbf{c}}_{d},\bm{\mathbf{z}}_{d})$ because $\bm{\mathbf{z}}_{d}$ is already optimized for $\bm{\mathbf{c}}_{d}$. In this case, the path assignments $\bm{\mathbf{c}}_{d}$ will be quickly trapped in a local optimum even if there are better path assignments. We also noticed that similar as in multinomial mixture models (Rigouste et al., 2007), the sampling of $\bm{\mathbf{c}}_{d}$ is almost deterministic, because $\log p(\bm{\mathbf{w}}_{d}|\bm{\mathbf{c}}_{d},\bm{\mathbf{z}}_{d})$ is a sum of log-likelihoods over all the words, and can differ by hundreds for different $\bm{\mathbf{c}}_{d}$’s. Therefore, it is difficult for a sampler to jump out of the local trap simply by its randomness.

We propose a remedy for this problem by sampling $\bm{\mathbf{c}}$ from $p(\bm{\mathbf{c}}|\bm{\mathbf{w}})$ directly instead of from $p(\bm{\mathbf{c}}|\bm{\mathbf{z}},\bm{\mathbf{w}})$ (Eq. 4) for the first $I$ iterations. In other words, we integrate $\bm{\mathbf{z}}$ out. In the first $I$ iterations, the sampler focuses on finding the optimal assignment $\bm{\mathbf{c}}$ for each document. Afterwards, the algorithm samples $p(\bm{\mathbf{c}}|\bm{\mathbf{z}},\bm{\mathbf{w}})$ to refine the model. Unfortunately, $p(\bm{\mathbf{c}}|\bm{\mathbf{w}})=\sum_{\bm{\mathbf{z}}}p(\bm{\mathbf{c}}|\bm{\mathbf{z}},\bm{\mathbf{w}})p(\bm{\mathbf{z}})$ has no closed-form representation. We approximate it with Monte-Carlo integration $p(\bm{\mathbf{c}}|\bm{\mathbf{w}})\approx\frac{1}{S}\sum_{\bm{\mathbf{z}}_{s}\sim p(\bm{\mathbf{z}})}p(\bm{\mathbf{c}}|\bm{\mathbf{z}}_{s},\bm{\mathbf{w}}),$ where $S$ is the number of samples, and $p(\bm{\mathbf{z}})=\int_{\bm{\mathbf{\theta}}}p(\bm{\mathbf{z}}|\bm{\mathbf{\theta}})p(\bm{\mathbf{\theta}})d\bm{\mathbf{\theta}}$ is a Polya distribution which is approximated with a uniform discrete distribution over levels.

4. System Implementation

Our distributed training system for hLDA consists of machine-level and thread-level parallelism, as shown in Fig. 2. On the machine level, we use MPI to synchronize the tree structure and the counts $(\bm{\mathbf{C}}_{\mathcal{I}},\bm{\mathbf{C}}_{\mathcal{C}},\bm{\mathbf{m}})$ across machines; and on the thread-level, a number of threads concurrently read and update the local counts.

In order to implement the system efficiently, several challenges must be addressed. Firstly, for each worker thread, the data layout and algorithm should be organized in a vectorization-friendly way for efficient memory access and computation. Moreover, expensive computation such as logarithms should be avoided as much as possible. Secondly, for efficient multi-thread parallelism, the shared data structure should be lock-free. Finally, the communication strategy need to be chosen carefully to minimize communication overhead and maximize the freshness of the counts. We now present solutions to address these challenges.

4.1. Vectorization and Pre-processing

We first discuss how to organize the computation for vectorization and pre-processing. The most time-consuming part of PCGS is the sampling of $\bm{\mathbf{c}}$ according to Eq. (4), or more concretely, computing $f_{I}(d,t)$ for each $t\in\mathcal{I}$ and $f_{C}(d,t)$ for each $t\in\mathcal{C}$. Because both $f_{I}(d,t)$ and $f_{C}(d,t)$ are very close to zero, we compute their logarithms. Rewrite Eq. (3):

[TABLE]

where $\mathcal{W}_{dl}$ is the set of all tokens in document $d$ that are assigned to level $l$, which can be computed by bucket sorting $\bm{\mathbf{z}}_{dn}$’s along with their positions $n$’s.

Eq. (5) can be vectorized as $\log\bm{\mathbf{f}}_{I}(d,l)=\sum_{v\in\mathcal{W}_{dl}}(\log\bm{\mathbf{\phi}}^{\top})_{v,\mathcal{I}_{l}}$, where $\mathcal{I}_{l}\subset\mathcal{I}$ is the subset of topic ids on level $l$, and $\bm{\mathbf{f}}_{I}(d,l)$ is the vector of $f_{I}(d,t)$’s for all topics in $\mathcal{I}_{l}$. The matrix $\log\bm{\mathbf{\phi}}^{\top}$ is the transpose of the element-wise logarithm of $\bm{\mathbf{\phi}}$, which is pre-processed to avoid the computation of logarithm. We store the matrix in row-major order, so that accessing some topics $\mathcal{I}_{l}$ for a certain word $v$ is continuous as long as the IDs $\mathcal{I}_{l}$ are continuous for each level, which is easy to achieve by sorting the topic ids according to their levels, since $\mathcal{I}$ do not change when sampling $\bm{\mathbf{c}}$. Therefore, computing $\log\bm{\mathbf{f}}_{I}(d,l)$ is fully vectorized by just adding the slices of the matrix $\log\bm{\mathbf{\phi}}^{\top}$ indexed by $\mathcal{W}_{dl}$ and $\mathcal{I}_{l}$.

Similarly, Eq. (2) can be rewritten as:

[TABLE]

where we convert the logarithm of multivariate beta function as the sum of logarithms (the derivation details can be found in the appendix). The term $h_{t}=\log\Gamma(C_{t}^{\neg d}+V\beta_{l_{t}})-\log\Gamma(C_{t}+V\beta_{l_{t}})$, and in $\mathcal{W}_{dl}^{\prime}$ we assign each token with an offset indicating which time does this word appear, e.g., if a word $v$ is in $\bm{\mathbf{w}}_{d}$ for three times, we put $(v,0)$, $(v,1)$ and $(v,2)$ into $\mathcal{W}_{dl}^{\prime}$.

Again, we vectorize the computation of $f_{C}(d,t)$ by computing $\log\bm{\mathbf{f}}_{C}(d,l)$ for each $l$, which is the vector of $\log f_{C}(d,t)$’s for $t\in\mathcal{C}_{l}$. $\mathcal{C}_{l}\subset\mathcal{C}$ is the subset of topic ids on level $l$, that can change during the sampling of $\bm{\mathbf{c}}$ due to the birth and death of topics. For efficient vectorization, the counts need to be stored such that $C_{tv}^{\neg d}$ is continuous for all $t\in\mathcal{C}_{l}$. To achieve this, we store separate count matrices for each level. For level $l$, $\bm{\mathbf{C}}_{\mathcal{C}_{l}}^{\top}$ is stored, which is made by concatenating all the columns $t\in\mathcal{C}_{l}$ of $\bm{\mathbf{C}}_{t}^{\top}$. When a new topic on level $l$ is created, we append it to $\bm{\mathbf{C}}_{\mathcal{C}_{l}}^{\top}$ as the rightmost column. The removal of columns is deferred after the sampling of $\bm{\mathbf{c}}$ finishes, and the result will not be affected since the dead topics correspond to zero columns in $\bm{\mathbf{C}}_{\mathcal{C}_{l}}^{\top}$.

Unlike computing $f_{I}(d,t)$, the logarithm in Eq. (6) cannot be pre-processed since the count $C_{tv}^{\neg d}$ changes during the sampling of $\bm{\mathbf{c}}$. Therefore, $f_{C}(d,t)$ is much more expensive to compute than $f_{I}(d,t)$, supporting our argument on the inefficiency of CGS in Sec. 3.1. Fortunately, for PCGS, the computation of logarithm is avoided as much as possible by keeping $\mathcal{C}$ a small set. To further accelerate the computation, we use the SIMD enabled Intel VML library for logarithms.

4.2. Concurrent Data Structures

In our system, the collapsed count matrices $\bm{\mathbf{C}}^{\top}_{\mathcal{C}}$ are concurrently read and updated by the worker threads, and the number of columns (topics) can grow over time. Since there are a lot of reads, the matrix must be read efficiently, i.e., lock free. Meanwhile, the consistency can be relaxed since a small deviation of the counts will not affect the result much. Therefore, we only ask the matrices to have eventual consistency, i.e., the values of the matrices should be eventually correct if no new updates are given. We adopt atomic writes to preserve eventual consistency, while the reads are relaxed as non-atomic operations, to maximize the reading performance.

The dynamic number of columns makes the implementation challenging. The straightforward implementation for growing the matrix involves allocating a new memory region, copying the original content to the new memory, and deallocating the original memory. However, this implementation cannot achieve eventual consistency because the updates during copying will not be incorporated.

Inspired by the lock-free design of a concurrent vector (Dechev et al., 2006), which is a one-dimensional version of our matrix, we provide an efficient implementation of the concurrent matrix. Internally, it holds a list of matrix blocks, where the $i$-th matrix block has the size $R\times 2^{c+i-1}$, while $c$ is a constant. The first matrix block represents the $[0,2^{c})$-th columns of the original matrix, the second matrix block represents the $[2^{c},3\times 2^{c})$-th columns of the original matrix, and so on. If there is a growing request that exceeds the current capacity, we allocates the next matrix block on the list. For every reading and updating request, the requested $(r,c)$ coordinate is converted to the $(r,b,c^{\prime})$ coordinate, where $b$ is the index of the matrix block on the list and $c^{\prime}$ is the column index within the matrix block. The coordinate conversion can be achieved with a BSR instruction in modern x86 systems (Dechev et al., 2006). Finally, to improve the locality, we defragment after each PCGS iteration, i.e., deallocating all the matrix blocks and concatenating their content to form a single larger matrix block.

4.3. Communication

For PCGS, we need to synchronize the instantiated count $\bm{\mathbf{C}}_{\mathcal{I}}$ across machines once per PCGS iteration, and the collapsed counts ($\bm{\mathbf{C}}_{\mathcal{C}}$, $\bm{\mathbf{m}}$) as frequently as possible. We now present an implementation of the synchronization by MPI.

Firstly, we synchronize $\bm{\mathbf{C}}_{\mathcal{I}}$ by the MPI_Allreduce operation. There are many approaches to synchronizing $\bm{\mathbf{C}}_{\mathcal{C}}$ and $\bm{\mathbf{m}}$. One possible solution is the parameter server as shown in Fig. 1(a). However, while parameter server typically assumes the amount of communication is high and the network bandwidth is the bottleneck, the amount of our PCGS communication is low and our main focus is on the latency, which determines how fresh the count is. The parameter server first merges the changes from individual worker machines at the server, and then pushes the new state to the workers. While the merging decreases the amount of communication, it increases the latency by sending the change through the server.

To optimize the latency, we design a decentralized communication strategy, in which all the worker nodes directly send their changes to all the other nodes, as illustrated in Fig. 1(c). There is a synchronization thread on each worker machine with a to_send buffer, a sending buffer and a receiving buffer. The worker threads write their changes to the to_send buffer, and the synchronization threads periodically exchange the content in the to_send buffer across machines, as follows: (1) Atomically exchange the to_send buffer and sending buffer, clear the new to_send buffer; (2) Gather the content of all sending buffers to the receiving buffer, by a MPI_Allgatherv operation; (3) Merge all the changes in the receiving buffer to the local copy of collapsed counts $\bm{\mathbf{C}}_{\mathcal{C}}$ and $\bm{\mathbf{m}}$.

5. Experiments

We evaluate the quality, efficiency and scalability of our algorithm and system on several datasets, including NIPS, NYTimes, PubMed from the UCI machine learning repository (Asuncion and Newman, 2007), and two subsets of the ClueWeb12 dataset (Project, 2013) (Table 1). The experiments are conducted on the Tianhe-2 supercomputer, which has two 12-core Xeon E5-2692v2 CPUs per node and an InfiniBand network. Our quantitative and qualitative results demonstrate the promise.

We quantitatively compare the quality of the inferred models by predictive log-likelihood using the document completion approach (Wallach et al., 2009). The corpus is divided as a training corpus $\bm{\mathbf{w}}^{\mbox{t}}$ and a testing corpus, and the testing corpus is further divided as an observed corpus $\bm{\mathbf{w}}^{\mbox{o}}$, which contains a random half of the tokens for each document in the testing corpus; and a heldout corpus $\bm{\mathbf{w}}^{\mbox{h}}$ of the other half of the tokens. The predictive log-likelihood is defined as $p(\bm{\mathbf{w}}^{\mbox{h}}|\bm{\mathbf{w}}^{\mbox{o}},\bm{\mathbf{\phi}})$, where the model $\bm{\mathbf{\phi}}$ is inferred from the training corpus $\bm{\mathbf{w}}^{\mbox{t}}$, and is approximated with a Monte-Carlo integration:

[TABLE]

where $\bm{\mathbf{c}}_{d}^{(s)}$ and $\bm{\mathbf{\theta}}_{d}^{(s)}$ are the samples from the posterior distribution $p(\bm{\mathbf{c}}_{d}^{(s)},\bm{\mathbf{\theta}}_{d}^{(s)}|\bm{\mathbf{w}}^{\mbox{o}},\bm{\mathbf{\phi}})$, which can be obtained with Gibbs sampling, and $S$ is the number of samples. Finally, we convert predictive log-likelihood to predictive perplexity

[TABLE]

where a lower perplexity score indicates a better model.

5.1. Quality of Inference

We first compare the model inferred by CGS and our proposed BGS and PCGS, and examine the effect of the initialization strategy (Sec. 3.4). We also include a comparison with the open source implementation for hLDA, hlda-c, which is a CGS algorithm but has a stick-breaking prior on $\bm{\mathbf{\theta}}$ instead of a Dirichlet prior (Blei et al., 2010).

Unlike parametric models, where the number of topics is fixed, nonparametric models such as hLDA produce different numbers of topics for different runs and various inference algorithms, even with the same hyper-parameter setting. It is not fair to directly compare the perplexity of two models with different numbers of topics. For a fair comparison, we choose a rich set of hyper-parameter configurations, run the algorithms for all these configurations, and plot the perplexity against the number of topics as in Fig. 3. In this experiment, we train a 4-layer model (i.e., $L=4$) on the NYTimes (subset) dataset and the NIPS dataset, and $\bm{\mathbf{\beta}}=(\beta_{0},0.5\beta_{0},0.25\beta_{0},0.25\beta_{0})$, where $\beta_{0}$ is chosen from $\{e^{-4.0},e^{-3.5},\dots,e^{2.0}\}$, $\gamma$ is chosen from $\{e^{-6.0},e^{-5.5},\dots,e^{0.0}\}$, and $\bm{\mathbf{\alpha}}=0.2\times\mathbf{1}$.

By comparing the perplexity produced by different algorithms, we have a number of observations:

• •

CGS and PCGS have similar quality, while BGS has worse results. This agrees with our previous analysis (Sec. 3.2) that BGS suffers from label switching and local optimum.

• •

Our initialization strategy helps obtain better results for both CGS and PCGS.

• •

Our result is not worse (actually better) than hlda-c. The discrepancy attributes to the different choice of prior on $\bm{\mathbf{\theta}}$.

5.2. Efficiency

We compare the efficiency of our algorithms against hlda-c. We run the serial version of all the algorithms for 70 iterations on the NYTimes (subset) dataset while setting $\bm{\mathbf{\beta}}=(1,0.5,0.25,0.1)$, and $\gamma$ is tuned to keep the number of topics around 300. The timing result is shown in Table 2. Our CGS implementation is 48 times faster than hlda-c. The significant gain of efficiency attributes to our vectorization and the conversion of the logarithm of gamma function to the sum of logarithms in Sec. 4.1 and the appendix. PCGS is 2.3 times faster than CGS, and BGS is 1.3 times faster than PCGS. These results match our analysis in Sec. 4.1 on that BGS and PCGS are more efficient than CGS. Overall, our PCGS implementation is 111 times faster than hlda-c.

Combining the results on inference quality and efficiency, we find PCGS to be a good trade-off between quality and efficiency by providing the CGS-level quality within BGS-level time consumption.

5.3. Sensitivity of Parameters

We now examine the impact of the hyper-parameters $M,I$ and $S$, which control the behavior of PCGS.

**Impact of $M$: ** $M$ is the threshold of the number of visits that we decide whether the topic distribution $\bm{\mathbf{\phi}}_{t}$ of a topic is rapid-changing or slow-changing. To investigate its effect, we run PCGS on the NYTimes dataset, setting $\bm{\mathbf{\beta}}=(1.0,0.5,0.25,0.125)$,and varying $M\in\{2^{0},\dots,2^{16},\infty\}$, while tuning $\gamma$ to keep the number of topics around 500. PCGS becomes CGS when $M=\infty$, and approaches BGS when $M\rightarrow 0$. As shown in Fig. 6, the perplexity goes down and the time consumption goes up as $M$ grows. We also find that $M\in[2^{6},2^{9}]$ provides a good trade-off between efficiency and quality. When $M=64$, there are 427 slow-changing topics which covers 99.7% documents, so the change of rapid-changing topic counts (amount of communication) is kept small.

**Impact of $I$: ** $I$ is the number of initializing iterations to sample from $p(\bm{\mathbf{c}}|\bm{\mathbf{w}})$. We run PCGS on the NYTimes (subset) dataset, setting $\beta_{0}=1$, and varying $I$ and $\gamma$. It can be seen from Fig. 6 that the perplexity steadily decreases for large $I$, which again shows that our initialization strategy is helpful. We select a moderate $I=32$ for all the experiments.

**Impact of $S$: ** The hyper-parameter $S$ is the number of Monte-Carlo samples to approximate $p(\bm{\mathbf{c}}|\bm{\mathbf{w}})$. When $S\rightarrow\infty$, we directly sample from $p(\bm{\mathbf{c}}|\bm{\mathbf{w}})$ in the first $I$ iterations. We run PCGS on the NYTimes (subset) dataset, with $\bm{\mathbf{\beta}}=(1.0,0.5,0.25,0.25)$ and $\gamma=10^{-40}$, and vary $S$ from 1 to 128. As shown in Fig. 6, $S$ has little impact on both the number of topics and the perplexity, implying that a small $S$, e.g., $S=5$, is adequate.

5.4. Scalability

Our experiments on scalability are in two folds: whether the quality of inference is affected by parallelization; and how good is the speedup. We first study the multi-thread setting, where the number of threads varies from 1 to 12 on the NYTimes corpus. The result is shown in Fig. 7(a), where we observe that there is no apparent increase of perplexity as the number of threads grows. The speedup with 12 threads is 8.56. The probable reasons of imperfect speedup include serial region, contention for atomic variables, and limited memory bandwidth.

For the multi-machine setting, there are two CPUs per machine. We run our implementation on the PubMed corpus on 1 to 10 CPUs as shown in Fig. 7(b), and on the larger ClueWeb12 (small) corpus for 10 to 100 CPUs as shown in Fig. 7(c). The speedup is 8.5 from 1 to 10 CPUs, and 7.15 from 10 to 100 CPUs. The perplexity is slightly affected by parallelization when the number of CPUs exceeds 7 and 80 on the two datasets, respectively, indicating that the dataset is not large enough to utilize that many CPUs.

Finally, to demonstrate the scalability, we learn a model with 1,722 topics of the 131-million-document ClueWeb12 (large) corpus with 50 machines, and the inference finishes in 7 hours. The results will be presented for qualitative evaluation in the next section.

5.5. Qualitative Analysis

We now demonstrate the topic hierarchy obtained from the ClueWeb12 (large) corpus, which is a crawl of web pages. The corpus is obtained by tokenizing the original ClueWeb12 dataset, randomly selecting about 30% documents, truncating the vocabulary size to 100,000 and keeping only the documents whose length is between $[50,500]$. We show the selected parts of the obtained tree in Fig. 8, where some topics whose number of occurrences does not pass a particular threshold are filtered out, and the font size of words is proportional to the 4th root of their frequency in the topic. The tree has 5 levels in total.111The visualization demo is available online at http://ml.cs.tsinghua.edu.cn/~jianfei/scalable-hlda.html.

Fig. 8(a) shows some selected topics on the first 3 levels. The root node contains the most commonly used words shared by all the documents. The second level contains a variety of general topics, such as “software”, “travel” and “city”, and the third level has more detailed concepts, e.g., the “city” topic on the second level splits as “shopping”, “city names”, and “locations”. We further show the topic subtrees of all the layers rooted at the highlighted nodes to examine the fine-grained concepts. For example, in Fig. 8(b) the “travel” topic is divided as “islands”, “India” and “vacations”, and the leaf level contains specific concepts, such as “ferry” and “diving”, which are correctly placed under the “islands” topic. In Fig. 8(c), the “computer” topic is divided as “website”, “windows”, “vps”, “programming”, “linux” and “forum”. To our knowledge, this is the first time that hLDA is applied to large-scale web data, and the results demonstrate our ability on automatically learning topic hierarchy from web data.

6. Conclusions and Discussions

We present a partially collapsed Gibbs sampling (PCGS) algorithm for the hierarchical latent Dirichlet allocation model, which is a combination of the collapsed weight algorithm and instantiated weight algorithm. The major feature of PCGS is that it is scalable and has high-quality inference. We also present an initialization strategy to further improve the model quality. To make PCGS scalable and efficient, we propose vectorization and pre-processing techniques, concurrent data structures, and an efficient communication strategy. The proposed algorithm and system are scalable to hundreds of millions of documents, thousands of topics, and thousands of CPU cores.

In the future, we plan to extend our method to the more sophisticated nested HDP model (Paisley et al., 2015; Ahmed et al., 2013). Developing sampling algorithms with sub-linear time complexity w.r.t. the number of topics via hashing is also an interesting direction.

Appendix A Derivation Details

A.1. Derivation of PCGS updates

Rewrite the joint distribution in Sec. 3.1 as: $p(\bm{\mathbf{w}},\bm{\mathbf{z}},\bm{\mathbf{\theta}},\bm{\mathbf{c}},\bm{\mathbf{\phi}}_{\mathcal{C}},\bm{\mathbf{\phi}}_{\mathcal{I}})$ =$\prod_{t\in\mathcal{C}}$ $p(\bm{\mathbf{\phi}}_{t})$ $\prod_{t\in\mathcal{I}}$ $p(\bm{\mathbf{\phi}}_{t})$ $\prod_{d=1}^{D}$ $p(\bm{\mathbf{c}}_{d}|\bm{\mathbf{c}}_{<d})$ $p(\bm{\mathbf{\theta}}_{d})$ $\prod_{d=1}^{D}$ $\prod_{n=1}^{N_{d}}$ $p(z_{dn}|\bm{\mathbf{\theta}})$ $p(w_{dn}|z_{dn},\bm{\mathbf{\phi}}_{\mathcal{C}},\bm{\mathbf{\phi}}_{\mathcal{I}}).$

Integrating out $\bm{\mathbf{\phi}}_{\mathcal{C}}$ and $\bm{\mathbf{\theta}}$, we have the marginal distribution $p(\bm{\mathbf{w}},\bm{\mathbf{z}},\bm{\mathbf{c}},\bm{\mathbf{\phi}}_{\mathcal{I}})$=\prod_{d=1}^{D}\big{[}p(\bm{\mathbf{c}}_{d}|\bm{\mathbf{c}}_{<d})\frac{B(\bm{\mathbf{C}}_{d}+\bm{\mathbf{\alpha}})}{B(\bm{\mathbf{\alpha}})}\big{]} $\prod_{t\in\mathcal{C}}\frac{B(\bm{\mathbf{C}}_{t}+\beta_{l_{t}})}{B(\beta_{l_{t}})}$ $\prod_{t\in\mathcal{I}}$ \big{[}\mbox{Dir}(\bm{\mathbf{\phi}}_{t};\beta_{l_{t}}) \prod_{v=1}^{V}(\phi_{tv})^{C_{tv}}\big{]}.

Utilizing the identity $\frac{B(\bm{\mathbf{\alpha}}+\bm{\mathbf{e}}_{k})}{B(\bm{\mathbf{\alpha}})}=\frac{\alpha_{k}}{\sum_{k}\alpha_{k}}$, where $\bm{\mathbf{e}}_{k}$ is a coordinate vector, we can derive the Gibbs sampling updates:

Sample $\bm{\mathbf{z}}$: Keeping only the terms relevant with $z_{dn}$, we have $p(z_{dn}=l|w_{dn}=v,\bm{\mathbf{w}},\bm{\mathbf{z}}_{\neg dn},\bm{\mathbf{c}},\bm{\mathbf{\phi}}^{\mathcal{I}})$ $\propto$ $\frac{B(\bm{\mathbf{C}}_{d}+\bm{\mathbf{\alpha}})}{B(\bm{\mathbf{\alpha}})}$ $\prod_{t\in\mathcal{C}}$ $\frac{B(\bm{\mathbf{C}}_{t}+\beta_{l_{t}})}{B(\beta_{l_{t}})}$ $\prod_{t\in\mathcal{I}}$ $\prod_{v=1}^{V}\phi_{tv}^{C_{tv}}$ $\propto$ $\frac{B(\bm{\mathbf{C}}_{d}^{\neg dn}+\bm{\mathbf{C}}_{d}^{dn}+\bm{\mathbf{\alpha}})}{B(\bm{\mathbf{C}}_{d}^{\neg dn}+\bm{\mathbf{\alpha}})}$ $\prod_{t\in\mathcal{C}}$ $\frac{B(\bm{\mathbf{C}}_{t}^{\neg dn}+\bm{\mathbf{C}}_{t}^{dn}+\beta_{l_{t}})}{B(\bm{\mathbf{C}}_{t}^{\neg dn}+\beta_{l_{t}})}$ $\prod_{t\in\mathcal{I}}$ $\prod_{v=1}^{V}$ $\phi_{tv}^{C_{tv}^{\neg dn}+C_{tv}^{dn}}$ $\propto(C_{dl}^{\neg dn}+\alpha_{l})$ $\begin{cases}\phi_{c_{dl},v}&c_{dl}\in\mathcal{I},\\ \frac{C_{c_{dl},v}^{\neg dn}+\beta_{l_{t}}}{C_{c_{dl}}^{\neg dn}+V\beta_{l_{t}}}&c_{dl}\in\mathcal{C}.\end{cases}$

Sample $\bm{\mathbf{c}}$: $p(\bm{\mathbf{c}}_{d}=\bm{\mathbf{c}}|\bm{\mathbf{w}},\bm{\mathbf{z}},\bm{\mathbf{c}}_{\neg d},\bm{\mathbf{\phi}}^{\mathcal{I}})$ $\propto$ $p(\bm{\mathbf{c}}_{d}|\bm{\mathbf{c}}_{\neg d})$ $\prod_{t\in\mathcal{C}}$ $\frac{B(\bm{\mathbf{C}}_{t}+\beta_{l_{t}})}{B(\beta_{l_{t}})}$ $\prod_{t\in\mathcal{I}}$ $\prod_{v=1}^{V}\phi_{tv}^{C_{tv}}$ $\propto$ $p(\bm{\mathbf{c}}_{d}|\bm{\mathbf{c}}_{\neg d})$ $\prod_{t\in\mathcal{C}}$ $\frac{B(\bm{\mathbf{C}}_{t}^{\neg d}+\bm{\mathbf{C}}_{t}^{d}+\beta_{l_{t}})}{B(\bm{\mathbf{C}}_{t}^{\neg d}+\beta_{l_{t}})}$ $\prod_{t\in\mathcal{I}}$ $\prod_{v=1}^{V}$ $\phi_{tv}^{C_{tv}^{d}}$ $\propto$ $\mbox{nCRP}(\bm{\mathbf{c}};\bm{\mathbf{\gamma}},\bm{\mathbf{c}}_{\neg d})$ $\prod_{l=1}^{L}$ $\begin{cases}f_{I}(d,c_{l})&c_{l}\in\mathcal{I},\\ f_{C}(d,c_{l})&c_{l}\in\mathcal{C}.\end{cases}$

Sample $\bm{\mathbf{\phi}}^{\mathcal{I}}$: For $t\in\mathcal{I}$, draw $\bm{\mathbf{\phi}}_{t}\sim\mbox{Dir}(\beta_{l_{t}}+\bm{\mathbf{C}}_{t})$.

A.2. Derivation of computing $\log f_{C}(d,t)$

We have $\log f_{C}(d,t)$ = $\log\frac{B(\bm{\mathbf{C}}_{t}^{\neg d}+\bm{\mathbf{C}}_{t}^{d}+\beta_{l_{t}})}{B(\bm{\mathbf{C}}_{t}^{\neg d}+\beta_{l_{t}})}$ =\big{[}\sum_{v=1}^{V}\sum_{i=0}^{C_{tv}^{d}-1}\log(C_{tv}^{\neg d}+i+\beta_{l_{t}})\big{]}+h_{t} =$\sum_{(v,o)\in\mathcal{W}^{\prime}_{dl}}\log(C_{tv}^{\neg d}+o+\beta_{l_{t}})+h_{t}$.

Acknowledgements.

The work was supported by the National Basic Research Program (973 Program) of China (No. 2013CB329403), National NSF of China (Nos. 61620106010, 61322308, 61332007), the Youth Top-notch Talent Support Program, Tsinghua Tiangong Institute for Intelligent Computing, and Special Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (the second phase). J Lu and S Liu are supported by National NSF of China (No. 61672308).