Scalable Training of Hierarchical Topic Models
Jianfei Chen† , Jun Zhu†∗, Jie Lu‡ , and Shixia Liu‡
†Dept. of Comp. Sci. & Tech., BNRist Center, State Key Lab for Intell. Tech. & Sys.
‡School of Software, BNRist Center, State Key Lab for Intell. Tech. & Sys.

Tsinghua University, Beijing, 100084, China{chenjian14,luj15}@mails.tsinghua.edu.cn; {dcszj,shixia}@tsinghua.edu.cn“computer”

ABSTRACT

file
files
serverdata
error

Large-scale topic models serve as basic tools for feature
extraction and dimensionality reduction in many practicalapplications. As a natural extension of flat topic models,hierarchical topic models (HTMs) are able to learn topics
of different levels of abstraction, which lead to deeper understanding and better generalization than their flat counterparts. However, existing scalable systems for flat topic
models cannot handle HTMs, due to their complicated data
structures such as trees and concurrent dynamically growing
matrices, as well as their susceptibility to local optima.
In this paper, we study the hierarchical latent Dirichlet allocation (hLDA) model which is a powerful nonparametric
Bayesian HTM. We propose an efficient partially collapsed
Gibbs sampling algorithm for hLDA, as well as an initialization strategy to deal with local optima introduced by treestructured models. We also identify new system challenges
in building scalable systems for HTMs, and propose efficient
data layout for vectorizing HTM as well as distributed data
structures including dynamic matrices and trees. Empirical studies show that our system is 87 times more efficientthan the previous open-source implementation for hLDA,
and can scale to thousands of CPU cores. We demonstrate
our scalability on a 131-million-document corpus with 28
billion tokens, which is 4-5 orders of magnitude larger thanpreviously used corpus. Our distributed implementation can
extract 1,722 topics from the corpus with 50 machines in just
7 hours.

“installation”

linux
packageinstallrpm
perl

driverwindows
drivers
download
bios“C++” “php/Java”“Linuxsoftware” “notebook”

const
virtualint
inlinedefinition

java
apacheproject
mysql
repository

notebook
driverdrivers
downloaddellphpmethodarray
objectjava

topic is a unigram distribution over words, and the highprobability words often present strong semantic correlation.
As a method for dimensionality reduction and feature extraction, topic models have been widely used in information retrieval [35], text analysis [7, 41], information visualization [33], online advertising [34], recommendation system [31], and many other application areas.
Practical applications call for scalable and efficient topic
modeling systems to handle big data. For example, online advertisement systems extract topics from billions ofsearch queries [34], and recommendation systems [1] needto handle millions of users and items. Various efforts has
been made to develop scalable topic modeling systems, including asynchronous distributed data parallel training [1,17], hybrid data-and-model-parallel training [37, 36], embarrassingly parallel BSP training [10, 38, 39], and GPUaccelerated training [40, 16]. These topic modeling systems
mainly handle the partition of the data and model and thesynchronization of the count matrix across machines.
These works are for flat topic models, such as latent Dirichlet allocation (LDA) [6], which assumes that all the topics
have the same level of abstraction. In reality, however, topics
are naturally organized in a hierarchy [21]. See Fig. 1 as an
example (a more complete version is in experiments). When

INTRODUCTION

Topic models are popular machine learning tools that extract a set of latent topics from an input text corpus. Each∗

“Linux”

class
stringpublicvoid
int

Figure 1: A subtree of the topic hierarchy learned from a
corpus with 131-million web pages (See experiments for details). Each node is a topic with its top-5 most frequentwords. The full tree has 1,722 topics. The text in the quote
is summarized by humans.

PVLDB Reference Format:
Jianfei Chen, Jun Zhu, Jie Lu and Shixia Liu. Scalable Trainingof Hierarchical Topic Models. PVLDB, 11(7): 826-839, 2018.
DOI: https://doi.org/10.14778/3192965.3192972

1.

“programming”

Corresponding author.

Permission to make digital or hard copies of all or part of this work forpersonal or classroom use is granted without fee provided that copies are
not made or distributed for profit or commercial advantage and that copies
bear this notice and the full citation on the first page. To copy otherwise, torepublish, to post on servers or to redistribute to lists, requires prior specificpermission and/or a fee. Articles from this volume were invited to present
their results at The 44th International Conference on Very Large Data Bases,
August 2018, Rio de Janeiro, Brazil.
Proceedings of the VLDB Endowment, Vol. 11, No. 7
Copyright 2018 VLDB Endowment 2150-8097/18/03... $ 10.00.
DOI: https://doi.org/10.14778/3192965.3192972

826

a topic about “programming” is observed in a document,
we are also likely to observe the more general topic suchas “computer” in the same document. By capturing such
relationships, hierarchical topic models can achieve deeper
understanding and better generalization [2, 21] of the corpus than the flat models. There are many approaches tolearn a topic hierarchy. For example, Google’s Rephil [19]puts a hierarchical noisy-or network on the documents; thesuper-topic approach learns topics of topics [18, 23]; andthe nested Chinese Restaurant Process (nCRP) [5, 2, 21]approach utilizes the nCRP as a prior on topic hierarchies,
with hierarchical latent Dirichlet allocation (hLDA) [5] as apopular example. These models have been successfully applied to document modeling [21], online advertising [19] andmicroblog location prediction [2], outperforming flat models.
Despite the success of hierarchical topic models (HTMs)
on small scale corpora, they have not been deployed on practical applications due to the lack of scalable systems as wellas efficient algorithms. Developing a scalable system for hierarchical topic models is challenging because of the complicated data structures, such as trees and dynamically growing matrices (e.g., in hLDA). We need to develop distributed
versions of these data structures while maintaining the efficiency and consistency. Moreover, the access patterns ofthese data structures are more complicated than the countmatrix in flat topic models, making efficient vectorization
more difficult. There are also algorithmic challenges, e.g.,dealing with the local optima during the learning.
In this paper, we make an initial attempt in scaling up
hierarchical topic models. We study the hierarchical latent
Dirichlet allocation (hLDA) [5] model, which can automatically learn the hierarchical topical structure with Gibbs sampling [14] by utilizing an nCRP prior. In hLDA, topics form
a tree with an nCRP prior, while each document is assignedwith 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. We propose a novel partially collapsed Gibbs
sampling (PCGS) algorithm that is more efficient than theoriginal collapsed Gibbs sampling (CGS) algorithm [5]. To
improve the model quality, we analyze the cause of local
optima and propose an effective initialization strategy. On
the system side, we propose a data layout for efficient vectorization, and design a concurrent dynamic matrix as well
as a distributed tree to support hLDA’s complicated data
structures for distributed training.
We extensively evaluate the model quality, efficiency andscalability of our algorithm and system. The experimental results show that our PCGS algorithm and the vectorization strategy lead to an 87 times single-thread speedup
than the existing open source implementation hlda-c [4],while having better model quality due to our initialization
strategy. Multi-threading and distributed computing further speed up our system, and we analyze the speedup andimpact to model quality as the amount of computational
resource grows. Finally, we demonstrate the scalablity ofour distributed system by extracting 1,722 topics from a
131-million-document corpus with 50 machines in 7 hours,
where the corpus is 4-5 orders of magnitude larger than thelargest corpus used for hLDA before. We believe that ouralgorithm and system can be extended to more hierarchical
topic models, such as the nested hierarchial Dirichlet process(nHDP) [21] model.

4

4

c1
3 c3c 1
4
c2
3
𝑝=
4.5

1
𝑝=
4.5

3
c1

0.5
𝑝=
4.5

1

c3
1

2c 1
2

3
0.5
×
4.5 3.53
2𝑝=×
4.5 3.5

𝑝=

𝑝=

3
1
×
4.5 3.5

(a)

c4

𝑝=

𝑝=

1
0.5
×
4.5 1.51
1
×4.5 1.5𝑝=

0.5 0.5×
4.5 0.5(b)(c)
Figure 2: The (a) CRP (b) nCRP distribution forp(c5 |c1 , . . . , c4 ) given previous 4 paths, with γ1 = γ2 = 0.5.
Numbers on the nodes are visited counts mt , and dashednodes are unvisited. p is the probability for each path. (c)Graphical model for hLDA.

2.

HIERARCHICAL LDA

We first review nCRP and hLDA for learning a topic hierarchy as well as the collapsed Gibbs sampling inference
algorithm proposed by Blei et al. [5].

2.1

Nested Chinese Restaurant Process

Nested Chinese Restaurant process (nCRP) [5] represents
a powerful nonparametric Bayesian method to learn a tree
structure, whose width and depth are unbounded. Following
previous work in hierarchical topic modeling [5], we present
the nCRP with (finite) L-levels, and the infinite-depth tree
can be approximated by setting L to a large number.
Suppose there is a tree with L levels, where each node
except the leaves has infinite children. A 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(c1 , c2 , . . . ) on the tree, where each path cd ∈ NL
+ consistsof L node IDs from the root to a certain leaf. nCRP defines
the joint distribution of paths asp(c1 , c2 , . . . ) =∞Yp(cd |c<d ),d=1where the subscript < d stands for all the possible indicesthat are smaller than d, i.e., c<d = {c1 , . . . , cd−1 }. Given
previous paths c<d , nCRP defines a probability distributionp(cd |c<d ) on the next path asp(cd |c<d ) =

L
Y
l=2

827

p(cdl |cd,l−1 , c<d ),which can be interpreted as a generative story: a person
starts from the root node 1 at level 1. Each time he walks toa child of the current node until he reaches a leaf node, where
the probability of going from cd,l−1 to cdl is p(cdl |cd,l−1 , c<d ).
The probability of the path he visited is then p(cd |c<d ).
The transition probability p(cdl |cd,l−1 , c<d ) follows the
Chinese Restaurant Process (CRP) [27]. We denote a visitof node t, if any previous path cd0 (d0 < d) passes through it.
The number of visits mt := #{(d0 , l)|cd0 = t, d0 < d}, where#{·} denotes the cardinality of a set. We assume that thepaths c<d have already visited T nodes in the infinite tree,with the node ID from 1 to T , and the current node cd,l−1
has K visited children with the IDs t1 , . . . , tK . This corresponds to a lazy strategy of assigning IDs to nodes: a node
does not have an ID until it is visited. In CRP, the probability of going to each child is proportional to its previous
times of visitmtip(cdl = ti |·) = γ +m,
k = 1, . . . , Ki
cd,l−1l(1)p(cdl = T + 1|·) = γ +mγl,clwhere Dir(·) is the Dirichlet distribution, and α and β are
Dirichlet hyper-parameters. See Fig. 2(c) for the probabilistic graphical model representation of hLDA.
There are two special cases of hLDA. When the tree degenerates to a chain, hLDA recovers the vanilla latent Dirichlet allocation (LDA) [6] with L topics, and when the tree has
two levels and the probability of assigning to the first levelθd1 is close to zero, hLDA recovers the Dirichlet Process
Mixture Model (DPMM) [20].

where γl is a hyper-parameter controlling the probabilityof going to an unvisited child. There are infinite unvisitedchildren and we may choose any of them to visit, but wealways assign the ID T +1 to the chosen one. This operationof assigning an ID to an unvisited node is referred as thegeneration of a new child in this paper.
Given c<d , the number of different possible paths for cd is
T , which is finite. There are Tleaf paths ending at existingleaves, and T − Tleaf paths forking from an internal visited
node. Fig. 2(a) and Fig. 2(b) illustrate the possible paths
and their probabilities for CRP and nCRP, respectively.

• The document-level count adl = #{n|zdn = l} is thenumber of tokens on level l in document d.

2.3

d,l−1

2.2

Notation of Counts

Before introducing the inference algorithm for hLDA, weneed to define some notations of counts. We will encounter
four types of counts. The first is the node visit count mt
defined in Sec. 2.1, which represents how many documentsthat visit node t. The remaining three counts are token
counts. Each token in the corpus has four attributes: its
document ID d, word ID wdn , level zdn , and the node ID
cd,zdn it is assigned to. The following counts represent the
number of tokens that have specific attributes• The topic-word count btv = #{(d, n)|wdn = v∧cd,zdn =t} is the number of tokens of the word v with topic t.
• The topic count st =with topic t.

P

v btv

is the number of tokens

Furthermore, bold uppercase letters denote the matrix formedby the counts, e.g., A is a D × L matrix composed by adn ’s.
Let ad denote the d-th row vector, a:,l be the l-th column
vector. Let m and s be T -dimensional vectors composed of
mt and st , respectively. Finally, d stands for document d,¬d represents all documents excluding document d, and ¬dnmeans all the tokens excluding token (d, n). Using theseas the superscript of counts means only counting the objects among the given documents / tokens. For instance,0
bdtv = #{n|wdn = v ∧cd,zdn = t} and b¬dtv = #{(d , n)|wd0 n =v ∧ cd0 ,zd0 n = t ∧ d0 6= d}.

Hierarchical Latent Dirichlet Allocation

hLDA is an nCRP-based topic model to learn a topic
hierarchy [5] from a corpus of D bag-of-words documents
Nd
W = {wd }D
d=1 , where each document wd = {wdn }n=1 has
Nd tokens, and each token is represented by its word ID
wdn ∈ {1, . . . , V } in the vocabulary of V unique words. In
hLDA, topics form an L-level tree, i.e., each tree node t
is a topic, and is associated with a distribution over words
φt ∈ ∆V −1 , where ∆V −1 is the (V −1)-simplex. Since nodes
and topics have one-to-one correspondence, we do not distinguish them in the sequel.
The hLDA model assumes that each document is assignedwith a path cd following the nCRP. The words of each document are modeled by a mixture of the topics in cd , withthe document-specific mixing proportion θ d . The generative
process for the corpus is:

2.4

Collapsed Gibbs Sampling (CGS)

The generative process of hLDA (Sec. 2.2) defines a joint
distributionp(W, z, θ, c, φ) =∞Yp(φt )

t=1• For each node t, draw φt ∼ Dir(βlt 1), where lt is thelevel of node t and 1 = (1, . . . , 1) is an all-one vector;

Nd
D Y
Y
d=1 n=1

DYp(cd |c<d )p(θ d )d=1

p(zdn |θ d )p(wdn |φcd,z),dn• For each document d:– Draw a path of topics cd ∼ p(cd |c<d ) according
to nCRP (Eq. 1);– Draw the topic mixing proportion θ d ∼ Dir(α);– For each position n, choose a node from the path
cd by drawing the level assignment zdn ∼ Mult(θ d ).
Then, draw a word from the chosen topic wdn ∼Mult(φcd,z ),
dn

828

D
where z = {zd }D
d=1 and c = {cd }d=1 . Given the training
corpus W, we want to infer the posterior distribution of
the underlying topics φt as well as the per-document latentvariables, including the paths cd , mixing proportions θ d andlayer assignments zd . We do not have closed-form solutions
for the posterior distribution p(z, θ, c, φ|W) [5]. Instead
we approximate the inference by drawing samples from it,
where Collapsed Gibbs sampling (CGS) [5] is one possible
sampler.

Based on the conjugacy between Dirichlet and multinomial distributions, CGS integrates out θ and φ to get thecollapsed posterior distribution:p(W, z, c) =

zd
wd

D
∞
Y

B(ad + α)  Y B(bt + βlt )p(cd |c<d ),B(α)
B(βlt 1)
t=1

zd
wd

(2)
where Γ(·) is the gamma function,
and B(·)
Q
P is the multivariate beta function, B(α) = k Γ(αk )/Γ( k αk ).
CGS samples from the collapsed posterior distribution in
Eq. (2) by Gibbs sampling [14], i.e., alternatively samplingz and c from their conditional distributions:

2

1

1

2

2

1

3

2

1

2

3

2

1

1

1

1

2

3

3

2

2

2

3

1

1

1

1

1

1

2

2

2

2

1

1

1

2

3

3

2

2

2

3

count

1

1

1

2

3

3

2

2

2

3

0

1

2

0

0

1

0

1

2

0

p(zdn = l|wdn = v, W, z¬dn , c, φI )φcdl ,v
cdl ∈ I,¬dn∝(adl + αl ) b¬dn
cdl ,v +βlt
 s¬dn +V β
cdl ∈ Ī.

fC (d, cl ),

l=1

d
B(b¬d
t + bt + βlt ).
B(b¬dt + βlt )Sample c: sample the path from(3)p(cd = c|W, z, c¬d , φI )(LYfI (d, cl ) cl ∈ I,∝p(c|c¬d )fC (d, cl ) cl ∈ Ī,
l=1

EFFICIENT TRAINING FOR HLDAThe existing algorithms for hLDA (e.g., CGS) are notefficient, which have limited the evaluation at a small-scale
and made hLDA not applicable in practice. To scale up
hLDA, there are two main algorithmic challenges:fI (d, t) =

3.2

t∈I

Time complexity analysis

log fI (d, t) = log

VYd(φtv )btv =v=1

X

log φtv ,(8)v∈Wd,lt

where lt is the level of node t, and Wdl = {wdn |zdn = l} isthe multiset of all tokens in document d that are assigned
to level l, which can be computed by bucket sorting the(zdn , wdn ) pairs.
Similarly, Eq. (3) can be rewritten as:X
log fC (d, t) =log(b¬d(9)
tv + o + βlt ) + ht ,θ,φĪB(ad + α)  Y B(bt + βlt )B(α)
B(βlt )YDir(φt ; βlt )(7)

We now discuss some implementation details of samplingc, and then we compare the time complexity of CGS and
PCGS. To sample c, we need to enumerate all the possible
T paths and compute their nCRP probability p(c|c¬d ) as
described in Sec. 2.1. Then we compute fI (d, t) for all t ∈ I
and fC (d, t) for all t ∈ Ī. Finally, we compute the posterior
probability p(cd = c|W, z, c¬d , φI ) for all the possible T
paths, and sample a path.
Because both fI (d, t) and fC (d, t) are very close to zero,
we compute their logarithms. Rewrite Eq. (6):Partially Collapsed Gibbs Sampling

t∈Ī

V
Y

(6)

βlt + bt.
V βlt + stφt ←We propose a partially collapsed Gibbs sampling (PCGS)
algorithm to avoid computing most of the logarithms, while
preserving the quality of inference.
We pick a subset I of the visited topic IDs {1, . . . , T }.
Unlike CGS which collapses out all the topic distributions
φt , PCGS only collapses out the topics that do not belong toI, while leaving the topics in I instantiated. Let Ī = N+ \I
be the complement of I. We allow matrices to use sets as
subscripts, e.g., φI . In this case, φI is a matrix composedof all the row vectors indexed by t ∈ I. PCGS samples fromthe marginal distributionZp(W, z, θ, c, φ)dθdφĪ
p(W, z, c, φI ) =d=1

d(φtv )btv .

Sample φI : For t ∈ I, draw φt ∼ Dir(βlt +bt ). To speed
up the mixing, we replace the sampling with expectation in
our implementation, i.e.,

To address these challenges, we propose a partially collapsed
Gibbs sampler (PCGS) with an effective initialization strategy.

p(cd |c<d )

V
Y

v=1

• Local optimum. The posterior distribution p(z, c|W)
is multi-modal, so the sampler can be trapped arounda local optimum, resulting in bad model quality.

=(5)

where fC is defined as Eq. (3), and• Efficiency. The main computational burden for CGS,
Eq. (3), needs lots of evaluation of the logarithms (aswe will show in Sec. 3.2), which is expensive.

D
Y


(4)

lt

cdl

fC (d, t) =

3.1

1

3

k)
= Pαkαk , where ek is a coUtilizing the identity B(α+eB(α)
kordinate vector, we can derive the Gibbs sampling updates(see the appendix for derivation details):
Sample z: draw the level assignment for each token from

where

3.

1

2

Figure 3: Illustration of Wd and Wd0b¬dncdl ,v + βl,p(zdn = l|wdn = v, W, z¬dn , c) ∝ (a¬dn+ αl ) ¬dn
dl
scdl + V βl
L
Y

2

1

Bucket sort

d=1

p(cd = c|W, z, c¬d ) ∝ p(c|c¬d )10(v,o)∈Wd,l(φtv )btv .

twhere we convert the logarithm of multivariate beta function as the sum of logarithms (the derivation details can be

v=1

829

Alg.
Both
CGS
PCGS

Step
Sample z
Sample c
Sample c
Sample c
Sample φITime Complexity
O(N L)
O(N T )
 × O(N T )(1 − ) × O(N T )O(T V )

Operation
Division
Logarithm
Logarithm
Addition
Divisiontaneously. More concretely, after a document d is assigned
to a certain path cd , its words are assigned to the levels zdof the current path. In the next iteration, even if there isanother path c0d such that p(wd |c0d ) is larger than p(wd |cd ),p(wd |c0d , zd ) is not likely to be larger than p(wd |cd , zd ) because zd is already optimized for cd . In this case, the pathassignments cd will be quickly trapped in a local optimum
even if there are better path. We also noticed that similar as multinomial mixture models [25], the sampling of cd
is almost deterministic, because p(wd |cd , zd ) can differ byhundreds orders of magnitudes for different cd ’s. Therefore,
it is difficult for a sampler to jump out of the local trap
simply by its randomness. In contrast, flat models such as
LDA do not have this type of local optima because their
topic assignments are directly chosen from all the T topics.
Since hLDA is sensitive to local optima, a proper initialization is crucial to get good results. We obtain a goodinitialization by sampling c from p(c|W) directly instead offrom p(c|z, W) (Eq. 5) for the first I iterations. In otherwords, we integrate z out. In the first I iterations, the sampler focuses on finding the optimal assignment c for eachdocument. Afterwards, the algorithm samples p(c|z, W) torefine the model.
P
The marginal distribution p(c|W) =z p(c|z, W)p(z)has no closed-form representation. WeP
approximate it with
Monte-Carlo integration p(c|W) ≈ S1 zs ∼p(z) p(c|zs , W),
R
where S is the number of samples, and p(z) = p(z|θ)p(θ)dθis a Polya distribution which is approximated with a uniform
discrete distribution over levels.
We also adopt the progressive online initialization strategy [5, 30], which starts with an empty collection of documents, and gradually adds documents by inferring the posterior of document-specific variables (cd , zd ) given all thepreviously observed documents. The documents are organized into mini-batches, and φI is sampled per mini-batch.

Table 1: Time complexity analysis of CGS and PCGS. N :corpus size, typically 106 − 1010 ; T : number of topics, typically 102 − 103 ; L: number of layers, typically 101 ; V : vocabulary size, typically 104 − 106 ;  = 0.05: proportion ofcollapsed topics.

found in the appendix). The term ht = log Γ(s¬dt + V βlt ) −0log Γ(st + V βlt ). The set Wdl
is similar as Wdl , but weassign each token with a count indicating which time does
this word appear, e.g., if a word v is in Wdl for three times,0
we put (v, 0), (v, 1) and (v, 2) into Wdl. See Fig. 3 as anillustration.
Theoretically, both CGS and PCGS can perform correct
inference since they both converge to the true posterior distribution. We also present some intuition on why PCGS
works: b¬d
tv + o + βlt in fC is a more fresh version of φtv
in fI , in the sense that b¬d
tv is updated instantly after eachsampling operation, but φtv is only updated once per iteration. However, if b¬d
tv is large, it will not change much during
the iteration, so PCGS should perform similarly with CGS.
Based on this intuition, we choose I to be the nodes t with
top 95% st , whose corresponding b¬d
tv should be large.
We now analyze the time complexity. By analyzing Eq. (8)and Eq. (9) we know that the time complexity of computing
both fI (d, t) and fC (d, t) for every document d and topic t
P
is at most O(N T ), where N = D
d=1 Nd is the size of the
corpus. Despite having the same time complexity, computing fI is much more efficient than computing fC . Because
log φtv only changes once per iteration according to Eq. (7),it can be pre-processed. Hence, computing fI only involves
N T floating-point additions. In contrast, the logarithm cannot be avoided during computing fC because the count b¬d
tv
is changing all the time. Since PCGS computes fC for only
5% of the topics, it is much more cheaper than CGS which
computes fI for 100% of the topics.
Sampling z requires enumerating the level for every token,
so its time complexity is O(N L), which is much cheaper
than O(N T ) of sampling c because the number of levels L
is much smaller than the number of topics T . Sampling
φI is O(T V ), which is again much smaller than O(N T )because the vocabulary size V , i.e., the number of uniquetokens, is much smaller than the total number of tokens
N . Therefore, we show that computing fC and fI is indeed
the main computational burden, and PCGS is much cheaper
than CGS because floating-point additions are much cheaper
than logarithms. Table 1 summarizes the time complexityof CGS and PCGS.

3.3

4.

SYSTEM IMPLEMENTATION

The tree structured model poses challenges on developing
the distributed system for training hLDA:
1. Efficient vectorization is difficult due to the complicated tree-related data structures.
2. The distributed count matrix has a dynamic numberof columns (topics).
3. The model needs a distributed tree which supports
adding and removing nodes as well as editing the counts.
4. Complicated consistency requirements for the tree andthe matrix. For example, we need to avoid adding
children to a non-existent node of the tree.
Previous flat topic modeling systems [10, 38, 17] that focus on the data / model partition and count matrix synchronization do not solve the aforementioned challenges ofhierarchical topic models. In this section, we propose a
distributed system for training hLDA. Our system follows
the data-parallel paradigm, where the documents are distributed across machines and the model is shared. Our system consists of machine-level and thread-level parallelism, as
shown in Fig. 4. On the machine level, we update the topic
distribution φ in a bulk-synchronous-parallel (BSP) fashion and update the tree structure and the counts (B, m)

Initialization Strategy

The tree structure of hLDA introduces more local optima
than flat topic models. As in Sec. 2.2, the topic assignment
for a particular token wdn is defined by a two-step allocation:
one first chooses a path cd from the infinite tree and thenchooses one topic cd,zdn from the path. If we want to change
the topic assignment, we need to change cd and zdn simul-

830

B

Other Worker Machines

B[1]

Phi[]

B[L]

WorkerWorkerWorkerThreadThreadThread

0

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

16

17

1(1)

B[1] B[2]
2(1)

3(1)

4(2)

5(2)

6(3)

B[3]

0

3

12

6

9

15

1

4

13

7

10

16

2

5

14

8

11

17

Figure 5: Data layout of our system. The number on thenode is “node ID (node rank)”. B is the counting matrix described in our algorithm and B[] is the actual stored format.

B[]
To summarize, the main data structures for hLDA are atree, and L matrices Phi[] and B[]. Phi[l] is a read-onlyand static matrix during the sampling of c since both log φtv
and I are static. However, B[l] is a read+write matrix withdynamically growing number of columns, which is challenging to implement. We will present the implementation ofB[l] in the next section.

Worker Machine
Async. communication
BSP communication
Read / write

Read only

Write only

Figure 4: An overview of the distributed inference system.

4.2
asynchronously; On the thread-level, a number of threads
concurrently read and update the local counts.

4.1

Concurrent Dynamic Matrices

The count matrix B[l] is concurrently read and updated
by the worker threads, and the number of columns (topics)
can grow over time. There are O(N T ) reads and O(N L)
updates per iteration, so the read operation must be efficient, i.e., lock free. On the other hand, the consistency canbe 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 updates to preserve eventual consistency,
while the reads are lock-free operations, to maximize thereading performance.
The dynamic number of columns makes the implementation challenging. The straightforward implementation forgrowing the matrix involves allocating a new memory space,copying the original content to the new space, and deallocating the original space. However, this implementation cannot achieve eventual consistency because the updates duringcopying will not be incorporated.
Inspired by the lock-free design of a concurrent vector [11],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 × 2c+i−1 , while c is a constant.
The first matrix block represents the [0, 2c )-th columns of
the original matrix, the second matrix block represents the[2c , 3 × 2c )-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, c0 ) coordinate, where b is theindex of the matrix block on the list and c0 is the column
index within the matrix block. According to Sec. 3.2, eachmatrix read is always accompanied with an evaluation of
the logarithm, so the overhead of the coordinate conversion
is acceptable. Finally, to improve the locality, we defragment after each PCGS iteration, i.e., deallocating all thematrix blocks and concatenating their content to form a single larger matrix block. We implement a parameter server
to synchronize the updates of the matrices across machines.

Vectorization and Data Layout

We first introduce the design of our data layout for efficient vectorization. As we already analyzed in Sec. 3.2, themost time-consuming part is the computation of fI and fC
for every t according to Eq. (8) and Eq. (9). We can seethat both fI and fC have similar access patterns. For eachdocument d and topic t, they read the tokens v in the multiset Wd,lt , and access the corresponding topic-word entry(t, v) of the topic distribution matrix φ or the count matrix
B. Hence, we only discuss the vectorization of fI . The computation of fI (d, t) for a single topic t cannot be vectorized
because the word IDs in Wd,lt is discontinuous and varies
by document. However, we can compute log fI (d, t) for all
the t’s in one level altogether. For vectorization, we store amatrix Phi[l] for each level l, whose v-th row Phi[l][v] isthe log φtv for all the topics t of level l. In this way, we cancompute all the log fI (d, t) for level l by vector summation:
summing up all the rows Phi[l][v] for v ∈ Wd,l .
Fig. 5 illustrates the data layout. For each node t, weassign a rank rt to it. On level l, the Il nodes in I are
assigned with the rank from 1 to Il , and the other nodes
have ranks from Il + 1. We store log φtv at Phi[lt ][v, rt ] andb¬dnat B[lt ][v, rt ], where lt is the level of topic t, and B[l]tv
has the same shape with Phi[l], but storing the counts.
We maintain a counter on each level to store the largest
rank on that level. If the sampler decides to generate a newnode on level l, the rank counter is increased by one and therank of the new node is set to be the rank counter. Then, wegrow B[l] to have one more column. After each iteration,
we reassign the rank by sorting the nodes on each level by
their topic count st ’s. We put t into I if st is among the top
95%, so that the rank of the nodes in I is still continuous.
We also need to permute the columns of Phi[l] and B[l]according to the new rank. The reassignment of the rank is
performed on a single worker machine and then broadcasted
to all the machines.

831

4.3

Distributed Tree

complexity is quadratic with respect to the number of topics,
which is not scalable to practical big data applications where
we often deal with thousands of topics.

To support the tree-structured model, the system needs
a distributed tree implementation, which supports adding
node, removing node, updating the counts m and querying
the current tree structure. The time complexity of the tree
operations is only O(DL) in total, which is much less thanthe O(N T ) complexity of sampling c. Therefore, we simply
use mutexes for a thread-safe version of the tree.
To share the tree across machines, one possible solution
is a master-client architecture, where the tree is only maintained on the master node, who answers all the queries via
RPC calls. This approach has O(DT ) communication complexity for the master. Moreover, the master cannot handle
requests simultaneously from multiple clients because it uses
mutexes. Therefore, this approach is not scalable.
We optimize this solution by maintaining the counts m ina parameter server. Each machine maintains its own copyof the tree, and the master is only responsible for adding
and removing nodes. The count updates are handled by theparameter server. There are O(L) count updates per document, so the cost of updating counts is O(DL). The numberof operations for adding or removing nodes is roughly O(T ),so the total amount of communication is O(T 2 + DL). Since
T  D and L  T , this complexity is much lower than theoriginal O(DT ). The consistency of the tree structure isguaranteed because the master handles all the node adding
and removing operations.

5.

5.2

• Style of parallelism: there is data parallel (DP) andhybrid data-and-model parallel (HDMP) approaches;• Style of synchronization: there are asynchronous updates by assuming that B changes slowly; and bulksynchronous-parallel (BSP) updates by formulating thealgorithm to be an embarrassingly parallel expectationmaximization (EM) or variational inference (VI) algorithm.
DP+async. Examples include Yahoo! LDA [1] and ParameterServer [17]. They use a parameter server, which is a distributed key-value store, to synchronize the counts. Thesesystems directly approximate the collapsed Gibbs sampling(CGS) algorithm for LDA. However, they assume that thecount matrix B fits in the main memory of a single machine,
thus cannot handle a large number of topics or a large vocabulary.
HDMP+sync. Examples include LightLDA, F+LDA andPeacock [37, 36, 34]. Based on a careful hybrid partitioning strategy, they implement the exact collapsed Gibbs sampling algorithm for LDA and thus have better model quality
than the DP+async. systems. They can also handle large
B by splitting it across machines. However, the throughput
is typically lower than pure DP systems due to the complicated scheduling and the large number of partitions, where
WarpLDA [10] is an exception who carefully organizes thescheduling.
DP+sync. Examples include ESCA [38] and Mr.LDA [39].
They use embarrassingly parallel EM or VI algorithms, where
the EM algorithm is shown to have similar model quality with CGS [38]. Due to the low system overhead, thethroughput can be very high [38], but again they cannot handle large B. This category also includes some recent GPUbased systems such as SaberLDA [16] and BIDMach [40].
Our system for hierarchical topic models is a data parallel one. We borrow the synchronization strategy from bothasynchronous and synchronous systems, where we asynchronously update B[] via a parameter server and synchronously
update Phi[] via an MPI Allreduce operation. Different
from existing flat topic modeling systems, we propose efficient data layouts and distributed data structures for hierarchical topic models. Exploiting better partition strategies
to support large B is one for our future work.

RELATED WORK

We discuss on the related work from both parallel inference for nonparametric Bayesian models and scalable systems for flat topic models.

5.1

Scalable Systems for Flat Models

There are various distributed systems developed for flat
topic models, particularly latent Dirichlet allocation (LDA).
These systems mainly focus on partitioning of the data andthe model, and the synchronization of the topic-word countmatrix B. We categorize existing topic modeling systems in
two dimensions

Parallel Inference for Nonparametric
Bayesian Models

Nonparametric Bayesian models [42] are flexible and they
can grow their structures as more data are observed. hLDA
is a nonparametric Bayesian model that is based on the
nCRP prior. Dirichlet process mixture models (DPMMs) [20]and hierarchical Dirichlet processes (HDPs) [28] are the flat
nonparametric Bayesian counterparts of hLDA which arealso based on CRP. There are many parallel algorithms for
DPMMs and HDPs, but none of them can be readily applied
to hLDAs. We now briefly review them.
Instantiated weight algorithms convert DPMMs and HDPs
to the equivalent stick-breaking process [26] formulation,
where an infinite-dimensional mixing proportion is explicitly instantiated. Variational inference [5, 32] or slice sampling [13] algorithms are then developed for the stick-breaking
process model variant. Instantiated weight algorithms are
efficient and embarrassingly parallel, but we find that theycannot obtain good results for hLDAs due to their susceptibility to local optima.
Split-merge algorithms rely on a split-merge Markov chain
Monte-Carlo transition kernel [15] that is able to jump outof local optima via non-local moves, such as splitting a topic
or merging two topics. Split-merge algorithms are developed
for parallel sampling for DPMMs [9, 8]. However, it is notclear how to split and merge topics on trees.
For hLDA, the only parallel algorithm we are aware ofis Wang and Blei’s instantiated weight variational inference
algorithm with split-merge moves [30]. However, its time

6.

EXPERIMENTS

We now present extensive results to evaluate the quality,
efficiency and scalability of our algorithm and system on
several datasets, including NIPS, 20NewsGroups, NYTimesand PubMed from the UCI machine learning repository [3],a subset of the NYTimes dataset, and two subsets of the
ClueWeb12 dataset [22] (Table 2). Our quantitative and

832

2800

Perplexity

Perplexity

3500
3000
2500

CGS
CGSi
PCGS
PCGSi
hlda-c

2600
2400
2200
2000

2000
(0, 33](33, 66](0, 33](66, 100](33, 66](66, 100]

Number of topics

Number of topics

Figure 6: Comparison of inference quality on sNYTimes (left) and NIPS (right). i denotes for our initialization strategy.
Table 2: Statistics of the datasets.
DatasetsNYTimes
NIPS
20NewsGroups
NYTimes
PubMedsClueWeb12
lClueWeb12

D
3 × 103
1.5 × 103
1.88 × 104
2.93 × 105
8.2 × 106
1.5 × 107
1.31 × 108# tokens
7.23 × 105
1.93 × 106
2.22 × 106
9.7 × 107
7.38 × 108
5.6 × 109
2.8 × 1010

Table 3: Text categorization results on 20NewsGroups.
V
101635
12375
60698
101635141043
100000100000hLDAalgorithm
None
CGS
CGSi
PCGS
PCGSi

qualitative results demonstrate the promise of our system.
The experiments are conducted on the Tianhe-2 supercomputer, which has two 12-core Xeon E5-2692v2 CPUs per
node and an InfiniBand network.
We quantitatively compare the quality of the inferred models by predictive log-likelihood, which indicates how well wemodel the testing corpus, using the document completion
approach [29]. To compute the predictive log-likelihood, thecorpus is divided as a training corpus Wt and a testing corpus, and the testing corpus is further divided as an observed
corpus Wo , which contains a random half of the tokens foreach document in the testing corpus; and a heldout corpus
Wh of the other half of the tokens. The predictive loglikelihood is defined as p(Wh |Wo , φ), where the model φ isinferred from the training corpus Wt , and is approximatedwith a Monte-Carlo integration:
D

S

Ld

L

Y 1 XYX
(s)
(s)
h
p(wh
p(Wh |wo , φ) ≈dn , zdn = l|cd , θ d , φ),S s=1 n=1
d=1

l=1

(s)
(s)
where cd and θ d are the samples from the posterior distri(s)(s)
bution p(cd , θ d |wo , φ), which can be obtained with Gibbs

8168
6500
8228
7382

Testing accuracy+tfidf +tfidf+lda
.818
.827
.828
.830
.834
.836
.812
.825
.828
.835

PCGS with an empty set as I. 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θ instead of a Dirichlet prior [5].
Unlike parametric models, whose numbers of topics are
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, wechoose a rich set of hyper-parameter configurations, run thealgorithms for all these configurations, and plot the perplexity against the number of topics as in Fig. 6. In this experiment, we train a 4-layer model (i.e., L = 4) on the sNYTimes dataset, which is a subset of the NYTimes dataset,and the NIPS dataset, and β = (β0 , 0.5β0 , 0.25β0 , 0.25β0 ),where β0 is chosen from {e−4.0 , e−3.5 , . . . , e2.0 }, γ is chosenfrom {e−6.0 , e−5.5 , . . . , e0.0 }, and α = 0.2 × 1.
By comparing the perplexity produced by different algorithms in Fig. 6, we have a number of observations:• CGS and PCGS have similar quality.
• Our initialization strategy helps obtain better results
for both CGS and PCGS.

sampling, and S is the number of samples. Finally, we convert predictive log-likelihood to predictive perplexity• Our result is not worse (actually better) than hlda-c.
The discrepancy attributes to the different choice ofprior on θ.

perplexity = exp(−log likelihood/number of tokens).
Intuitively, a perplexity of 1,000 means that the model is
as confused as the test data as if it chooses uniformly from
1,000 possibilities per word. The perplexity is between 1 andthe vocabulary size V , and a lower perplexity score indicatesa better model.

6.1

Perplexity

Besides perplexity, we also present accuracy results of thesupervised text categorization task on the 20NewsGroups
dataset, whose each document is labeled as one of the 20
categories. We utilize topic models for feature extraction,T
where the feature
P for document d is the count vector (fdt )t=1 ,where fdt = v bdtv indicates the numbers of occurrences of
topic t in document d. We also incorporate a popular flat
topic model, latent Dirichlet allocation (LDA) for comparison. The result is in Table 3. For each algorithm, we report
the perplexity, and the testing accuracy of a linear support

Quality of Inference

We first compare the model inferred by our implementation of CGS and PCGS and examine the effect of the initialization strategy (Sec. 3.3), where CGS is a special case of833

4.40 x102

CGS-L CGS PCGS
NIPS

105

2.74 xSpeedup

5.63 xTime (s)Time (s)7
3.58 x4.09 x102

CGS-L CGS PCGSNYTimesC

CGS-L CGS PCGS
PubMed
Z

phi

other

Time (s)Time (s)tures. The first feature set “+tfidf” contains standard tf-idf
word features [24] and hLDA topic features; and the second feature set “+tfidf+lda” further has LDA features. We
report the accuracy for different hLDA training algorithms,
and the baseline “None” does not utilize any hLDA features.
CGSi +tfidf features give 0.834 accuracy, outperforming
the 0.827 accuracy by the None+tfidf+LDA baseline. Utilizing the features produced by both hLDA and LDA, the
CGSi +tfidf+LDA approach further improves the accuracy
to 0.836. The comparison between different training algorithms for hLDA is similar with our previous perplexity
study of Fig. 6. PCGSi +tfidf+LDA gives similar 0.835 accuracy with CGSi +tfidf+LDA, and the algorithms without
our initialization strategies, CGS and PCGS, are considerably worse in terms of both perplexity and testing accuracy.

6.2

300

1000

(c) Number of topics
C
Z

1
0.8
0.6
0.4
0.2
0
CGS
PCGS
(b) Time decomposition
1
0.8
0.6
0.4
0.2
0
CGS
PCGS
(d) Time decomposition
phi
other

Figure 8: Speedup of PCGS over CGS w.r.t. the number of
documents. The three
bars
3.58
x in are 3000, 30,000, and 300,000
2.74 xin (d).
documents in (b) ; and 100, 300, and 1,000 topics

4

CGS
PCGS
PubMedvector machine (liblinear [12])with two different 10
set0 of fea-

5

total

10
Figure7: Running time comparison.
The numbers
3.58 x104 on the2.74 xplot are the speedups over the previous implementation, e.g.,2
the10speedupof CGS over CGS-L. CGS-L is 4.52 times fasterthan hlda-c on sNYTimes.
102
CGS-L

6

4
100

100

100

30000 300000(a) Number of documents

104

total

3
2
3000

100
CGS-L CGS PCGSsNYTimes

4

Proportion of execution time100

7.68 xSpeedup2.76 xTime (s)Time (s)6.97 x1025Proportion of execution time104

4.52 xsampling operations, such as maintaining the tree and thesystem overhead, are denoted as “other”. PCGS gives 2.745.63 times speedup over CGS, by avoiding the computing of
CGS-L
CGS
PCGS
logarithms and cheaper
read-only data structures.
The datalayout and vectorization
proposed
in
Sec.
4.1
give
3.58-7.68
PubMedtimes speedup. Furthermore, CGS-L is 4.52 times fasterthan hlda-c due to the implementation, mostly because wecast difference of logarithm of gamma function as sums oflogarithms as Eq. 3. Overall, our PCGS implementation is
87 times faster than hlda-c on the sNYTimes dataset.

6.3

Impact of the Problem SizeWe further examine the impact factors to the speedup of
PCGS over CGS. PCGS mainly reduces the cost of samplingc. By Table 1, the time complexity of sampling c is proportional to both the corpus size N and the number of topics
T . Therefore, we expect the speedup of PCGS over CGS
is large when the cost of sampling c dominates, i.e., when
N and T are large. We report the speedup results on the
NYTimes dataset, varying N and T , while keeping all theother factors fixed. Fig. 8(a) shows the speedup of PCGS
over CGS increased from 2.1 to 4.6 times when the numberof documents increases from 3,000 to 300,000. The increaseof speedup is due to the increased proportion of time consumption for sampling c, as shown in Fig. 8(b). Similarly,
Fig. 8(c) shows that the speedup of PCGS over CGS increased from 4.1 to 6.1 times when the number of topics
increases from 100 to 1,000; and the time proportion of eachstep is shown in Fig. 8(d).

Efficiency

We now study the impact of our PCGS algorithm and thedata layout to the efficiency of system. We compare the running time of three implementations, PCGS, CGS, and CGS
without the proposed data layout and vectorization (“CGSL”) on sNYTimes, NIPS, NYTimes and PubMed datasets.
The number of topics is kept around 300 by tuning thehyper-parameters. To keep the running time reasonable, we
utilize a single thread on sNYTimes and NIPS, 12 threads
on NYTimes, and ten 12-core CPUs on PubMed. On sNYTimes, we also report the running time of hlda-c, whose
running time is longer than 2 hours on all the other datasets.
The running time is shown in Fig. 7. We report the total
running time as well as the running time for each sampling
step, i.e., sampling C, Z and φ as shown in Table 1. Non-

6.4

Impact of the Machine SizeWe also study whether of the number of threads and CPUs
affects the speedup of PCGS over CGS. We run PCGS andCGS from 1 to 12 threads on the NYTimes corpus, and834

PCGS vs CGS
PCGS speedup
CGS speedup

10

SpeedupSpeedup

8
6
4

PCGS vs CGS
PCGS speedup
CGS speedup

6
4

130

2.2

Perplexity

8

Number of topics

12

120
110

1.8
2

2
3 4 5 6 7 8 9 10(b) Number of machines

Perplexity

6.6

2200
2000

2

8

32 128S

Scalability

Our experiments on scalability are in two folds: whether
the quality of inference is affected by parallelization; andhow good is the speedup. We repeat each experiment for 10
times. For each experiment, we choose the largest dataset
such that the data fit in the main memory and the running time is less than two hours. Based on this criteria, wechoose NYTimes for the multi-thread setting, PubMed for1-10 machines, and sClueWeb12 for 10-100 machines, which
is a subset of the ClueWeb12 corpus.
For the multi-thread setting, where the number of threads
varies from 1 to 12 on the NYTimes corpus. The result isshown in Fig. 12(a), where we observe that there is no apparent increase of the 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 permachine. We run our implementation on the PubMed corpus on 1 to 10 CPUs as shown in Fig. 12(b), and on thelarger sClueWeb12 corpus for 10 to 100 CPUs as shown in
Fig. 12(c). The speedup is 8.5 from 1 to 10 CPUs, and 7.15
from 10 to 100 CPUs. The perplexity is slightly affected byparallelization 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 modelwith 1,722 topics of the 131-million-document lClueWeb12
corpus with 50 machines, and the inference finishes in 7
hours. The results will be presented for qualitative evaluation in the next section.

I=1
I=2
I=4
I=8
I=16
I=32
I=64

2400

(50, 100] (100, 150]

Number of topics

Figure 10: Impact of I to perplexity.

from 3 to 10 CPUs on the PubMed corpus, because CGS
exceeds the main memory capacity on 1-2 machines. CGS
needs to access the concurrent dynamic matrix more often
than PCGS. The scalability of the concurrent dynamic matrix is worse than a read-only matrix which PCGS mostly
accesses, due to the contention of the atomic update operations. Therefore, the speedup of PCGS over CGS increases
as the number of threads grows, as shown in Fig. 9(a).
Meanwhile, the updates of the distributed matrices and trees
are done asynchronously, so the speedup of PCGS over CGS
remains relatively stable with respect to the number of machines, as shown in Fig. 9(b).

6.5

32 128Figure 11: Impact of the number of Monte-Carlo samples
S.

Figure 9: Speedup of PCGS over CGS w.r.t. the numberof threads and machines. The two lines are the speedup of
PCGS / CGS as the number of threads / machines grows,
and the bar is the speedup of PCGS over CGS.

(0, 50]8

S

0
1 2 3 4 5 6 7 8 9 101112(a) Number of threads

2

100

2
0

#103

Sensitivity of Parameters

6.7

We now examine the impact of the hyper-parameters I
and S of our initialization strategy.
Impact of I:
I is the number of initializing iterations
to sample from p(c|w). We run PCGS on the sNYTimes
dataset, setting β0 = 1, and varying I and γ. It can beseen from Fig. 10 that the perplexity steadily decreases forlarge I, which again shows that our initialization strategy is
helpful. We select a moderate I = 32 for all the experiments.
Impact of S: S is the number of Monte-Carlo samples toapproximate p(c|w). When S → ∞, we directly sample
from p(c|w) in the first I iterations. We run PCGS onsNYTimes, with β = (1.0, 0.5, 0.25, 0.25) and γ = 10−40 ,and vary S from 1 to 128. As shown in Fig. 11, S has little
impact on both the number of topics and the perplexity,
implying that a small S, e.g., S = 5, is adequate.

Qualitative Analysis

We now demonstrate the topic hierarchy obtained fromthe lClueWeb12 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 are between [50, 500]. We show theselected parts of the obtained tree in Fig. 13, 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.1
Fig. 13(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
1
The visualization demo is available online at http://ml.
cs.tsinghua.edu.cn/~jianfei/scalable-hlda.html.

835

#103
10

4

SpeedUp

Perplexity

4.1

3.9

The complicated data structures of hierarchical topic models pose challenges for building scalable training systems.
To solve these challenges, we design a distributed training
system with an efficient data layout for vectorization, concurrent dynamic matrices and distributed trees. The proposed system can handle hundreds of millions of documents,
thousands of topics, with thousands of CPU cores.
In the future, we plan to extend our method to the moresophisticated nested hierarchical Dirichlet process model [21,2]. Developing sampling algorithms with a sub-linear timecomplexity with respect to the number of topics is also a
promising direction.

Ideal

5

3.8
3.7

0
1

12

1

Number of threads

12

Number of threads(a) NYTimes corpus
3.4

#103

10

SpeedUp

Perplexity

Ideal

3.2

8.

5

3
0
1

10

1

Number of CPUs

ACKNOWLEDGMENTS

We thank Arnab Bhadury for his help in proofreading
the paper. The work was supported by the National NSF of
China (Nos. 61620106010, 61621136008, 61332007, 61761136020),
Beijing Natural Science Foundation (No. L172037), and aSpecial Program for Applied Research on Super Computation of the NSFC-Guangdong Joint Fund (the 2nd phase).

10

Number of CPUs(b) PubMed corpus

Appendix#10310

Derivation of PCGS updates: Rewrite the joint distribution in Sec. 2.4 as:YYp(W, z, θ, c, φĪ , φI ) =p(φt )p(φt )

Ideal

SpeedUp

Perplexity

3.6
3.4
3.2

5

3
0
10

100

Number of CPUs

10

100

p(cd |c<d )p(θ d )Number of CPUs

Nd
D Y
Y

p(zdn |θ)p(wdn |zdn , φĪ , φI ).

d=1 n=1d=1Integrating out φC and θ, we have the marginal distribution(c) sClueWeb12 corpus
Figure 12: Perplexity and speedup as the amount of computational resource increases.

p(W, z, c, φI ) =DYB(ad + α)  Y B(bt + βlt )p(cd |c<d )B(α)
B(βlt )

d=1

Yvariety 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. 13(b)the “travel” topic is divided as “islands”, “India” and “vacations”, and the leaf level contains specific concepts, suchas “ferry” and “diving”, which are correctly placed underthe “islands” topic. In Fig. 13(c), the “computer” topic is
divided as “website”, “windows”, “vps”, “programming”,“linux” and “forum”. To our knowledge, this is the firsttime that hLDA is applied to large-scale web data, and theresults demonstrate our ability on automatically learning
topic hierarchy from web data.

7.

t∈I

t∈Ī
D
Y

t∈I

t∈Ī

Dir(φt ; βlt )VY(φtv )btv .

v=1

k)Utilizing the identity B(α+e= Pαkαk , where ek is a coB(α)
kordinate vector, we can derive the Gibbs sampling updates:
Sample z: Keeping only the terms relevant with zdn , weQ
d +α)
have p(zdn = l|wdn = v, w, z¬dn , c, φI ) ∝ B(C
t∈C
B(α)
QV
B(Ct +βlt ) Q
B(C¬dn
+Cdn
Ctv
d
d +α)
∝
t∈I
v=1 φtv
B(βlt )
B(C¬dn
+α)
d
dn
¬dn
dn
Q
Q
Q
B(C¬dn
+C
+β
)
Ctv
+Ctv
V
l
t
t
t
.
t∈C
t∈I
v=1 φtv
B(C¬dn+βlt )t
Q
Sample c: p(cd = c|w, z, c¬d , φI ) ∝ p(cd |c¬d ) t∈C
¬d
d
QV
Q
B(Ct +Ct +βlt )B(Ct +βlt ) Q
Ctv
t∈I
v=1 φtv ∝ p(cd |c¬d )t∈CB(βlt )B(C¬dt +βlt )(d
Q
QVQ
fI (d, cl ) cl ∈ I,
Ctv
L
t∈I
v=1 φtv ∝ nCRP(c; γ, c¬d )l=1
fC (d, cl ) cl ∈ C.
Sample φI : For t ∈ I, draw φt ∼ Dir(βlt + bt ).
Derivation of log fC (d, t):

CONCLUSIONS AND DISCUSSIONS

We address the problem of scalable training of hierarchical topic models, taking hierarchical latent Dirichlet allocation (hLDA) as a concrete example. We present a partiallycollapsed Gibbs sampling algorithm for the hierarchical latent Dirichlet allocation model. PCGS is efficient because
it avoids most computation of logarithms. To deal with thelocal minima introduced by the tree structure, we present
an initialization strategy that improves the model quality.

B(b¬d +bd +βlttt
We have log fC (d, t) = log B(b¬d +β )
ltt
 PV Pbdtv −1=log(b¬dtv + i + βlt ) + ht
i=0
v=1
P
= (v,o)∈W 0 log(b¬d
tv + o + βlt ) + ht .
dl

836

)

news

email
free
onlinebusiness
hotelscrackslomalka
pro

hotel

area
citiesphone
park

travelcrack

city
airportserials

order
free

city

price

student

series

php

power

code

model documentation programs
dataversioneducationsystemssourceschoolprice

shipping
cart

girl

accessories
steel
products

staff

girls
black
lady
fashionshipping

price notes
local west organic prescription john pages flowers energy
car power file board cracks college parts macromedia electric product martinbuild research islands france india holiday
price
shirtsgas products jean
onlinefishing talkprologin optimal campus charger
auto products files
serverplant electron
country ans kerala apartments shipping san north naturalforesttumblrphase
race special flower
loswoodpaul accessories tags
food
kit
italy
image treesearch comment skyrock goa
east
honda
cable server password lomalka graduatebuy
pricescrack
china
travelmagnetic mercedes
england wiki
suite
seed
serial international republiccarangeles south proteinmhz
generic
vacation kitchen
cracksinch
pro academic
data member
archivedelhi
aismachine michael
structuresize
william log
table
pacific
seeds
river
web
panasonic
fat
bmw
supply error
arab
forum
pharmacy
real
hotel
perfect
david
love
price
wall
bag
coursesshirt(a) Top 3 levels. The highlighted nodes are expanded as below.

hotelshotel

travelcity
airportislands
countrykeralacomment

loveweddings
wedding

costa
rica
surf
rican
estate

spa
spas

ferry

ferries

thailandpattayaaustralia
japan
phuketwritingsrichinmoyquotespoetry diving
poems
poem
poets

greecepoet

italy
cosmicsanbrazil
toll

tour

scuba

united
apr

padi

diver

japan

dictionary

flightticketjaipur

lacs
plots

cusco
spain

taj

mahal

kathmandu
bhutanflights
airlines

tourism
music

airvideosconnectingindia

airways
airlinesaiholidaysrentlondonparis
rooms

trainsagra

flights

hotel
hotelstickets
cheaprajasthan
tours

hotel japanese travel gurgaon tour
nepal
booking
tours trekking
quikr
peru
kanji
feb
noida
tokyo
tours
wildlife
trek

dive

italy
vacation
hotel

goa

traveldelhi

republicarab

belizeislandladiesholiday
apartments

india

train

stationavailability
china
railcard

italy

breakfast
bed
stars
rome

holiday
airlinesinn

airport vacation
rentalsmaui
beachparkingairlinetickets
flights

hawaiiwarsaw
hotelessuitestravelape

hoteleaccommodations
chicagopara
polandvillas
villa
bali
ski

thailand

car

hire

rentalairportcars(b) Travel subtree

file

files

server
data
error
php

driverwindowsnet

mailing
mirrordrivers

download
bios

list

java
file
apache extension ncurses
includesextensions deprecatedorg
line
catalina
core

stats
rand

ereg

includephp

functionline
deprecatedmaxdb
mysqli
ibase
oci
odbc

downloaddell

unmeteredmbit
servers

hosting

packageinstall
rpm
perl

int

trac

years
svn

changeset
edgewalldatemessagetrac
threadlinux

stringpublicvoid

oracle
fader
quick
error
pro
dedicatedtns
automation invalid
vps
recordserversql

notebook windows
driverdrivers

class

vps

windows

virtualphp
methodinlinedefinition

object
java

const
int

array

dll

exe

windows

java
apacheproject

download
mysql
errorsrepository

wiki

text freebsd
aaestring linux thread
backup fpml
sharpdevelop
int
date
tgz
netbsdbytes

comparevoid
object(c) Computer subtree
Figure 13: Selected subtrees on the topic hierarchy extracted from ClueWeb12 (large).

837

git
rfc

png

9.

REFERENCESIn ICML, 2006.
[19] K. P. Murphy. Machine learning: a probabilistic
perspective. MIT press, 2012.
[20] R. M. Neal. Markov chain sampling methods fordirichlet process mixture models. Journal ofcomputational and graphical statistics, 9(2):249–265,
2000.
[21] J. Paisley, C. Wang, D. M. Blei, and M. I. Jordan.
Nested hierarchical dirichlet processes. IEEE
Transactions on Pattern Analysis and MachineIntelligence, 37(2):256–270, 2015.
[22] T. L. Project. The clueweb12 dataset.
http://lemurproject.org/clueweb12/, 2013.
[23] J. Pujara and P. Skomoroch. Large-scale hierarchical
topic models. In NIPS Workshop on Big Learning,
2012.
[24] J. Ramos et al. Using tf-idf to determine word
relevance in document queries. In Proceedings of thefirst instructional conference on machine learning,volume 242, pages 133–142, 2003.
[25] L. Rigouste, O. Cappé, and F. Yvon. Inference andevaluation of the multinomial mixture model for text
clustering. Information processing & management,
43(5):1260–1280, 2007.
[26] J. Sethuraman. A constructive definition of dirichletpriors. Statistica sinica, pages 639–650, 1994.
[27] Y. W. Teh. Dirichlet process. In Encyclopedia ofmachine learning, pages 280–287. Springer, 2011.
[28] Y. W. Teh, M. I. Jordan, M. J. Beal, and D. M. Blei.
Sharing clusters among related groups: Hierarchicaldirichlet processes. In NIPS, 2004.
[29] H. M. Wallach, I. Murray, R. Salakhutdinov, and
D. Mimno. Evaluation methods for topic models. In
ICML, 2009.
[30] C. Wang and D. M. Blei. Variational inference for thenested chinese restaurant process. In NIPS, 2009.
[31] C. Wang and D. M. Blei. Collaborative topic modelingfor recommending scientific articles. In Proceedings ofthe 17th ACM SIGKDD international conference on
Knowledge discovery and data mining, pages 448–456.
ACM, 2011.
[32] C. Wang and D. M. Blei. Truncation-free online
variational inference for bayesian nonparametricmodels. In Advances in neural information processing
systems, pages 413–421, 2012.
[33] X. Wang, S. Liu, J. Liu, J. Chen, J. J. H. Zhu, and
B. Guo. Topicpanorama: A full picture of relevant
topics. TVCG, 22(12):2508–2521, 2016.
[34] Y. Wang, X. Zhao, Z. Sun, H. Yan, L. Wang, Z. Jin,
L. Wang, Y. Gao, J. Zeng, Q. Yang, et al. Towards
topic modeling for big data. ACM Transactions onIntelligent Systems and Technology, 9(4), 2014.
[35] X. Wei and W. B. Croft. Lda-based document models
for ad-hoc retrieval. In SIGIR, 2006.
[36] H.-F. Yu, C.-J. Hsieh, H. Yun, S. Vishwanathan, andI. S. Dhillon. A scalable asynchronous distributed
algorithm for topic modeling. In Proceedings of the24th International Conference on World Wide Web,
pages 1340–1350. International World Wide Web
Conferences Steering Committee, 2015.
[37] J. Yuan, F. Gao, Q. Ho, W. Dai, J. Wei, X. Zheng,[1] A. Ahmed, M. Aly, J. Gonzalez, S. Narayanamurthy,
and A. Smola. Scalable inference in latent variable
models. In WSDM, 2012.
[2] A. Ahmed, L. Hong, and A. J. Smola. Nested chinese
restaurant franchise process: Applications to user
tracking and document modeling. In ICML, 2013.
[3] A. Asuncion and D. Newman. Uci machine learning
repository, 2007.
[4] D. Blei. hlda-c. http://www.cs.columbia.edu/~blei/
downloads/hlda-c.tgz, 2009.
[5] D. M. Blei, T. L. Griffiths, and M. I. Jordan. The
nested chinese restaurant process and bayesiannonparametric inference of topic hierarchies. Journalof the ACM (JACM), 57(2):7, 2010.
[6] D. M. Blei, A. Y. Ng, and M. I. Jordan. Latent
dirichlet allocation. JMLR, 3:993–1022, 2003.
[7] J. L. Boyd-Graber, D. M. Blei, and X. Zhu. A topic
model for word sense disambiguation. In
EMNLP-CoNLL, 2007.
[8] T. Campbell, J. Straub, J. W. Fisher III, and J. P.
How. Streaming, distributed variational inference forbayesian nonparametrics. In NIPS, 2015.
[9] J. Chang and J. W. Fisher III. Parallel sampling of dp
mixture models using sub-cluster splits. In Advancesin Neural Information Processing Systems, pages
620–628, 2013.
[10] J. Chen, K. Li, J. Zhu, and W. Chen. Warplda: acache efficient o(1) algorithm for latent dirichlet
allocation. PVLDB, 9(10):744–755, 2016.
[11] D. Dechev, P. Pirkelbauer, and B. Stroustrup.
Lock-free dynamically resizable arrays. In
International Conference On Principles Of Distributed
Systems, 2006.
[12] R.-E. Fan, K.-W. Chang, C.-J. Hsieh, X.-R. Wang,
and C.-J. Lin. Liblinear: A library for large linear
classification. Journal of machine learning research,9(Aug):1871–1874, 2008.
[13] H. Ge, Y. Chen, E. CAM, M. Wan, andZ. Ghahramani. Distributed inference for dirichlet
process mixture models. In ICML, 2015.
[14] S. Geman and D. Geman. Stochastic relaxation, gibbs
distributions, and the bayesian restoration of images.
IEEE Transactions on pattern analysis and machine
intelligence, (6):721–741, 1984.
[15] S. Jain and R. M. Neal. A split-merge markov chainmonte carlo procedure for the dirichlet process
mixture model. Journal of Computational andGraphical Statistics, 13(1):158–182, 2004.
[16] K. Li, J. Chen, W. Chen, and J. Zhu. Saberlda:
Sparsity-aware learning of topic models on gpus. In
Proceedings of the Twenty-Second International
Conference on Architectural Support for Programming
Languages and Operating Systems, pages 497–509.
ACM, 2017.
[17] M. Li, D. G. Andersen, J. W. Park, A. J. Smola,
A. Ahmed, V. Josifovski, J. Long, E. J. Shekita, and
B.-Y. Su. Scaling distributed machine learning withthe parameter server. In OSDI, volume 1, page 3,
2014.
[18] W. Li and A. McCallum. Pachinko allocation:
Dag-structured mixture models of topic correlations.

838

E. P. Xing, T.-Y. Liu, and W.-Y. Ma. Lightlda: Big
topic models on modest computer clusters. In
Proceedings of the 24th International Conference on
World Wide Web, pages 1351–1361. International
World Wide Web Conferences Steering Committee,
2015.
[38] M. Zaheer, M. Wick, J.-B. Tristan, A. Smola, and
G. L. Steele Jr. Exponential stochastic cellular
automata for massively parallel inference. In
AISTATS, 2016.
[39] K. Zhai, J. Boyd-Graber, N. Asadi, and M. L.
Alkhouja. Mr. lda: A flexible large scale topic
modeling package using variational inference inmapreduce. In WWW, 2012.
[40] H. Zhao, B. Jiang, J. F. Canny, and B. Jaros. Samebut different: Fast and high quality gibbs parameter
estimation. In Proceedings of the 21th ACM SIGKDDInternational Conference on Knowledge Discovery andData Mining, pages 1495–1502. ACM, 2015.
[41] J. Zhu, A. Ahmed, and E. P. Xing. MedLDA:
Maximum margin supervised topic models. JMLR,
13:2237–2278, 2012.
[42] J. Zhu, J. Chen, W. Hu, and B. Zhang. Big learningwith bayesian methods. National Science Review, pagenwx044, 2017.

839