\section*{Appendix B: MCMC via Auxiliary Variable and Gibbs Sampling}
\label{app:aux}

\section[Estimating hyperparameters eta and gamma]
{Estimating hyperparameters $\eta$ and $\gamma$}
\label{sec:est-eta-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 \cite[Appendix 
B]{minka2000estimating} in our development. 

\begin{eqnarray}
\frac{\Gamma (x)}{\Gamma(n + x)} &\geq& 
\frac{\Gamma (\hat{x}) \exp\left((\hat{x} - x)b\right)}
{\Gamma(n + \hat{x})} \label{eq:minka-bound1}\\ 
b &=& \Psi(n + \hat{x}) - \Psi(\hat{x}) \nonumber \\[1em] 
\frac{\Gamma(n + x)}{\Gamma (x)} &\geq& cx^a \: \:
\text{ if } n \geq 1, x \geq 1 \label{eq:minka-bound2}\\ 
a &=& \left(\Psi(n + \hat{x}) - \Psi(\hat{x})\right) \hat{x} 
\nonumber\\
c &=&  \frac{\Gamma(n + \hat{x})}{\Gamma (\hat{x})} \hat{x}^{-a}
\nonumber
\end{eqnarray}

From the cLDA hierarchical model 
\eqref{eq:cLDA-e}--\eqref{eq:cLDA-a}, 
after integrating out $\bm \beta$'s, we get the marginal posterior of 
$(\bm \pi, {\bm z})$, given ${\bm w}$ and $\eta$ as:  
\begin{equation}
p_{\eta, {\bm w}}(\bm \pi, {\bm z}) \propto \prod_{k = 1}^K 
\left [ \frac{\Gamma(V\eta)}{\Gamma(m_k + V\eta)} 
\prod_{v = 1}^V \frac{\Gamma(m_{kv} + \eta)} {\Gamma(\eta)}
\right]
\end{equation}
Using \eqref{eq:minka-bound1} and \eqref{eq:minka-bound2}, we write 
it as: 
\begin{eqnarray}
p_{\eta, {\bm w}}(\bm \pi, {\bm z}) &\geq& \prod_{k = 1}^K 
\left [ 
\frac{\Gamma(V\eta_0) \exp\left(V(\eta_0 - \eta)b_k\right)}
{\Gamma(m_k + V\eta_0)} 
\prod_{v = 1}^V c_{kv} \eta^{a_{kv}}
\right] \\ 
a_{kv} &=& \eta_0 \Big( \Psi(m_{kv} + \eta_0) - \Psi (\eta_0) 
\Big) \\ 
b_k &=& \Psi(m_k + V\eta_0) - \Psi(V\eta_0)\\ 
c_{kv} &=& \frac{\Gamma(m_{kv} + \eta_0)}{\Gamma(\eta_0)} 
\eta_0^{-a_{kv}} 
\end{eqnarray}
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   
\begin{equation}
{\cal L}_\eta = \sum_{k=1}^K \sum_{v = 1}^V 
\Big( \log c_{kv} + a_{kv} \log \eta \Big) + 
\sum_{k=1}^K \Big( \log \Gamma (V\eta_0) - 
\log \Gamma(m_k + V\eta_0) + V (\eta_0 - \eta) b_k \Big).
\end{equation}
Taking derivatives with respect to $\eta$, 
\begin{equation}
\frac{\partial {\cal L}_\eta}{\partial \eta} = 
\frac{1}{\eta} \sum_{k=1}^K \sum_{v = 1}^V a_{kv} - 
V \sum_{k=1}^K b_k 
\end{equation}
and setting it to zero, we get 
\begin{equation}
\eta = \frac{\sum_{k=1}^K \sum_{v = 1}^V a_{kv}}{V\sum_{k=1}^K b_k}
= \frac{\eta_0}{V} 
\frac{\sum_{k=1}^K \sum_{v = 1}^V 
	\Big[ \Psi(m_{kv} + \eta_0) - \Psi(\eta_0)\Big]} 
{\sum_{k=1}^K \Big[ \Psi(m_k + V\eta_0) - \Psi(V\eta_0) \Big]}
\end{equation}
We can compute the maximum via the fixed point iteration 
\citep{Minka:2000}. 

Similarly, from the cLDA hierarchical model  
\eqref{eq:cLDA-e}--\eqref{eq:cLDA-a}, 
after integrating out $\bm \theta$'s, we get the marginal posterior of 
$(\bm \pi, {\bm z})$, given ${\bm w}$ and $\gamma$ as:  
\begin{equation}
p_{\gamma, {\bm w}}(\bm \pi, {\bm z}) \propto \prod_{j = 1}^J \prod_{d = 
1}^{D_j}   
\left [ \frac{\Gamma(\gamma)}{\Gamma(n_{jd} + \gamma)} 
\prod_{k = 1}^K 
\frac{\Gamma(n_{jdk} + \gamma \pi_{jk})} {\Gamma(\gamma \pi_{jk})}
\right]
\end{equation}
Using \eqref{eq:minka-bound1} and \eqref{eq:minka-bound2}, we  
write it as: 
\begin{eqnarray}
p_{\gamma, {\bm w}}(\bm \pi, {\bm z}) &\geq& \prod_{j = 1}^{J} \prod_{d = 
1}^{D_j} 
\left [ 
\frac{\Gamma(\gamma_0) \exp\left((\gamma_0 - \gamma)b_{jd}\right)}
{\Gamma(n_{jd} + \gamma_0)} 
\prod_{k = 1}^K c_{jdk} (\gamma \pi_{jk})^{a_{jdk}}
\right] \\ 
a_{jdk} &=& \gamma_0 \pi_{jk} \Big( \Psi(n_{jdk} + \gamma_0 
\pi_{jk}) - \Psi (\gamma_0 \pi_{jk}) \Big) \\ 
b_k &=& \Psi(n_{jdk} + \gamma_0) - \Psi(\gamma_0) \\ 
c_{jdk} &=& \frac{\Gamma(n_{jdk} + \gamma_0 \pi_{jk})}
{\Gamma(\gamma_0 \pi_{jk})}{(\gamma_0 \pi_{jk})}^{-a_{jdk}} 
\end{eqnarray}
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   
\begin{eqnarray}
{\cal L}_\gamma &=& \sum_{j = 1}^J \sum_{d = 1}^{D_j} \sum_{k = 1}^K
\Big[ \log c_{jdk} + a_{jdk} \log \gamma + a_{jdk} \log \pi_{jk} 
\Big] + \\ 
&\hspace{1cm}& \sum_{j = 1}^J \sum_{d = 1}^{D_j} 
\Big[ \log \Gamma (\gamma_0) - \log \Gamma(n_{jd} + \gamma_0) 
+ (\gamma_0 - \gamma) b_{jd} \Big] \nonumber.
\end{eqnarray}
Taking derivatives with respect to $\gamma$, 
\begin{equation}
\frac{\partial {\cal L}_\gamma}{\partial \gamma} = 
-\sum_{j = 1}^{J} \sum_{d = 1}^{D_j} b_{jd} + 
\frac{1}{\gamma} \sum_{j = 1}^{J} \sum_{d = 1}^{D_j} 
\sum_{k = 1}^{K} a_{jdk} 
\end{equation} 
and setting it to zero, we get 
\begin{eqnarray}
\gamma &=& \frac{\sum_{j = 1}^{J} \sum_{d = 1}^{D_j} \sum_{k = 
1}^{K} a_{jdk}}{\sum_{j = 1}^{J} \sum_{d = 1}^{D_j} b_{jd}} 
\nonumber \\
&=& \gamma_0 \frac{\sum_{j = 1}^{J} \sum_{d = 1}^{D_j} 
\sum_{k = 1}^{K} \pi_{jk} \Big[ 
\Psi (n_{jdk} + \gamma_0 \pi_{jk}) - \Psi (\gamma_0 \pi_{jk}) \Big]}
{\sum_{j = 1}^{J} \sum_{d = 1}^{D_j} \Big[ \Psi (n_{jd} + \gamma_0) 
- \Psi (\gamma_0) \Big]} 
\end{eqnarray}
We can compute the maximum via the fixed point iteration 
\citep{Minka:2000}. 

\section{Langevin Monte Carlo}
\label{app:mmala}

The Metropolis Adjusted Langevin Algorithm 
~\cite[MALA]{Girolami:2011}, as described in Section~\ref{sec:mgs}, 
is given the following steps 
\begin{enumerate}
	\item Propose $\bm \pi^{(*)}$ via Langevin dynamics 
	\eqref{eq:langevin-discrete-eq}
	\item Calculate the MH acceptance ratio 
	\begin{equation}
	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)}
	\end{equation}
	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. 
\end{enumerate}

\eat{
In this section, we first develop the Metropolis Adjusted Langevin 
Algorithm (MALA)~\cite{Girolami:2011} to sample collection level 
random variables $\bm \pi$\footnote{We ignore the subscript $j$ and 
other dependencies on ${\bm z}$ and ${\bm w}$ to keep the notation 
uncluttered.} with density $p(\bm \pi)$. We then discuss a few 
shortcomings of MALA and propose potential improvements.   

Briefly, MALA induces a Markov chain with invariant density 
$p(\bm \pi)$ via the Metropolis-Hastings (MH) algorithm  
\cite{MetropolisEtal:1953, Hastings:1970} with Langevin updates. 
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 
\begin{equation}
a (\bm \pi^{(t)}, \bm \pi^{(*)}) = \text{min} \left(1, 
\frac{\tilde{p}(\bm \pi^{(*)}) q(\bm \pi^{(t)} \leftarrow \bm \pi^{(*)})} 
{\tilde{p}(\bm \pi^{(t)}) q(\bm \pi^{(*)} \leftarrow \bm \pi^{(t)}) } 
\right)
\label{eq:mh-accept-ration}
\nonumber 
\end{equation}
where $\tilde{p}(\bm \pi)$ denotes the unnormalized density 
of $\bm \pi$. The accept-reject step ensures that the proposed 
Markov chain is reversible with respect to the stationary target 
density and satisfies detailed balance. The goal of proposal 
distribution is to simulate 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~\cite{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^{(t)})$ is defined by the stochastic 
differential equation (SDE) 
\begin{equation}
\text{d}\bm \pi(t) = \frac{1}{2} \nabla_{\bm \pi} {\cal L}(\bm \pi^{(t)}) \text{d}t + \text{d}{\bm b}(t)   
\label{eq:langevin-sde}
\end{equation} 
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~\eqref{eq:langevin-sde} as
\begin{equation}
\bm \pi^{(*)} = {\bm \mu}(\bm \pi^{(t)}, \varepsilon) + \varepsilon {\bm \xi}^{(t)},
\label{eq:langevin-eq}
\end{equation}
where $\varepsilon$ is a user defined step-size for the 
discretization and 
\begin{equation}
{\bm \mu}(\bm \pi^{(t)}, \varepsilon) = \bm \pi^{(t)} + 
\frac{\varepsilon^2}{2} \nabla_{\bm \pi} {\cal L}(\bm \pi^{(t)}), \:\:
{\bm \xi}^{(t)} \stackrel{\text{iid}}{\sim} \mathcal{N}_K(0, \mathds{1}_K).
\end{equation}
The random variable ${\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~\eqref{eq:langevin-eq} may be 
transient and is no longer reversible with respect to the 
stationary density $p(\bm \pi)$~\cite{RobertsTweedie:1996}. 
However, one can ensure convergence to
the target density $p(\bm \pi)$ via a MH accept-reject scheme 
\cite[p. 591]{Besag:1994}, which we denote by MALA, as follows. 
\begin{enumerate}
	\item Propose $\bm \pi^{(*)}$ via Langevin dynamics 
	\eqref{eq:langevin-eq}
	\item Calculate the MH acceptance ratio 



					
		
					

	\begin{equation}
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)}
\end{equation}
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. 
\end{enumerate}
 

\begin{notes-to-stay}
To run this update, we need a closed form expression for the  
gradient $\nabla_{\bm \pi} \log p(\bm \pi_j \,|\, {\bm z}_j, {\bm w}_j) $. Let 
$\varPsi (x)$ be $\frac{\partial \log\Gamma (x)}{\partial x}$, 
i.e., the digamma function. From 
\eqref{eq:posterior-pi}, for $k = 1, \ldots, K$, we can write 
\begin{eqnarray}
\frac{\partial }{\partial \pi_{jk}} \log p(\bm \pi_j \,|\, {\bm z}_j, {\bm w}_j) 
& = & \sum_{d=1}^{D_j} \gamma \, \Bigl [
\varPsi(\gamma \pi_{jk} + n_{jdk}) - 
\varPsi(\gamma \pi_{jk}) \Bigr ] + 
\frac{1}{\pi_{jk}} \Bigl (\alpha - 1 \Bigr ),
\end{eqnarray} 
which can be easily computed by any implementation of 
\textsl{Polygamma} functions. 
\end{notes-to-stay}










One may notice this MALA update as a special case of 
Hamiltonian Monte Carlo~\cite[Section 5.5.2]{Neal:2010}. However, 
the proposed MALA update for $\bm \pi^{(t+1)}$ has some 
shortcomings. First, $\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. Second, typical Dirichlet priors over the probability 
simplex put most of their probability mass on the edges and corners 
of the simplex---e.g. in linguistic models such LDA 
\cite{Patterson:2013}. Computing gradients for these models become 
unstable when the probabilities are close to zero and causes issues 
for MALA updates. Third, 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~\cite{Girolami:2011}. 
\begin{notes-to-stay}
Note: Since we use an unnormalized quantity given by 
\eqref{eq:posterior-pi} to compute gradients, large collections may 
give very large gradients, and vice versa. One natural option is to 
normalize gradients during the implementation.  
\end{notes-to-stay}

}

\subsection{Riemann Manifold Metropolis Adjusted Langevin 
Algorithm}


Girolami et al. \cite[Section 5]{Girolami:2011} 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 
\begin{equation}
\text{d}\bm \pi(t) = \frac{1}{2} \widetilde{\nabla}_{\bm \pi} {\cal L} ( \bm \pi^{(t)} ) \text{d}t + \text{d}\tilde{{\bm b}}(t)   
\label{eq:reimann-langevin}
\end{equation} 
where the natural gradient is 
\begin{equation}
\widetilde{\nabla}_{\bm \pi} {\cal L}( \bm \pi^ {(t)} ) = {\bm G}\{ \bm \pi (t) \}^{-1} \nabla_{\bm \pi} {\cal L}(\bm \pi^{(t)})
\end{equation} 
and the Brownian motion is 
\begin{equation}
\text{d}\tilde{{\bm b}}_k(t) = {| {\bm G}\{ \bm \pi (t) \} |}^{-1/2} \sum_{k^\prime = 1}^{K} \frac{\partial }{\partial \pi_{k^\prime}} \Bigl [ {\bm G}\{ \bm \pi (t) \}^{-1}_{kk^\prime} \, {| {\bm G}\{ \bm \pi (t) \} |}^{1/2} \Bigl ] \text{d}t 
+ {\Bigl [ {{\bm G}\{ \bm \pi (t) \}}^{-1/2} \text{d}{\bm b}(t)  \Bigl ]}_k
\label{eq:reimann-langevin-brownian}
\end{equation} 
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 \eqref{eq:reimann-langevin} will be 
reduced to the standard Langevin SDE.  

By expanding the gradients in~\eqref{eq:reimann-langevin-brownian} 
and discretizing \eqref{eq:reimann-langevin} via the 
first-order Euler integration, we get the proposal for 
the Riemann Manifold Metropolis Adjusted Langevin Algorithm (MMALA):
\begin{equation}
\bm \pi^{(*)} = {\bm \mu}(\bm \pi^{(t)}, \varepsilon) + \varepsilon  {{\bm G}(\bm \pi^{(t)})}^{-1/2} {\bm \xi}^{(t)},
\label{eq:riemann-langevin-euler}
\end{equation}
where $k^\text{th}$ element in vector ${\bm \mu}(\bm \pi^{(t)}, \varepsilon)$ is 
given by 
\begin{equation}
\begin{split}
{\bm \mu}(\bm \pi^{(t)}, \varepsilon)_k &= \pi_k^{(t)} 
+ \frac{\varepsilon^2}{2} \Bigl \{ {{\bm G} ( \bm \pi^{(t)} )}^{-1} 
\nabla_{\bm \pi} {\cal L}(\bm \pi^{(t)}) \Bigr \}_k \\ 
&\hspace{5mm}- \varepsilon^2 \sum_{k^\prime = 1}^{K} \Bigl \{ 
{\bm G}(\bm \pi^{(t)} )^{-1} \frac{\partial {\bm G}(\bm \pi^{(t)} )}{\partial 
\pi_{k^\prime}} {\bm G}(\bm \pi^{(t)} )^{-1} \Bigr \}_{kk^\prime} \\ 
&\hspace{5mm}+ \frac{\varepsilon^2}{2} \sum_{k^\prime = 1}^{K} \bigl \{{\bm G}(\bm \pi^{(t)} )^{-1} \bigl \}_{kk^\prime} \text{tr} \Bigl \{ {\bm G}(\bm \pi^{(t)} )^{-1} \frac{\partial {\bm G}(\bm \pi^{(t)} )}{\partial \pi_{k^\prime}}  \Bigr \} 
\end{split}
\end{equation}
The corresponding proposal density is given by 
\begin{equation}
q(\bm \pi^{(*)} \leftarrow \bm \pi^{(t)}) = 
\mathcal{N}_K(\bm \pi^{(*)} \,|\, {\bm \mu}(\bm \pi^{(t)}, \varepsilon), \, 
\varepsilon^2 {\bm G}(\bm \pi^{(t)} )^{-1})
\end{equation}

\subsection{Langevin Updates on Probability Simplices: Boundary 
Considerations}

\eat{
To perform valid moves (e.g. via MALA or MMALA) of $\bm \pi$ that 
lies on the probability simplex, one needs to consider boundary 
conditions. One natural solution to handle boundaries is to 
re-parameterize $\bm \pi$~\cite{Patterson:2013}. Here, we revisit the 
original parameter $\bm \pi_j \in \mathbb{S}_K$, $j = 1, 2, \ldots, J$ 
of interest and the unnormalized conditional density 
\eqref{eq:posterior-pi}.  Let ${\bm \varphi}_j 
= (\varphi_{j1}, \ldots, \varphi_{jK}) \in \mathbb{R}^K$. We take the 
prior on ${\bm \varphi}_j$ as a product of i.i.d. Gamma random variables 
as 
\begin{equation}
\begin{split}
|\varphi_{jk}| &\stackrel{\text{iid}}{\sim} \text{Gamma} \Bigl (\alpha, 1 \Bigr ), \\
p({\bm \varphi}_j) &\propto \prod_{k = 1}^{K} {|\varphi_{jk}|}^{\alpha - 1} e^{-|\varphi_{jk}|}. 
\end{split}
\end{equation}
We define $|\varphi_{j.}| := \sum_{k=1}^K |\varphi_{jk}|$. Let  
$\pi_{jk}$ be $|\varphi_{jk}| / |\varphi_{j.}|$, for each $k = 1, 
2, \ldots, K$. This choice keeps the prior on $\bm \pi_j$ a    
Dirichlet density. We can now re-write the unnormalized conditional 
density~\eqref{eq:posterior-pi} on $\bm \pi_j$ as 
}

We first re-write the unnormalized conditional 
density~\eqref{eq:posterior-pi} on $\bm \pi_j$ as 
\begin{equation}
\tilde{p}({\bm \varphi}_j \,|\, {\bm z}_j, {\bm w}_j) \propto  
\prod_{d=1}^{D_j} \prod_{k=1}^K
\Biggl(\frac{ \Gamma(\gamma \frac{|\varphi_{jk}|}{|\varphi_{j.}|} + n_{jdk})}
						{\Gamma(\gamma \frac{|\varphi_{jk}|}{|\varphi_{j.}|})}
\Biggr)
\prod_{k=1}^K {|\varphi_{jk}|}^{\alpha - 1} e^{-|\varphi_{jk}|}
\label{eq:posterior-pi-re-parametrization}
\end{equation}
The MMALA update for the new parametrization is given by 
\begin{equation}
{\bm \varphi}_j^{(*)} = {\bm \mu}({\bm \varphi}_j^{(t)}, \varepsilon) + \varepsilon  {{\bm G}({\bm \varphi}_j^{(t)})}^{-1/2} {\bm \xi}^{(t)}
\label{eq:riemann-langevin-euler-reparam}
\end{equation}
where for $k = 1, \ldots, K$, we have  
\begin{equation}
\begin{split}
{\bm \mu}({\bm \varphi}_j^{(t)}, \varepsilon)_k &= \varphi_{jk}^{(t)}  
+ \frac{\varepsilon^2}{2} \Bigl \{ {{\bm G} ( {\bm \varphi}_j^{(t)} )}^{-1} \nabla_{{\bm \varphi}_j} {\cal L}({\bm \varphi}_j^{(t)}) \Bigr \}_k \\ 
&\hspace{5mm}- \varepsilon^2 \sum_{k^\prime = 1}^{K} \Bigl \{ {\bm G}({\bm \varphi}_j^{(t)} )^{-1} \frac{\partial {\bm G}({\bm \varphi}_j^{(t)} )}{\partial \varphi_{jk^\prime}} {\bm G}({\bm \varphi}_j^{(t)} )^{-1} \Bigr \}_{kk^\prime} \\ 
&\hspace{5mm}+ \frac{\varepsilon^2}{2} \sum_{k^\prime = 1}^{K} \bigl \{{\bm G}({\bm \varphi}_j^{(t)} )^{-1} \bigl \}_{kk^\prime} \text{tr} \Bigl \{ {\bm G}({\bm \varphi}_j^{(t)} )^{-1} \frac{\partial {\bm G}({\bm \varphi}_j^{(t)} )}{\partial \varphi_{jk^\prime}}  \Bigr \}
\end{split} 
\end{equation}
The proposal density of this diffusion process is given by 
\begin{equation}
\begin{split}
q({\bm \varphi}_j^{(*)} \leftarrow {\bm \varphi}_j^{(t)}) &= 
\mathcal{N}_K({\bm \varphi}_j^{(*)} \,|\, {\bm \mu}({\bm \varphi}_j^{(t)}, \varepsilon), \, \varepsilon^2 {\bm G}({\bm \varphi}_j^{(t)} )^{-1}) \\ 
&\propto {| {\bm G}({\bm \varphi}_j^{(t)})^{-1} |}^{-1/2} 
\exp \Bigl \{ -\frac{1}{2 \varepsilon^2} {({\bm \varphi}_j^{(*)} - {\bm \mu}({\bm \varphi}_j^{(t)}, \varepsilon))}^\text{T} 
{\bm G}({\bm \varphi}_j^{(t)}) {({\bm \varphi}_j^{(*)} - {\bm \mu}({\bm \varphi}_j^{(t)}, \varepsilon))}
\Bigr\} \\
\end{split}
\label{eq:pi-transition-density-re-parameterization}
\end{equation}
We take the metric tensor ${\bm G} ( {\bm \varphi}^{(t)}_j )$ 
as $\diag(|{\bm \varphi}^{(t)}_j|)^{-1}$ 
(\cite{Patterson:2013} suggested this choice in a different 
context.), as it gives simplified expressions for  
\begin{equation}
{\bm \mu}({\bm \varphi}_j^{(t)}, \varepsilon)_k = \varphi_{jk}^{(t)} 
+ \frac{\varepsilon^2}{2} \, \Bigl \{ \diag(|{\bm \varphi}_j^{(t)}|) \, \nabla_{{\bm \varphi}_j} {\cal L}({\bm \varphi}_j^{(t)}) \Bigr \}_k
+  \frac{\varepsilon^2}{2} \, \sign(\varphi^{(t)}_{jk}) 
\end{equation}
and
\begin{equation}
q({\bm \varphi}_j^{(*)} \leftarrow {\bm \varphi}_j^{(t)}) 
\propto \Bigl [ \prod_{k=1}^{K} {|\varphi_{jk}^{(t)}|}^{-1/2} \Bigr ]
\exp \Bigl \{ -\frac{1}{2 \varepsilon^2} {({\bm \varphi}_j^{(*)} - {\bm \mu}({\bm \varphi}_j^{(t)}, \varepsilon))}^\text{T} 
\diag({|{\bm \varphi}_j^{(t)}|})^{-1} {({\bm \varphi}_j^{(*)} - {\bm \mu}({\bm \varphi}_j^{(t)}, \varepsilon))}
\Bigr\} 
\end{equation}
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: 
\begin{equation}
\log \tilde{p}({\bm \varphi}_j \,|\, {\bm z}_j, {\bm w}_j) \propto  \sum_{k = 1}^{K} \Bigl [ \Bigl \{ (\alpha - 1) \log |\varphi_{jk}| \Bigr \} - |\varphi_{jk}| \Bigr ]
+ \sum_{d=1}^{D_j} \sum_{k=1}^K \Bigl [ \log\Gamma(\gamma \frac{|\varphi_{jk}|}{|\varphi_{j.}|} + n_{jdk}) - \log\Gamma(\gamma \frac{|\varphi_{jk}|}{|\varphi_{j.}|})  \Bigr ] \nonumber 
\end{equation}
We then define $\nabla_{{\bm \varphi}_j} {\cal L}({\bm \varphi}_j^{(t)})$ 
by the partial derivatives  
\begin{equation}
\begin{split}
\frac{\partial \log \tilde{p}({\bm \varphi}_j \,|\, {\bm z}_j, {\bm w}_j)}{\partial \varphi_{jk^*}}  
&= \sign (\varphi_{jk^*}) \Bigl \{ 
\frac{\gamma}{{|\varphi_{j.}|}^2}  \sum_{d=1}^{D_j} \sum_{k = 1}^K  {|\varphi_{jk}|} 
\Bigl [  \varPsi(\gamma \frac{{|\varphi_{jk}|}}{{|\varphi_{j.}|}}) - \varPsi (\gamma \frac{{|\varphi_{jk}|}}{{|\varphi_{j.}|}} + n_{jdk}) \Bigr ] \\
&\hspace{4mm}- \frac{\gamma}{{|\varphi_{j.}|}} \sum_{d=1}^{D_j} 
\Bigl [  \varPsi(\gamma \frac{{|\varphi_{jk^*}|}}{{|\varphi_{j.}|}}) - \varPsi (\gamma \frac{{|\varphi_{jk^*}|}}{{|\varphi_{j.}|}} + n_{jdk^*}) \Bigr ] 
+ \Bigl [ \frac{(\alpha - 1)}{{|\varphi_{jk^*}|}} - 1 \Bigr ] \Bigr \} \\
\end{split}
\end{equation}

\noindent
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 \ref{alg:gibbs-mmala-pi-z}. 





\section{Additional Experimental Results}
\label{app:nips-results}
This section gives additional results of experiments on real world  
corpora that are discussed in Section~\ref{sec:experiments}. 

\begin{figure}[htbp!]
	\centering
	\subfloat[cLDA topics: \texttt{hclust} dendrogram]
	{\includegraphics[width=.55\textwidth]
		{clda-topics-hclust-dendrogram}} 
	\subfloat[cLDA topics' size (log scale)]
	{\includegraphics[width=.44\textwidth]
		{clda-topics-word-assignments}} \\ 
	\subfloat[LDA topics: \texttt{hclust} dendrogram]
	{\includegraphics[width=.55\textwidth]
		{lda-topics-hclust-dendrogram}}
	\subfloat[LDA topics' size (log scale)]
	{\includegraphics[width=.44\textwidth]
		{lda-topics-word-assignments}} \\ 
	\subfloat[HDP topics: \texttt{hclust} dendrogram]
	{\includegraphics[width=.55\textwidth]
		{hdp-topics-hclust-dendrogram}}
	\subfloat[HDP topics' size (log scale)]
	{\includegraphics[width=.44\textwidth]
		{hdp-topics-word-assignments}} 
	\caption{Comparing topics learned via cLDA-AGS, LDA-CGS, and 
		HDP-CRF algorithms. Left-column shows \texttt{hclust} 
		dendrograms built 
		from topic-to-topic similarity matrices calculated based on 
		the 
		\emph{manhattan} distance for all three algorithms. For 
		each 
		method, right-column shows barplots of topic-sizes (i.e. 
		the 
		total number of 
		words assigned to a topic). Topics in the barplot $x$-axis  
		were ordered based on the order of topics in the 
		corresponding dendrogram leaf-nodes to ease comparison. 
		HDP found redundant set of topics and many topics with 
		too low topic-size---e.g. the children of the first/left 
		child in the dendrogram in Plot (e). See discussion in 
		Section \ref{sec:eval-performance}.        
	}
	\label{fig:topics-hclust-dendrogram-size}
\end{figure}

\begin{figure}[htbp!]
	\centering
	\includegraphics[width=.7\textwidth]
	{clda-lda-hdp-doc-hclust-silhouette}
	\caption{Boxplot statistics (median, lower hinge, and upper 
		hinge) of silhouette widths computed on \emph{documents}' 
		\texttt{hclust} clusters with various values of the number 
		of clusters, for corpus NIPS $00$-$18$. Algorithm 
		\texttt{hclust} was applied on documents' topic 
		distributions estimated via cLDA, LDA, and HDP sampling 
		algorithms.   
		Silhouette widths of clusters based on HDP are relatively 
		constant with different values of the number of clusters. 
		HDP estimated large set of minute topics ($K = 204$), which 
		may have helped clustering documents. cLDA performs better 
		than LDA, and is comparable with HDP or better than HDP 
		with the right number of \texttt{hclust} clusters (e.g. 
		from $50$ to $150$). See discussion in Section 
		\ref{sec:eval-performance}.  
	}
	\label{fig:clda-lda-hdp-doc-hclust-silhouette}
\end{figure}

\begin{figure}[htbp!]
	\centering
	\includegraphics[width=.7\textwidth]
		{lda-clda-hdp-topic-coherences}
	\caption{Estimated topic \emph{coherences} for models 
		LDA, cLDA, and HDP  for corpus NIPS $00$-$18$. See 
		discussion in Section \ref{sec:eval-performance}.   }
	\label{fig:lda-clda-hdp-topic-coherences} 
\end{figure}

\begin{figure}[t!] 
	\centering
	\includegraphics[width=.98\linewidth]
	{19nipsproceedings-clda-pi-samples}
	\caption{Estimates of topic distributions for four collections 
		(defined on timespans $1988$-$1992$, $1993$-$1997$, $1998$-$2002$, 
		and $2003$-$2005$) of the NIPS $00$-$18$ corpus via the cLDA AGS 
		algorithm.
	
	}
	\label{fig:nips-clda-pi-samples}
\end{figure}

\eat{
\begin{figure}[t!] 
\centering
\includegraphics[width=.7\linewidth]
{perp-19nipsproceedings-J4-K90-D2741-V9156-seed1983-e0d25-a0d5-g1-2016Jun18162135-clda-pi-samples}
\caption{Estimates of topic distributions for four collections 
(defined on timespans $1988$-$1992$, $1993$-$1997$, $1998$-$2002$, 
and $2003$-$2005$) of the NIPS $00$-$18$ corpus via the cLDA AGS 
algorithm. See each topic's description in 
Table~\ref{tab:nips-topic-words}.}
\label{fig:nips-clda-pi-samples-old}
\end{figure}
}
\eat{
\begin{table*}[t!] 
	\centering
	\begin{threeparttable}[b] 
	\caption{Most probable words from a $90$-topic cLDA 
	model for corpus NIPS $00$-$18$ }
\label{tab:nips-topic-words}
	\begin{tabular}{c c l}
		\toprule
		Topic ID\tnote{a} & 
		Description & 
		Most Probable Words\tnote{b} \\ 
		\midrule
		1 & neuro-science & cells, visual, receptive, spatial, 
		cortex, cortical \\ 
		3 & Bayesian inference & bayesian, prior, posterior, 
		distribution, 
		data, evidence \\ 
		6 & time series & time, sequence, state, 
		recurrent, temporal, hidden, markov \\
		8 & pattern recognition & object, feature, model, 
		recognition, image, parts, matching \\
		9 & image processing & image, pixel, color, 
		texture, shape, local, vision \\
		11 & neural networks & network, neural, neuron, 
		state, model, hopfield, equation \\
		15 & optimization & problem, solution, optimization, 
		constraints, linear, optimal, cost \\		
		19 & PCA & pca, data, analysis, 
		principal, component, feature, linear \\
		20 & mixture models & model, data, parameters, 
		mixture, em, likelihood, distribution \\
		22 & stop-words & set, algorithm, number, 
		used, problem, results, method \\
		25 & SVM & kernel, svm, vector, 
		support, training, data, space \\
		28 & neural networks & memory, patterns, capacity, 
		associative, stored, input, recall \\ 
		29 & stop-words & figure, function, value, case, rate, 
		number, results\\ 
		33 & information theory & information, entropy, mutual, 
		distribution, probability, fisher, log \\ 
		37 & matrix factorization & matrix, vector, eigenvalues, 
		linear, diagonal, eigenvectors, operator \\ 
		44 & experimentals & data, error, test, set, 
		training, prediction, estimate, performance \\ 
		48 & classification & classification, class, classifier, 
		training, feature, decision, data \\ 
		54 & neural networks & units, hidden, weights, learning, 
		input, weight, patterns \\ 
		55 & neural networks & learning, training, examples, task, 
		target, active \\ 
		58 & neural networks & network, input, neural, output, 
		training, layer\\ 
		64 & optimization & learning, gradient, algorithm, error, 
		weight, descent\\
		74 & signal processing & noise, signal, snr, filter, power, 
		frequency\\ 
		81 & theory & let, theorem, function, probability, bound, 
		given, case\\ 
		82 & decision trees & node, tree, decision, leaf, time, root, 
		pruning\\ 
		88 & image processing & points, space, figure, regions, 
		surface, local\\ 
		\bottomrule
	\end{tabular}
	\begin{tablenotes} 
	\item [a] The model is learned via the AGS algorithm. 
	\item [b] Here, we discarded repeated words with a common stem 
	for brevity.	
\end{tablenotes} 
\end{threeparttable} 
\end{table*}
}
\begin{figure*}[ht!] 
	\centering 
	\subfloat[topic $3$ (Science)]{
		\includegraphics[width=.3\linewidth]
		{ags-cgs-16Newsgroups-J4-K30-D10764-V9208-seed1983-e0d05-a1-g2d07-2018Jan31125617-clda-topic-3}}\hfill%
	\subfloat[topic $10$ (Uninformative)]{
		\includegraphics[width=.3\linewidth]
		{ags-cgs-16Newsgroups-J4-K30-D10764-V9208-seed1983-e0d05-a1-g2d07-2018Jan31125617-clda-topic-10}}\hfill%
	\subfloat[topic $21$ (Recreation)]{
		\includegraphics[width=.3\linewidth]
		{ags-cgs-16Newsgroups-J4-K30-D10764-V9208-seed1983-e0d05-a1-g2d07-2018Jan31125617-clda-topic-12}}\\
	\subfloat[topic $13$ (Computers)]{
		\includegraphics[width=.3\linewidth]
		{ags-cgs-16Newsgroups-J4-K30-D10764-V9208-seed1983-e0d05-a1-g2d07-2018Jan31125617-clda-topic-13}}\hfill%
	\subfloat[topic $28$ (Politics, gun trade)]{
		\includegraphics[width=.3\linewidth]
		{ags-cgs-16Newsgroups-J4-K30-D10764-V9208-seed1983-e0d05-a1-g2d07-2018Jan31125617-clda-topic-28}}\hfill%
	\subfloat[topic $29$ (Computers)]{
		\includegraphics[width=.3\linewidth]
		{ags-cgs-16Newsgroups-J4-K30-D10764-V9208-seed1983-e0d05-a1-g2d07-2018Jan31125617-clda-topic-29}}
	\caption{$20$ most probable words for three selected topics from a 
		$30$-topic cLDA model trained on corpus $16$newsgroups.	The 
		$x$-axis gives the corresponding (estimated) probabilities of words 
		given a topic.}
	\label{fig:16news-topics-new}
\end{figure*}




\begin{figure}[t!] 
	\centering
	\includegraphics[width=.8\linewidth]
	{clda-alpha-search-perplexity-boxplots}
	\caption{Boxplots of perplexity scores of cLDA models with 
	various values of $\alpha$ keeping $\eta$ and $\gamma$ fixed, 
	for corpus NIPS $00$-$18$.}
	\label{fig:nips-clda-alpha-sensitivity}
\end{figure}
\section{Perplexity Calculation for cLDA and LDA}
\label{app:perplexity}

In this section, we derive expressions for Perplexity (defined in 
Section~\ref{sec:eval-performance}) for both cLDA and LDA 
models. To compute the perplexity score~\eqref{eq:c-lda-perplexity}, 
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 
\begin{equation}
\ell_{w_{jdi}} (\bm \beta, \bm \theta_{jd}, z_{jdi}) 
= \prod_{k = 1}^K \prod_{v = 1}^V \beta_{kv}^{z_{jdik} w_{jdiv}} 
= \sum_{k = 1}^K z_{jdik} \prod_{v = 1}^V \beta_{kv}^{w_{jdiv}}
\end{equation} 
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 \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 
\begin{equation}
\begin{split}
\ell_{w_{jdi}} (\bm \beta, \bm \theta_{jd}) 
&= E \Bigl [ \sum_{k = 1}^K z_{jdik} \prod_{v = 1}^V 
\beta_{kv}^{w_{jdiv}} \Bigl ] 
= \sum_{k = 1}^K \Bigl [ E(z_{jdik}) \prod_{v = 1}^V 
\beta_{kv}^{w_{jdiv}} \Bigl ] \\ 
&= \sum_{k = 1}^K \Bigl [ \theta_{jdk} \prod_{v = 1}^V 
\beta_{kv}^{w_{jdiv}} \Bigl ] \\ 
\end{split}
\label{eq:c-lda-perplexity-pred-ll}
\end{equation}      
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~\eqref{eq:c-lda-perplexity} by the Monte 
Carlo average 
\begin{equation}
\frac{1}{S} \sum_{s = 1}^S \sum_{k = 1}^K \Bigl [ \theta^{[s]}_{jdk} \prod_{v = 1}^V {\beta^{[s]}}_{kv}^{w_{jdiv}} \Bigl ]
\label{eq:c-lda-perplexity-est}
\end{equation}   
We can use any of the augmented chains on $(\bm \beta, \bm \pi, \bm \theta, 
{\bm z})$ described in Section~\ref{sec:mcmc} to compute this 
average.   
An elegant alternative is to substitute the following estimates of 
$\theta^{[s]}_{jdk}$ and ${\beta^{[s]}}_{kv}$ in 
\eqref{eq:c-lda-perplexity-est}  
\begin{equation}
\begin{split}
\hat{\theta}^{[s]}_{jdk} 
&= E_{\bm \pi_j, {\bm z}_d, {\bm w}_d} \left (\theta^{[s]}_{jdk} \,|\, {\bm z}_d, 
\bm \pi_j, {\bm w}_d \right) = \frac{n_{jdk}^{\text{train}} + \gamma 
\bm \pi^{[s]}_{j}}{n_{jd.}^{\text{train}} + \gamma} \\ 
\hat{\beta}^{[s]}_{kv} 
&= E_{{\bm z}, {\bm w}} \left (\beta^{[s]}_{kv} \,|\, {\bm z}, {\bm w} \right) =  
\frac{m_{..kv}^{\text{train}} + \eta}{m_{..k.}^{\text{train}} + V 
\eta} \\ 
\end{split}
\end{equation}
which can be computed for every sample in a chain on $(\bm \pi, {\bm z})$. 
The resulting likelihood estimate dominates the original estimate 
\eqref{eq:c-lda-perplexity-est} in terms of variance. This approach 
is sometimes called as Rao-Blackwellization, see, e.g., 
\cite[Chapter 4]{RobertCasella:2005}. 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~\ref{app:vem}) into~\eqref{eq:c-lda-perplexity-pred-ll} 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]}$ \citep{GriffithsSteyvers:2004}, we have 
\begin{equation}
\hat{p} (w_{di} \,|\, {\bm w}^{\text{train}})
= \frac{1}{S} \sum_{s = 1}^S \sum_{k = 1}^K \Bigl [ 
\hat{\theta}^{[s]}_{dk} \prod_{v = 1}^V 
{\hat{\beta}^{[s]}}_{kv}^{w_{div}} \Bigl ], 
\label{eq:lda-perplexity}
\end{equation} 
where 
\begin{equation}
\begin{split}
\hat{\theta}^{[s]}_{dk} 
&= \frac{n_{dk}^{\text{train}} + \alpha}{n_{d.}^{\text{train}} + 
K\alpha} \\ 
\hat{\beta}^{[s]}_{kv} 
&= \frac{m_{.kv}^{\text{train}} + \eta}{m_{.k.}^{\text{train}} + 
	V\eta}. 
\end{split}
\end{equation} 
     

\section{Expressions for the Prior, Likelihood, and Posterior}
\label{app: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 
\eqref{eq:cLDA-e}--\eqref{eq:cLDA-a} given in Section 
\ref{sec:hierarchical-model}, we can write the prior $p_h(\bm \psi)$ 
as 

\begin{equation*}
p_h ({\bm z} \,|\, \bm \theta, \bm \pi, \bm \beta) 
p_h (\bm \theta \,|\, \bm \pi) p_h (\bm \pi) p_h (\bm \beta), 
\label{eq:prior-psi}
\end{equation*}
where the individual density functions are given by 
\eqref{eq:cLDA-e}--\eqref{eq:cLDA-b} 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
\eqref{eq:cLDA-e}--\eqref{eq:cLDA-b}, we obtain

\begin{eqnarray}
p_h(\bm \psi) &=& 
  \prod_{j=1}^J 
  \left(\frac{\Gamma\bigl( K \alpha \bigr)}
  {\Gamma(\alpha)^K} 
  \prod_{k=1}^K \pi_{jk}^{\alpha - 1} 
  \right)
	\Biggl[ 
	\prod_{j=1}^J \prod_{d=1}^{D_j} \prod_{k=1}^K 
	\theta_{jdk}^{n_{jdk}} 
	\Biggr] 
	\nonumber \\ 
& &	\Biggl[ 
		\prod_{j=1}^J \prod_{d=1}^{D_j} 
		\Biggl(\frac{\Gamma\bigl( \gamma \bigr)}
								{\prod_{k=1}^K {\Gamma(\gamma \pi_{jk})}} 
					 \prod_{k=1}^K \theta_{jdk}^{\gamma \pi_{jk} - 1} 
		\Biggr)
  \Biggr] \nonumber \\ 
	& & \Biggl[ \prod_{k=1}^K \biggl( 
	\frac{\Gamma(V\eta)}{\Gamma(\eta)^V} \prod_{v=1}^V \beta^{\eta 
	- 1}_{kv} \biggr)
  \Biggr].
\label{eq:dens-of-btz}
\end{eqnarray}

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 
\eqref{eq:cLDA-b} induces the likelihood function 
$\ell_{{\bm w}} (\bm \psi) := p({\bm w} \,|\, {\bm z}, \bm \theta, \bm \pi, \bm \beta)$  
\begin{equation}
  \label{eq:lda-likelihood}
  \begin{split}
 
    p ({\bm w} \,|\, {\bm z}, \bm \theta, \bm \pi, \bm \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{split}
\end{equation}
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~\eqref{eq:dens-of-btz}, 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~\eqref{eq:lda-likelihood} and the
prior~\eqref{eq:dens-of-btz} into the expression for the Bayes rule 
\eqref{eq:posterior-psi}, and 
absorbing constants in Dirichlet normalizing constants into an 
overall constant of proportionality, we have the posterior density 
\begin{eqnarray} 
	\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] \\ 
	& \hspace{4.7mm}
	\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} \nonumber 
\end{eqnarray}

\section{Conditional posterior for $z_{jdi}$}
\label{app:posterior-z}

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  
\begin{eqnarray}
\begin{split}
p_{h, {\bm w}} \left(z_{jdik} = 1 \,|\, {\bm z}^{(-jdi)}, \bm \pi \right) 
& \propto p \left(z_{jdik} = 1, w_{jdiv} = 1 \,|\, {\bm z}^{(-jdi)}, 
{\bm w}^{(-jdi)}, \bm \pi\right)  \\ 
& = \frac{p \left(z_{jdik} = 1, w_{jdiv} = 1, {\bm z}^{(-jdi)} \,|\, 
	{\bm w}^{(-jdi)}, \bm \pi \right)}
{p \left({\bm z}^{(-jdi)} \,|\, {\bm w}^{(-jdi)}, \bm \pi \right)} \\ 
& \propto  
\frac{\gamma \pi_{jk} + n^{(-jdi)}_{jdk}}
{\gamma + n^{(-jdi)}_{jd.}} \,
\frac{\eta  + m^{(-jdi)}_{..kv}}
{V\eta + m^{(-jdi)}_{..k.}}	
\end{split}
\label{eq:posterior-z-derivation}
\end{eqnarray}
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 \eqref{eq:pi-z-posterior}
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~\citep[CGS]{GriffithsSteyvers:2004}.
\section{Variational Inference}
\label{app:vem}

We now develop variational methods \citep{Jordan: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 
\cite{BleiNgJordan: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~\citep{Jordan:1999, Bishop: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 
\citep{Bishop:2006}
\begin{equation}
	\log m(h) = {\cal L}(q, p_{h}) + \text{KL}(q, p_{h, {\bm w}})
	\label{eq:log-marginal-prob}
\end{equation}
where\footnote{The 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.}
\begin{equation}
{\cal L}(q, p_{h}) = \int \sum_{{\bm z}}
q(\bm \psi) \log \left \{ \frac{p_{h}(\bm \psi, {\bm w})}{q(\bm \psi)} \right \}
d\bm \beta d\bm \pi d\bm \theta  
\label{eq:variational-lowerbound-original}
\end{equation}
and
\begin{equation}
	\text{KL}(q, p_{h, {\bm w}}) = -\int
	\sum_{{\bm z}} q(\bm \psi) \log \left \{
		\frac{p_{h, {\bm w}} (\bm \psi)}{q(\bm \psi)} \right \} d\bm \beta d\bm \pi 
d\bm \theta.
	\label{eq:kl-divergence}
\end{equation}
Note that ${\cal L}(q, p_{h})$ 
in~\eqref{eq:variational-lowerbound-original} is a 
functional of the distribution $q(\bm \psi)$ and a function of the 
hyperparameters $h$. The Kullback-Leibler (KL) divergence specified
in~\eqref{eq:kl-divergence} 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~\eqref{eq:log-marginal-prob}, ${\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    
\begin{equation}
q(\bm \psi \,|\, \bm \psi^{\prime}) = 
\Bigl [ \prod_{k = 1}^{K} q(\bm \beta_k \,|\, {\bm \lambda}_k) \Bigr ]  
\Bigl [ \prod_{j = 1}^{J} q(\bm \pi_j \,|\, {\bm \tau}_j) 
\Bigl [ \prod_{d = 1}^{D_j} q(\bm \theta_d \,|\, {\bm \rho}_d) 
\Bigl (
\prod_{i = 1}^{n_{jd}} q({\bm z}_{di} \,|\, \bm \phi_{jdi}  ) 
\Bigr )
\Bigr ] 
\Bigr ]  
\label{eq:variational-posterior}
\end{equation}
We take its independent component distributions 
\cite{BleiNgJordan:2003} are 
\begin{equation}
\begin{split}
\bm \beta_k &\sim \Dir_{V}({\bm \lambda}_k) \\ 
\bm \pi_j &\sim \Dir_{K}({\bm \tau}_j) \\ 
\bm \theta_{jd} &\sim \Dir_{K}({\bm \rho}_{jd}) \\ 
{\bm z}_{jdi} &\sim \Mult_{K}(\bm \phi_{jdi}) \\ 
\end{split}
\label{eq:variational-distributions}
\end{equation}
Note that the lower-bound 
\eqref{eq:variational-lowerbound-original}, which is an expectation 
with respect to \eqref{eq:variational-posterior}, 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  
\cite{BleiNgJordan:2003}'s scheme, but updating 
parameter ${\bm \tau}$ does not have a closed form expression. 
\cite{Kim: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 \ref{alg:vem} and is 
denoted by the acronym VEM. 



\begin{algorithm}[t]
	\KwData{Observed words ${\bm w}$ and document metadata}
	\KwResult{Optimal variational parameters $({\bm \lambda}, {\bm \tau}, 
	{\bm \rho}, 
		\bm \phi)$} \vspace{.3cm}
	
	initialize $({\bm \lambda}^{(0)}, {\bm \tau}^{(0)})$\;
	
	\While{not converged}{ 
		
		\tcp{Step $1$: Expectation} 
		
		\vspace{.1cm}
		
		\For{document $d = 1, \ldots, D_j$, $j = 1, \ldots, J$}{
			
			initialize $\rho_{jdk}^{(0)} = \frac{\gamma \tau_{jk}}
			{\tau_{j.}} + \frac{n_{jd}}{K}, \, k = 1, \ldots, K$\;
			
			\tcp{Variational updates for each document}
			
			\While{not converged}{ 
				
				\For{word $w_{jdi}$, $i = 1, \ldots, n_{jd}$}{
					variational Multinomial update for 
					$\bm \phi_{jdi}$ 
					via~\eqref{eq:var-multinomial-phi}\; 
				}
				variational Dirichlet update for ${\bm \rho}_{jd}$ via 
				\eqref{eq:var-dirichlet-rho}\; 
				
			}
			
		}
		
		\tcp{Variational updates for each topic}
		
		variational Dirichlet update for ${\bm \lambda}_{k}, \, k = 1, 
		\ldots, K$ via~\eqref{eq:var-dirichlet-lambda}\; 
		
		\vspace{.2cm}
		
		\tcp{Step $2$: Maximization}
		
		\vspace{.1cm}
		
		\tcp{Constraint Newton updates for collection-level topic 
			mixtures}
		
		\For{collection $j = 1, \ldots, J$}{
			initialize $(a^{(0)}_j, {\bm \omega}^{(0)}_j)$ based on the 
			current ${\bm \tau}_j$\; 
			\While{not converged}{ 
				constraint Newton update for ${\bm \omega}_{j}$ via 
				\eqref{eq:constraint-newton-step}\; 
				Newton update for $a_j$ 
				via~\eqref{eq:newton-step-a-j}\; 
			}
			set ${\bm \tau}_j = a^{(\text{final})}_j * 
			{\bm \omega}^{(\text{final})}_j$ 
		}
		optimize hyperparameter $h = (\alpha, \gamma, \eta)$\; 
		
	}
	\caption{Variational expectation maximization (VEM)}
	\label{alg:vem}
\end{algorithm}

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~\eqref{eq:variational-lowerbound-original}  as 
follows.  
\begin{equation}
\begin{split}
{\cal L}(q, p_{h}) &= 
\sum_{k = 1}^{K} \mathbb{E}_{q_k}[ \log p_{\eta}(\bm \beta_k) ] 
+ \sum_{j = 1}^{J} \mathbb{E}_{q_j}[ \log p_{\alpha}(\bm \pi) ] 
+ \sum_{j = 1}^{J} \sum_{d = 1}^{D_j} \mathbb{E}_{q_{jd}}[ \log 
p_{\gamma}(\bm \theta_{jd} \,|\, \bm \pi_{j}) ] \\ 
&+ \sum_{j = 1}^{J} \sum_{d = 1}^{D_j} \sum_{i = 1}^{n_{dj}} 
\mathbb{E}_{q_{jdi}}[ \log p(z_{jdi} \,|\, \bm \theta_{jd}) ] 
+ \sum_{j = 1}^{J} \sum_{d = 1}^{D_j} \sum_{i = 1}^{n_{dj}} 
\mathbb{E}_{q_{jdi}}[ \log p(w_{jdi} \,|\, z_{jdi}, \bm \beta) ] \\ 
&- \sum_{k = 1}^{K} \mathbb{E}_{q_k}[ \log q(\bm \beta_k \,|\, {\bm \lambda}_k) 
] 
- \sum_{j = 1}^{J} \mathbb{E}_{q_j}[ \log q(\bm \pi_j \,|\, {\bm \tau}_j)  ] \\
&- \sum_{j = 1}^{J} \sum_{d = 1}^{D_j} \mathbb{E}_{q_{jd}}[ \log 
q(\bm \theta_{jd} \,|\, {\bm \rho}_{jd}) ]
- \sum_{j = 1}^{J} \sum_{d = 1}^{D_j} \sum_{i = 1}^{n_{dj}} 
\mathbb{E}_{q_{jdi}}[ \log q({\bm z}_{jdi} \,|\, \bm \phi_{jdi})  ] \\ 
\end{split}
\label{eq:variational-lowerbound-full}
\end{equation}

To ease notation, we denote $\beta_{k.} = \sum_{v = 1}^{V} 
\beta_{kv}$. Let $\bm \theta \sim \Dir_L ({\bm \rho})$. We evaluate 
the expectation of the log of a single probability component 
$\theta_k$ analytically, using \cite[Appendix 
A.1]{BleiNgJordan:2003}: 
\[\mathbb{E}_{q}[ \log \theta_k \,|\, \rho_k ] = \varPsi(\rho_{k}) - \varPsi(\rho_{.})\] 
Similarly, we evaluate the following expectations analytically:   
\begin{equation}
\begin{split}
\mathbb{E}_{q}[ \log \beta_{kv} \,|\, \lambda_{kv} ] &= \varPsi(\lambda_{kv}) - \varPsi(\lambda_{k.}) \\ 
\mathbb{E}_{q}[ \log \pi_{jk} \,|\, \tau_{jk} ] &= \varPsi(\tau_{jk}) - \varPsi(\tau_{j.}) \\ 
\mathbb{E}_q[\log \theta_{jdk} \,|\, \rho_{jdk}  ] &= \varPsi(\rho_{jdk}) - \varPsi(\rho_{jd.}) \\ 
\end{split}
\end{equation}

\noindent
Let $\bm \theta \sim \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~\citep[Theorem 3.1]{Kim:2013} 
\begin{equation}
\mathbb{E}[ \log\Gamma (\alpha \theta_k)] \leq 
\log\Gamma (\alpha \mathbb{E}[\theta_k]) + 
\frac{\alpha}{\rho_{.}}(1 - \mathbb{E}[\theta_k]) + 
(1 - \alpha\mathbb{E}[\theta_k]) \Bigl[ \log \mathbb{E}[\theta_k] + \varPsi(\rho_{.}) - \varPsi(\rho_{k}) \Bigr] 
\label{eq:kim-lemma}
\end{equation}    

\noindent
We then write the individual expectations as follows:  
\begin{equation}
\begin{split}
\mathbb{E}_{q_k}[ \log p_{\eta}(\bm \beta_k) ] &=
\log\Gamma (V \eta) - 
V \log\Gamma (\eta) + 
\sum_{v = 1}^{V} (\eta - 1) \bigl [ \varPsi(\lambda_{kv}) - \varPsi(\lambda_{k.}) \bigr ] \\
\mathbb{E}_{q_j}[ \log p_{\alpha}(\bm \pi) ] &=
\log\Gamma (K \alpha) - K \log\Gamma (\alpha) + 
\sum_{k = 1}^{K} (\alpha - 1) \bigl [ \varPsi(\tau_{jk}) - \varPsi(\tau_{j.}) \bigr ] \\
\mathbb{E}_{q_{jd}}[ \log p_{\gamma}(\bm \theta_{jd} \,|\, \bm \pi_{j}) ] &= 
\mathbb{E}_q[ \log\Gamma (\gamma)] - \sum_{k - 1}^{K}  \mathbb{E}_q[\log\Gamma (\gamma \pi_{jk})]  
+ \sum_{k = 1}^{K} \mathbb{E}_q[ (\gamma \pi_{jk} - 1) \log \theta_{jdk} ] \\
&\geq \log\Gamma (\gamma) - \sum_{k = 1}^{K} 
\Bigl [ 
\log\Gamma (\gamma \mathbb{E}_q[\pi_{jk}] ) + \frac{\gamma}{\tau_{j.}} (1 - \mathbb{E}_q[\pi_{jk}]) \\ 
&\hspace{3.2cm} + (1 - \gamma \mathbb{E}_q[\pi_{jk}]) \bigl [ \log \mathbb{E}_q[\pi_{jk}] + \varPsi(\tau_{j.}) - \varPsi(\tau_{jk}) \bigr ] 
\Bigr ] \\ 
&\hspace{2cm} + \sum_{k = 1}^{K} 
\Bigl [ 
\gamma \mathbb{E}_q[\pi_{jk}] \mathbb{E}_q[\log \theta_{jdk}] - \mathbb{E}_q[\log \theta_{jdk}] 
\Bigr ] \\
&\geq \log\Gamma (\gamma) - \frac{\gamma}{\tau_{j.}} (K - 1) 
- (\gamma - K) \Bigl [
\log \tau_{j.} - \varPsi(\tau_{j.}) + \varPsi(\rho_{jd.})  
\Bigr ] \\
&\hspace{2cm} - \sum_{k = 1}^{K} 
\Bigl [  
\log\Gamma (\frac{\gamma \tau_{jk}}{\tau_{j.}}) + (1 - \frac{\gamma \tau_{jk}}{\tau_{j.}}) \bigl [ \log (\tau_{jk}) - \varPsi(\tau_{jk}) + \varPsi(\rho_{jdk}) \bigr ] 
\Bigr ] \\
\mathbb{E}_{q_{jdi}}[ \log p(z_{jdi} \,|\, \bm \theta_{jd}) ] &=
\sum_{k = 1}^{K} \phi_{jdik} \bigl [ \varPsi(\rho_{jdk}) - \varPsi(\rho_{jd.}) \bigr ] \\
\mathbb{E}_{q_{jdi}}[ \log p(w_{jdi} \,|\, z_{jdi}, \bm \beta) ] &=
\sum_{k = 1}^{K} \sum_{v = 1}^{V} \phi_{jdik} w_{jdiv} 
\bigl [ \varPsi(\lambda_{kv}) - \varPsi(\lambda_{k.}) \bigr ] \\
\mathbb{E}_{q_k}[ \log q(\bm \beta_k \,|\, {\bm \lambda}_k) ] &=
\log\Gamma (\lambda_{k.}) - \sum_{v = 1}^{V} \log\Gamma (\lambda_{kv}) + 
\sum_{v = 1}^{V} (\lambda_{kv} - 1) \bigl [ \varPsi(\lambda_{kv}) - \varPsi(\lambda_{k.}) \bigr ] \\
\mathbb{E}_{q_j}[ \log q(\bm \pi_j \,|\, {\bm \tau}_j)  ]  &=
\log\Gamma (\tau_{j.}) - \sum_{k = 1}^{K} \log\Gamma (\tau_{jk}) + 
\sum_{k = 1}^{K} (\tau_{jk} - 1) \bigl [ \varPsi(\tau_{jk}) - \varPsi(\tau_{j.}) \bigr ] \\
\mathbb{E}_{q_{jd}}[ \log q(\bm \theta_{jd} \,|\, {\bm \rho}_{jd})  ] &= 
\log\Gamma (\rho_{jd.}) - \sum_{k = 1}^{K} \log\Gamma (\rho_{jdk}) + 
\sum_{k = 1}^{K} (\rho_{jdk} - 1) \bigl [ \varPsi(\rho_{jdk}) - \varPsi(\rho_{jd.}) \bigr ] \\
\mathbb{E}_{q_{jdi}}[ \log q({\bm z}_{jdi} \,|\, \bm \phi_{jdi}) ] &=  
\sum_{k = 1}^{K} \phi_{jdik} \log \phi_{jdik} \\
\end{split}
\label{eq:var-exp}
\end{equation} 
For $\mathbb{E}_{q_{jd}}[ \log p_{\gamma}(\bm \theta_{jd} \,|\, \bm \pi_{jd})]$, 
the second step uses~\eqref{eq:kim-lemma} 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~\citep[Appendix]{BleiNgJordan:2003}. 


\subsection{Variational Dirichlet Update for Collections} 
\label{sec:var-dir-tau}

Grouping the expectations that contain ${\bm \tau}_j$ from the 
lower-bound ${\cal L}$\footnote{We ignore the arguments of ${\cal 
L}$ to ease notation.}, we get: 
\begin{equation}\nonumber 
{\cal L}_{[\tau_j]} = \mathbb{E}_{q_j}[ \log p_{\alpha}(\bm \pi) ] 
- \mathbb{E}_{q_j}[ \log q(\bm \pi_j \,|\, {\bm \tau}_j)  ] 
+ \sum_{d = 1}^{D_j} \mathbb{E}_{q_{jd}}[ \log p_{\gamma}(\bm \theta_{jd} \,|\, \bm \pi_{j}) ] 
\end{equation}
This lower-bound does not produce a closed form expression for   
updating $\tau_j$'s. \cite{Kim: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 \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
\begin{equation} 
\begin{split}
{\cal L}_{[a_j {\bm \omega}_j]} &= 
\sum_{k = 1}^{K} \Bigl[ 
(\alpha - a_j \omega_{jk})  \bigl[ \varPsi(a_j \omega_{jk}) - \varPsi(a_j) \bigr]  + 
\log\Gamma (a_j \omega_{jk}) 
\Bigr] - 
\log\Gamma (a_j ) \\ 
&\hspace{.4cm} - \sum_{d = 1}^{D_j} \Bigl[ 
\frac{\gamma}{a_j} (K - 1) 
+ (\gamma - K) \bigl [
\log a_j - \varPsi(a_j) + \varPsi(\rho_{jd.})  
\bigr ] \Bigr]\\
&\hspace{.4cm} - \sum_{d = 1}^{D_j} \sum_{k = 1}^{K} \Bigl[ 
\log\Gamma (\gamma \omega_{jk}) + (1 - \gamma \omega_{jk}) \bigl [ \log (a_j \omega_{jk}) - \varPsi(a_j \omega_{jk}) + \varPsi(\rho_{jdk}) \bigr ] 
\Bigr]
\end{split}
\end{equation}
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 
\begin{equation} 
\begin{split}
{\cal L}_{[\omega_{jk}]} &= 
\varPsi(a_j \omega_{jk}) \Bigl[ \alpha + D_j - a_j \omega_{jk} - \gamma D_j \omega_{jk} \Bigr]  
+ \gamma \omega_{jk} \sum_{d = 1}^{D_j} \varPsi(\rho_{jdk}) \\ 
&\hspace{.4cm} + \log\Gamma (a_{j} \omega_{jk}) 
- D_j \Bigl[ \log\Gamma (\gamma \omega_{jk}) + (1 - \gamma \omega_{jk}) \log (a_{j} \omega_{jk}) \Bigr]  
\end{split}
\end{equation}
Its first and second derivatives, denoted by $g_{jk}$ and $h_{jk}$, 
are: 
\begin{eqnarray} 
\frac{\partial}{\partial \omega_{jk}} {\cal L}_{[\omega_{jk}]} 
&=& a_j\varPsi^{\prime}(a_j \omega_{jk}) \Bigl[ \alpha + D_j - a_j \omega_{jk} - \gamma D_j \omega_{jk} \Bigr]  
-  \gamma D_j \varPsi(a_j \omega_{jk}) + \gamma \sum_{d = 1}^{D_j} \varPsi(\rho_{jdk}) \nonumber \\ 
&& - \: D_j \Bigl[ \gamma \varPsi(\gamma \omega_{jk}) + \frac{1}{\omega_{jk}}  - \gamma - \gamma \log (a_{j} \omega_{jk}) \Bigr]  
\label{eq:omega-grad} \\ 
\frac{\partial^{2}}{\partial \omega_{jk}} {\cal L}_{[\omega_{jk}]} 
&=& a_j^2\varPsi^{\prime\prime}(a_j \omega_{jk}) \Bigl[ \alpha + D_j - a_j \omega_{jk} - \gamma D_j \omega_{jk} \Bigr] 
- a_{j}\varPsi^{\prime}(a_j \omega_{jk}) \Bigl[ a_j + 2 \gamma D_j \Bigr] \nonumber \\ 
&& - \: D_j \Bigl[ \gamma^2 \varPsi^{\prime}(\gamma \omega_{jk}) - \frac{1}{\omega^2_{jk}} - \frac{\gamma}{\omega_{jk}}  \Bigr]  
\label{eq:omega-hess}
\end{eqnarray}
We can see that the hessian given by 
\eqref{eq:omega-hess} 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   
\begin{equation} 
\begin{bmatrix} \text{diag}(\bm{h}) & \mathbf{1} \\ {\mathbf{1}}^{\text{T}} & 0 \end{bmatrix}
\left[ \begin{array}{c} \Delta \omega_{jk} \\ u \end{array} \right] =  
\left[ \begin{array}{c} -\bm{g} \\ 0 \end{array} \right],
\label{eq:constraint-newton}
\end{equation}
that yields  
\begin{equation} 
\Delta \omega_{jk} = \Bigl\{ \frac{ \sum_{k=1}^{K} 
\frac{g_{jk}}{h_{jk}} }{ \sum_{k=1}^{K} \frac{1}{h_{jk}} } \Bigr\}   
\left[ 
\begin{array}{c} 
\frac{1}{h_{j1}} \\ 
\cdots \\ 
\frac{1}{h_{jK}} 
\end{array} 
\right] - 
\left[ 
\begin{array}{c} 
\frac{g_{j1}}{h_{j1}} \\ 
\cdots \\ 
\frac{g_{jK}}{h_{jK}}
\end{array} 
\right]
\label{eq:constraint-newton-step}
\end{equation}
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}$~\citep{Kim: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 
\begin{equation} 
\begin{split}
{\cal L}_{[a_j]} &= 
\sum_{k = 1}^{K} \bigl[ \varPsi(a_j \omega_{jk}) - \varPsi(a_j) \bigr]
\Bigl( \alpha + D_j - a_j \omega_{jk} - \gamma D_j \omega_{jk} \Bigr)  
+ \sum_{k = 1}^{K} \log\Gamma (a_j \omega_{jk}) \\ 
&\hspace{.4cm}  - \log\Gamma (a_j) - \frac{\gamma D_j (K - 1)}{a_j}
\end{split}
\end{equation}
Its first and second derivatives, denoted by $g_{j}$ and $h_{j}$, 
are: 
\begin{eqnarray} 
\frac{\partial}{\partial a_j} {\cal L}_{[a_j]} 
&=& \sum_{k = 1}^{K} \bigl[ \omega_{jk} \varPsi^{\prime}(a_j \omega_{jk}) - \varPsi^{\prime}(a_j) \bigr]
\Bigl( \alpha + D_j - a_j \omega_{jk} - \gamma D_j \omega_{jk} \Bigr) \nonumber \\ 
&& + \: (K - 1) \gamma D_j a_j^{-2} 
\label{eq:a-grad} \\ 
\frac{\partial^{2}}{\partial a_j} {\cal L}_{[a_j]} 
&=& \sum_{k = 1}^{K} \bigl[ \omega^2_{jk} \varPsi^{\prime\prime}(a_j \omega_{jk}) - \varPsi^{\prime\prime}(a_j) \bigr]
\Bigl( \alpha + D_j - a_j \omega_{jk} - \gamma D_j \omega_{jk} \Bigr) \nonumber \\ 
&& - \: \sum_{k = 1}^{K}  \omega_{jk} \bigl[ \omega_{jk} \varPsi^{\prime}(a_j \omega_{jk}) - \varPsi^{\prime}(a_j) \bigr]
- 2(K - 1) \gamma D_j a_j^{-3}
\label{eq:a-hess}
\end{eqnarray}
The Newton update step for $a_j$ is then given by 
\begin{equation}
\Delta a_j = - h^{-1}_{j} g_{j}. 
\label{eq:newton-step-a-j}
\end{equation}
We alternately maximize ${\cal L}_{[a_j {\bm \omega}_j]}$ with 
respect to $a_j$ and ${\bm \omega}_j$ until convergence. 
 


\subsection{Variational Multinomial Update for Words}
\label{sec:var-dir-phi}

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~\eqref{eq:var-exp}, and applying 
Lagrangian multipliers as
\begin{equation} 
\begin{split}
{\cal L}_{[\phi_{jdik}]} 
&= \phi_{jdik} \Bigl[ [\varPsi(\rho_{jdk}) - \varPsi(\rho_{jd.})] + 
 [\varPsi(\lambda_{kv}) - \varPsi(\lambda_{k.})] - \log \phi_{jdik} \Bigr] \\ 
&+ \mu_{jdi} \Bigl[ \sum_{k=1}^{K} \phi_{jdik} - 1 \Bigr]
\end{split}
\end{equation}
Taking derivatives with respect to $\phi_{jdik}$, we get 
\begin{equation} 
\frac{\partial}{\partial \phi_{jdik}}{\cal L}_{[\phi_{jdik}]} 
= [\varPsi(\rho_{jdk}) - \varPsi(\rho_{jd.})] + 
 [\varPsi(\lambda_{kv}) - \varPsi(\lambda_{k.})] - \log \phi_{jdik} 
- 1 + \mu_{jdi} 
\end{equation}
Setting this to zero yields the maximum value for $\phi_{jdik}$ 
\begin{equation} 
\phi_{jdik} 
\propto \exp \Bigl( [\varPsi(\rho_{jdk}) - \varPsi(\rho_{jd.})] + 
 [\varPsi(\lambda_{kv}) - \varPsi(\lambda_{k.})] \Bigr)
\label{eq:var-multinomial-phi}
\end{equation}
  
\subsection{Variational Dirichlet Update for Documents} 
\label{sec:var-dir-rho}

We maximize the lower-bound ${\cal L}$ with respect to $\rho_{jdk}$. 
Collecting terms that contain $\rho_{jd}$ from~\eqref{eq:var-exp}, 
we get 
\begin{equation} 
\begin{split}
{\cal L}_{[\rho_{jd}]} 
&= \sum_{k = 1}^{K} \Bigl( \frac{\gamma \tau_{jk}}{\tau_{j.}} + \sum_{i = 1}^{n_{jd}}\phi_{jdik} - \rho_{jdk} \Bigr)
 \bigl [ \varPsi(\rho_{jdk}) - \varPsi(\rho_{jd.}) \bigr ] \\ 
&- \log\Gamma (\rho_{jd.}) + \sum_{k = 1}^{K} \log\Gamma (\rho_{jdk})
\end{split}
\end{equation}
Taking derivatives with respect to $\rho_{jdk}$, we get 
\begin{equation} 
\frac{\partial}{\partial \rho_{jdk}}{\cal L}_{[\rho_{jd}]}
=  \bigl [ \varPsi^{\prime}(\rho_{jdk}) - \varPsi^{\prime}(\rho_{jd.}) \bigr ] \Bigl( \frac{\gamma \tau_{jk}}{\tau_{j.}} + \sum_{i = 1}^{n_{jd}}\phi_{jdik} - \rho_{jdk} \Bigr)
\end{equation}
Setting this to zero yields the maximum value for $\rho_{jdk}$
\begin{equation} 
\rho_{jdk} = \frac{\gamma \tau_{jk}}{\tau_{j.}} + \sum_{i = 1}^{n_{jd}}\phi_{jdik}
\label{eq:var-dirichlet-rho}
\end{equation}

\subsection{Variational Dirichlet Update for Topics}
\label{sec:var-dir-lambda}

We maximize the lower-bound ${\cal L}$ with respect to 
$\lambda_{kv}$. Collecting terms that contain ${\bm \lambda}_{k}$ from 
\eqref{eq:var-exp}, we get 
\begin{equation} 
\begin{split}
{\cal L}_{[{\bm \lambda}_{k}]} 
&= \sum_{k = 1}^{K} \sum_{v = 1}^{V} \Bigl( \eta - \lambda_{kv} + \sum_{j = 1}^{J} \sum_{d = 1}^{D_j} \sum_{i = 1}^{n_{jd}}\phi_{jdik} w_{jdiv} \Bigr)
 \bigl [ \varPsi(\lambda_{kv}) - \varPsi(\lambda_{k.}) \bigr ] \\ 
&- \sum_{k = 1}^{K} \bigl [ \log\Gamma (\lambda_{k.}) - \sum_{v = 1}^{V} \log\Gamma (\lambda_{kv}) \bigr ]
\end{split}
\end{equation}
Taking derivatives with respect to $\lambda_{kv}$, we get 
\begin{equation} 
\frac{\partial}{\partial \lambda_{kv}}{\cal L}_{[{\bm \lambda}_{k}]}
=  \bigl [ \varPsi^{\prime}(\lambda_{kv}) - \varPsi^{\prime}(\lambda_{k.}) \bigr ] \Bigl( \eta - \lambda_{kv} + \sum_{j = 1}^{J} \sum_{d = 1}^{D_j} \sum_{i = 1}^{n_{jd}}\phi_{jdik} \Bigr)
\end{equation}
Setting this to zero yields the maximum value for $\lambda_{kv}$
\begin{equation} 
\lambda_{kv} = \eta + \sum_{j = 1}^{J} \sum_{d = 1}^{D_j} \sum_{i = 1}^{n_{jd}}\phi_{jdik} 
\label{eq:var-dirichlet-lambda}
\end{equation}

\subsection{Optimize Hyperparameters} 
\label{sec:opt-h}

As in the variational EM algorithm of 
LDA~\citep{BleiNgJordan:2003}, 
the E-step of the cLDA VEM algorithm updates the variational 
parameters based on the expressions provided above (see 
Algorithm~\ref{alg:vem}). 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: 
\begin{eqnarray} 
{\cal L}_{[\alpha]} &=& J \log\Gamma (K \alpha) - JK \log\Gamma (\alpha) + 
\sum_{j = 1}^{J} \sum_{k = 1}^{K} \alpha \bigl [ \varPsi(\tau_{jk}) - 
\varPsi(\tau_{j.}) \bigr ] \\ 
{\cal L}_{[\eta]} &=& K\log\Gamma (V \eta) - KV \log\Gamma (\eta) + 
\sum_{k = 1}^{K} \sum_{v = 1}^{V} \eta \bigl [ \varPsi(\lambda_{kv}) - 
\varPsi(\lambda_{k.}) \bigr ]  \\ 
{\cal L}_{[\gamma]} &=& \sum_{j = 1}^{J} D_j [\log\Gamma (\gamma) - \frac{\gamma}{\tau_{j.}} 
(K - 1)] - \gamma \sum_{j = 1}^{J} \sum_{d = 1}^{D_j} \Bigl [
\log \tau_{j.} - \varPsi(\tau_{j.}) + \varPsi(\rho_{jd.})  
\Bigr ] \nonumber \\
&& - \sum_{j = 1}^{J} \sum_{d = 1}^{D_j} \sum_{k = 1}^{K} 
\Bigl [  
\log\Gamma (\frac{\gamma \tau_{jk}}{\tau_{j.}}) - 
\frac{\gamma \tau_{jk}}{\tau_{j.}} \bigl [ \log (\tau_{jk}) - 
\varPsi(\tau_{jk}) + \varPsi(\rho_{jdk}) \bigr ] 
\Bigr ]
\end{eqnarray}
The first and second derivatives are the following: 
\begin{eqnarray} 
\frac{\partial}{\partial \alpha}{\cal L}_{[\alpha]} &=& JK \bigl [ 
\varPsi (K \alpha) - \varPsi (\alpha) \bigr ] + 
\sum_{j = 1}^{J} \sum_{k = 1}^{K} \bigl [ \varPsi(\tau_{jk}) - 
\varPsi(\tau_{j.}) \bigr ] \\ 
\frac{\partial^2}{\partial \alpha}{\cal L}_{[\alpha]} &=& JK^2 
\varPsi^{\prime} (K \alpha) - JK \varPsi^{\prime} (\alpha) \\ 
\frac{\partial}{\partial \eta}{\cal L}_{[\eta]} &=& KV \bigl [ 
\varPsi (V \eta) - \varPsi (\eta) \bigr ] + 
\sum_{k = 1}^{K} \sum_{v = 1}^{V} \bigl [ \varPsi(\lambda_{kv}) - 
\varPsi(\lambda_{k.}) \bigr ]  \\ 
\frac{\partial^2}{\partial \eta}{\cal L}_{[\eta]} &=& KV^2
\varPsi^{\prime} (V \eta) - KV \varPsi^{\prime} (\eta) \\ 
\frac{\partial}{\partial \gamma}{\cal L}_{[\gamma]} &=& 
\sum_{j = 1}^{J} D_j [\varPsi (\gamma) - \frac{1}{\tau_{j.}} 
(K - 1)] - \sum_{j = 1}^{J} \sum_{d = 1}^{D_j} \Bigl [
\log \tau_{j.} - \varPsi(\tau_{j.}) + \varPsi(\rho_{jd.})  
\Bigr ] \nonumber \\
&& - \sum_{j = 1}^{J} \sum_{d = 1}^{D_j} \sum_{k = 1}^{K} 
\frac{\tau_{jk}}{\tau_{j.}}
\Bigl [  
\varPsi (\frac{\gamma \tau_{jk}}{\tau_{j.}}) - 
\bigl [ \log (\tau_{jk}) - 
\varPsi(\tau_{jk}) + \varPsi(\rho_{jdk}) \bigr ] 
\Bigr ] \\ 
\frac{\partial^2}{\partial \gamma}{\cal L}_{[\gamma]} &=& 
\sum_{j = 1}^{J} D_j \varPsi^{\prime} (\gamma) 
 - \sum_{j = 1}^{J} \sum_{d = 1}^{D_j} \sum_{k = 1}^{K} 
\frac{\tau^2_{jk}}{\tau^2_{j.}}
\varPsi^{\prime} (\frac{\gamma \tau_{jk}}{\tau_{j.}})
\end{eqnarray}  
Using these derivatives, one can update $\alpha$, $\gamma$, $\eta$ 
via a Newton method~\citep{BleiNgJordan:2003, Minka:2000}.  
\section{Summary}
\label{sec:discussion}

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 \citep{TehEtal: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.
 
\begin{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.     
\end{notes-to-stay}



\section{Experimental Analysis}
\label{sec:experiments}

This section compares the performance of the proposed cLDA model 
with alternatives, LDA and HDP, on real-world corpora (Section 
\ref{sec:datasets}). 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~\ref{sec:eval-performance}), and (b) the applicability of 
cLDA model to real-world corpora 
(Section~\ref{sec:usability-study}). We also briefly discuss 
guidelines for selecting the number of topics $K$ in cLDA 
(Section~\ref{sec:eval-performance}). 


\subsection{Datasets}
\label{sec:datasets}

We use three document corpora based on: (i) the NIPS $00$-$18$ 
dataset~\citep{Chechik:2007}, (ii) customer reviews from 
Yelp\footnote{\url{http://uilab.kaist.ac.kr/research/WSDM11/}}, and 
(iii) the $20$Newsgroups  
dataset\footnote{\url{http://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 
$0$ 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.



\subsection{Model Evaluation and Model Selection}
\label{sec:eval-performance}

\subsubsection{Comparisons with Models in the Prior Work} 
\cite[Section 3]{WallachMimnoMcCallum:2009} proposed 
\emph{non-uniform} base measures instead of the typical 
\emph{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{\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{\text{iid}}{\sim} G_d, \, i = 1, \ldots, n_d$, and that for $d = 1, 
\ldots, D$, $G_d \stackrel{\text{iid}}{\sim} {\cal D}_{G_0,\alpha}$, the Dirichlet process 
(DP) with base probability measure $G_0$ and precision parameter 
$\alpha > 0$ \citep{Ferguson:1973, Ferguson: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 \emph{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~\citep{TehEtal: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 DP\footnote{Code by Wang and Blei 
(2010): \url{https://github.com/blei-lab/hdp}}. We use the popular Chinese Restaurant Franchise~\citep[CRF]{TehEtal:2006} sampling scheme for inference. 


\subsubsection{Criteria for Evaluation} 
We use the perplexity 
computation scheme that is specified as in \cite{Patterson:2013, 
WallachMurrayMimmo:2009}. 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  
\begin{equation}
{\cal S}({\bm w}^{\text{test}} 
\,|\, 
{\bm w}^{\text{train}}) = \exp \frac{-\sum_{ 
w_{jdi} 
		\in  {\bm w}^{\text{test}} } \log p (w_{jdi} \,|\, 
		{\bm w}^{\text{train}}
	)}{| {\bm w}^{\text{test}} |}   
\label{eq:c-lda-perplexity}
\end{equation} 
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~\ref{app:perplexity}). 


\cite{MimnoEtAl: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  
\textsl{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 
\citep{GriffithsSteyvers:2004}, AGS, and  CRF). 

Another option is the topic \textsl{coherence score} 
\citep{MimnoEtAl: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 \textsl{coherence score}
\citep{MimnoEtAl:2011}, for topic $k = 1, 2, \ldots, K$, as 
\begin{equation}
\text{topic-coherence}_k = \sum_{i=2}^{m} \sum_{j=1}^{i} \log 
\frac{\text{df}(v^{(k)}_i, v^{(k)}_j) + 
	1}{\text{df}(v^{(k)}_j)} 
\label{eq:mimno-topic-coherence}
\end{equation}       
The intuition behind this score is that group of words 
belonging to a topic possibly co-occur with in a 
document.


\subsubsection{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. 
\figurename~\ref{fig:lda-clda-K-search-perplexity-lineplots} 
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. 

\begin{figure}[htbp!]
	\centering
	\subfloat[$16$Newsgroups]
	{\includegraphics[width=.66\textwidth]
		{16newsgroups-lda-clda-K-search-perplexity-lineplots}} \\ 
	\subfloat[Yelp]
	{\includegraphics[width=.65\textwidth]
		{yelp-lda-clda-K-search-perplexity-lineplots}}\\
	\subfloat[NIPS $00$-$18$]
	{\includegraphics[width=.65\textwidth]
		{19nipsproceedings-lda-clda-K-search-perplexity-lineplots}} 
	\caption{Estimated (average) per-word perplexities for models 
		LDA and cLDA with single collection ($J = 1$),  
		multi-collection  ($J > 1$) and various 
		configurations of the number of topics $K$ using corpora 
		$16$Newsgroups, Yelp, and NIPS $00$-$18$. }
	\label{fig:lda-clda-K-search-perplexity-lineplots}
\end{figure}

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.  

\textsl{Selecting $K$ and Hyperparameters $\alpha, \gamma, \eta$:} 
\figurename~\ref{fig:lda-clda-K-search-perplexity-lineplots} 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 \ref{sec:clda-estimate-eta-gamma}) 
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 
\eqref{eq:gem-est-eta} and \cite{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.              

\eat{
We now describe the experimental setup for evaluating the 
performance of cLDA models using corpora $16$Newsgroups, NIPS 
$00$-$18$, and Yelp. It should be noted that the number of 
topics and hyperparameters are fixed for both LDA and cLDA models. 
One option is to select them via a model selection procedure (e.g., 
based on perplexity).
\begin{notes-to-stay}
\figurename~\ref{fig:select-k-perplexity} 
gives plots of per-word 
perplexities of the held-out words for corpora $16$Newsgroups and 
NIPS $00$-$18$ via the chains CGS AGS, for various choices of $K$. 
As we can see, perplexities of cLDA models with $K = 64, 80$, and 
$128$  are reasonably close for corpus $16$Newsgroups. Similarly, 
perplexities of cLDA models with $K = 90$ and $128$ are relatively 
similar for corpus NIPS $00$-$18$. For cLDA and LDA, we set $K = 
64$ for the $16$Newsgroups corpus and $K = 90$ for the NIPS $00$-$18$ 
corpus as they yield parsimonious models. 
\begin{figure}[ht!] 
\centering 
\subfloat[AGS: $16$Newsgroups]{
\includegraphics[width=.48\linewidth]
{select-k-15Newsgroups-J4-D10764-V9208-seed1983-e0d25-a0d5-g1-2016Jun17152718-clda-perplexities}}
\subfloat[CGS: $16$Newsgroups]{
\includegraphics[width=.48\linewidth]
{select-k-15Newsgroups-J4-D10764-V9208-seed1983-e0d25-a0d5-g1-2016Jun17152718-lda-perplexities}}
\subfloat[AGS: NIPS $00$-$18$]{
\includegraphics[width=.48\linewidth]
{select-k-19nipsproceedings-J4-D2741-V9156-seed1983-e0d25-a0d5-g1-2016Jun17152643-clda-perplexities}}
\subfloat[CGS: NIPS $00$-$18$]{
\includegraphics[width=.48\linewidth]
{select-k-19nipsproceedings-J4-D2741-V9156-seed1983-e0d25-a0d5-g1-2016Jun17152643-lda-perplexities}}
\caption{Per-word perplexities of the held-out words via the chains 
CGS and AGS with various choices of $K$, for corpora $16$Newsgroups 
and NIPS $00$-$18$.}
\label{fig:select-k-perplexity}
\end{figure}
\end{notes-to-stay}
We ran the Markov chains AGS and CGS for $1000$ iterations. 
\figurename~\ref{fig:perplexity-20news:b} gives plots of per-word 
perplexities for the test words using algorithms CGS and AGS for 
corpus $16$Newsgroups. From the plots, we see that the cLDA model   
outperforms the LDA model reasonably well in terms of predictive 
power.    


\begin{figure}[t!] 
\centering 
\includegraphics[width=1\linewidth]
{perp-15Newsgroups-J4-K64-D10764-V9208-seed1983-e0d25-a0d5-g1-2016Jun18162245-perplexity2}
\caption{Per-word perplexities of the held-out words for corpus 
$16$Newsgroups via the LDA CGS algorithm (dotted line) and the cLDA 
AGS algorithm with single collection $J = 1$ (dashed line) and 
multiple collections $J = 4$. }
\label{fig:perplexity-20news:b}
\end{figure}


In \figurename~\ref{fig:perplexity-20news:b}, the dashed 
line represents the per-word perplexities for the test words for a 
cLDA model with one collection---i.e., we assumed that all 
documents belong to a single collection. 
\figurename~\ref{fig:perplexity-20news:b}
shows that the proposed cLDA models are superior to 
the LDA model, 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. 


\figurename~\ref{fig:nips-perplexity-aux-vem-cgs} gives plots of 
per-word perplexities for the test words via algorithms AGS and CGS 
for corpus NIPS $00$-$18$. We took the number of collections $J = 
4$ and the number of topics $K = 90$ for this corpus. The plots  
shows that both AGS and CGS chains perform almost equal, but the 
AGS chain has a slight edge. We believe this is partly due to the
nature of partitions defined for the corpus; collections in this 
corpus share many common topics, compared to the 16Newsgroups 
corpus discussed above. It also suggests that the gain in using 
cLDA is the more when we model corpora with separable collections. 
Lastly, \figurename~\ref{fig:yelp-perplexity-aux-vem-cgs}
gives plots of per-word perplexities via algorithms AGS and CGS for 
corpus Yelp. We took the number of collections $J = 5$ and the 
number of topics $K = 60$ for this corpus. This plot also confirms 
the superiority of the cLDA model over the LDA model.

\begin{figure}[t!] 
\includegraphics[width=1\linewidth]
{perp-19nipsproceedings-J4-K90-D2741-V9156-seed1983-e0d25-a0d5-g1-2016Jun18162135-ags-vs-cgs}
\caption{Per-word perplexities of the held-out words for 
the NIPS $00$-$18$ corpus via the LDA CGS algorithm and the cLDA 
AGS algorithm with the 
number of collections $J = 4$.}
\label{fig:nips-perplexity-aux-vem-cgs}
\end{figure}

\begin{figure}[t!] 
\includegraphics[width=1\linewidth]
{perp-yelp-J5-K60-D24310-V9517-seed1983-e0d25-a0d5-g1-2016Jul06194443-ags-vs-cgs}
\caption{Per-word perplexities of the held-out words for the  
Yelp corpus via the LDA CGS algorithm and the cLDA AGS algorithm 
with the number of collections $J = 5$.}
\label{fig:yelp-perplexity-aux-vem-cgs}
\end{figure}
}


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. 
\figurename~\ref{fig:lda-clda-topic-coherences-sizes} 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.       

\begin{figure}[htbp!]
	\centering
	\subfloat[Topic size]
	{\includegraphics[width=.65\textwidth]
		{lda-clda-topic-sizes}} \\ 
	\subfloat[Topic coherence]
	{\includegraphics[width=.65\textwidth]
		{lda-clda-topic-coherences}}
	\caption{Estimated topic sizes and coherences for models 
		LDA and cLDA  for corpus $16$Newgroups, NIPS $00$-$18$, and Yelp. }
	\label{fig:lda-clda-topic-coherences-sizes}
\end{figure}


\textsl{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~\ref{tab:computation-costs-extended}). 
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  
\texttt{R}-\texttt{C++} (using the efficient \texttt{Rcpp} and 
\texttt{RcppArmadillo} libraries) programming environment. 

\begin{table*}%
	\centering
	\begin{threeparttable}[b] 
		\caption{Average per iteration execution time (in seconds) for  
			various topic modeling algorithms. } 
		\label{tab:computation-costs-extended}
		\begin{tabular}{c c r r r r}
			\toprule 
			Corpus & Number of Topics & CGS & AGS & 
			MGS\tnote{a} & VEM\tnote{b} \\  
			\midrule
			$16$Newsgroups & $64$ & $0.6101$ & $0.6879$ & $1.8616$ & 
			$52.3518$\\ 
			NIPS $00$-$18$ & $90$ & $3.9568$ & $4.2588$ & $12.5580$ & 
			$320.9931$ 
			\\ 	
			Yelp & $60$ & $1.2997$ & $1.4743$ & $3.9524$ & $164.0285$ \\
			\bottomrule 
		\end{tabular}
		\begin{tablenotes} 
			\item [a] MGS chains took approximately twice or more the 
			execution time of CGS chains for all three corpora.
			\item [b] We can see that per iteration running cost for VEM is 
			relatively high for these corpora. We believe that this due to 
			the single-threaded vanilla implementation for VEM; it may be 
			improved by an efficient parallel implementation, which we 
			leave to future research. 	
		\end{tablenotes} 
	\end{threeparttable} 
\end{table*}


\subsubsection{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 \texttt{hclust} 
algorithm on topic distributions. \texttt{hclust} first computes a 
topic-to-topic similarity matrix based on the \emph{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 
\ref{app:nips-results}). HDP, on the other hand, found too many 
redundant topics, as shown in 
\figurename~\ref{fig:lda-clda-hdp-topic-sizes}. To evaluate 
clusters induced by \texttt{hclust} quantitatively, we look at 
silhouette widths computed on topics' \texttt{hclust} clusters. We 
favor methods with high silhouette widths. 
\figurename~\ref{fig:clda-lda-hdp-topic-hclust-silhouette} gives 
boxplot statistics (i.e. median, lower hinge, and upper hinge) of 
silhouette widths computed on topics' \texttt{hclust} clusters with 
various values of the number of clusters (a user specified value in 
\texttt{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 
(\figurename~\ref{fig:clda-lda-hdp-doc-hclust-silhouette})
\begin{figure}[htbp!]
	\centering
	\includegraphics[width=.65\textwidth]
	{clda-lda-hdp-topic-hclust-silhouette}
	\caption{Boxplot statistics (median, lower hinge, and upper 
		hinge) of silhouette widths computed on \emph{topics}' 
		\texttt{hclust} clusters (i.e. hierarchical clustering of 
		topic distributions) with various values of the number 
		of clusters,  for corpus NIPS $00$-$18$ 
	}
	\label{fig:clda-lda-hdp-topic-hclust-silhouette}
\end{figure}
\begin{figure}[htbp!]
	\centering
	\includegraphics[width=.65\textwidth]
		{lda-clda-hdp-topic-sizes}
	\caption{Estimated topic sizes for models 
		LDA, cLDA, and HDP  for corpus NIPS $00$-$18$. }
	\label{fig:lda-clda-hdp-topic-sizes}
\end{figure}


\subsection{Usability Study}
\label{sec: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 \figurename~\ref{fig:nips-clda-lda-topics} 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. 


\begin{figure*}[t!] 
	\centering
	\subfloat[cLDA: t$8$ (ANN)\label{fig:nips-clda-lda-topics:a}]{
		\includegraphics[width=.3\linewidth]
		{ags-cgs-19nipsproceedings-J4-K90-D2741-V9156-seed1983-e0d026-a1-g6d72-2018Jan31125540-clda-topic-8}}\hfill%
	\subfloat[cLDA: t$57$ (MCMC)\label{fig:nips-clda-lda-topics:b}]{
		\includegraphics[width=.3\linewidth]
		{ags-cgs-19nipsproceedings-J4-K90-D2741-V9156-seed1983-e0d026-a1-g6d72-2018Jan31125540-clda-topic-57}}\hfill%
	\subfloat[cLDA: t$85$ (Bayesian models)\label{fig:nips-clda-lda-topics:c}]{
		\includegraphics[width=.3\linewidth]
		{ags-cgs-19nipsproceedings-J4-K90-D2741-V9156-seed1983-e0d026-a1-g6d72-2018Jan31125540-clda-topic-85}}\\
	\subfloat[LDA: t$8$ (ANN)\label{fig:nips-clda-lda-topics:d}]{
		\includegraphics[width=.3\linewidth]
		{ags-cgs-19nipsproceedings-J4-K90-D2741-V9156-seed1983-e0d026-a1-g6d72-2018Jan31125540-lda-topic-8}}\hfill%
	\subfloat[LDA: t$57$ (MCMC)\label{fig:nips-clda-lda-topics:e}]{
		\includegraphics[width=.3\linewidth]
		{ags-cgs-19nipsproceedings-J4-K90-D2741-V9156-seed1983-e0d026-a1-g6d72-2018Jan31125540-lda-topic-57}}\hfill%
	\subfloat[LDA: t$85$ (Bayesian models)\label{fig:nips-clda-lda-topics:f}]{
		\includegraphics[width=.3\linewidth]
		{ags-cgs-19nipsproceedings-J4-K90-D2741-V9156-seed1983-e0d026-a1-g6d72-2018Jan31125540-lda-topic-85}}
	\caption{$20$ most probable words for topics from $90$-topic cLDA (top row) and LDA (bottom row) 
	models trained on corpus NIPS $00$-$18$. Bars represent the corresponding (estimated) 
	probabilities of words given a topic.}
	\label{fig:nips-clda-lda-topics}
\end{figure*}


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~\ref{sec:introduction}. 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). 

\figurename~\ref{fig:nips-sel-pi} 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 
\figurename~\ref{fig:nips-sel-pi:a}, 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, 
\figurename~\ref{fig:nips-clda-pi-samples}.)
\figurename~\ref{fig:nips-sel-pi:b} 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 
\figurename~\ref{fig:nips-sel-pi:c}, 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 
\figurename~\ref{fig:nips-sel-pi:a}. 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~\ref{app:nips-results}.) 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.   


\begin{figure}[t!] 
\centering 
\subfloat[M$1$: AGS run on the corpus with four collections 
\label{fig:nips-sel-pi:a}]{
	\includegraphics[width=.7\linewidth]
	{19nipsproceedings-clda-sel-pi-samples}}\\% 
\subfloat[M$2$: Single CGS run on the whole 
corpus\label{fig:nips-sel-pi:b}]{
	\includegraphics[width=.7\linewidth]
	{19nipsproceedings-lda-sel-pi-est}}\\% 
\subfloat[M$3$: Single CGS run on each of the four partitions of 
the 
corpus\label{fig:nips-sel-pi:c}]{
	\includegraphics[width=.7\linewidth]
	{19nipsproceedings-split-lda-sel-pi-est}}
\caption{Estimates of topic distributions (of five randomly 
	selected topics) for four timespans $1988$-$1992$, 
	$1993$-$1997$, 
	$1998$-$2002$, and $2003$-$2005$ of the NIPS conference 
	proceedings 
	via three topic modeling approaches. Clearly, cLDA identifies 
	interpretable topic patterns compared to other two LDA 
	approaches. 
	See discussion in the text.
}
\label{fig:nips-sel-pi}
\end{figure}

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. 
\figurename~\ref{fig:yelp-clda-pi} 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~\ref{tab:yelp-topic-desc} 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 \figurename~\ref{fig:yelp-lda-pi}.    


\begin{figure}[t!] 
	\centering 
	\subfloat[M$1$: AGS run on the corpus with 
	five collections 
	\label{fig:yelp-clda-pi}]{
		\includegraphics[width=.8\linewidth]
		{ags-cgs-yelp-J5-K40-D24310-V9517-seed1983-e0d034-a1-g3d44-2018Jan31125648-clda-sel-pi-samples}}
	\\%
	\subfloat[M$2$: Single CGS run on the whole corpus
	\label{fig:yelp-lda-pi}]{
		\includegraphics[width=.8\linewidth]
		{ags-cgs-yelp-J5-K40-D24310-V9517-seed1983-e0d034-a1-g3d44-2018Jan31125648-lda-sel-pi-est}}
	\caption{Estimates of topic proportions for five customer 
		ratings $1$, $2$, $3$, $4$, and $5$ via algorithms AGS and 
		CGS 
		using the Yelp corpus. Table~\ref{tab:yelp-topic-desc} 
		gives 
		additional details about the listed topics. Clearly, cLDA 
		provides meaningful and interpretable patterns of topics in 
		the corpus for different customer ratings, compared to its  
		non-hierarchical precursor, LDA. In both of these models, 
		topics are aligned to ease comparison, based on most 
		probable 
		words given a topic. See description in the text.     
	}
	\label{fig:yelp-sel-pi}	
\end{figure} 


\begin{table}[t!] 
	\centering
	\begin{threeparttable}[b] 
		\caption{A subset of topics from a 
			$40$-topic cLDA model for corpus Yelp}
		\label{tab:yelp-topic-desc}
		\begin{tabular}{c c l}
			\toprule
			Topic ID\tnote{a} & Sentiment 
			& 
			Short Description \\ 
			\midrule
			3 & positive & food \\  
			11 & positive & food, place \\ 
			26 & positive & food, ambiance, group dining \\ 
			33 & positive & food, atmosphere \\ 
			5 & moderately positive & food, service \\ 
			6 & moderately positive & food \\ 
			13 & moderately positive & food, waiting time \\ 
			35 & moderately positive & food, sandwich, salad \\ 
			39 & moderately positive & place, ambiance \\ 
			8 & neutral & Mexican food\\ 
			9 & neutral & place, food \\ 
			16 & neutral & bar, food, atmosphere \\ 
			20 & neutral & location, parking\\ 
			34 & neutral & location, ambiance \\ 
			1 & moderately negative &  ambiance, night-life \\ 
			7 & moderately negative &  food, service \\ 
			31 & moderately negative & food, meal quantity \\  
			14 & negative & food, service \\ 
			18 & negative & place, food \\   
			\bottomrule
		\end{tabular}
		\begin{tablenotes} 
			\item [a] The topics are ordered in the 
			non-decreasing order of the positive 
			sentiment or polarity of the most probable 
			words given a topic.
		\end{tablenotes} 
	\end{threeparttable} 
\end{table}


  
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. 
\figurename~\ref{fig:16news-clda-pi-samples-new} 
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. \figurename~\ref{fig:16news-topics-new} 
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 
\figurename~\ref{fig:16news-clda-pi-samples-new}: this conveys the 
fact that corpus $16$newsgroups may be surely separable via the 
subject matter of newsgroups. 

\begin{figure}[t!] 
	\centering
	\includegraphics[width=.8\linewidth]
	{ags-cgs-16Newsgroups-J4-K30-D10764-V9208-seed1983-e0d05-a1-g2d07-2018Jan31125617-clda-pi-samples}
	\caption{Estimates of topic distributions for collections Computers, 
		Politics, Recreation, and Science in corpus $16$newsgroups.}
	\label{fig:16news-clda-pi-samples-new}
\end{figure}


  
\section{A Compound Hierarchical Model} 
\label{sec: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~\ref{sec:experiments}.) 
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})$. 

\begin{figure}[ht!]
\centering 
\begin{tikzpicture}[x=4.5cm,y=2cm]

\node[obs]                      (w)       {$w_{jdi}$} ; %
\node[latent, left=.4 of w]     (z)       {$z_{jdi}$} ; %
\node[latent, left=.4 of z]     (theta)   {$\bm \theta_{jd}$}; %
\node[const, above=.7 of theta](gamma)   {$\gamma$};
\node[latent, left=.4 of theta, style=dashed] (pi)      {$\bm \pi_j$}; 
\node[const, above=.7 of pi]   (alpha)   {$\alpha$};

\node[latent, above=.6 of w]        (beta)     {$\bm \beta_k$}; %
\node[const, left=.55 of beta]    (eta)     {$\eta$}; %

\edge[semithick,style=dashed]    {alpha}   {pi} ;
\edge[semithick,style=dashed]    {pi}      {theta} ;
\edge[semithick]                 {gamma}   {theta} ;
\edge[semithick]                 {theta}   {z} ;
\edge[semithick]                 {eta}     {beta} ;
\edge[semithick]                 {z, beta}  {w} ;

\plate {word-plate} { %
	(w) %
	(z) %
} {$i = 1, \ldots, n_{jd}$}; %

\plate {doc-plate} { %
	(word-plate) %
	(theta) %
} {$d = 1, \ldots, D_j$} ; %

\plate [style=dashed]{collection-plate} { %
	(doc-plate) %
	(pi) %
} {$j = 1, \ldots, J$} ;


\plate {topic-plate} { %
	(beta) %
} {$k = 1, \ldots, K$} ; %

\end{tikzpicture}
\caption{Graphical model of the latent Dirichlet allocation (LDA) 
model (the inner structure) and the compound latent Dirichlet 
allocation (cLDA) model (the outer dashed-structure is the extension 
to LDA): Nodes denote random variables, shaded nodes denote observed 
variables, edges denote conditional dependencies, and plates denote  
replicated processes.\label{fig:gm}}
\end{figure}

We use $\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 
$\Dir_L(a, \ldots, a)$ denotes an $L$-dimensional Dirichlet 
distribution with the constant base measure $(1, \ldots, 1)$. 
$\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 \figurename~\ref{fig:gm}).
\begin{eqnarray}
\bm \beta_k &\stackrel{\text{iid}}{\sim}& \Dir_V(\eta, \ldots, \eta), 
\text{ for topic } k = 1, \ldots, K \label{eq:cLDA-e} \\ 
\bm \pi_j &\stackrel{\text{iid}}{\sim}& \Dir_K(\alpha, \ldots, \alpha), 
\text{ for collection } j = 1, \ldots, J \label{eq:cLDA-d} \\
\bm \theta_{jd} &\stackrel{\text{iid}}{\sim}& \Dir_K(\gamma \pi_{j1}, \ldots, \gamma 
\pi_{jK}), \text{ for document } jd 
\label{eq:cLDA-c}\\ 
z_{jdi} &\stackrel{\text{iid}}{\sim}& \Mult_K(\bm \theta_{jd}), 
\text{ for each word } w_{jdi} \label{eq:cLDA-b}\\ 
w_{jdi} &\stackrel{\text{ind}}{\sim}& \Mult_V(\bm \beta_{z_{jdi}}) \label{eq:cLDA-a}
\end{eqnarray}

 

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~\citep{BleiNgJordan: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.

\subsection*{Related Work}

\cite{Zhai: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~\citep{TehEtal:2006} and Nested Chinese Restaurant 
Process~\citep{Blei: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 
\ref{sec:eval-performance}). 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~\citep{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~\ref{sec:inference}).




In light of adding supervision, several modifications to LDA model 
have been proposed in the literature. Supervised LDA 
\citep[sLDA]{McauliffeBlei:2008} 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.\eat{For regression, variable $y_d$ is normally 
distributed with the mean as the empirical distribution. For 
classification, one can assume a discrete likelihood for the class 
labels $y_d$.} 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~\citep{ZhuXing:2014}.  cLDA, on the other hand, proposes 
a generative framework incorporating the collection-level 
characteristics, without much computational burden (details,  
Section~\ref{sec:experiments}). Also, the objectives of these 
related models are different from the focus of this paper.


\begin{notes-to-stay}
One can view the hierarchical model of cLDA as a special instance 
of the model proposed by \citet{Kim:2013}. \citet{Kim: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 
\citet{Kim: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. 
\end{notes-to-stay} 


\section{Introduction}
\label{sec: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 
\citep{SaltonWongYang:1975}, latent semantic indexing  
\citep[LSI]{Deerwester:1990}, and probabilistic topic models, e.g. 
probabilistic LSI~\citep[pLSI]{Hofmann:1999}, latent Dirichlet 
allocation~\citep[LDA]{BleiNgJordan:2003}, 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 (\figurename~\ref{fig:gm}) 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 \citep{Zhai:2004}. Suppose, we have 
\textsl{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.

\begin{figure}[t!] 
	\centering 
	\includegraphics[width=.65\linewidth]
	{19nipsproceedings-clda-sel-pi-samples-labels}
	\caption{
		Estimated topic proportions of five topics for the four 
		time spans 
		of the NIPS conference proceedings. 
	}
	\label{fig:nips-clda-sel-pi}
\end{figure}

For example, \figurename~\ref{fig:nips-clda-sel-pi} 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~\ref{sec:experiments}). 
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 \ref{sec:eval-performance}). 

In this paper, we introduce the compound latent Dirichlet 
allocation (cLDA, Section~\ref{sec:hierarchical-model}) 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 \figurename~\ref{fig:nips-clda-sel-pi}. 
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~\citep[CGS]{GriffithsSteyvers:2004} (details in 
Sections~\ref{sec:pi-sampling-illustration}, 
\ref{sec:experiments}, \ref{sec:algorithm-illustration}, and 
\ref{app:nips-results}), 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 
\ref{sec:pi-sampling-illustration} and Section 
\ref{sec:algorithm-illustration}) 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~\ref{sec:hierarchical-model} formally defines the cLDA
hierarchical model. Section~\ref{sec:inference} describes 
algorithms for posterior inference and evaluates correctness of the 
algorithms using a synthetically corpus. In 
Section~\ref{sec:experiments}, 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 
\citep{WallachEtal:2009}, and external 
measures such as topic coherence \citep{MimnoEtAl: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~\ref{sec:discussion} concludes this paper with a  
summary of our work.  

\eat{
\IEEEPARstart{N}{ewspapers}, 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. Among many others (e.g. part-of-speech 
tagging, entity recognition, and sentiment analysis) core text 
mining tasks include text categorization, summarization, and 
semantic analysis, which are known to be challenging tasks for 
computers. To address these tasks in an unsupervised manner, 
numerous statistical methods such as TF-IDF 
\citep{SaltonWongYang:1975}, latent semantic indexing  
\citep[LSI]{Deerwester:1990}, and topic models have emerged in the 
literature. With minimal human effort, they can be generalized to 
arbitrary text in any natural language.  

Text analysis tools such as topic models are developed to capture 
the underlying semantic structure of a corpus (i.e. a document 
collection) based on co-occurrences of words. They can be used for 
efficient browsing as well as large-scale analyses of text. Both   
probabilistic \citep{Hofmann:1999, BleiNgJordan:2003} and 
nonprobabilistic (e.g. TF-IDF, LSI) schemes of topic modeling have 
been developed in the literature: the former gives parsimonious and 
more general solutions, but the latter is more computationally 
efficient. However, traditional topic models center on a single 
collection of text documents. This assumption is insufficient for 
comparative analyses of text, especially, for partitionable text 
corpora. 

Suppose, we deal with corpora of (i) News articles tagged into 
categories based on subject matters Politics, Business, Opinion, or 
Sports, (ii) articles accepted in various workshops or consecutive 
proceedings of a conference, or (iii) 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 \citep{Zhai:2004}. Moreover, the category 
labels, workshop names, article timestamps, or geotags can provide 
significant prior knowledge regarding the unique structure of those 
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.   

Challenges particular to multi-collection corpora topic modeling 
include: (i) Blending multi-collection information to the model and 
exploring hidden topics across different collections are non 
trivial. (ii) The discriminative element is more challenging, 
i.e., for each topic identified, it is of interest to find the 
information specific to each collection. It also expects an  
``alignment'' of various collections based on the corpus-level 
topics identified. Zhai et al. \citep{Zhai:2004} refer to this 
problem of  multi-collection corpora modeling as comparative text 
mining (CTM). They use probabilistic latent semantic 
indexing~\citep[pLSI]{Hofmann:1999} as the building blocks of CTM 
and their framework has several user configurable parameters. 

Latent Dirichlet Allocation~\citep[LDA]{BleiNgJordan:2003}, which 
extends over pLSI, is a popular probabilistic topic model. LDA, 
pLSI, and similar topic models consider documents as bags of words, 
considering the order of words in a document uninformative. 
LDA (\figurename~\ref{fig:gm}) 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. 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) Any predefined structure or organization of 
collections is ignored, and as a result, we 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. 



In this paper, we introduce the compound latent Dirichlet  
allocation (cLDA) model that incorporates any organization 
hierarchy of documents in a corpus into a unified modeling 
framework, assuming 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 (definition, 
Section~\ref{sec:hierarchical-model}). All these collection-level 
and document-level variables share a common set of topics at the 
corpus level, which enables us to perform interesting inferences, 
as shown in \figurename~\ref{fig:nips-clda-sel-pi}. cLDA exhibits a 
certain degree of supervision, since we assume that the collection 
membership for each word (and document) in a corpus is known 
implicitly.  On social bookmarking websites, it is natural for the 
readers to annotate or tag written text. cLDA can produce a 
one-to-one correspondence between possible topics and these 
user-defined tags. cLDA can aid visual thematic analyses of large 
sets of documents, including corpus, collection, and document 
specific views. 
 

The parameters of interest in cLDA are hidden and can be inferred 
via posterior inference. However, exact inference is intractable in 
cLDA; and, non-conjugate relationships in the model hierarchy 
further makes approximate inference challenging. Popular methods 
for approximate  posterior inference in topic models are Markov 
chain Monte Carlo (MCMC) and variational methods.  MCMC methods 
enable us to sample hidden variables of interest from the 
intractable posterior and are quite reliable regarding their 
convergence properties. We consider two MCMC methods for cLDA: (a) 
one uses the traditional auxiliary variable updates within Gibbs 
sampling scheme, and (b) the other uses Langevin dynamics within 
Gibbs sampling, a method that received recent attention. The former 
method gives superior performance without any computational 
overhead (Section~\ref{sec:pi-sampling-illustration}; 
Sections~\ref{sec:algorithm-illustration} and 
\ref{app:nips-results} in the Supplement), and is the main focus 
in this paper. We will also see that the former method has only a 
little computational overhead compared to the popular collapsed 
Gibbs sampling algorithm for LDA~\citep{GriffithsSteyvers:2004} 
(details, Section~\ref{sec:experiments}). Variational methods are 
often used in the topic modeling community as they give fast, 
parallel implementations by construction. Even though they converge 
rather quickly, our studies show that (which we discuss only in the 
Supplement) their solutions are suboptimal compared to 
the results of MCMC schemes we proposed.

The contributions of this paper are three-fold: (i) we propose a 
probabilistic, hierarchical model cLDA that can capture the topic 
structure of a corpus including organization hierarchy of 
documents, (ii) we discuss an efficient MCMC sampling scheme 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. 

The remainder of the paper is organized as follows. In 
Section~\ref{sec:hierarchical-model}, we formally define the 
proposed hierarchical model of cLDA. 
Section~\ref{sec:inference} develop the MCMC algorithm for 
posterior inference and evaluates correctness of the algorithms 
using a synthetically corpus. In Section~\ref{sec:experiments}, we 
assess the performance of the cLDA model and conclude that it 
exhibits superior predictive performance, quantitatively, using 
perplexity---a popular scheme for evaluating the performance of 
topic models. We will also provide a usability study for the cLDA 
model in this section. Section~\ref{sec:discussion} concludes this 
paper with a discussion of future research directions.  
}



\section*{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.



\bibliographystyle{plainnat}

\section{Posterior Sampling}
\label{sec:inference}

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$, \eqref{eq:cLDA-e}--\eqref{eq:cLDA-a} 
in the hierarchical model induce a 
prior distribution $p_h(\bm \psi)$ on $\bm \psi$. 
Equation~\eqref{eq:cLDA-a} 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
\begin{equation}
p_{h,{\bm w}}(\bm \psi) \propto \ell_{{\bm w}}(\bm \psi) p_h(\bm \psi) \\
\label{eq:posterior-psi}
\end{equation}
Using \eqref{eq:cLDA-e}--\eqref{eq:cLDA-a} of the 
hierarchical model, we can write 
\eqref{eq:posterior-psi} as (Section \ref{app:posterior} gives 
additional details)
\begin{equation}
\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] \\ 
	& \hspace{.5cm} \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}
\label{eq:posterior-psi-exp}
\end{equation}
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  
\cite{GriffithsSteyvers:2004}) and variational methods (e.g. 
see~\cite{BleiNgJordan: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 
\cite{TehNewmanWelling: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~\ref{app:vem}, to stay our discussion 
focused. 

\subsection{Inference via Markov chain Monte Carlo Methods}
\label{sec:mcmc}


According to the hierarchical model 
\eqref{eq:cLDA-e}--\eqref{eq:cLDA-a}, $\bm \theta_{jd}$'s and 
$\bm \beta_{k}$'s are independent, and by 
inspecting the posterior~\eqref{eq:posterior-psi}, given 
$(\bm \pi, {\bm z})$, we get: 
\begin{equation}
\begin{split}
\bm \theta_{jd} & \sim \Dir_K \biggl( n_{jd1} + \gamma \pi_{j1}, 
\ldots, n_{jdK} + \gamma \pi_{jK} \biggr),  \\ 
\bm \beta_{k} & \sim \Dir_V \biggl( m_{k1} + \eta, \ldots, m_{kV} 
+ \eta \biggr), 
\end{split}
\label{eq:theta-beta-posterior}
\end{equation}
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 \eqref{eq:theta-beta-posterior}
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  
\begin{eqnarray}
  \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]  
  \\ & \hspace{4.6mm}
	\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}
	\label{eq:pi-z-posterior}
\end{eqnarray}

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~\eqref{eq:pi-z-posterior}, we 
obtain a closed form expression for the conditional posterior 
distribution for $z_{jdi}$, given ${\bm z}^{(-jdi)}$ and $\pi_{jk}$, 
as 
\begin{equation}
p_{h,{\bm w}} \left(z_{jdik} = 1 \,|\, . \right) 
\propto  
\frac{\gamma \pi_{jk} + n^{(-jdi)}_{jdk}}
{\gamma + n^{(-jdi)}_{jd}} 
\frac{\eta  + m^{(-jdi)}_{kv}}
{V\eta + m^{(-jdi)}_{k}}	
\label{eq:posterior-z}
\end{equation}
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 \ref{app:posterior-z}) 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  
\begin{equation}
\tilde{p}_{{\bm w}_j}(\bm \pi_j \,|\, {\bm 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}
\label{eq:posterior-pi}
\end{equation} 
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. 

\subsubsection{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~\citep{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~\citep{Abramowitz:1974}
\begin{equation}
\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,
\label{eq:gammafn-sterling-exp}
\end{equation}
into the marginal posterior density 
\eqref{eq:posterior-pi} as  
\begin{equation}
\tilde{p}_{{\bm w}_j}(\bm \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}
\label{eq:posterior-pi-augmented}
\end{equation}
(\cite{Newman:2009} and \cite{TehEtal:2006}  used a 
similar idea in a different hierarchical 
model.) The 
expression~\eqref{eq:posterior-pi-augmented} 
introduces auxiliary variable $s_{jdk}$. 
By 
inspecting~\eqref{eq:posterior-pi-augmented}, 
we get a closed form 
expression for sampling $\bm \pi_j$, for $j = 1, 
\ldots, J$:   
\begin{equation}
\bm \pi_j \sim \Dir_{K} \left (\sum_{d = 1}^{D_j} s_{jd1} + \alpha, 
\ldots, \sum_{d = 1}^{D_j} s_{jdK} + \alpha 
\right),
\label{eq:aux-pi}
\end{equation}
and we update the auxiliary variable $s_{jdk}$ 
by the Antoniak sampling scheme 
\citep[Appendix A]{Newman:2009}, i.e., the 
Chinese restaurant process~\citep[CRP]{Aldous:1985} 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
\begin{equation}
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)
\end{equation}
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 
\ref{alg:gibbs-aux-pi-z}. We use AGS to denote this chain. 

\begin{algorithm}[t]
 \KwData{Observed words ${\bm w}$ and document metadata}
 \KwResult{A Markov chain on $(\bm \pi, {\bm z})$} \vspace{.3cm}
 initialize $(\bm \pi^{(0)}, {\bm z}^{(0)})$\;
 \For{Gibbs iteration $t$}{ 
 
	\tcp{Sampling word topic indices} 
	
	\For{word $w_{jdi}$, $i = 1, \ldots, n_{jd}$, $d = 1, \ldots, D_j$, $j = 1, \ldots, J$}{
		given $\bm \pi^{(t)}_j$, sample $z^{(t + 1)}_{jdi}$ via $p_{{\bm w}}(z_{jdi} 
		\,|\, {\bm z}^{(-jdi)}, \bm \pi^{(t)}_j)$ given by
		\eqref{eq:posterior-z}\;
		update count statistics $n_{jdk}$, $n_{jd}$, $m_{kv}$, 
		and $m_{k}$, according to $z^{(t + 1)}_{jdi}$\;
	}

	\tcp{Auxiliary variable sampling for collection-level topic mixtures}
	\For{collection $j = 1, \ldots, J$}{
		given $(\bm \pi^{(t)}_j, {\bm z}^{(t + 1)})$, update $s_{jdk}$ via 
		the Antoniak sampling scheme\;
		given $s_{jdk}$, sample $\bm \pi_j$ via~\eqref{eq:aux-pi}\;
		discard $s_{jdk}$ and update $\bm \pi^{(t + 1)}_j$ as $\bm \pi_j$\;
	}

 }
 \caption{Augmented Gibbs sampler (AGS)}
\label{alg:gibbs-aux-pi-z}
\end{algorithm}

\subsubsection{Metropolis Adjusted Langevin Monte Carlo} 
\label{sec:mgs}
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 \citep{MetropolisEtal:1953, 
Hastings:1970} with Langevin dynamics~\citep{Girolami:2011}.
This scheme has become popular for sampling on the simplex 
recently~\citep{Patterson: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 
\begin{equation}
a (\bm \pi^{(t)}, \bm \pi^{(*)}) = \text{min} \left(1, 
\frac{\tilde{p}(\bm \pi^{(*)}) q(\bm \pi^{(t)} \leftarrow \bm \pi^{(*)})} 
{\tilde{p}(\bm \pi^{(t)}) q(\bm \pi^{(*)} \leftarrow \bm \pi^{(t)}) } 
\right)
\label{eq:mh-accept-ration}
\nonumber 
\end{equation}
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~\citep{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) 
\begin{equation}
\text{d}\bm \pi(t) = \frac{1}{2} \nabla_{\bm \pi} {\cal L}(\bm \pi^{(t)}) 
\text{d}t + \text{d}{\bm b}(t)   
\label{eq:langevin-sde}
\end{equation} 
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~\eqref{eq:langevin-sde} as 
\begin{eqnarray}
\bm \pi^{(*)} &=& {\bm \mu}(\bm \pi^{(t)}, \varepsilon) + \varepsilon 
{\bm \xi}^{(t)}, \label{eq:langevin-discrete-eq} \\
{\bm \mu}(\bm \pi^{(t)}, \varepsilon) &=& \bm \pi^{(t)} + 
\frac{\varepsilon^2}{2} \nabla_{\bm \pi} {\cal L}(\bm \pi^{(t)}) 
\nonumber\\
{\bm \xi}^{(t)} &\stackrel{\text{iid}}{\sim}& \mathcal{N}_K(0, \mathds{1}_K) \nonumber
\end{eqnarray}
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~\eqref{eq:langevin-discrete-eq} 
may be 
transient and is no longer reversible with respect to the 
stationary density $p(\bm \pi)$~\citep{RobertsTweedie:1996}. 
However, one can ensure convergence to
the target density $p(\bm \pi)$ via a MH accept-reject scheme 
\citep[p. 591]{Besag:1994}, which we denote by MALA (Section 
\ref{app:mmala}). One may notice this MALA update as a 
special case of Hamiltonian Monte Carlo~\citep[Section 
5.5.2]{Neal:2010}. 

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~\citep{Girolami: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 
\citep{Patterson: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, \cite[Section 
5]{Girolami:2011} 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 \eqref{eq:langevin-sde} 
(details in Section \ref{app:mmala}). We thus use 
the corresponding Riemannian Manifold Metropolis Adjusted Langevin 
Algorithm \citep[MMALA]{Girolami:2011} 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$~\citep{Patterson: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 
\begin{equation}
\begin{split}
|\varphi_{k}| &\stackrel{\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}
\end{equation}
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~\eqref{eq:posterior-pi} on $\bm \pi$ in terms of 
${\bm \varphi}$, and derive the MMALA updates, accordingly 
(Section \ref{app:mmala}). The Markov 
chain on $(\bm \pi, {\bm z})$ induced by this scheme is denoted by MGS 
(Algorithm \ref{alg:gibbs-mmala-pi-z}).  

\begin{algorithm}[t!]
	\KwData{Observed words ${\bm w}$ and document metadata}
	\KwResult{A Markov chain on $(\bm \pi, {\bm z})$} \vspace{.3cm}
	initialize $(\bm \pi^{(0)}, {\bm z}^{(0)})$\;
	initialize ${\bm \varphi}^{(0)}$, i.e., the re-parametrization for 
	$\bm \pi^{(0)}$\;
	\For{Gibbs iteration $t$}{ 
		
		\tcp{Sampling word topic indices} 
		
		\For{word $w_{jdi}$, $i = 1, \ldots, n_{jd}$, $d = 1, 
		\ldots, 
			D_j$, $j = 1, \ldots, J$ }{
			given $\bm \pi^{(t)}_j$, sample $z^{(t + 1)}_{jdi}$ via 
			$p(z_{jdi} 
			\,|\, {\bm z}^{(-jdi)}, {\bm w}, \bm \pi^{(t)}_j)$ given by
			\eqref{eq:posterior-z}\;
			update count statistics $n_{jdk}$, $n_{jd.}$, 
			$m_{..kv}$, 
			and $m_{..k.}$, according to $z^{(t + 1)}_{jdi}$\;
		}
		
		\tcp{Metropolis Hastings updates for collection-level 
		topic mixtures}
		
		\For{collection $j = 1, \ldots, J$}{
			propose the MMALA update ${\bm \varphi}_j^{(*)}$ 
			via~\eqref{eq:riemann-langevin-euler-reparam}\;
			calculate the acceptance ratio 
			$\text{a}({\bm \varphi}_j^{(t)}, 
			{\bm \varphi}_j^{(*)})$, based on the unnormalized density 
			\eqref{eq:posterior-pi-re-parametrization} and 
			the transition density 
			\eqref{eq:pi-transition-density-re-parameterization} \;
			
			\eIf{$\text{Uniform}(0, 1) < \text{min}(1, 
			a({\bm \varphi}_j^{(t)}, {\bm \varphi}_j^{(*)}))$ } { 
				set ${\bm \varphi}_j^{(t+1)} = {\bm \varphi}_j^{(*)}$ \tcp{MH 
				accept}
			}{
				set ${\bm \varphi}_j^{(t+1)} = {\bm \varphi}_j^{(t)}$ \tcp{MH 
				reject}
			}
			
			set $\pi_{jk} = 
			\frac{\left|\varphi^{(t+1)}_{jk}\right|}{\sum_{k=1}^{K}{\left|\varphi^{(t+1)}_{jk}
			 \right|}}, \, k = 1, 2, \ldots, K$\;
		}
	
		
	}
	\caption{MMALA updates within Gibbs sampler (MGS)}
	\label{alg:gibbs-mmala-pi-z}
\end{algorithm}

 

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~\eqref{eq:theta-beta-posterior}. 
The augmented chain will hold the  
convergence properties of the chain on $(\bm \pi, {\bm z})$.


\subsection{Empirical Evaluation of Samples $\bm \pi$}
\label{sec:pi-sampling-illustration}

We consider a synthetic corpus by simulating  
\eqref{eq:cLDA-e}--\eqref{eq:cLDA-a} 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 
\eqref{eq:cLDA-d} 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~\ref{alg:gibbs-mmala-pi-z} and 
~\ref{alg:vem}) 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 iterations\footnote{An implementation of all 
algorithms 
and datasets discussed in this paper is available as an R package 
at \url{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~\ref{tab:ags-mgs-vem} 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. 
\figurename~\ref{fig:pi-2-random-walk} 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~\ref{sec:algorithm-illustration}.

\begin{figure*}[t!] 
	\centering 
	\subfloat[AGS: $\bm \pi^{(2000)}_{2} = (.347, .015, .636)$]{
		\includegraphics[trim={0 1.2cm 0 
		1.2cm},clip,width=.45\linewidth]
		{perp-J2-K3-D200-V40-cs100-ds200-seed1983-a0d1-g1-e0d25-2016Jun06105546-aux-pi2-path}}
	\subfloat[MGS: $\bm \pi^{(2000)}_{2} = (.379, .005, .615)$]{
		\includegraphics[trim={0 1.2cm 0 
		1.2cm},clip,width=.45\linewidth]
		{perp-J2-K3-D200-V40-cs100-ds200-seed1983-a0d1-g1-e0d25-2016Jun06105546-mmala-pi2-path}}\\
	\subfloat[VEM: $\bm \pi^{(45)}_{2} = (.258, .155, .585)$]{
		\includegraphics[trim={0 1.2cm 0 
			1.2cm},clip,width=.45\linewidth]
		{perp-J2-K3-D200-V40-cs100-ds200-seed1983-a0d1-g1-e0d25-2016Jun06105546-vem-pi2-path}}
	\caption{Plots of values of $\bm \pi_{2}$ via algorithms AGS, MGS, 
	and VEM. With approximately $23$ iterations the AGS chain 
	reached the optimal region, i.e., $.07$ from the true value 
	$\bm \pi^{\text{true}}_{2} = (.584, .386, .030)$, but the MGS 
	chain took $139$ iterations to reach there. Algorithm VEM never 
	reached the optimal regions: at convergence, VEM hit points 
	that are $0.18$ far from $\bm \pi^{\text{true}}_{2}$.}
	\label{fig:pi-2-random-walk}
\end{figure*}

\begin{table}%
	\centering
	\begin{threeparttable}[b] 
		\caption{Estimated values of $\bm \pi$ via   	
			algorithms AGS, MGS, and VEM} 
		\label{tab:ags-mgs-vem}
		\begin{tabular}{c c c c}
			\toprule 
			Method & 
			$\hat{\bm \pi}_{1}$\tnote{$\dagger$} &
			$\hat{\bm \pi}_{2}$\tnote{$\dagger$} &
			Iterations \\  
			\midrule
			AGS & 
			$(\epsilon, .996, .003)$ & 
			$(.347, .015, .636)$ & 	
			$2{,}000$ \\ 
			MGS & 
			$(.001, .997, \epsilon)$ & 
			$(.379, .005, .615)$ & 	
			$2{,}000$ \\ 
			VEM & 
			$(.057, .935, .006)$ & 
			$(.258, .155, .585)$ & 	
			$45$ \\ 
			\bottomrule 
		\end{tabular}
		\begin{tablenotes} 
			\item [$\dagger$] $\bm \pi^{\text{true}}_{1} = (.002, 
			\epsilon, .997)$, $\bm \pi^{\text{true}}_{2} = (.584, 
			.386, .030)$
		\end{tablenotes} 
	\end{threeparttable} 
\end{table}

\subsection{Selecting Hyperparameters in cLDA} 
\label{sec:clda-estimate-eta-gamma}
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} = \argmax_h m_{{\bm w}}(h)$. However, for models such as LDA 
and cLDA, the function $m_{{\bm w}}(h)$ is analytically intractable 
(Section \ref{sec:inference}). Works in the literature  
\citep{BleiNgJordan:2003,Wallach:2006} suggest an approximate EM 
algorithm to estimate $\argmax_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 
\citep{DempsterLairdRubin:1977} is a natural candidate to estimate 
$\argmax_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 
\citep{BleiNgJordan:2003}) or Markov chain Monte Carlo (e.g. 
Gibbs-EM \citep{Wallach:2006}). We follow the latter approach for 
cLDA, which includes: 
\begin{enumerate}
	\item Initialize $h_0$ and $\bm \psi_0$
	\item Until convergence  
	\begin{enumerate}
		\item[\textbf{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 
		\item[\textbf{M-step}] $h_{(t+1)} = \argmax_h \sum_{s = 
		1}^S \log p_{h_t}(\bm \psi_t^{(s)}, {\bm w})$ \label{e-step}
	\end{enumerate}
\end{enumerate}    

One solution to M-step is via the fixed point 
iteration \citep{minka2000estimating}. We derive the expressions for 
maximizing hyperparameter  $\gamma$ and $\eta$ as follows (see 
Section \ref{sec:est-eta-gamma}): 
\begin{equation}
\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]}
\label{eq:gem-est-eta} 
\end{equation} 
\begin{equation}
\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]}
\label{eq:gem-est-gamma} 
\end{equation}
where $\sum_{x, y, \ldots}$ represents $\sum_{x = 1}^X \sum_{y = 
1}^Y \ldots$. To illustrate the performance of the proposed method, 
\figurename~\ref{fig:clda-gem-eta-gamma-estimates-synth} 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. 



\begin{figure}[htbp!]
	\centering
	\includegraphics[width=.48\textwidth]
	{clda-ags-em-J2-K3-D200-V40-cs100-ds200-seed1983-a1-g0d8-e0d5-2017Dec22152501-gem-gamma-eta}
	\caption{Independent Gibbs-EM estimates ($\circ$) of 
	hyperparameters $(\eta, \gamma)$ in cLDA with constant 
	hyperparmeter $\alpha = 1$, for a synthetic corpora with 
	hyperparameters $\alpha_\text{true} = 1, 
	\gamma_\text{true} = .8, \eta_\text{true} 
	= .5$ (\textcolor{blue}{$\times$}). }
	\label{fig:clda-gem-eta-gamma-estimates-synth}
\end{figure}    
\section{Comparison of AGS, MGS, and VEM on a Synthetic Corpus}
\label{sec:algorithm-illustration}

This section gives additional details for the comparative study 
given in Section~\ref{sec:pi-sampling-illustration} of the main 
paper. 
\begin{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 
iterations\footnote{An implementation of all algorithms and 
datasets 
discussed in this paper is available as an R package at 
\url{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. 
\end{notes-to-stay} 
\figurename~\ref{fig:pi-1-random-walk} 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~\eqref{eq:theta-beta-posterior} 
(see Section \ref{sec:inference}). 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., \cite{GriffithsSteyvers: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 \ref{tab:ags-mgs-vem}).


\begin{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}$.
\end{notes-to-stay} 

\begin{figure*}[t!] 
\centering 
\subfloat[AGS: $\bm \pi^{(2000)}_{1} = (\epsilon, .996, .003)$]{
\includegraphics[width=.45\linewidth]
{perp-J2-K3-D200-V40-cs100-ds200-seed1983-a0d1-g1-e0d25-2016Jun06105546-aux-pi1-path}}
\subfloat[MGS: $\bm \pi^{(2000)}_{1} = (.001, .997, \epsilon)$]{
\includegraphics[width=.45\linewidth]
{perp-J2-K3-D200-V40-cs100-ds200-seed1983-a0d1-g1-e0d25-2016Jun06105546-mmala-pi1-path}}
\\ \subfloat[VEM: $\bm \pi^{(45)}_{1} = (.057, .935, .006)$]{
\includegraphics[width=.45\linewidth]
{perp-J2-K3-D200-V40-cs100-ds200-seed1983-a0d1-g1-e0d25-2016Jun06105546-vem-pi1-path}}
 
\caption{Plots of values of $\bm \pi_{1}$ via algorithms AGS, MGS, and 
VEM. Here, the variable $\epsilon$ denotes a small number. With 
approximately $42$ iterations the AGS chain reached the optimal 
region, i.e., $.003$ from the true value $\bm \pi^{\text{true}}_{1} = 
(.002, \epsilon, .997)$, but the MGS chain took $248$ iterations to 
reach there. Algorithm VEM never reached the optimal 
regions: at convergence, VEM hit points that are $0.08$ far 
from $\bm \pi^{\text{true}}_{1}$. See discussion in the text.}
\label{fig:pi-1-random-walk}
\end{figure*}



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  \figurename~\ref{fig:pi-1-random-walk} 
and~\ref{fig:pi-2-random-walk}. \figurename~\ref{fig:aux-mmala-acf} 
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.  


\begin{figure*}[t!] 
  \centering 
  \subfloat[$\pi_{11}$: AGS]{
  \includegraphics[width=.45\linewidth]
  {perp-J2-K3-D200-V40-cs100-ds200-seed1983-a0d1-g1-e0d25-2016May19132346-aux-acf-pi11}}
  \subfloat[$\pi_{22}$: AGS]{
  \includegraphics[width=.45\linewidth]
  {perp-J2-K3-D200-V40-cs100-ds200-seed1983-a0d1-g1-e0d25-2016May19132346-aux-acf-pi22}}
  \\ \subfloat[$\pi_{11}$: MGS]{
  \includegraphics[width=.45\linewidth]
  {perp-J2-K3-D200-V40-cs100-ds200-seed1983-a0d1-g1-e0d25-2016May19132346-mmala-acf-pi11}}
  \subfloat[$\pi_{22}$: MGS]{
  \includegraphics[width=.45\linewidth]
  {perp-J2-K3-D200-V40-cs100-ds200-seed1983-a0d1-g1-e0d25-2016May19132346-mmala-acf-pi22}}
  	
  \caption{Plots of ACF for samples $\pi_{11}$ and $\pi_{22}$  
	from the Markov chains MGS and AGS.}
	\label{fig:aux-mmala-acf}
\end{figure*}

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



