Analyses of Multi-collection Corpora via Compound Topic Modeling

TL;DR

This paper introduces the compound latent Dirichlet allocation (cLDA) model for analyzing multiple text collections simultaneously, capturing shared topics and their variations across subcollections, with improved inference methods and demonstrated effectiveness.

Contribution

The paper proposes the cLDA model that extends traditional topic models to multiple collections, incorporating prior knowledge and offering efficient inference techniques.

Findings

cLDA effectively captures shared and varying topics across collections

The model outperforms existing methods in qualitative and quantitative evaluations

Efficient MCMC method improves parameter estimation

Abstract

As electronically stored data grow in daily life, obtaining novel and relevant information becomes challenging in text mining. Thus people have sought statistical methods based on term frequency, matrix algebra, or topic modeling for text mining. Popular topic models have centered on one single text collection, which is deficient for comparative text analyses. We consider a setting where one can partition the corpus into subcollections. Each subcollection shares a common set of topics, but there exists relative variation in topic proportions among collections. Including any prior knowledge about the corpus (e.g. organization structure), we propose the compound latent Dirichlet allocation (cLDA) model, improving on previous work, encouraging generalizability, and depending less on user-input parameters. To identify the parameters of interest in cLDA, we study Markov chain Monte Carlo (MCMC) and variational inference approaches extensively, and suggest an efficient MCMC method. We evaluate cLDA qualitatively and quantitatively using both synthetic and real-world corpora. The usability study on some real-world corpora illustrates the superiority of cLDA to explore the underlying topics automatically but also model their connections and variations across multiple collections.

Figures

Equations

$${\mathbf{\beta }}_k$$

$${\mathbf{\pi }}_j$$

$${\mathbf{\theta }}_{jd}$$

$$z_{jdi}$$

$$w_{jdi}$$

$$p_{h,\mathbf{w}}(\mathbf{\psi })\propto {\ell }_{\mathbf{w}}(\mathbf{\psi })p_h(\mathbf{\psi })$$

$$\begin{split}p_{h,{\bm{w}}}(\bm{\psi})&\propto\Biggl{[}\prod_{j=1}^{J}\prod_{d=1}^{D_{j}}\frac{\prod_{k=1}^{K}\theta_{jdk}^{n_{jdk}+\gamma\pi_{jk}-1}}{\prod_{k=1}^{K}{\Gamma(\gamma\pi_{jk})}}\Biggr{]}\Biggl{[}\prod_{j=1}^{J}\prod_{k=1}^{K}\pi_{jk}^{\alpha-1}\Biggr{]}\\ &\hskip 14.22636pt\Biggl{[}\prod_{k=1}^{K}\prod_{v=1}^{V}\beta_{kv}^{\sum_{j=1}^{J}\sum_{d=1}^{D_{j}}m_{jdkv}+\eta-1}\Biggr{]}\end{split}$$

$$\begin{split}\bm{\theta}_{jd}&\sim\operatorname{Dir}_{K}\biggl{(}n_{jd1}+\gamma\pi_{j1},\ldots,n_{jdK}+\gamma\pi_{jK}\biggr{)},\\ \bm{\beta}_{k}&\sim\operatorname{Dir}_{V}\biggl{(}m_{k1}+\eta,\ldots,m_{kV}+\eta\biggr{)},\end{split}$$

$$\displaystyle\begin{split}p_{h,{\bm{w}}}(\bm{\pi},{\bm{z}})&\propto\Biggl{[}\prod_{j=1}^{J}\prod_{d=1}^{D_{j}}\prod_{k=1}^{K}\frac{\Gamma(\gamma\pi_{jk}+n_{jdk})}{\prod_{k=1}^{K}{\Gamma(\gamma\pi_{jk})}}\Biggr{]}\\ &\hskip 13.08827pt\Biggl{[}\prod_{j=1}^{J}\prod_{k=1}^{K}\pi_{jk}^{\alpha-1}\Biggr{]}\Biggl{[}\prod_{k=1}^{K}\frac{\prod_{v=1}^{V}\Gamma(m_{..kv}+\eta)}{\Gamma(m_{..k.}+V\eta)}\Biggr{]}\end{split}$$

$$p_{h,\mathbf{w}}\left(z_{jdik}=1|.\right)\propto \frac{\gamma {\pi }_{jk}+n_{jdk}^{(-jdi)}}{\gamma +n_{jd}^{(-jdi)}}\frac{\eta +m_{kv}^{(-jdi)}}{V\eta +m_k^{(-jdi)}}$$

$${\tilde{p}}_{{\mathbf{w}}_j}({\mathbf{\pi }}_j|{\mathbf{z}}_j)\propto \prod _{d=1}^{D_j}\prod _{k=1}^K\frac{\Gamma (\gamma {\pi }_{jk}+n_{jdk})}{\Gamma (\gamma {\pi }_{jk})}\prod _{k=1}^K{\pi }_{jk}^{\alpha -1}$$

$$\frac{\Gamma (\gamma {\pi }_{jk}+n_{jdk})}{\Gamma (\gamma {\pi }_{jk})}=\sum _{s=0}^{n_{jdk}}\text{S}(n_{jdk},s){(\gamma {\pi }_{jk})}^s,$$

$${\tilde{p}}_{{\mathbf{w}}_j}({\mathbf{\pi }}_j|.)\propto \prod _{d=1}^{D_j}\prod _{k=1}^K\text{S}(n_{jdk},s_{jdk}){(\gamma {\pi }_{jk})}^{s_{jdk}}\prod _{k=1}^K{\pi }_{jk}^{\alpha -1}$$

$${\mathbf{\pi }}_j\sim {\mathrm{Dir}}_K\left(\sum _{d=1}^{D_j}{s}_{jd1}+\alpha ,\dots ,\sum _{d=1}^{D_j}{s}_{jdK}+\alpha \right),$$

$$s_{jdk}=\sum _{l=1}^{n_{jdk}}{s}^{(l)},s^{(l)}\stackrel{\text{iid}}{\sim }\text{Bernoulli}\left(\frac{\gamma {\pi }_{jk}}{\gamma {\pi }_{jk}+l-1}\right)$$

$$a({\mathbf{\pi }}^{(t)},{\mathbf{\pi }}^{(\ast )})=\text{min}\left(1,\frac{\tilde{p}({\mathbf{\pi }}^{(\ast )})q({\mathbf{\pi }}^{(t)}\leftarrow {\mathbf{\pi }}^{(\ast )})}{\tilde{p}({\mathbf{\pi }}^{(t)})q({\mathbf{\pi }}^{(\ast )}\leftarrow {\mathbf{\pi }}^{(t)})}\right)$$

$$\text{d}\mathbf{\pi }(t)=\frac{1}{2}{\nabla }_{\mathbf{\pi }}\mathcal{L}({\mathbf{\pi }}^{(t)})\text{d}t+\text{d}\mathbf{b}(t)$$

$${\mathbf{\pi }}^{(\ast )}$$

$$\mathbf{\mu }({\mathbf{\pi }}^{(t)},\epsilon )$$

$${\mathbf{\xi }}^{(t)}$$

$$\begin{split}|\varphi_{k}|&\stackrel{{\scriptstyle\text{iid}}}{{\sim}}\text{Gamma}\Bigl{(}\alpha,1\Bigr{)},\\ p({\bm{\varphi}})&\propto\prod_{k=1}^{K}{|\varphi_{k}|}^{\alpha-1}e^{-|\varphi_{k}|}.\end{split}$$

$$\eta_{(t+1)}=\frac{\eta_{t}}{V}\frac{\sum_{s,k,v}\Big{[}\Psi(m^{(s)}_{kv}+\eta_{t})-\Psi(\eta_{t})\Big{]}}{\sum_{s,k}\Big{[}\Psi(m^{(s)}_{k}+V\eta_{t})-\Psi(V\eta_{t})\Big{]}}$$

$$\gamma_{(t+1)}=\gamma_{t}\frac{\sum_{s,j,d,k}\pi^{(s)}_{jk}\Big{[}\Psi(n^{(s)}_{jdk}+\gamma_{t}\pi^{(s)}_{jk})-\Psi(\gamma_{t}\pi^{(s)}_{jk})\Big{]}}{\sum_{s,j,d}\Big{[}\Psi(n^{(s)}_{jd}+\gamma_{t})-\Psi(\gamma_{t})\Big{]}}$$

$$\mathcal{S}({\mathbf{w}}^{\text{test}}|{\mathbf{w}}^{\text{train}})=\exp \frac{-{\sum }_{w_{jdi}\in {\mathbf{w}}^{\text{test}}}\log p(w_{jdi}|{\mathbf{w}}^{\text{train}})}{|{\mathbf{w}}^{\text{test}}|}$$

$${\text{topic-coherence}}_k=\sum _{i=2}^m\sum _{j=1}^i\log \frac{\text{df}(v_i^{(k)},v_j^{(k)})+1}{\text{df}(v_j^{(k)})}$$

$$p_h(\mathbf{z}|\mathbf{\theta },\mathbf{\pi },\mathbf{\beta })p_h(\mathbf{\theta }|\mathbf{\pi })p_h(\mathbf{\pi })p_h(\mathbf{\beta }),$$

$$p_h(\mathbf{\psi })$$

$$\begin{array}{rl}p(\mathbf{w}|\mathbf{z},\mathbf{\theta },\mathbf{\pi },\mathbf{\beta })& =\prod _{j=1}^J\prod _{d=1}^{D_j}\prod _{i=1}^{n_{jd}}\prod _{k:z_{jdik}=1}\prod _{v=1}^V{\beta }_{kv}^{w_{jdiv}}\\ & =\prod _{j=1}^J\prod _{d=1}^{D_j}\prod _{k=1}^K\prod _{v=1}^V\prod _{i\in S_{jdk}}{\beta }_{kv}^{w_{jdiv}}\\ & =\prod _{j=1}^J\prod _{d=1}^{D_j}\prod _{k=1}^K\prod _{v=1}^V{\beta }_{kv}^{{\sum }_{i\in S_{jdk}}{w}_{jdiv}}\\ & =\prod _{j=1}^J\prod _{d=1}^{D_j}\prod _{k=1}^K\prod _{v=1}^V{\beta }_{kv}^{m_{jdkv}},\end{array}$$

$$\displaystyle\begin{split}p_{h,{\bm{w}}}(\bm{\psi})&\propto\Biggl{[}\prod_{k=1}^{K}\prod_{v=1}^{V}\beta_{kv}^{\sum_{j=1}^{J}\sum_{d=1}^{D_{j}}m_{jdkv}+\eta-1}\Biggr{]}\\ &\hskip 13.37277pt\Biggl{[}\prod_{j=1}^{J}\prod_{d=1}^{D_{j}}\frac{\prod_{k=1}^{K}\theta_{jdk}^{n_{jdk}+\gamma\pi_{jk}-1}}{\prod_{k=1}^{K}{\Gamma(\gamma\pi_{jk})}}\Biggr{]}\Biggl{[}\prod_{j=1}^{J}\prod_{k=1}^{K}\pi_{jk}^{\alpha-1}\Biggr{]}\end{split}$$

$$\begin{array}{rl}{p}_{h,\mathbf{w}}\left(z_{jdik}=1|{\mathbf{z}}^{(-jdi)},\mathbf{\pi }\right)& \propto p\left(z_{jdik}=1,w_{jdiv}=1|{\mathbf{z}}^{(-jdi)},{\mathbf{w}}^{(-jdi)},\mathbf{\pi }\right)\\ & =\frac{p\left(z_{jdik}=1,w_{jdiv}=1,{\mathbf{z}}^{(-jdi)}|{\mathbf{w}}^{(-jdi)},\mathbf{\pi }\right)}{p\left({\mathbf{z}}^{(-jdi)}|{\mathbf{w}}^{(-jdi)},\mathbf{\pi }\right)}\\ & \propto \frac{\gamma {\pi }_{jk}+n_{jdk}^{(-jdi)}}{\gamma +n_{jd.}^{(-jdi)}}\frac{\eta +m_{..kv}^{(-jdi)}}{V\eta +m_{..k.}^{(-jdi)}}\end{array}$$

$$a(\bm{\pi}^{(t)},\bm{\pi}^{(*)})=\frac{p(\bm{\pi}^{(*)})}{p(\bm{\pi}^{(t)})}\frac{\exp\bigl{(}-\frac{1}{2\varepsilon^{2}}{\|\bm{\pi}^{(t)}-{\bm{\mu}}(\bm{\pi}^{(*)},\varepsilon)\|}^{2}\bigr{)}}{\exp\bigl{(}-\frac{1}{2\varepsilon^{2}}{\|\bm{\pi}^{(*)}-{\bm{\mu}}(\bm{\pi}^{(t)},\varepsilon)\|}^{2}\bigr{)}}$$

$$\text{d}\mathbf{\pi }(t)=\frac{1}{2}{\tilde{\nabla }}_{\mathbf{\pi }}\mathcal{L}({\mathbf{\pi }}^{(t)})\text{d}t+\text{d}\tilde{\mathbf{b}}(t)$$

$${\tilde{\nabla }}_{\mathbf{\pi }}\mathcal{L}({\mathbf{\pi }}^{(t)})=\mathbf{G}\{\mathbf{\pi }(t){\}}^{-1}{\nabla }_{\mathbf{\pi }}\mathcal{L}({\mathbf{\pi }}^{(t)})$$

Full text

Keywords— Statistical learning, Unsupervised learning, Text analysis, Topic models

1 Introduction

Newspapers, magazines, scientific journals, and social media messages being composed in daily living produce routinely an enormous volume of text data. The corresponding content comes from diverse backgrounds and represent distinct themes or ideas; modeling and analyzing such heterogeneity in large-scale is crucial in any text mining frameworks. Typically, text mining aims to extract relevant and interesting information from the text by the process of structuring the written text (e.g. via semantic parsing, stemming, lemmatization), inferring hidden patterns within the structured data, and finally, deciphering the results. To address these tasks in an unsupervised manner, numerous statistical methods such as TF-IDF (Salton et al., 1975), latent semantic indexing (Deerwester et al., 1990, LSI), and probabilistic topic models, e.g. probabilistic LSI (Hofmann, 1999, pLSI), latent Dirichlet allocation (Blei et al., 2003, LDA), have emerged in the literature. With minimal human effort, they can be generalized to arbitrary text in any natural language.

Topic models such as LDA are developed to capture the underlying semantic structure of a corpus (i.e. a document collection) based on co-occurrences of words. They consider documents as bags of words, considering the order of words in a document uninformative. LDA (Figure 2) assumes that (a) a topic is a latent (i.e., unobserved) distribution on the corpus vocabulary, (b) each document in the corpus is described by a latent mixture of topics, and (c) each observed word in a document, there is a latent variable representing a topic from which the word is drawn. This model is suitable for documents coming from a single collection. This assumption is insufficient for comparative analyses of text, especially, for partition-able text corpora, as we describe next.

Suppose, we deal with corpora of (i) articles accepted in various workshops or consecutive proceedings of a conference or (ii) blogs/forums from people different countries. It may be of interest to explore research topics across multiple/consecutive proceedings or workshops of a conference or cultural differences in blogs and forums from different countries (Zhai et al., 2004). Suppose, we have political news articles from different sources such as New York Times, Washington Post, The Wall Street Journal, and Reuters. Although, all talk about topic politics, modeling articles from each news source as a collection may help to explore each sources’ article style, policies, region influences, culture, etc. Moreover, the category labels, workshop names, article timestamps, or geotags can provide significant prior knowledge regarding the unique structure of these corpora. We aim to include these prior structures and characteristics of the corpus into a probabilistic model adding less computational burden, for expert analyses of the corpus, collection, and document-level characteristics.

For example, Figure 1 gives the results of an experiment on a real document corpus employing the proposed model in this paper. The corpus consists of articles accepted in the proceedings of the NIPS conference for the years from $1988$ to $2005$. We wish to analyze topics that evolve over this timespan. We thus partitioned the corpus into four collections based on time (details in Section 4). The plot shows estimates of topic proportions for all four collections. Evidently, some topics got increased attention (e.g. Markov chain Monte Carlo (MCMC), generative models) and some topics got decreased popularity (e.g. neural networks) over the years of the conference. For some topics, the popularity is relatively constant (e.g. signal processing) from the beginning of NIPS. The algorithm used is unsupervised, and only takes documents, words, and their collection labels as input.

Assuming a flat structure to all documents in a corpus, LDA or its variants is not well suited to the multi-collection corpora setting: (1) Ignoring predefined structure or organization of collections may end up having topics that describe only some, not all of the collections. (2) There is no direct way to define which set of topics describe the common information across collections and which topic describes information specific to a particular collection. A crude solution is to consider each collection as a separate corpus and fit an LDA model for each corpus. There are reasons why one may want to avoid it: (i) one needs to solve the nontrivial alignment of topics in each model, for any useful comparison of topics among collections (Moreover, the topics inferred from individual corpus partitions and the whole corpora can themselves be different) and (ii) information loss due to modeling collections separately, especially for small datasets (details, Section 4.2).

In this paper, we introduce the compound latent Dirichlet allocation (cLDA, Section 2) model that incorporates any prior knowledge on the organization of documents in a corpus into a unified modeling framework. cLDA assumes a shared, collection-level, latent mixture of topics for documents in each collection. This collection-level mixture is used as the base measure to derive the (latent) topic mixture for each document. All collection-level and document-level variables share a common set of topics at the corpus level, which enables us to perform exciting inferences, as shown in Figure 1. cLDA exhibits a certain degree of supervision incorporating the collection membership for each word (and document) in a corpus implicitly in the modeling framework. cLDA can thus aid visual thematic analyses of large sets of documents, and include corpus, collection, and document specific views.

The parameters of interest in cLDA are hidden and are inferred via posterior inference. However, exact inference is intractable in cLDA; and, non-conjugate relationships in the model further make approximate inference challenging. Popular approximate posterior inference methods in topic models are MCMC and variational methods. MCMC methods enable us to sample hidden variables of interest from the intractable posterior with convergence guarantees. We consider two MCMC methods for cLDA: (a) one uses the traditional auxiliary variable updates within Gibbs sampling, and (b) the other uses Langevin dynamics within Gibbs sampling, a method that received recent attention. Our experimental evidence suggests the former method, which gives superior performance with only a little computational overhead compared to the collapsed Gibbs sampling algorithm for LDA (Griffiths and Steyvers, 2004, CGS) (details in Sections 3.2, 4, F, and H), and is the main focus in this paper. Variational methods are often used in topic modeling as they give fast, parallel implementations by construction. Although they converge rather quickly, our studies show that (Section 3.2 and Section F) their solutions are suboptimal compared to the results of the other two MCMC schemes.

The contributions of this paper are three-fold: (i) we propose a probabilistic model cLDA that can capture the topic structure of a corpus including organization hierarchy of documents, (ii) we study efficient methods for posterior inference in cLDA, and (iii) we perform an empirical study of the real-world applicability of the cLDA model—for example, (a) analyzing topics that evolves overtime, (b) analyzing patterns of topics on customer reviews, and (a) summarizing topic structure of document collections in a corpus—via three text corpora used in the research community. Also, note that the inference about collection-level topic mixtures may be of interest to the general perspective of posterior sampling on the probability simplex in statistics.

The remainder of the paper is organized as follows. Section 2 formally defines the cLDA hierarchical model. Section 3 describes algorithms for posterior inference and evaluates correctness of the algorithms using a synthetically corpus. In Section 4, we assess the performance of the cLDA model and conclude that it exhibits superior performance, both quantitatively (e.g. via perplexity—a popular scheme for evaluating the predictive performance of topic models (Wallach et al., 2009b), and external measures such as topic coherence (Mimno et al., 2011)) and qualitatively. We also compare cLDA with other popular models in the literature, and provide a usability study for cLDA in this section. Section 5 concludes this paper with a summary of our work.

2 A Compound Hierarchical Model

We first set up some terminology and notation. Vectors are denoted by bold, lower case alphabets (e.g. $\bm{\pi}$) and scalar values are denoted by normal, lowercase letters (e.g. $\pi_{jk}$). Matrices or tensors are denoted by bold, upper case Latin alphabets (e.g. ${\bm{G}}$) or bold Greek alphabets without subscripts (e.g. $\bm{\beta}$). There is a vocabulary $\mathcal{V}$ of $V$ terms in the corpus; in general, $\mathcal{V}$ is considered as the union of all the word tokens in all the documents of the corpus, after removing stop-words and normalizing tokens (e.g. stemming). The number of topics $K$ is assumed to be known. (Discussion of how to handle this issue in practice is Section 4.) By definition, a topic is a distribution over ${\cal V}$, i.e., a point in the $V$-$1$ dimensional simplex $\mathbb{S}_{V}$. We will form a $K\times V$ matrix $\bm{\beta}$, whose $k^{\text{th}}$ row is the $k^{\text{th}}$ topic (how $\bm{\beta}$ is formed will be described shortly). Thus, the rows of $\bm{\beta}$ are vectors $\bm{\beta}_{1},\ldots,\bm{\beta}_{K}$, all lying in $\mathbb{S}_{V}$. There are $J$ collections in the corpus and for $j=1,2,\ldots,J$, collection $j$ has $D_{j}$ documents. For $d=1,\ldots,D_{j}$, document $d$ in collection $j$ (i.e. document $jd$) has $n_{jd}$ words, $w_{jd1},\ldots,w_{jdn_{jd}}$. Each word is represented by the index or id of the corresponding term from the vocabulary. We represent document $jd$ by the vector ${\bm{w}}_{jd}=(w_{jd1},\ldots,w_{jdn_{jd}})$, collection $j$ by the concatenated vector ${\bm{w}}_{j}=({\bm{w}}_{j1},\ldots,{\bm{w}}_{jD_{j}})$, and the corpus by the concatenated vector ${\bm{w}}=({\bm{w}}_{1},\ldots,{\bm{w}}_{J})$.

We use $\operatorname{Dir}_{L}(a\omega_{1},\ldots,a\omega_{L})$ to denote the finite-dimensional Dirichlet distribution on the ($L$-$1$)-dimensional simplex. This has two parameters, a scale (concentration) parameter $a$ and a base measure $(\omega_{1},\ldots,\omega_{L})$ on the ($L$-$1$)-dimensional simplex. Thus $\operatorname{Dir}_{L}(a,\ldots,a)$ denotes an $L$-dimensional Dirichlet distribution with the constant base measure $(1,\ldots,1)$. $\operatorname{Mult}_{L}(b_{1},\ldots,b_{L})$ to represents the multinomial distribution with number of trials equal to $1$ and probability vector $(b_{1},\ldots,b_{L})$. Let $h=(\eta,\alpha,\gamma)\in(0,\infty)^{3}$ be the hyperparameters in the model. We formally define the cLDA model as (see Figure 2).

[TABLE]

The distribution of $z_{jd1},\ldots,z_{jdn_{jd}}$ will depend on the document-level variable $\bm{\theta}_{jd}$ that represents a distribution on the topics for document $jd$. A single $\bm{\theta}_{jd}$ for document $jd$ encourages different words in document $jd$ to share the same document level characteristics. The distribution of $\bm{\theta}_{j1},\ldots,\bm{\theta}_{jD_{j}}$ will depend on the collection-level variable $\bm{\pi}_{j}$ which indicates a distribution on the topics for collection $j$. A single $\bm{\pi}_{j}$ for collection $j$ encourages different documents in collection $j$ to have the same collection level properties. A single $\bm{\beta}$ is shared among all documents, which encourages documents in various collections in the corpus to share the same set of topics. Note that the standard LDA model (Blei et al., 2003) is a special case of the proposed cLDA model, where there exist single parameters $\bm{\pi}_{j}$ and $\bm{\beta}$ that fail to capture potential heterogeneity amongst the predefined collections, an objective that cLDA is designed for.

Related Work

Zhai et al. (2004) refer to the problem of multi-collection corpora modeling as comparative text mining (CTM), which uses pLSI as a building block. Comparing with CTM model, cLDA employs an efficient and generalizable LDA-based framework that has several advantages over pLSI—for example, LDA incorporates Dirichlet priors for document topic structures in a natural way to deal with newly encountered documents. Furthermore, cLDA combines collection-specific characteristics in a natural probabilistic framework enabling efficient posterior sampling, depending less on user-defined parameters as in CTM.

The cLDA model shares some similarities, but also exhibits differences from: Hierarchical Dirichlet Process (Teh et al., 2006) and Nested Chinese Restaurant Process (Blei et al., 2004), which are introduced to learn document and topic hierarchies from the data non-parametrically. In nonparametric models, as we observe more and more data, the data representations grow structurally. Instead of manually specifying the number of topics K, these nonparametric topic models infer K from the data by assuming a Dirichlet process prior for document topic distributions and topics. A major hurdle in these frameworks is that inference can be computationally challenging for large datasets. Our experimental evidence also shows that HDP produces too many fragmented topics, which may lead the practitioner to bear the additional burden of post processing (Section 4.2). Here, similar to LDA, cLDA is a parametric model, i.e., the data representational structure is fixed and does not grow as more data are observed. cLDA assumes that K is fixed and can be inferred from the data directly, e.g., empirical Bayes methods (George, 2015), or by cross-validation. We thus have simple Dirichlet priors in cLDA without adding much burden to the model and inference (details appear in Section 3).

In light of adding supervision, several modifications to LDA model have been proposed in the literature. Supervised LDA (Mcauliffe and Blei, 2008, sLDA) is an example; for each document $d$, sLDA introduces a response variable $y_{d}$ that is assumed to be generated from document $d$’s empirical topic mixture distribution. In practice, the posterior inference in such a setting can be inefficient due to the high non-linearity of the discrete distribution on the empirical parameters (Zhu and Xing, 2014). cLDA, on the other hand, proposes a generative framework incorporating the collection-level characteristics, without much computational burden (details, Section 4). Also, the objectives of these related models are different from the focus of this paper.

{notes-to-stay}

One can view the hierarchical model of cLDA as a special instance of the model proposed by Kim et al. (2013). Kim et al. (2013) extended the LDA model to have multiple layers in the hierarchical model, for corpora with many categories and subcategories. In cLDA, we are interested in corpora with just one layer of top level categories. This model is adequate for many real-world corpora such as books with various chapters and conference proceedings from different years. Additionally, we will show that Kim et al. (2013)’s variational methods for the multilayer hierarchical model is suboptimal compared to the proposed MCMC methods for two layered data.

A popular scheme for evaluating the performance of topic models is perplexity (or the predictive probability of held-out documents), which we use here to assess cLDA’s performance quantitatively. In the text domain, one can use trained cLDA models for tasks such as classifying text documents and summarizing document collections and corpora. cLDA also gives excellent options to visualize and browse documents. The model and algorithm we describe can be used to form a hierarchy of documents in a collection including relevant topic information. One can then use the document hierarchy to enlighten one’s perception of the contents of the collection.

3 Posterior Sampling

The parameters of interest in the cLDA model, i.e., (a) corpus-level topics, (b) collection-level mixture of topics, (c) document-level mixture of topics, and (d) topic indices of words are hidden. We identify these hidden variables given the observed word statistics and document organization hierarchy in the corpus via posterior inference.

Let $\bm{\theta}=(\bm{\theta}_{11},\ldots,\bm{\theta}_{1D_{1}},\ldots,\bm{\theta}_{J1},\ldots,\bm{\theta}_{JD_{J}})$, ${\bm{z}}_{jd}=(z_{jd1},\ldots,z_{jdn_{jd}})$ for $d=1,\ldots,D_{j}$, $j=1,\ldots,J$, ${\bm{z}}=({\bm{z}}_{11},\ldots,{\bm{z}}_{1D_{1}},\ldots,{\bm{z}}_{J1},\ldots,{\bm{z}}_{JD_{J}})$, and $\bm{\pi}=(\bm{\pi}_{1},\ldots,\bm{\pi}_{J})$. We will then use $\bm{\psi}$ to denote the latent variables $(\bm{\beta},\bm{\pi},\bm{\theta},{\bm{z}})$ in the cLDA model. For any given $h$, (1)–(5) in the hierarchical model induce a prior distribution $p_{h}(\bm{\psi})$ on $\bm{\psi}$. Equation (5) gives the likelihood $\ell_{{\bm{w}}}(\bm{\psi})$. The words ${\bm{w}}$ and their document and collection labels are observed. We are interested in $p_{h,{\bm{w}}}(\bm{\psi})$, the posterior distribution of $\bm{\psi}$ given ${\bm{w}}$ corresponding to the prior $p_{h}(\bm{\psi})$. Applying Bayes rule, we can write the posterior distribution $p_{h,{\bm{w}}}$ of $\bm{\psi}$ as

[TABLE]

Using (1)–(5) of the hierarchical model, we can write (6) as (Section A gives additional details)

[TABLE]

where $n_{jdk}$ is the number of words in document $d$ in collection $j$ that are assigned to topic $k$, and $m_{jdkv}$ is the number of words in document $d$ in collection $j$ for which the latent topic is $k$ and the index of the word in the vocabulary is $v$. These count statistics depend on both ${\bm{z}}$ and ${\bm{w}}$. Note that the constants in the Dirichlet normalizing constants are absorbed into the overall constant of proportionality. Unfortunately, the normalizing constant of the posterior $p_{h,{\bm{w}}}(\bm{\psi})$, is the likelihood of the data with all latent variables integrated out, is a non-trivial integral. This makes exact inference difficult in cLDA.

Popular methods for approximate posterior inference in topic models are Markov chain Monte Carlo (e.g. see Griffiths and Steyvers (2004)) and variational methods (e.g. see Blei et al. (2003)). Although variational methods may give a fast and scalable approximation for the posterior, due to optimizing the proxy lower-bound, it may not produce optimal solutions as the MCMC methods in practice (e.g. see Teh et al. (2007)). That is the case in this model setting, as our experimental evidence suggests. Hence, we leave the details of our development of variational methods (VEM) for cLDA in Section E, to stay our discussion focused.

3.1 Inference via Markov chain Monte Carlo Methods

According to the hierarchical model (1)–(5), $\bm{\theta}_{jd}$’s and $\bm{\beta}_{k}$’s are independent, and by inspecting the posterior (6), given $(\bm{\pi},{\bm{z}})$, we get:

[TABLE]

where $m_{kv}=\sum_{j=1}^{J}\sum_{d=1}^{D_{j}}m_{jdkv}$ and $d=1,\ldots,D_{j},j=1,\ldots,J,k=1,\ldots,K$. Note that (8) implicitly dependent on the observed data ${\bm{w}}$.

We can integrate out $\bm{\theta}$’s and $\bm{\beta}$’s to get the marginal posterior distribution of $(\bm{\pi},{\bm{z}})$ (up to a normalizing constant) as

[TABLE]

Let the vector ${\bm{z}}^{(-jdi)}$ be the topic assignments of all words in the corpus except for word $w_{jdi}$. And, we define $n_{jd}:=\sum_{k=1}^{K}n_{jdk}$ and $m_{k}:=\sum_{v=1}^{V}m_{kv}$. By inspecting (9), we obtain a closed form expression for the conditional posterior distribution for $z_{jdi}$, given ${\bm{z}}^{(-jdi)}$ and $\pi_{jk}$, as

[TABLE]

where the superscript $(-jdi)$ for the count statistics $n_{jdk}$, $n_{jd}$, $m_{kv}$, and $m_{k}$ means that we discard the contribution of word $w_{jdi}$ for counting. (see Section B) This enables us to build a Gibbs sampling chain on ${\bm{z}}$, by sampling $z_{jdi}$, given $\pi_{jk}$ and ${\bm{z}}^{(-jdi)}$.

Given ${\bm{z}}_{j}$ and the observed data ${\bm{w}}_{j}$, we have the unnormalized posterior probability density function for $\bm{\pi}_{j}$ as

[TABLE]

Here we use the fact that $\bm{\pi}_{j}$’s are independent of $\bm{\beta}_{k}$’s. We wish to sample $\bm{\pi}_{j}$’s from this distribution, however the normalized density function $p_{{\bm{w}}_{j}}(\bm{\pi}_{j}\,|\,{\bm{z}}_{j})$ is computationally intractable: the density $p_{{\bm{w}}_{j}}(\bm{\pi}_{j}\,|\,{\bm{z}}_{j})$ has a non-conjugate relationship between $\bm{\pi}_{j}$’s and $\bm{\theta}_{dj}$’s, which makes its normalizer intractable. Next, we describe two Markov chain Monte Carlo schemes that enable us to sample from this posterior density. We shall only focus on the first scheme due to its superior performance in our numerical experiments.

3.1.1 Auxiliary Variable Sampling

The first scheme is based on the traditional auxiliary variable sampling scheme. The idea of auxiliary variable sampling is that one can sample from a distribution $f(x)$ for variable $x$ by sampling from some augmented distribution $f(x,s)$ for variable $x$ and auxiliary variable $s$, such that the marginal distribution of $x$ is $f(x)$ under $f(x,s)$. One can build a Markov chain using this idea in which, auxiliary variable $s$ is introduced temporarily and discarded, only leaving the value of $x$. Since $f(x)$ is the marginal distribution of $x$ under $f(x,s)$, this update for $x$ will leave $f(x)$ invariant (Neal, 2000). Suppose $\text{S}(n_{jdk},s)$ denotes the unsigned Stirling number of the first kind. We can then get the expression for augmented sampling by plugging in the factorial expansion (Abramowitz, 1974)

[TABLE]

into the marginal posterior density (11) as

[TABLE]

(Newman et al. (2009) and Teh et al. (2006) used a similar idea in a different hierarchical model.) The expression (13) introduces auxiliary variable $s_{jdk}$. By inspecting (13), we get a closed form expression for sampling $\bm{\pi}_{j}$, for $j=1,\ldots,J$:

[TABLE]

and we update the auxiliary variable $s_{jdk}$ by the Antoniak sampling scheme (Newman et al., 2009, Appendix A), i.e., the Chinese restaurant process (Aldous, 1985, CRP) with concentration parameter $\gamma\pi_{jk}$ and the number of customers $n_{jdk}$. The auxiliary variable $s_{jdk}$ is typically updated by drawing $n_{jdk}$ Bernoulli variables as

[TABLE]

In our experience, this augmented update for collection level parameters will add a low computational overhead to the collapsed Gibbs sampling chain on ${\bm{z}}$, and is easy to implement in practice. The Markov chain on $(\bm{\pi},{\bm{z}})$ based on the auxiliary variable update within Gibbs sampling is given by Algorithm 1. We use AGS to denote this chain.

3.1.2 Metropolis Adjusted Langevin Monte Carlo

Another option to define a Markov chain with invariant density $p_{{\bm{w}}_{j}}(\bm{\pi}_{j}\,|\,{\bm{z}}_{j})$ is to employ the Metropolis-Hastings (MH) algorithm (Metropolis et al., 1953; Hastings, 1970) with Langevin dynamics (Girolami and Calderhead, 2011). This scheme has become popular for sampling on the simplex recently (Patterson and Teh, 2013).

Typically, MH algorithm proposes a transition $\bm{\pi}^{(t)}\rightarrow\bm{\pi}^{(*)}$ with density $q(\bm{\pi}^{(*)}\leftarrow\bm{\pi}^{(t)})$—i.e. the proposal—for the current step $t$, and then accept it with probability

[TABLE]

where $\tilde{p}(\bm{\pi})$ denotes the unnormalized density of $\bm{\pi}$. Here, we ignore the subscript $j$ and other dependencies for brevity. The accept-reject step ensures that the proposed Markov chain is reversible with respect to the stationary target density and satisfies detailed balance. The proposal distribution simulates random-walks—e.g. $q(\bm{\pi}^{(*)}\leftarrow\bm{\pi}^{(t)})=\mathcal{N}_{K}(\bm{\pi}^{(*)}\,|\,\bm{\pi}^{(t)},\bm{\Sigma})$, a $K$-dimensional normal distribution with mean $\bm{\pi}^{(t)}$ and covariance matrix $\bm{\Sigma}$. A key challenge of MH in practice is to find a proposal with reasonable acceptance rate, especially when $K$ is large.

Recent developments show that Langevin dynamics (Kennedy, 1990) is an ideal option to define a proposal distribution. Langevin dynamics proposes random walks by a combination of gradient updates and Gaussian noise as follows. We denote the log-density at state $t$ by ${\cal L}(\bm{\pi}^{(t)}):=\log p(\bm{\pi}^{(t)})$. The Langevin diffusion with stationary distribution $p(\bm{\pi})$ is defined by the stochastic differential equation (SDE)

[TABLE]

where ${\bm{b}}$ denotes a $K$-dimensional Brownian motion. Given the current state $t$, we then define a proposal based on the first-order Euler discretization of (16) as

[TABLE]

where $\varepsilon$ is a user defined step-size for the discretization and ${\bm{\xi}}^{(t)}$ is distributed according to a zero mean $K$-dimensional multivariate normal distribution with the identity covariance matrix $\mathds{1}_{K}$. This induces a proposal density $q(\bm{\pi}^{(*)}\leftarrow\bm{\pi}^{(t)})=\mathcal{N}_{K}(\bm{\pi}^{(*)}\,|\,{\bm{\mu}}(\bm{\pi}^{(t)},\varepsilon),\varepsilon^{2}\mathds{1}_{K})$. Note that the discretized process (17) may be transient and is no longer reversible with respect to the stationary density $p(\bm{\pi})$ (Roberts and Tweedie, 1996). However, one can ensure convergence to the target density $p(\bm{\pi})$ via a MH accept-reject scheme (Besag, 1994, p. 591), which we denote by MALA (Section C). One may notice this MALA update as a special case of Hamiltonian Monte Carlo (Neal, 2010, Section 5.5.2).

However, the proposed MALA update for $\bm{\pi}^{(t+1)}$ has some shortcomings. First, the drift term ${\bm{\xi}}^{(t)}$ in the MALA proposal is based on an isotropic diffusion and may be inefficient for strongly correlated variables with widely differing variances, forcing the step size $\varepsilon$ to accommodate variates with smallest variance (Girolami and Calderhead, 2011). Second, $\bm{\pi}^{(t+1)}\in\mathbb{S}_{K}$, the probability simplex $\mathbb{S}_{K}$ is compact, and it needs to handle the cases when MALA proposes a path that’s outside the simplex. Third, typical Dirichlet priors over the probability simplex put most of their probability mass on the edges and corners of the simplex—e.g. in models such LDA (Patterson and Teh, 2013). Computing gradients for these models become unstable when the probabilities are close to zero and causes issues for MALA updates. Our approaches to handle these issues are discussed next.

To handle the first issue, (Girolami and Calderhead, 2011, Section 5) suggested using a preconditioning matrix ${{\bm{G}}}(\bm{\pi})$. They defined ${{\bm{G}}}(\bm{\pi})$ as an arbitrary metric tensor on a Riemannian manifold induced by the parameter space of a statistical model. This requires us to update the natural gradient $\nabla_{\bm{\pi}}{\cal L}(\bm{\pi}^{(t)})$ and Brownian motion $\text{d}{\bm{b}}(t)$ in (16) (details in Section C). We thus use the corresponding Riemannian Manifold Metropolis Adjusted Langevin Algorithm (Girolami and Calderhead, 2011, MMALA) in this paper.

To perform valid moves of $\bm{\pi}$ that lies on the probability simplex (the second and third issues), one needs to consider boundary conditions. A natural solution to handle boundaries is to re-parameterize $\bm{\pi}$ (Patterson and Teh, 2013). Let ${\bm{\varphi}}=(\varphi_{1},\ldots,\varphi_{K})\in\mathbb{R}^{K}$. We take the prior on ${\bm{\varphi}}$ as a product of i.i.d. Gamma random variables as

[TABLE]

We define $|\varphi_{.}|:=\sum_{k=1}^{K}|\varphi_{k}|$. Let $\pi_{k}$ be $|\varphi_{k}|/|\varphi_{.}|$, for each $k=1,2,\ldots,K$. This choice keeps the prior on $\bm{\pi}$ a Dirichlet density. We can then re-write the unnormalized conditional density (11) on $\bm{\pi}$ in terms of ${\bm{\varphi}}$, and derive the MMALA updates, accordingly (Section C). The Markov chain on $(\bm{\pi},{\bm{z}})$ induced by this scheme is denoted by MGS (Algorithm 2).

Note that we can easily augment any of these two chains on $(\bm{\pi},{\bm{z}})$: AGS and MGS to a Markov chain on $\bm{\psi}$ with invariant distribution $p_{h,{\bm{w}}}(\bm{\psi})$ by using the conditional distribution of $(\bm{\beta},\bm{\theta})$, given by (8). The augmented chain will hold the convergence properties of the chain on $(\bm{\pi},{\bm{z}})$.

3.2 Empirical Evaluation of Samples $\bm{\pi}$

We consider a synthetic corpus by simulating (1)–(5) of the cLDA hierarchical model with the number of collections $J=2$ and the number of topics $K=3$. We did this solely so that we can visualize the results of the algorithms. We also took the vocabulary size $V=40$, the number of documents in each collection $D_{j}=100$, and the hyperparameters $h_{\text{true}}=(\alpha,\gamma,\eta)=(.1,1,.25)$. Collection-level Dirichlet sampling via (2) with $\alpha$ produced two topic distributions $\bm{\pi}^{\text{true}}_{1}=(.002,\epsilon,.997)$ and $\bm{\pi}^{\text{true}}_{2}=(.584,.386,.030)$, where $\epsilon$ denotes a small number.

We study the ability of the proposed algorithms AGS, MGS, and VEM (details, Algorithms 2 and  3) for cLDA to recover parameters $\bm{\pi}^{\text{true}}_{1}$ and $\bm{\pi}^{\text{true}}_{2}$. We do this by comparing samples of $\bm{\pi}_{j}$ from the AGS and MGS chains on $(\bm{\pi},{\bm{z}})$ with variational estimates of $\bm{\pi}_{j}$ from VEM iterations111An implementation of all algorithms and datasets discussed in this paper is available as an R package at https://github.com/clintpgeorge/clda. We initialized $\bm{\pi}^{(0)}_{1}=\bm{\pi}^{(0)}_{2}=(.33,.33,.33)$ and use hyperparameters $h_{\text{true}}$ for all three algorithms. Using the data ${\bm{w}}$, we ran both chains AGS and MGS for $2000$ iterations, and algorithm VEM converged after $45$ EM iterations.

Table 1 gives the values of $\bm{\pi}_{1}$ and $\bm{\pi}_{2}$ on the $2$-simplex from the last iteration of all three algorithms. We can see that both AGS and MGS chains were able to recover the values $\bm{\pi}^{\text{true}}_{1}$ and $\bm{\pi}^{\text{true}}_{2}$ reasonably well, although AGS chain has an edge. To converge to optimal regions $\|\bm{\pi}^{\text{true}}_{1}-\bm{\pi}^{(s)}_{1}\|\leq.003$ and $\|\bm{\pi}^{\text{true}}_{2}-\bm{\pi}^{(s)}_{2}\|\leq.07$, AGS chain took cycles $42$ and $23$ only; but, MGS chain required cycles $248$ and $139$, respectively. ($\|.\|$ denotes L-$1$ norm on the simplex $\mathbb{S}_{K}$.) Lastly, algorithm VEM never reached the optimal regions: at convergence, VEM hit points that are $0.08$ far from $\bm{\pi}^{\text{true}}_{1}$ and $0.18$ far from $\bm{\pi}^{\text{true}}_{2}$. Here, the reported values are after the necessary alignment of topics with the true set of topics, for all the three methods. Figure 3 gives trace plots of values of $\bm{\pi}_{2}$ on the $2$-simplex for all three algorithms. Note that due to the superior performance of AGS chain over MGS chain, we use it in our experimental analyses. Additional details of this experiment is provided in Section F.

3.3 Selecting Hyperparameters in cLDA

The hyperparameters $h=(\alpha,\gamma,\eta)$ in cLDA can affect the prior and posterior distributions of parameters $\psi$ and should be selected carefully. A natural solution to select them is via maximum likelihood: for example, we let $m_{{\bm{w}}}(h)$ denote the marginal likelihood of the data as a function of $h$, and use $\hat{h}=\operatorname*{arg\,max}_{h}m_{{\bm{w}}}(h)$. However, for models such as LDA and cLDA, the function $m_{{\bm{w}}}(h)$ is analytically intractable (Section 3). Works in the literature (Blei et al., 2003; Wallach, 2006) suggest an approximate EM algorithm to estimate $\operatorname*{arg\,max}_{h}m(h)$ for LDA, which can be described as follows. We consider ${\bm{w}}$ “observed data,” and $\bm{\psi}$ “missing data.” We have the “complete data likelihood” $p_{h}(\bm{\psi},{\bm{w}})$ available, then the EM algorithm (Dempster et al., 1977) is a natural candidate to estimate $\operatorname*{arg\,max}_{h}m(h)$, since $m(h)$ is the “incomplete data likelihood.” However, the E-step in EM involves calculating an expectation with respect to the intractable posterior $p_{h,{\bm{w}}}$ of the model. One solution is to approximate this expectation via variational inference (e.g. Variational-EM (Blei et al., 2003)) or Markov chain Monte Carlo (e.g. Gibbs-EM (Wallach, 2006)). We follow the latter approach for cLDA, which includes:

Initialize $h_{0}$ and $\bm{\psi}_{0}$ 2. 2.

Until convergence

1. E-step

Sample $\bm{\psi}_{t}^{(1)},\ldots,\bm{\psi}_{t}^{(S)}$ from $p_{h_{t},{\bm{w}}}(\bm{\psi})$ via chain AGS or MGS 2. M-step

$h_{(t+1)}=\operatorname*{arg\,max}_{h}\sum_{s=1}^{S}\log p_{h_{t}}(\bm{\psi}_{t}^{(s)},{\bm{w}})$

One solution to M-step is via the fixed point iteration (Minka, 2000a). We derive the expressions for maximizing hyperparameter $\gamma$ and $\eta$ as follows (see Section D):

[TABLE]

[TABLE]

where $\sum_{x,y,\ldots}$ represents $\sum_{x=1}^{X}\sum_{y=1}^{Y}\ldots$. To illustrate the performance of the proposed method, Figure 4 shows $20$ independent Gibbs-EM estimates of cLDA hyperparameters $(\eta,\gamma)$ with default $\alpha=1$. We used a synthetic corpus with number of collections $J=2$, number of topics $K=3$, vocabulary size $V=40$, collection size = $100$, and document size = $200$ and the “true” hyperparameters $\alpha_{\text{true}}=1,\gamma_{\text{true}}=.8,\eta_{\text{true}}=.5$. We initialized the algorithm with $(\eta_{0},\gamma_{0})=(1,1)$. We notice that Gibb-EM recovered the true hyperparameters $\eta$ and $\gamma$ with a margin of error in all cases.

One approach to deal with hyperparameter $\alpha$ is to put a Gamma prior on $\alpha$ and include a posterior sampling scheme for $\alpha$. But, our sensitivity analyses on $\alpha$ shows no significant impact for $\alpha$ on the predictive power of cLDA, once we select $\eta$ and $\gamma$. We thus use a default value for $\alpha$ in our experiments.

4 Experimental Analysis

This section compares the performance of the proposed cLDA model with alternatives, LDA and HDP, on real-world corpora (Section 4.1). We explore two performance metrics: (a) the quality of the inferred model, which we measure by the predictive power of the model (e.g. via perplexity) and external topic evaluation scores (e.g. topic coherence) (Section 4.2), and (b) the applicability of cLDA model to real-world corpora (Section 4.3). We also briefly discuss guidelines for selecting the number of topics $K$ in cLDA (Section 4.2).

4.1 Datasets

We use three document corpora based on: (i) the NIPS $00$-$18$ dataset (Globerson et al., 2007), (ii) customer reviews from Yelp222http://uilab.kaist.ac.kr/research/WSDM11/, and (iii) the $20$Newsgroups dataset333http://qwone.com/~jason/20Newsgroups. Standard corpus preprocessing involve tokenizing text and discarding standard stop-words, numbers, dates, and infrequent words from the corpus. The NIPS $00$-$18$ dataset is a popular dataset used in the topic modeling community. It consists of papers published in proceedings [math] to $18$ of the Neural Information Processing Systems (NIPS) conference (i.e. years from $1988$ to $2005$). Along with standard preprocessing, we discarded words with length less than two or with length two and not in (“ml”, “ai”, “kl”, “bp”, “em”, “ir”, “eb”) from the corpus. Finally, this corpus consists of $2{,}741$ articles and $9{,}156$ unique words. A question of interest is to see how various research topics evolve over the years $1988$-$2005$ of NIPS proceedings. Typically, new topics do not emerge in consecutive years. So, we partition the NIPS $00$-$18$ corpus into four collections based on time periods $1988$-$1992$, $1993$-$1997$, $1998$-$2002$, and $2003$-$2005$ in our analyses.

Yelp is a crowd-sourced platform for user reviews and recommendations of restaurants, shopping, nightlife, entertainment, etc. In this paper, we used a collection of restaurant reviews from Yelp. Each text review (i.e. a document) in the collection is associated with a customer rating on the scale $1$ to $5$, where $5$ being the best and $1$ being the worst. We discarded reviews with less than $50$ words. After standard preprocessing, this corpus consists of $24{,}310$ restaurant reviews and $9{,}517$ unique words.

The $20$Newsgroups dataset consists of approximately $20{,}000$ news articles that are distributed across $20$ different newsgroups. According to the subject matter of articles, these $20$ newsgroups are further partitioned into six groups. We took a subset of this dataset with the subject groups computers, recreation, science, and politics. After standard preprocessing, this corpus consists of $10{,}764$ articles and $9{,}208$ unique words. We consider the subject group of a document as its collection label to fit cLDA models. We denote this corpus by $16$Newsgroups.

4.2 Model Evaluation and Model Selection

4.2.1 Comparisons with Models in the Prior Work

(Wallach et al., 2009a, Section 3) proposed non-uniform base measures instead of the typical uniform base measures used in the Dirichlet priors over document-level topic distributions $\bm{\theta}_{d}$’s in LDA. The purpose was to get a better asymmetric Dirichlet prior for LDA, which is generally used for corpora with a flat organization hierarchy for documents. We can consider one such model as a particular case of the cLDA model proposed here: their model is closely related to a cLDA model with a single collection. To illustrate, we perform a comparative study of perplexity scores of cLDA with single and multi-collection in our experiments.

Very briefly, HDP can be described as follows. Suppose we have $D$ populations, and that for population $d,\,d=1,\ldots,D$, there are observations $w_{di}\stackrel{{\scriptstyle\text{ind}}}{{\sim}}F_{z_{di},\sigma_{di}},\,i=1,\ldots,n_{d}$. Here, $F_{z_{di},\sigma_{di}}$ is a distribution depending on some latent variable $z_{di}$ and possibly also on some other known parameter $\sigma_{di}$ particular to individual $di$. We assume that $z_{di}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}G_{d},\,i=1,\ldots,n_{d}$, and that for $d=1,\ldots,D$, $G_{d}\stackrel{{\scriptstyle\text{iid}}}{{\sim}}{\cal D}_{G_{0},\alpha}$, the Dirichlet process (DP) with base probability measure $G_{0}$ and precision parameter $\alpha>0$ (Ferguson, 1973, 1974). In applications such as topic modeling, it is desirable to model the distributions of the $w_{di}$’s as mixtures, and to have mixture components shared among the distributions of the $w_{di}$’s in different populations. We can obtain this by taking $G_{0}$ itself to have a Dirichlet process prior, $G_{0}\sim{\cal D}_{{\cal K},\gamma}$, where ${\cal K}$ is a probability distribution and $\gamma>0$. In the topic modeling setting, for word $w_{di}$ of document $d$, we imagine that there exists a topic $z_{di}$—a distribution on ${\cal V}$, from which the word is drawn. Typically, the distribution ${\cal K}$ is a member of a known parametric family such as Dirichlet with parameter $\omega$. The hierarchical model described here is a two-level DP (Teh et al., 2006), which is (conceptually) similar to a one-collection cLDA model. We thus compare the results of one-collection cLDA and LDA with the results of two-level DP444Code by Wang and Blei (2010): https://github.com/blei-lab/hdp. We use the popular Chinese Restaurant Franchise (Teh et al., 2006, CRF) sampling scheme for inference.

4.2.2 Criteria for Evaluation

We use the perplexity computation scheme that is specified as in Patterson and Teh (2013); Wallach et al. (2009c). We first divide documents in the corpus into a training set and a held-out set. Second, we partition words in every document ${\bm{w}}_{jd}$ in the held-out set to two sets of words, ${\bm{w}}_{jd}^{\text{train}}$ and ${\bm{w}}_{jd}^{\text{test}}$. We also use ${\bm{w}}_{jd}^{\text{train}}$ to denote words in a training document $jd$. We define the vector ${\bm{w}}^{\text{train}}$ for the training corpus combining ${\bm{w}}_{jd}^{\text{train}}$, $j=1,\ldots,J$, $d=1,\ldots,D_{j}$. Similarly, we define the vector ${\bm{w}}^{\text{test}}$. We compute the per-word perplexity for the held-out words ${\bm{w}}^{\text{test}}$ (uses the fact that $z_{jdi},\,i=1,\ldots,n_{jd}$ are conditionally independent given $\bm{\theta}_{jd}$) as

[TABLE]

where $|{\bm{w}}^{\text{test}}|$ is the length of vector ${\bm{w}}^{\text{test}}$. Exact computation of this score is intractable. However, we can estimate this score via MCMC or variational methods for both cLDA and LDA models (see Section G).

Mimno et al. (2011) suggested some generic evaluation scores based on human coherence judgments of estimated topics via topic models such as LDA. One such score is the topic size, i.e., the number of words assigned to a topic in the corpus. One can estimate it by samples from the posterior of the topic latent variable ${\bm{z}}$, given observed words ${\bm{w}}$ (e.g. via Markov chains: collapsed Gibbs sampling (CGS) for LDA (Griffiths and Steyvers, 2004), AGS, and CRF).

Another option is the topic coherence score (Mimno et al., 2011), which is computed based on the most probable words in an estimated topic for a corpus. We find the $m$ most probable words $v_{1}^{(k)},v_{2}^{(k)},\ldots,v_{m}^{(k)}$ for topic $k$ by sorting the vocabulary words in topic $k$ in the descending order of topic specific probabilities (i.e. $\beta_{kt},t=1,\ldots,V$, estimated via Markov chains such as CGS, AGS, and CRF) and picking the top $m$ words. Let $\text{df}(v_{t})$ be the document frequency of term $v_{t}$, i.e., the number of documents in the corpus which have the term $v_{t}$. Let $\text{df}(v_{i},v_{j})$ be the co-document frequency of the terms $v_{i}$ and $v_{j}$, i.e., the number of documents in the corpus which have both of the terms $v_{i}$ and $v_{j}$. We then define coherence score (Mimno et al., 2011), for topic $k=1,2,\ldots,K$, as

[TABLE]

The intuition behind this score is that group of words belonging to a topic possibly co-occur with in a document.

4.2.3 Comparing Performance of cLDA and LDA Models

We first look at the criterion perplexity scores for both LDA and cLDA models with various values of the number of topics keeping each model hyperparameter fixed (and default). We ran the Markov chains AGS and CGS for $1000$ iterations. Figure 5 gives (average) per-word perplexities for the held-out (test) words using algorithms CGS and AGS (with single collection, i.e., $J=1$ and multi-collection, $J>1$) for corpora NIPS $00$-$18$, $16$Newsgroups, and yelp. From the plots, we see that cLDA outperforms LDA well in terms of predictive power, except for corpus NIPS $00$-$18$, for which cLDA has a slight edge over LDA. We believe the marginal performance for corpus NIPS $00$-$18$ is partly due to the nature of partitions defined in the corpus; collections in this corpus share many common topics, compared to corpus $16$Newsgroups. It also suggests that the gain in using cLDA is larger when we model corpora with separable collections.

In terms of perplexity, we see a small improvement for cLDA models with multi-collection over single collection. But, the advantage of cLDA over LDA is clear, even with a single collection. Also, note that cLDA has a better selection of priors as well as the ability to incorporate document hierarchy in a corpus into the modeling framework.

Selecting $K$ and Hyperparameters $\alpha,\gamma,\eta$: Figure 5 also gives some insights on the number of topics $K$ should be used for each corpus. Perplexities of cLDA models go down quickly with increase in the number of topics in the beginning of the curves, and then the rate of decrease go steady after reasonable number of topics (e.g. $K=90$ corpus NIPS $00$-$18$). LDA, on the other hand, has a “U” curve with increase in perplexity after $K=30$ or $K=40$, for corpora $16$Newsgroups and Yelp. However, for corpus NIPS $00$-$18$, LDA shows a behavior similar to cLDA. We thus use $K=90$ for corpus NIPS $00$-$18$, $K=30$ for corpus $16$Newsgroups, and $K=40$ for corpus Yelp, in our comparative analyses in this section and case studies.

To estimate hyperparameters $\eta$ and $\gamma$ in cLDA, we use the Gibbs-EM algorithm (Section 3.3) for all three corpora with the selected $K$ and constant $\alpha=1$. Estimated $(\hat{\eta},\hat{\gamma})$’s for cLDA models at convergence for corpora $16$Newsgroups, NIPS $00$-$18$, and Yelp are $(.05,2.07)$, $(.026,6.72)$, and $(.034,3.44)$, respectively. Similarly, one can estimate hyperparameters $\alpha$ and $\eta$ in LDA based on a similar Gibbs-EM algorithm (see (19) and Wallach (2006)). Estimated $(\hat{\eta},\hat{\alpha})$’s for LDA models at convergence for corpora $16$Newsgroups, NIPS $00$-$18$, and Yelp are $(.053,.048)$, $(.027,.06)$, and $(.029,.12)$, respectively. Note that estimates $\hat{\eta}$’s for LDA and cLDA are quite similar for all three corpora.

We now evaluate distributions of topics learned by cLDA and LDA models for corpora $16$Newsgroups, NIPS $00$-$18$, and Yelp with chosen $K$ and hyperparameters $h$, quantitatively. Figure 6 shows boxplots of estimated topic sizes and coherences for models LDA and cLDA for all three corpora. Note that LDA models have uniform topic sizes compared to cLDA models. Coherence scores of cLDA topics are better than LDA topics except for corpus $16$Newsgroups, which we think, is partially due to having relatively easily separable topics compared to the other two complex corpora.

Execution time: In our experience, both AGS (with a reasonable number of collections) and CGS chains have comparable computational cost for all three corpora (details in Table 2). On average, the combined execution-time of multiple CGS runs for a corpus that is partitioned on collection labels was fairly equivalent to the execution-time of a single CGS run on the whole corpus. We implemented all these algorithms single-threaded on the R-C++ (using the efficient Rcpp and RcppArmadillo libraries) programming environment.

4.2.4 Comparing Performance of cLDA with LDA and HDP

We now compare the performance of one-collection cLDA and LDA with two-level DP using corpus NIPS $00$-$18$. We ran all the three chains AGS, CGS, and CRF for $1{,}000$ iterations and used the last sample from each chain for our analysis. We used fixed $K=90$ for both algorithms AGS and CGS, but algorithm CRF inferred $204$ topics from the corpus. We first evaluate the quality of topics learned for each method, by applying the hclust algorithm on topic distributions. hclust first computes a topic-to-topic similarity matrix based on the manhattan distance. Then, it builds a hierarchical tree in a bottom-up approach. We noticed that the hierarchies formed by cLDA and LDA are comparable. By looking closely at topic word distributions, we notice that cLDA produces better topics, exploiting the asymmetric hierarchical Dirichlet prior on document $\bm{\theta}_{d}$’s that which is absent in LDA. (details, Section H). HDP, on the other hand, found too many redundant topics, as shown in Figure 8. To evaluate clusters induced by hclust quantitatively, we look at silhouette widths computed on topics’ hclust clusters. We favor methods with high silhouette widths. Figure 7 gives boxplot statistics (i.e. median, lower hinge, and upper hinge) of silhouette widths computed on topics’ hclust clusters with various values of the number of clusters (a user specified value in hclust). Overall, cLDA topics outperform HDP topics, and cLDA topics are comparable or better than LDA topics. Our analysis on clustering on learned document topic proportions, i.e., $\bm{\theta}_{d}$s, also show similar results (Figure 16)

4.3 Usability Study

In text mining, one can use trained cLDA models for tasks such as classifying text documents and summarizing document collections and corpora. cLDA also gives excellent options to visualize and browse documents. To illustrate some of these, we provide a usability study for the proposed cLDA model next. We are interested in two aspects of a learned cLDA model: (a) interpretability of the learned topics and topic structures, and (b) visualizing corpora in terms of learned topics.

cLDA infers collection-level and document-level topic distributions based on a common set of corpus-level topics. So it is natural to check whether cLDA’s topics are as interpretable as LDA’s topics. We present three examples of LDA and cLDA topics t$8$, t$57$, and t$85$ in Figure 9 to furnish a qualitative comparison of these models, for corpus NIPS $00$-$18$. For each topic $k$, we show $20$ most probable words based on the non-decreasing order of the $\beta_{kv}$ values over words $v$’s (denoted by bars in each plot). Note that we did any necessary alignments for the set of topics from each of these models to ease comparison. We labeled topics based on their most probable words. (We do these post processing steps for this section and all the following.) Note that cLDA topics are meaningful, easily interpretable, and comparable with LDA topics.

A typical question of interest in topic modeling is to identify topics that evolve over time. cLDA provides a way to model time evolving corpus. On the other hand, LDA has some limitations due to several issues such as the required alignment of topics and information loss of segmentation, as mentioned in Section 1. To illustrate some of these issues involved, we study corpus NIPS $00$-$18$—consists of four collections based on document timestamps—using both cLDA and LDA models. Once we fit a cLDA model (e.g. via AGS chain) for the corpus, we can use each estimated collection-level topic distribution for a time period as its natural topic allocation. Recall that LDA does not have collection-level topic mixtures by the model construction. However, one can estimate them via each collection’s word topic vector ${\bm{z}}_{j}$ (e.g. sampled from the CGS chain).

Figure 10 shows estimates of four collection-level topic distributions $\bm{\pi}_{1},\bm{\pi}_{2},\bm{\pi}_{3}$, and $\bm{\pi}_{4}$ for corpus NIPS $00$-$18$ based on three different topic modeling approaches, which we will describe below. In Figure 10a, we applied the cLDA AGS algorithm on the whole corpus with four collections, i.e., $J$ = 4 (denoted by M$1$). Each barplot represents the value of a topic element in the last sample $\bm{\pi}^{(1000)}$ from the AGS chain. (We show values for a few selected topics here; see an extended list of topics in the Appendix, Figure 18.) Figure 10b is based on the LDA CGS algorithm on the whole corpus (denoted by M$2$). To estimate collection $j$’s topic mixture, we used the last sample ${\bm{z}}_{j}^{(1000)}$ from the CGS chain. In Figure 10c, we considered each collection in the corpus as a separate corpus. We then ran the LDA CGS algorithm on each sub-corpus (denoted by M$3$) one by one. We used the last sample ${\bm{z}}_{j}^{(1000)}$ from collection $j$’s CGS chain to estimate its topic mixture.

One can infer interesting patterns from the estimates of collection-level topic mixtures via method M$1$, as shown in Figure 10a. Several research topics got increased and decreased interest, and some topics were relatively constant in popularity over the time span $1988$-$2005$ of NIPS. As we would have expected, topic t$57$, a topic on MCMC and inference, and topic t$85$, a topic on generative models got increased popularity over the time span, which is interesting to watch. One topic that lost popularity is t$8$, which is about neural networks and related topics. (Additional details are given in Section H.) We can also see that the estimates via method M$1$ are superior to the estimates based on other two methods M$2$ and M$3$. Note that for method M$3$, one must solve the non-trivial topic alignment problem for each LDA model learned for every sub-corpus. Additionally, we notice that method M$2$ gives better estimates than method M$3$. We believe this is due to the lack of information sharing among the partitions in method M$3$ that might have affected the quality of the learned LDA models. On the other hand, cLDA considers the corpus as a whole to incorporate the organization hierarchy of documents into the model, eliminating the issues discussed.

Modeling customer behaviors is the core interest in operations management community and market research. Analyzing patterns of topics that pervade through customer reviews can give interesting insights about customer needs. We wish to study topic patterns for different customer ratings, and we thus use corpus Yelp. A natural choice is to partition the corpus into collections based on customer review ratings on the scale $1$-$5$, and fit a cLDA model. Note that the two other corpora used in this paper have collections that occurred more naturally; but, for this corpus, we follow this partitioning scheme just for experimentation. Figure 11a shows estimated collection-level topic distributions of a $40$-topic cLDA model trained on corpus Yelp via AGS chain. Specifically, each barplot represents the value of a topic element for a collection from the last sample $\bm{\pi}^{(1000)}$ of the AGS chain. We can see a gradual increase of probabilities for some topics (e.g. t$3$, t$11$, t$26$) and gradual decrease of probabilities for some topics (e.g. t$7$, t$14$, t$18$), for rating from $1$ to $5$. We noticed that the topics with gradual increase in probabilities correspond to positive reviews, and the topics with gradual decrease in probabilities correspond to negative reviews.

Table 3 gives additional details of these topics and their manually assigned sentiment. In our experience, people do not give detail reviews about the place, ambiance, food, or customer service, if they really like them (e.g. reviews with rating $5$). If they partially like a place (e.g. reviews with rating $3$ and $4$), they provide comments with all the details. Additionally, topics that show little or relatively small variability for different ratings are general topics about food styles. For completeness, we include the collection-level topic distributions estimated via a single CGS chain (i.e. via method M$2$) for the same Yelp corpus with $K=40$. Note that LDA is inferior in giving any interpretable analyses from the results, as shown in Figure 11b.

A key feature of the cLDA model, compared to its non-hierarchical alternative, LDA, is that it provides collection level topics, which summarize themes or topics of individual collections in a corpus. We perform a study on the applicability of this feature using corpus $16$newsgroups. Recall that corpus $16$newsgroups is partitioned based on four subject groups. It is natural to check whether each collection’s topic distribution exhibits proper weights to the topics that are related to the subject matter of the collection. Figure 12 shows estimated collection-level topic distributions, i.e., $\bm{\pi}_{j},\,j=1,\ldots,4$, from a learned $30$-topic cLDA model via the AGS algorithm. (We found that topic $10$, which is dominant in all collections, is a potential collection of stop-words in the corpus. Figure 19 provides examples of topics for reference.) The results show that cLDA enables us to summarize topics of individual collections and perform a meaningful comparison of topics among collections, an aspect that is not well exploited in a flat modeling framework such as LDA. We also notice sparse allocations for topics in each collection and no major sharing of topics among collections in Figure 12: this conveys the fact that corpus $16$newsgroups may be surely separable via the subject matter of newsgroups.

5 Summary

In this paper, we developed a topic model, compound latent Dirichlet allocation (cLDA), that can incorporate several characteristics of a corpus including the organization hierarchy of documents. One can employ cLDA model to (a) explore research topics across multiple/consecutive proceedings or workshops of a conference, (b) discover cultural differences in blogs and forums from different countries, or (c) compare nuanced review summarization on various customer ratings (for customer relations management). We proposed posterior inference methods for cLDA (e.g. AGS, MGS, VEM), which can be used for analyzing collections of documents and we recommend algorithm AGS. We also discussed guidelines for selecting $K$ and hyperparameters in cLDA, which perform well empirically. Additionally, we have seen that proposed parametric hierarchy in cLDA adds only a little computational burden compared to the non-hierarchical alternative, LDA. cLDA model makes several natural assumptions about the data generating framework, which we believe helped us getting a much simpler inference scheme, compared to nonparametric alternatives such as Hierarchical Dirichlet Process (Teh et al., 2006). We have also provided an empirical comparison of cLDA with two-level DP.

Lastly, we studied the applicability and performance of cLDA model in both synthetic and real-world corpora. cLDA is quite useful for analyzing topic distributions of collections in a corpus, and it can better understand the underlying thematic structure of the corpus. cLDA provides relevant insights about topics that evolve over time, as noted in the experiments on NIPS conference proceedings. This model also presents some ideas to employ metadata such as customer ratings for “soft” segmentation while analyzing the thematic structure of customer reviews. Although such partitions do not emerge naturally, together with cLDA, they may provide diverse perspectives of customer behaviors.

{notes-to-stay}

We also wish to study the applicability of this model for a corpus where documents exhibit multiple collection memberships. For example, a corpus built from the Wikipedia articles, where an article can be assigned to one or more Wikipedia categories.

Appendix A Expressions for the Prior, Likelihood, and Posterior

This section derives expressions of the prior $p_{h}(\bm{\psi})$, likelihood $\ell_{{\bm{w}}}(\bm{\psi})$, and the posterior $p_{h,{\bm{w}}}(\bm{\psi})$. From the hierarchical model (1)–(5) given in Section 2, we can write the prior $p_{h}(\bm{\psi})$ as

[TABLE]

where the individual density functions are given by (1)–(4) of the model. Let $n_{jdk}=\sum_{i=1}^{n_{jd}}z_{jdik}$, i.e. $n_{jdk}$ is the number of words in document $d$ in collection $j$ that are assigned to topic $k$, and let $n_{j.k}=\sum_{d=1}^{D_{j}}\sum_{i=1}^{n_{jd}}z_{jdik}$, i.e. $n_{j.k}$ is the number of words in all documents in collection $j$ that are assigned to topic $k$. Using the Dirichlet and multinomial distributions specified in (1)–(4), we obtain

[TABLE]

For $j=1,\ldots,J$, $d=1,\ldots,D_{j}$ and $k=1,\ldots,K$, let $S_{jdk}=\{i:1\leq i\leq n_{jd}\text{ and }z_{jdik}=1\}$, which is the set of indices of all words in document $d$ in collection $j$ whose latent topic variable is $k$. With this notation, Equation (4) induces the likelihood function $\ell_{{\bm{w}}}(\bm{\psi}):=p({\bm{w}}\,|\,{\bm{z}},\bm{\theta},\bm{\pi},\bm{\beta})$

[TABLE]

where $m_{jdkv}=\sum_{i\in S_{jdk}}w_{jdiv}$ counts the number of words in document $d$ in collection $j$ for which the latent topic is $k$ and the index of the word in the vocabulary is $v$. Recalling the definition of $n_{jdk}$ given just before (23), and noting that $\sum_{i\in S_{jdk}}w_{jdiv}=\sum_{i=1}^{n_{jd}}z_{jdik}w_{jdiv}$, we see that $m_{jdkv}=\sum_{i=1}^{n_{jd}}z_{jdik}w_{jdiv}$ and $\sum_{v=1}^{V}m_{jdkv}=n_{jdk}$

Plugging the likelihood (24) and the prior (23) into the expression for the Bayes rule (6), and absorbing constants in Dirichlet normalizing constants into an overall constant of proportionality, we have the posterior density

[TABLE]

Appendix B Conditional posterior for $z_{jdi}$

This section derives a closed form expression for the conditional posterior $p_{h,{\bm{w}}}\left(z_{jdik}=1\,|\,{\bm{z}}^{(-jdi)},\bm{\pi}\right)$ for $k=1,\ldots,K$. The vector ${\bm{z}}^{(-jdi)}$ contains topic assignments of all words in the corpus except for word $w_{jdi}$. We can then write

[TABLE]

Here, we used the fact that $p\left({\bm{z}}^{(-jdi)}\,|\,{\bm{w}}^{(-jdi)},\bm{\pi}\right)\propto p_{h,{\bm{w}}}(\bm{\pi},{\bm{z}})^{(-jdi)}$ given by (9) and $\Gamma(x+1)=x\Gamma(x)$. The superscript $(-jdi)$ means that we discard the contribution of word $w_{jdi}$ in count statistics $n_{jdk}$, $n_{jd.}$, $m_{..kv}$, and $m_{..k.}$. This development is motivated by the LDA collapsed Gibbs sampling algorithm (Griffiths and Steyvers, 2004, CGS).

Appendix C Langevin Monte Carlo

The Metropolis Adjusted Langevin Algorithm  (Girolami and Calderhead, 2011, MALA), as described in Section 3.1.2, is given the following steps

Propose $\bm{\pi}^{(*)}$ via Langevin dynamics (17) 2. 2.

Calculate the MH acceptance ratio

[TABLE]

and set $\bm{\pi}^{(t+1)}=\bm{\pi}^{(*)}$ with probability min$(1,a(\bm{\pi}^{(t)},\bm{\pi}^{(*)}))$, and set $\bm{\pi}^{(t+1)}=\bm{\pi}^{(t)}$ with the remaining probability.

C.1 Riemann Manifold Metropolis Adjusted Langevin

Algorithm

Girolami et al. (Girolami and Calderhead, 2011, Section 5) suggests the preconditioning matrix ${{\bm{G}}}(\bm{\pi})$ as an arbitrary metric tensor on a Riemannian manifold induced by the parameter space of a statistical model. We write the stochastic differential equation for the Langevin diffusion on the Riemannian manifold as

[TABLE]

where the natural gradient is

[TABLE]

and the Brownian motion is

[TABLE]

where the subscript $k$ indicates $k^{\text{th}}$ element in vector $\text{d}\tilde{{\bm{b}}}(t)$. Note that in the Euclidean space the metric tensor ${{\bm{G}}}(\bm{\pi})$ is an identity matrix, and thus (27) will be reduced to the standard Langevin SDE.

By expanding the gradients in (29) and discretizing (27) via the first-order Euler integration, we get the proposal for the Riemann Manifold Metropolis Adjusted Langevin Algorithm (MMALA):

[TABLE]

where $k^{\text{th}}$ element in vector ${\bm{\mu}}(\bm{\pi}^{(t)},\varepsilon)$ is given by

[TABLE]

The corresponding proposal density is given by

[TABLE]

C.2 Langevin Updates on Probability Simplices: Boundary

Considerations

We first re-write the unnormalized conditional density (11) on $\bm{\pi}_{j}$ as

[TABLE]

The MMALA update for the new parametrization is given by

[TABLE]

where for $k=1,\ldots,K$, we have

[TABLE]

The proposal density of this diffusion process is given by

[TABLE]

We take the metric tensor ${\bm{G}}({\bm{\varphi}}^{(t)}_{j})$ as $\operatorname{diag}(|{\bm{\varphi}}^{(t)}_{j}|)^{-1}$ (Patterson and Teh (2013) suggested this choice in a different context.), as it gives simplified expressions for

[TABLE]

and

[TABLE]

To derive $\nabla_{{\bm{\varphi}}_{j}}{\cal L}({\bm{\varphi}}_{j}^{(t)})$, we first write $\tilde{p}({\bm{\varphi}}_{j}\,|\,{\bm{z}}_{j},{\bm{w}}_{j})$ in a convenient logarithmic form:

[TABLE]

We then define $\nabla_{{\bm{\varphi}}_{j}}{\cal L}({\bm{\varphi}}_{j}^{(t)})$ by the partial derivatives

[TABLE]

Similar to the AGS chain, in this scheme, we implement a Markov chain on $(\bm{\pi},{\bm{z}})$ via MMALA updates within Gibbs sampling (MGS), as in Algorithm 2.

Appendix D Estimating hyperparameters $\eta$ and $\gamma$

In this section, we derive expressions for the fixed point iterations of hyperparameters $\eta$ and $\gamma$ in the cLDA model. We use the following two bounds by (Minka, 2000a, Appendix B) in our development.

[TABLE]

From the cLDA hierarchical model (1)–(5), after integrating out $\bm{\beta}$’s, we get the marginal posterior of $(\bm{\pi},{\bm{z}})$, given ${\bm{w}}$ and $\eta$ as:

[TABLE]

Using (40) and (41), we write it as:

[TABLE]

Here, $\eta_{0}$ represents the value of $\eta$ at the current state. We now denote the lower bound on $\log p_{\eta,{\bm{w}}}(\bm{\pi},{\bm{z}})$ as ${\cal L}_{\eta}$, which we write

[TABLE]

Taking derivatives with respect to $\eta$,

[TABLE]

and setting it to zero, we get

[TABLE]

We can compute the maximum via the fixed point iteration (Minka, 2000b).

Similarly, from the cLDA hierarchical model (1)–(5), after integrating out $\bm{\theta}$’s, we get the marginal posterior of $(\bm{\pi},{\bm{z}})$, given ${\bm{w}}$ and $\gamma$ as:

[TABLE]

Using (40) and (41), we write it as:

[TABLE]

Here, $\gamma_{0}$ represents the value of $\gamma$ at the current state. We now denote the lower bound on $\log p_{\gamma,{\bm{w}}}(\bm{\pi},{\bm{z}})$ as ${\cal L}_{\gamma}$, which we write

[TABLE]

Taking derivatives with respect to $\gamma$,

[TABLE]

and setting it to zero, we get

[TABLE]

We can compute the maximum via the fixed point iteration (Minka, 2000b).

Appendix E Variational Inference

We now develop variational methods (Jordan et al., 1999) to approximate the intractable posterior $p_{h,{\bm{w}}}(\bm{\psi})$ in the cLDA model. Our approach can be viewed as an extension of Blei et al. (2003)’s inference scheme in the LDA model. Briefly, in variational methods, one considers a restricted family of distributions instead of working on the intractable posterior, and then seeks the member of the family that is “closest” to the posterior (Jordan et al., 1999; Bishop et al., 2006). One way to restrict the family of approximating distributions is to use a parametric distribution (i.e., variational distribution) that is governed by a set of parameters (i.e., variational parameters). Typically, this parametric distribution is much simpler to work with than the original posterior by assuming independence between respective variables. The goal is then to identify the parameters which give the tightest lower-bound with in the family.

Let $p_{h}(\bm{\psi},{\bm{w}})$ be the joint probability of $\bm{\psi}=(\bm{\beta},\bm{\pi},\bm{\theta},{\bm{z}})$ and ${\bm{w}}$ based on the cLDA model. Suppose $q(\bm{\psi})$ is any parametric distribution over latent variables $\bm{\psi}$. We can then write the log marginal probability of the data ${\bm{w}}$ as (Bishop et al., 2006)

[TABLE]

where555The summation $\sum_{{\bm{z}}}$ represents the summation over all $z_{jdi}$s. We use summation instead of an integral because $z_{jdi}$s are discrete.

[TABLE]

and

[TABLE]

Note that ${\cal L}(q,p_{h})$ in (59) is a functional of the distribution $q(\bm{\psi})$ and a function of the hyperparameters $h$. The Kullback-Leibler (KL) divergence specified in (60) satisfies $\text{KL}(q,p_{h,{\bm{w}}})\geq 0$—by the positivity of the KL divergence, with equality if, and only if, $q(\bm{\psi})$ equals the posterior $p_{h,{\bm{w}}}(\bm{\psi})$. Following (58), ${\cal L}(q,p_{h})$ is a lower-bound for the log marginal probability. We can maximize the lower-bound ${\cal L}(q,p_{h})$ with respect to $q(\bm{\psi})$, which is also equivalent to minimizing $\text{KL}(q,p_{h,{\bm{w}}})$. The tightest lower-bound occurs when the KL divergence vanishes, i.e., when $q(\bm{\psi})$ equals the posterior distribution (but it is intractable to work with). Thus, in variational methods, one considers a restricted family of distributions $q(\bm{\psi})$ instead of working on the intractable posterior, and then seeks the member of the family for which the lower-bound ${\cal L}(q,p_{h})$ is maximized.

For cLDA, we define a fully factorized variational distribution with the variational parameters $\bm{\psi}^{\prime}=({\bm{\lambda}},{\bm{\tau}},{\bm{\rho}},\bm{\phi})$ as

[TABLE]

We take its independent component distributions Blei et al. (2003) are

[TABLE]

Note that the lower-bound (59), which is an expectation with respect to (61), is intractable due to the non-conjugate relationships between $\bm{\pi}_{j}$’s and $\bm{\theta}_{jd}$’s. We will also see that estimation of variational parameters $({\bm{\lambda}},{\bm{\rho}},\bm{\phi})$ follows closely to Blei et al. (2003)’s scheme, but updating parameter ${\bm{\tau}}$ does not have a closed form expression. Kim et al. (2013) proposed a solution to an expectation of similar form in a different context, which we employ here. The corresponding variational expectation maximization scheme for cLDA is described in Algorithm 3 and is denoted by the acronym VEM.

We now describe a way to estimate the variational parameters $({\bm{\lambda}},{\bm{\tau}},{\bm{\rho}},\bm{\phi})$ via minimizing the KL divergence between the posterior $p_{h,{\bm{w}}}(\bm{\psi})$ and the variational distribution $q(\bm{\psi}\,|\,\bm{\psi}^{\prime})$. We first write down the variational lower-bound (59) as follows.

[TABLE]

To ease notation, we denote $\beta_{k.}=\sum_{v=1}^{V}\beta_{kv}$. Let $\bm{\theta}\sim\operatorname{Dir}_{L}({\bm{\rho}})$. We evaluate the expectation of the log of a single probability component $\theta_{k}$ analytically, using (Blei et al., 2003, Appendix A.1):

[TABLE]

Similarly, we evaluate the following expectations analytically:

[TABLE]

Let $\bm{\theta}\sim\operatorname{Dir}_{L}({\bm{\rho}})$, $\mathbb{E}[\theta_{k}]=\rho_{k}/\rho_{.}$, and $\alpha\in[0,\infty)$. We can expand the intractable expectation $\mathbb{E}[\log\Gamma(\alpha\theta_{k})]$ as (Kim et al., 2013, Theorem 3.1)

[TABLE]

We then write the individual expectations as follows:

[TABLE]

For $\mathbb{E}_{q_{jd}}[\log p_{\gamma}(\bm{\theta}_{jd}\,|\,\bm{\pi}_{jd})]$, the second step uses (65) and the independence assumption of the variational distribution, $\mathbb{E}_{q}[\pi_{jk}\log\theta_{jdk}]=\mathbb{E}_{q}[\pi_{jk}]\mathbb{E}_{q}[\log\theta_{jdk}]$, and the third step uses the result $\mathbb{E}_{q}[\pi_{jk}]=\tau_{jk}/\tau_{j.}$.

We have the variational parameters $({\bm{\lambda}},{\bm{\tau}},{\bm{\rho}},\bm{\phi})$ for the latent variables $(\bm{\beta},\bm{\pi},\bm{\theta},{\bm{z}})$ in the cLDA model. We will see in the following subsections that the updates for the variational parameters, except for ${\bm{\tau}}_{j}$’s, follow closely that of the variational Dirichlet and Multinomial updates of the LDA model (Blei et al., 2003, Appendix).

E.1 Variational Dirichlet Update for Collections

Grouping the expectations that contain ${\bm{\tau}}_{j}$ from the lower-bound ${\cal L}$666We ignore the arguments of ${\cal L}$ to ease notation., we get:

[TABLE]

This lower-bound does not produce a closed form expression for updating $\tau_{j}$’s. Kim et al. (2013) suggested to use a Newton’s update with equality constraints for update in a similar hierarchical modeling context. A similar approach is followed here. We first break down each Dirichlet parameter ${\bm{\tau}}_{j}$ into a scale parameter $a_{j}$ and a base measure ${\bm{\omega}}_{j}$ that satisfies the equality constraint $\sum_{k=1}^{K}\omega_{jk}=1$. The corresponding variational distribution for $\bm{\pi}_{j}$ is redefined as $\bm{\pi}_{j}\sim\operatorname{Dir}_{K}(a_{j}{\bm{\omega}}_{j})$. We will see that this decomposition will enable us to perform a Newton’s update with equality constraints. Utilizing the equality constraint and collecting terms that contain $a_{j}$ and ${\bm{\omega}}_{j}$, we get

[TABLE]

To maximize ${\cal L}_{[a_{j}{\bm{\omega}}_{j}]}$ with respect to $\omega_{jk}$, we first form the objective function for $\omega_{jk}$ by collecting terms as

[TABLE]

Its first and second derivatives, denoted by $g_{jk}$ and $h_{jk}$, are:

[TABLE]

We can see that the hessian given by (70) is diagonal. We use $u$ to denote the dual variable for the sums to one constraint. We then form the constraint Newton step $\Delta\omega_{jk}$ by solving the set of linear equations

[TABLE]

that yields

[TABLE]

By construction, this update satisfies $\sum_{k=1}^{K}\Delta\omega_{jk}=0$ and preserves the sums to one constraint of the variational parameter $\omega_{jk}$ (Kim et al., 2013).

To maximize ${\cal L}_{[a_{j}{\bm{\omega}}_{j}]}$ with respect to $a_{j}$, we first form the objective function for $a_{j}$, by collecting terms as

[TABLE]

Its first and second derivatives, denoted by $g_{j}$ and $h_{j}$, are:

[TABLE]

The Newton update step for $a_{j}$ is then given by

[TABLE]

We alternately maximize ${\cal L}_{[a_{j}{\bm{\omega}}_{j}]}$ with respect to $a_{j}$ and ${\bm{\omega}}_{j}$ until convergence.

E.2 Variational Multinomial Update for Words

We derive the expression for updating the variational parameter $\phi_{jdik}$—the probability that $jdi^{\text{th}}$ word is generated by topic $k$—via maximizing the lower-bound ${\cal L}$ with respect to the constraint $\sum_{k=1}^{K}\phi_{jdik}=1$. We first form the Lagrangian by collecting the terms that contain $\phi_{jdik}$ from (66), and applying Lagrangian multipliers as

[TABLE]

Taking derivatives with respect to $\phi_{jdik}$, we get

[TABLE]

Setting this to zero yields the maximum value for $\phi_{jdik}$

[TABLE]

E.3 Variational Dirichlet Update for Documents

We maximize the lower-bound ${\cal L}$ with respect to $\rho_{jdk}$. Collecting terms that contain $\rho_{jd}$ from (66), we get

[TABLE]

Taking derivatives with respect to $\rho_{jdk}$, we get

[TABLE]

Setting this to zero yields the maximum value for $\rho_{jdk}$

[TABLE]

E.4 Variational Dirichlet Update for Topics

We maximize the lower-bound ${\cal L}$ with respect to $\lambda_{kv}$. Collecting terms that contain ${\bm{\lambda}}_{k}$ from (66), we get

[TABLE]

Taking derivatives with respect to $\lambda_{kv}$, we get

[TABLE]

Setting this to zero yields the maximum value for $\lambda_{kv}$

[TABLE]

E.5 Optimize Hyperparameters

As in the variational EM algorithm of LDA (Blei et al., 2003), the E-step of the cLDA VEM algorithm updates the variational parameters based on the expressions provided above (see Algorithm 3). We can use the optimal lower-bound ${\cal L}(q^{*},p_{h})$ as the tractable approximation for the log marginal likelihood $\log m(h)$. In the M-step of VEM, we can then update the hyperparameters $h=(\alpha,\gamma,\eta)$ by maximizing the optimal lower-bound with respect to $h$. We collect terms that contain each hyperparameter, separately, as follows:

[TABLE]

The first and second derivatives are the following:

[TABLE]

Using these derivatives, one can update $\alpha$, $\gamma$, $\eta$ via a Newton method (Blei et al., 2003; Minka, 2000b).

Appendix F Comparison of AGS, MGS, and VEM on a Synthetic Corpus

This section gives additional details for the comparative study given in Section 3.2 of the main paper. {notes-to-stay} Here we compare the performance of three approximate inference methods, AGS, MGS, and VEM, for the cLDA model. For this purpose, we created a synthetic corpus by simulating lines $1$-$5$ of the cLDA hierarchical model with the following configurations. We took the number of collections $J=2$ and the number of topics $K=3$. We did this solely so that we can visualize the results of the algorithms. Besides, we took the vocabulary size $V=40$, the number of documents in each collection $D_{j}=100$, and the hyperparameters $h_{\text{true}}=(\alpha,\gamma,\eta)=(.1,1,.25)$. The collection-level Dirichlet sampling with hyperparameter $\alpha$ (via line $2$ of the hierarchical model) produced two topic distributions $\bm{\pi}^{\text{true}}_{1}=(.002,.997,\epsilon)$ and $\bm{\pi}^{\text{true}}_{2}=(.584,.030,.386)$, where $\epsilon$ denotes a small number. The data ${\bm{w}}$ is generated by simulating line $5$.

We first study the ability of these algorithms to recover parameters $\bm{\pi}^{\text{true}}_{1}$ and $\bm{\pi}^{\text{true}}_{2}$. We do this by comparing samples of $\bm{\pi}_{j}$ from the AGS and MGS chains on $(\bm{\pi},{\bm{z}})$ with variational estimates of $\bm{\pi}_{j}$ from the VEM algorithm iterations777An implementation of all algorithms and datasets discussed in this paper is available as an R package at https://github.com/clintpgeorge/clda. We initialized $\bm{\pi}^{(0)}_{1}=\bm{\pi}^{(0)}_{2}=(.33,.33,.33)$ for all three algorithms. Using the data ${\bm{w}}$, we ran both AGS and MGS chains for $2000$ iterations, and the VEM algorithm converged after $45$ EM iterations.

Figure 13 gives trace plots of values of $\bm{\pi}_{1}$ and $\bm{\pi}_{2}$ on the $2$-simplex for all three algorithms. Consider any of the choices of hyperparameters $h$ and the number of topics $K$. Markov chains AGS and MGS induces a chain on $\bm{\psi}$ with invariant distribution $p_{h,{\bm{w}}}(\bm{\psi})$ by using the conditional distribution of $(\bm{\beta},\bm{\theta})$, given by (8) (see Section 3). They essentially give us a sequence $(\bm{\beta}^{(1)},\bm{\pi}^{(1)},\bm{\theta}^{(1)},{\bm{z}}^{(1)}),\ldots,(\bm{\beta}^{(S)},\bm{\pi}^{(S)},\bm{\theta}^{(S)},{\bm{z}}^{(S)})$. Consider the component $\bm{\pi}_{j}^{(s)}$ of $\bm{\pi}^{(s)}$. While both $\bm{\pi}_{j}^{(s)}$ and $\bm{\pi}_{j}^{\text{true}}$ are points in the $K$-$1$ simplex, their interpretations are different: $\bm{\pi}_{j}^{(s)}$ is a distribution on the $K$ topics $\beta_{1}^{(s)},\ldots,\beta_{K}^{(s)}$, while $\bm{\pi}_{j}^{\text{true}}$ is a distribution on the $K$ topics $\beta_{1}^{\text{true}},\ldots,\beta_{K}^{\text{true}}$, and these are different sets of topics. This is also the case with $K$ components of $\bm{\theta}_{d}^{(s)}$ and $\bm{\theta}_{d}^{\text{true}}$ (see, e.g., Griffiths and Steyvers (2004); George (2015)). The results of VEM described in the paper hold a similar case. Thus, before making any comparison between the values of $\bm{\pi}_{j}^{(s)}$ and $\bm{\pi}_{j}^{\text{true}}$, one needs to align the corresponding sets of topics. For simple corpora such as the one used in our empirical evaluation, one can trivially re-align topics to compare $\bm{\pi}_{j}^{\text{true}}$ and $\bm{\pi}_{2}^{(s)}$s (see Table 1).

{notes-to-stay}

From the plots, we can see that both AGS and MGS chains were able to recover the values $\bm{\pi}^{\text{true}}_{1}$ and $\bm{\pi}^{\text{true}}_{2}$ reasonably well, even though the AGS chain has an edge. Converging to optimal regions $\|\bm{\pi}^{\text{true}}_{1}-\bm{\pi}^{(s)}_{1}\|<10^{-2}$ and $\|\bm{\pi}^{\text{true}}_{2}-\bm{\pi}^{(s)}_{2}\|<10^{-2}$, the AGS chain took cycles $42$ and $23$ only; but, the MGS chain required cycles $248$ and $139$, respectively. ($\|.\|$ denotes L-$1$ norm on the simplex $\mathbb{S}_{K}$.) Lastly, the VEM algorithm started off nicely, but never reached the optimal regions; at convergence, VEM hit points that are $0.1432$ far from $\bm{\pi}^{\text{true}}_{1}$ and $0.2507$ far from $\bm{\pi}^{\text{true}}_{2}$.

We now compare the mixing rates of the two chains AGS and MGS. People often use diagnostics such as trace plots and auto-correlation function (ACF) plots for this purpose. Although both chains appear to converge in reasonable cycles, the AGS chain mixes faster as shown in Figure 13 and 3. Figure 14 further supports this fact: it gives plots of ACF’s for two random elements $\pi_{11}$ and $\pi_{22}$ of the $\bm{\pi}$ matrix, from each iteration of the chains AGS and MGS. These plots suggests that the AGS chain mixes faster, iterations separated by a lag of $25$ or $30$ are essentially uncorrelated. Note that even though the VEM algorithm did not reach the optimal regions here, in our experience, it converges relatively quickly. For modeling corpora with large document collections, we still recommend the reader to use the VEM algorithm as a practical alternative, considering its speed gains and parallelization capabilities. On the other hand, AGS gives more accurate results; hence, we use AGS for our future analysis.

Note that sampling ${\bm{z}}$’s in both AGS and MGS chains is quite similar to sampling ${\bm{z}}$’s in the LDA CGS (Griffiths and Steyvers, 2004) chain. The CGS chain is a well-studied chain, see, e.g., George (2015); so, we do not report diagnostics for samples of ${\bm{z}}$ from the chains AGS and MGS here.

Appendix G Perplexity Calculation for cLDA and LDA

In this section, we derive expressions for Perplexity (defined in Section 4.2) for both cLDA and LDA models. To compute the perplexity score (21), we first need an expression for the predictive likelihood $p(w_{jdi}\,|\,{\bm{w}}^{\text{train}})$. We obtain this by marginalizing the likelihood over $\bm{\beta},\bm{\pi},\bm{\theta}$ and $z_{jdi}$. There is a slight abuse of notation here: the variables $\bm{\beta},\bm{\pi},\bm{\theta}$ are all dependent on the training data ${\bm{w}}^{\text{train}}$ only, but we ignore this in the notation. From the hierarchical model, we can write the predictive likelihood for word $w_{jdi}$ as

[TABLE]

where we used the property that $\prod_{k=1}^{K}\prod_{v=1}^{V}\beta_{kv}^{z_{jdik}w_{jdiv}}=\prod_{v=1}^{V}\beta_{k^{\prime}v}^{w_{jdiv}}$ for some $k^{\prime}$ such that $z_{jdik^{\prime}}=1$. The dependence of the likelihood on $\bm{\pi}_{j}$’s is ignored in the notation as it’s included in $\bm{\theta}_{jd}$’s. Since $z_{jdi}\sim\operatorname{Mult}_{K}(\bm{\theta}_{jd})$, $E(z_{jdi})=\theta_{jdk}$, for $j=1,\ldots,J$, $d=1,\ldots,D_{j}$, $i=1,\ldots,n_{jd}$ and $k=1,\ldots,K$. We then have

[TABLE]

Let $(\bm{\beta}^{[1]},\bm{\theta}^{[1]}),\ldots,(\bm{\beta}^{[S]},\bm{\theta}^{[S]})$ be a Markov chain with invariant distribution $\nu_{h,{\bm{w}}^{\text{train}}}(\bm{\beta},\bm{\theta})$. One can then estimate the marginal likelihood $p(w_{jdi}\,|\,{\bm{w}}^{\text{train}})$ in (21) by the Monte Carlo average

[TABLE]

We can use any of the augmented chains on $(\bm{\beta},\bm{\pi},\bm{\theta},{\bm{z}})$ described in Section 3.1 to compute this average. An elegant alternative is to substitute the following estimates of $\theta^{[s]}_{jdk}$ and ${\beta^{[s]}}_{kv}$ in (97)

[TABLE]

which can be computed for every sample in a chain on $(\bm{\pi},{\bm{z}})$. The resulting likelihood estimate dominates the original estimate (97) in terms of variance. This approach is sometimes called as Rao-Blackwellization, see, e.g., (Robert and Casella, 2005, Chapter 4). To compute these estimates, we only need a Markov chain on $(\bm{\pi},{\bm{z}})$, which has reduced computational cost compared to a Markov chain on $(\bm{\beta},\bm{\pi},\bm{\theta},{\bm{z}})$. Note that one can plug in variational estimates of $\theta_{jdk}$ and $\beta_{kv}$ via the cLDA variational EM algorithm (see Section E) into (96) to estimate the marginal likelihood $p(w_{jdi}\,|\,{\bm{w}}^{\text{train}})$.

Using similar arguments, we can derive an expression for the estimate of the marginal likelihood $p(w_{di}\,|\,{\bm{w}}^{\text{train}})$ for the LDA model. Given the CGS chain $z^{[1]},\ldots,z^{[S]}$ (Griffiths and Steyvers, 2004), we have

[TABLE]

where

[TABLE]

Appendix H Additional Experimental Results

This section gives additional results of experiments on real world corpora that are discussed in Section 4.

Acknowledgments

This work is supported by grants NIH #$7$ R$21$ GM$101719$-$03$ to George Michailidis and Institute of Education Sciences #R$305160004$, and the University of Florida Informatics Institute.