---
abstract: 'We present a new algorithm for identifying the transition and emission probabilities of a hidden Markov model (HMM) from the emitted data. Expectation-maximization becomes computationally prohibitive for long observation records, which are often required for identification. The new algorithm is particularly suitable for cases where the available sample size is large enough to accurately estimate second-order output probabilities, but not higher-order ones. We show that if one is only able to obtain a reliable estimate of the pairwise co-occurrence probabilities of the emissions, it is still possible to uniquely identify the HMM if the emission probability is *sufficiently scattered*. We apply our method to hidden topic Markov modeling, and demonstrate that we can learn topics with higher quality if documents are modeled as observations of HMMs sharing the same emission (topic) probability, compared to the simple but widely used bag-of-words model.'
bibliography:
- 'hmm\_refs.bib'
title: |
    **Learning Hidden Markov Models from Pairwise Co-occurrences\
    with Application to Topic Modeling**
---

Introduction
============

Hidden Markov models (HMMs) are widely used in machine learning when the data samples are time *dependent*, for example in speech recognition, language processing, and video analysis. The graphical model of a HMM is shown in Figure \[fig:hmm\]. HMM models a (time-dependent) sequence of data $\{Y_t\}_{t=0}^T$ as indirect observations of an underlying Markov chain $\{X_t\}_{t=0}^T$ which is not available to us. Homogeneous HMMs are parsimonious models, in the sense that they are fully characterized by the transition probability $\Pr[X_{t+1}|X_{t}]$ and the emission probability $\Pr[Y_t|X_t]$ even though the size of the given data $\{Y_t\}_{t=0}^T$ can be very large.

Consider a homogeneous HMM such that:

-   a latent variable $X_t$ can take $K$ possible outcomes $x_1,...,x_K$;

-   an ambient variable $Y_t$ can take $N$ possible outcomes $y_1,...,y_N$.

Recall that [@rabiner1986introduction; @Ghahramani2001]:

-   Given both $\{X_t\}_{t=0}^T$ and $\{Y_t\}_{t=0}^T$, the complete joint probability factors, and we can easily estimate the transition probability $\Pr[X_{t+1}|X_{t}]$ and the emission probability $\Pr[Y_t|X_t]$.

-   Given only $\{Y_t\}_{t=0}^T$, but assuming we know the underlying transition and emission probabilities, we can calculate the observation likelihood using the forward algorithm, estimate the most likely hidden sequence using the Viterbi algorithm, and compute the posterior probability of the hidden states using the forward-backward algorithm.

The most natural problem setting, however, is when neither the hidden state sequence nor the underlying probabilities are known to us—we only have access to a sequence of observations, and our job is to reveal the HMM structure, characterized by the transition matrix $\Pr[X_{t+1}|X_{t}]$ and the emission probability $\Pr[Y_t|X_t]$ from the set of observations $\{Y_t\}_{t=0}^T$.

Related work {#sec:related}
------------

The traditional way of learning a HMM from $\{Y_t\}_{t=0}^T$ is via expectation-maximization (EM) [@rabiner1986introduction], in which the expectation step is performed by calling the forward-backward algorithm. This specific instance of EM is also called the Baum-Welch algorithm [@baum1970maximization; @Ghahramani2001]. However, the complexity of Baum-Welch is prohibitive when $T$ is relatively large—the complexity of the forward-backward algorithm is ${\mathcal{O}}(K^2T)$, but EM converges slowly, so the forward-backward algorithm must be called many times. This is a critical issue, because a HMM can only be learned with high accuracy when the number of observation samples $T$ is large enough.

One way of designing scalable algorithms for learning HMMs is to work with sufficient statistics—a summary of the given observation sequence, whose size does not grow with $T$. Throughout this paper we assume that the HMM process is stationary (time-invariant), which is true almost surely if the underlying Markov process is ergodic and the process has been going on for a reasonable amount of time. With $T$ large enough, we can accurately estimate the co-occurrence probability between two consecutive emissions $\Pr[Y_t,Y_{t+1}]$. According to the graphical model shown in Figure \[fig:hmm\], it is easy to see that given the value of $X_t$, $Y_t$ is conditionally independent of all the other variables, leading to the factorization $$\begin{aligned}
\label{eq:hmm_fac2}
\Pr[Y_t,Y_{t+1}] 
= \sum_{k,j=1}^{K}&\Pr[Y_t|X_t=x_k]\Pr[Y_{t+1}|X_{t+1}=x_j]\Pr[X_t=x_k,X_{t+1}=x_j]\end{aligned}$$ Let ${\boldsymbol{\varOmega}}\in{\mathbf{R}}^{N\times N}$, ${\boldsymbol{M}}\in{\mathbf{R}}^{N\times K}$, and ${\boldsymbol{\varTheta}}\in{\mathbf{R}}^{K\times K}$, with their elements defined as $$\begin{aligned}
\varOmega_{n\ell} &= \Pr[Y_t=y_n,Y_{t+1}=y_\ell],\\
M_{nk} &= \Pr[Y_t=y_n|X_t=x_k],\\
\varTheta_{kj} &= \Pr[X_t=x_k,X_{t+1}=x_j].\end{aligned}$$ Then, equations can be written compactly as $$\begin{aligned}
\label{eq:MTM}
{\boldsymbol{\varOmega}}= {\boldsymbol{M}}{\boldsymbol{\varTheta}}{\boldsymbol{M}}^{{\!\top\!}}.\end{aligned}$$ Noticing that $\eqref{eq:MTM}$ is a nonnegative matrix tri-factorization with a number of inconsequential constraints for ${\boldsymbol{M}}$ and ${\boldsymbol{\varTheta}}$ to properly represent probabilities, @Vanluyten2008 [@Lakshminarayanan2010; @Cybenko2011] proposed using nonnegative matrix factorization (NMF) to estimate the HMM probabilities. However, NMF-based methods have a serious shortcoming in this context: the tri-factorization  is in general not unique, because it is fairly easy to find a nonsingular matrix ${\boldsymbol{Q}}$ such that both ${\boldsymbol{M}}{\boldsymbol{Q}}\geq0$ and ${\boldsymbol{Q}}^{-1}{\boldsymbol{\varTheta}}{\boldsymbol{Q}}^{-{{\!\top\!}}}\geq0$, and then $\widetilde{{\boldsymbol{M}}}={\boldsymbol{M}}{\boldsymbol{Q}}$ and $\widetilde{{\boldsymbol{\varTheta}}}={\boldsymbol{Q}}^{-1}{\boldsymbol{\varTheta}}{\boldsymbol{Q}}^{-{{\!\top\!}}}$ are equally good solutions in terms of reconstructing the co-occurrence matrix ${\boldsymbol{\varOmega}}$. When we use $({\boldsymbol{M}},{\boldsymbol{\varTheta}})$ and $(\widetilde{{\boldsymbol{M}}},\widetilde{{\boldsymbol{\varTheta}}})$ to perform HMM inference, such as estimating hidden states or predicting new emissions, the two models often yield completely different results, unless ${\boldsymbol{Q}}$ is a permutation matrix.

A number of works propose to use *tensor* methods to overcome the identifiability issue. Instead of working with the pairwise co-occurrence probabilities, they start by estimating the joint probabilities of three consecutive observations $\Pr[Y_{t-1},Y_t,Y_{t+1}]$. Noticing that these three random variables are conditionally independent given $X_t$, the triple-occurrence probability factors into $$\begin{aligned}
\Pr[Y_{t-1},Y_t,Y_{t+1}] = \sum_{k=1}^{K}\Pr[X_t=x_k]\Pr[Y_{t-1}|X_t=x_k]\Pr[Y_{t}|X_t=x_k]\Pr[Y_{t+1}|X_t=x_k],\end{aligned}$$ which admits a tensor *canonical polyadic decomposition* (CPD) model [@Hsu2009; @Anandkumar2012; @Anandkumar2014]. Assuming $K\leq N$, the CPD is essentially unique if two of the three factor matrices have full column rank, and the other one is not rank one [@harshman1970foundations]; in the context of HMMs, this is equivalent to assuming ${\boldsymbol{M}}$ and ${\boldsymbol{\varTheta}}$ both have linearly independent columns, which is a relatively mild condition. The CPD is known to be unique under much more relaxed conditions [@sidiropoulos2017tensor], but in order to uniquely retrieve the transition probability using the relationship $$\Pr[Y_{t+1}|X_{t}] = \sum_{j=1}^{K}
\Pr[Y_{t+1}|X_{t+1}\!=\!x_j]\Pr[X_{t+1}\!=\!x_j|X_t],$$ $K\leq N$ is actually the best we can achieve using triple-occurrences without making further assumptions. [^1] A salient feature in this case is that if the triple-occurrence probability $\Pr[Y_{t-1},Y_t,Y_{t+1}]$ is exactly given (meaning the rank of the triple-occurrence tensor is indeed smaller than $N$), the CPD can be efficiently calculated using generalized eigendecomposition and related algebraic methods [@sanchez1990tensorial; @leurgans1993decomposition; @DeLathauwer2004a]. These methods do not work well, however, when the low-rank tensor is perturbed; e.g., due to insufficient mixing / sample averaging of the triple occurrence probabilities.

It is also possible to handle cases where $K>N$. The key observation is that, given $X_t$, $Y_t$ is conditionally independent of $Y_{t-1},...,Y_{t-\tau}$ and $Y_{t+1},...,Y_{t+\tau}$. Then, grouping $Y_{t-1},...,Y_{t-\tau}$ into a single categorical variable taking $N^\tau$ possible outcomes, and $Y_{t+1},...,Y_{t+\tau}$ into another one, we can construct a much bigger tensor of size $N^\tau\times N^\tau\times N$, and then uniquely identify the underlying HMM structure with $K\gg N$ as long as certain linear independence requirements are satisfied for the conditional distribution of the *grouped* variables [@Allman2009; @bhaskara14a; @Huang2016; @Sharan2017]. It is intuitively clear that for fixed $N$, we need a much larger realization length $T$ in order to accurately estimate $(2\tau+1)$-occurrence probabilities as $\tau$ grows, which is the price we need to pay for learning a HMM with a larger number of hidden states.

This paper
----------

The focus of this paper is on cases where $K \leq N$, and $T$ is large enough to obtain accurate estimate of $\Pr[Y_t,Y_{t+1}]$, but not large enough to accurately estimate triple or higher-order occurrence probabilities. We [*prove*]{} that it is actually possible to recover the latent structure of an HMM only from pairwise co-occurrence probabilities $\Pr[Y_t,Y_{t+1}]$, provided that the underlying emission probability $\Pr[Y_t|X_t]$ is *sufficiently scattered*. Compared to the existing NMF-based HMM learning approaches, our formulation employs a different (determinant-based) criterion to ensure identifiability of the HMM parameters. Our matrix factorization approach resolves cases that cannot be handled by tensor methods, namely when $T$ is insufficient to estimate third-order probabilities, under an additional condition that is quite mild: that the emission probability matrix ${\boldsymbol{M}}$ must be *sufficiently scattered*, rather than simply full column-rank.

We apply our method to hidden topic Markov modeling (HTMM) [@gruber2007hidden], in which case the number of hidden states (topics) is indeed much smaller than the number of ambient states (words). HTMM goes beyond the simple and widely used bag-of-words model by assuming that (ordered) words in a document are emitted from a hidden topic sequence that evolves according to a Markov model. We show improved performance on real data when using this simple and intuitive model to take word ordering into account when learning topics, which also benefits from our identifiability guaranteed matrix factorization method.

As an illustrative example, we showcase the inferred topic of each word in a news article (removing stop words) in Figure \[fig:topic\], taken from the Reuters21578 data set obtained at [@reuters21578]. As we can see, HTMM gets much more consistent and smooth inferred topics compared to that obtained from a bag-of-words model (cf. supplementary material for details). This result agrees with human understanding.

Second-order vs. Third-order Learning
=====================================

We start by arguing that for the same observation data $\{Y_t\}_{t=0}^T$, the estimate of the pairwise co-occurrence probability $\Pr[Y_t,Y_{t+1}]$ is always more accurate than that of the triple co-occurrence probability $\Pr[Y_{t-1},Y_t,Y_{t+1}]$.

Let us first explicitly describe the estimator we use for these probabilities. For each observation $Y_t$, we define a coordinate vector ${\boldsymbol{\psi}}_t\in{\mathbf{R}}^K$, and ${\boldsymbol{\psi}}_t=\bm{e}_k$ if $Y_t=y_k$. The natural estimator for the pairwise co-occurrence probability matrix ${\boldsymbol{\varOmega}}$ is $$\label{eq:omega2}
\widehat{{\boldsymbol{\varOmega}}} = \frac{1}{T}\sum_{t=0}^{T-1}{\boldsymbol{\psi}}_t{\boldsymbol{\psi}}_{t+1}^{{\!\top\!}},$$ and similarly for the triple co-occurrence probability $\underline{{\boldsymbol{\varOmega}}_3}$ $$\label{eq:omega3}
\widehat{\underline{{\boldsymbol{\varOmega}}_3}} = \frac{1}{T-1}\sum_{t=1}^{T-1}{\boldsymbol{\psi}}_{t-1}\circ{\boldsymbol{\psi}}_t\circ{\boldsymbol{\psi}}_{t+1},$$ where $\circ$ denotes vector outer-product. [^2]

The first observation is that both $\widehat{{\boldsymbol{\varOmega}}}$ and $\widehat{\underline{{\boldsymbol{\varOmega}}_3}}$ are unbiased estimators: Obviously $\operatorname{E}({\boldsymbol{\psi}}_t{\boldsymbol{\psi}}_{t+1}^{{\!\top\!}})={\boldsymbol{\varOmega}}$ and likewise for the triple-occurrences, and taking their averages does not change the expectation. However, the individual terms in the summation are not independent of each other, making it hard to determine how fast estimates converge to their expectation. The state-of-the-art concentration result for HMMs [@kontorovich2006measure] states that for any 1-Lipschitz function $f$ $$\Pr[|f(\{Y_t\})-\operatorname{E}f(\{Y_t\})|>\epsilon]
\leq 2\exp\left(-T\epsilon^2/c\right),$$ where $c$ is a constant that only depends on the specific HMM structure but not on the function $f$ (cf. [@kontorovich2006measure] for details). Taking $f$ as any entry in $\widehat{{\boldsymbol{\varOmega}}}$ or $\widehat{\underline{{\boldsymbol{\varOmega}}_3}}$, we can check that indeed it is 1-Lipschitz, meaning as $T$ goes to infinity, both estimators converge to their expectation with negligible fluctuations.

We now prove that for a given set of observations $\{Y_t\}_{t=0}^T$, $\widehat{{\boldsymbol{\varOmega}}}$ is always going to be more accurate than $\widehat{\underline{{\boldsymbol{\varOmega}}_3}}$. Since both of them represent probabilities, we use two common metrics to measure the differences between the estimators and their expectations, the Kullback-Leibler divergence $D_\text{KL}(\cdot)$ and the total-variation difference $D_\text{TV}(\cdot)$.

\[ppst:2&gt;3\] Let $\widehat{{\boldsymbol{\varOmega}}}$ and $\widehat{\underline{{\boldsymbol{\varOmega}}_3}}$ be obtained from the same set of observations $\{Y_t\}_{t=0}^T$, we have that $$\begin{aligned}
D_\textup{KL}(\widehat{{\boldsymbol{\varOmega}}}\|{\boldsymbol{\varOmega}}) &\leq 
D_\textup{KL}(\widehat{\underline{{\boldsymbol{\varOmega}}_3}}\|\underline{{\boldsymbol{\varOmega}}_3})
\qquad\text{and}\\
D_\textup{TV}(\widehat{{\boldsymbol{\varOmega}}}\|{\boldsymbol{\varOmega}}) &\leq 
D_\textup{TV}(\widehat{\underline{{\boldsymbol{\varOmega}}_3}}\|\underline{{\boldsymbol{\varOmega}}_3}).\end{aligned}$$

The proof of Proposition \[ppst:2&gt;3\] is relegated to the supplementary material.

Identifiability of HMMs from Pairwise Co-occurrence Probabilities
=================================================================

The arguments made in the previous section motivate going back to matrix factorization methods for learning a HMM when the realization length $T$ is not large enough to obtain accurate estimates of triple co-occurrence probabilities. As we have explained in §\[sec:related\], the co-occurrence probability matrix ${\boldsymbol{\varOmega}}$ admits a nonnegative matrix tri-factorization model . There are a number of additional equality constraints. Columns of ${\boldsymbol{M}}$ represent conditional distributions, so ${\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{M}}={\boldsymbol{\mathit{1}}}^{{\!\top\!}}$. Matrix ${\boldsymbol{\varTheta}}$ represents the joint distribution between two consecutive Markovian variables, therefore ${\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{\varTheta}}{\boldsymbol{\mathit{1}}}=1$. Furthermore, we have that ${\boldsymbol{\varTheta}}{\boldsymbol{\mathit{1}}}$ and ${\boldsymbol{\varTheta}}^{{\!\top\!}}{\boldsymbol{\mathit{1}}}$ represent $\Pr[X_t]$ and $\Pr[X_{t+1}]$ respectively, and since we assume that the Markov chain is stationary, they are the same, i.e., ${\boldsymbol{\varTheta}}{\boldsymbol{\mathit{1}}}={\boldsymbol{\varTheta}}^{{\!\top\!}}{\boldsymbol{\mathit{1}}}$. Notice that this does not imply that ${\boldsymbol{\varTheta}}$ is symmetric, and in fact it is often not symmetric.

@huang2016nips considered a factorization model similar to in a different context, and showed that identifiability can be achieved under a reasonable assumption called *sufficiently scattered*, defined as follows.

\[def:suf\_scat\] Let $\operatorname{cone}({\boldsymbol{M}}^{{\!\top\!}})^*$ denote the polyhedral cone $\{{\boldsymbol{x}}:{\boldsymbol{M}}{\boldsymbol{x}}\geq 0\}$, and ${\mathcal{C}}$ denote the elliptical cone $\{{\boldsymbol{x}}:\|{\boldsymbol{x}}\|\leq{\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{x}}\}$. Matrix ${\boldsymbol{M}}$ is called **sufficiently scattered** if it satisfies that: (i) $\operatorname{cone}({\boldsymbol{M}}^{{\!\top\!}})^*\subseteq{\mathcal{C}}$, and (ii) $\operatorname{cone}({{\boldsymbol{M}}}^{{\!\top\!}})^\ast\cap{\rm bd}{\mathcal{C}}=\{\lambda {\bm e}_k:\lambda\geq 0,k=1,...,K\}$, where ${\rm bd}{\mathcal{C}}$ denotes the boundary of ${\mathcal{C}}$, $\{{\boldsymbol{x}}:\|{\boldsymbol{x}}\|={\boldsymbol{\mathit{1}}}^{{{\!\top\!}}}{\boldsymbol{x}}\}$.

The sufficiently scattered condition was first proposed in [@huang2014tsp] to establish uniqueness conditions for the widely used *nonnegative matrix factorization* (NMF). For the NMF model ${\boldsymbol{\varOmega}}=\bm{WH}^{{\!\top\!}}$, if both $\bm{W}$ and $\bm{H}$ are sufficiently scattered, then the nonnegative decomposition is unique up to column permutation and scaling. Follow up work strengthened and extended the identifiability results based on this geometry inspired condition [@fu2015bss; @huang2016nips; @fu2017spl]. A similar tri-factorization model was considered in [@huang2016nips] in the context of bag-of-words topic modeling, and it was shown that among all feasible solutions of , if we find one that minimizes $|\det{\boldsymbol{\varTheta}}|$, then it recovers the ground-truth latent factors ${\boldsymbol{M}}$ and ${\boldsymbol{\varTheta}}$, assuming the ground-truth ${\boldsymbol{M}}$ is sufficiently scattered. In our present context, we therefore propose the following problem formulation:

\[prob:main\] $$\begin{aligned}
\operatorname*{\textrm{minimize}}_{{\boldsymbol{\varTheta}},{\boldsymbol{M}}}~~~ & |\det{\boldsymbol{\varTheta}}| \\
\textrm{subject to}~~~ & {\boldsymbol{\varOmega}}={\boldsymbol{M}}{\boldsymbol{\varTheta}}{\boldsymbol{M}}^{{\!\top\!}}, \label{eq:Omega}\\
& {\boldsymbol{\varTheta}}\geq0, {\boldsymbol{\varTheta}}{\boldsymbol{\mathit{1}}}={\boldsymbol{\varTheta}}^{{\!\top\!}}{\boldsymbol{\mathit{1}}}, {\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{\varTheta}}{\boldsymbol{\mathit{1}}}=1, \label{eq:Theta}\\
& {\boldsymbol{M}}\geq0, {\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{M}}={\boldsymbol{\mathit{1}}}^{{\!\top\!}}. \end{aligned}$$

Regarding Problem , we have the following identifiability result.

\[thm:unique\] [@huang2016nips] Suppose ${\boldsymbol{\varOmega}}$ is constructed as ${\boldsymbol{\varOmega}}={\boldsymbol{M}}_\natural{\boldsymbol{\varTheta}}_\natural{\boldsymbol{M}}_\natural^{{\!\top\!}}$, where ${\boldsymbol{M}}_\natural$ and ${\boldsymbol{\varTheta}}_\natural$ satisfy the constraints in , and in addition (i) $\operatorname{rank}({\boldsymbol{\varTheta}}_\natural)=K$ and (ii) ${\boldsymbol{M}}_\natural$ is . Let $({\boldsymbol{M}}_\star,{\boldsymbol{\varTheta}}_\star)$ be an optimal solution for , then there must exist a permutation matrix ${\boldsymbol{\varPi}}\in{\mathbf{R}}^{K\times K}$ such that $${\boldsymbol{M}}_\natural = {\boldsymbol{M}}_\star{\boldsymbol{\varPi}}, \qquad
{\boldsymbol{\varTheta}}_\natural = {\boldsymbol{\varPi}}^{{\!\top\!}}{\boldsymbol{\varTheta}}_\star{\boldsymbol{\varPi}}.$$

One may notice that in [@huang2016nips], there are no constraints on the core matrix ${\boldsymbol{\varTheta}}$ as we do in . In terms of identifiability, it is easy to see that if the ground-truth ${\boldsymbol{\varTheta}}_\natural$ satisfies , solving  even without  will produce a solution ${\boldsymbol{\varTheta}}_\star$ that satisfies , thanks to uniqueness. In practice when we are given a less accurate ${\boldsymbol{\varOmega}}$, such “redundant” information will help us improve the estimation error, but that goes beyond identifiability consederations.

[.3]{}

[.3]{}

[.3]{}

The proof of Theorem \[thm:unique\] is referred to [@huang2016nips]. Here we provide some insights on this geometry-inspired sufficiently scattered condition, and discuss why it is a reasonable (and thus practical) assumption. The notation $\operatorname{cone}({\boldsymbol{M}}^{{\!\top\!}})^*=\{{\boldsymbol{x}}:{\boldsymbol{M}}{\boldsymbol{x}}\geq0\}$ comes from the convention in convex analysis that it is the *dual cone* of the conical hull of the row vectors of ${\boldsymbol{M}}$, i.e., $\operatorname{cone}({\boldsymbol{M}}^{{\!\top\!}})=\{{\boldsymbol{M}}^{{\!\top\!}}\bm{\alpha}:\bm{\alpha}\geq0\}$. Similarly, we can derive that the dual cone of ${\mathcal{C}}$ is ${\mathcal{C}}^*=\{{\boldsymbol{x}}:\|{\boldsymbol{x}}\|\leq{\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{x}}/\sqrt{K-1}\}$. A useful property of the dual cone is that for two convex cones $\mathcal{A}$ and $\mathcal{B}$, $\mathcal{A}\subseteq\mathcal{B}$ iff $\mathcal{B}^*\subseteq\mathcal{A}^*$. Therefore, the first requirement of sufficiently scattered in Definition \[def:suf\_scat\] equivalently means $${\mathcal{C}}^*\subseteq\operatorname{cone}({\boldsymbol{M}}^{{\!\top\!}}).$$ We give a geometric illustration of the sufficiently scattered condition in Figure \[fig:geo-scattered\] for $K=3$, and we focus on the 2-dimensional plane ${\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{x}}=1$. The intersection between this plane and the nonnegative orthant is the probability simplex, which is the triangle in Figure \[fig:geo-scattered\]. The outer circle represents ${\mathcal{C}}$, and the inner circle represents ${\mathcal{C}}^*$, again intersecting with the plane, respectively. The rows of ${\boldsymbol{M}}$ are scaled to sum up to one, and they are represented by black dots in Figure \[fig:geo-scattered\]. Their conical hull is represented by the shaded region. The polygon with dashed lines represents the dual of $\operatorname{cone}({\boldsymbol{M}}^{{\!\top\!}})$, which is indeed a subset of ${\mathcal{C}}$, and touches the boundary of ${\mathcal{C}}$ only at the coordinate vectors.

Figure \[fig:geo-separable\] shows a special case of sufficiently scattered called *separability*, which first appeared in [@donoho2004does] also to establish uniqueness of NMF. In this case, all the coordinate vectors appear in rows of ${\boldsymbol{M}}$, therefore $\operatorname{cone}({\boldsymbol{M}})$ equals the nonnegative orthant. It makes sense that this condition makes the identification problem easier, but it is also a very restrictive assumption. The sufficiently scattered condition, on the other hand, only requires that the shaded region contains the inner circle, as shown in Figure \[fig:geo-scattered\]. Intuitively this requires that the rows of ${\boldsymbol{M}}$ be “well scattered” in the probability simplex, but not to the extent of “separable”. Separability-based HMM identification has been considered in [@Barlier2015; @Glaude2015]. However, the way they construct second-order statistics is very different from ours. Figure \[fig:geo-no\_id\] shows a case where ${\boldsymbol{M}}$ is not sufficiently scattered, and it also happens to be a case where ${\boldsymbol{M}}$ is not identifiable.

As we can see, the elliptical cone ${\mathcal{C}}^*$ is tangent to all the facets of the nonnegative orthant. As a result, for ${\boldsymbol{M}}$ to be sufficiently scattered, it is necessary that there are enough rows of ${\boldsymbol{M}}$ lie on the boundary of the nonnegative orthant, i.e., ${\boldsymbol{M}}$ is relatively sparse. Specifically, if ${\boldsymbol{M}}$ is sufficiently scattered, then each column of ${\boldsymbol{M}}$ contains at least $K-1$ zeros [@huang2014tsp]. This is a very important insight, as exactly checking whether a matrix is sufficiently scattered may be computationally hard. In the present paper we further show the following result.

\[ppst:volume\] The ratio between the volume of the hyperball obtained by intersecting ${\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{x}}=1$ and ${\mathcal{C}}^*$ and the probability simplex is $$\label{eq:ratio}
\frac{1}{\sqrt{\pi K}}\left(\frac{4\pi}{K(K-1)}\right)^{\frac{K-1}{2}}
\Gamma\left(\frac{K}{2}\right).$$

The proof is given in the supplementary material. As $K$ grows larger, the volume ratio goes to zero at a super-exponential decay rate. This implies that the volume of the inner sphere quickly becomes negligible compared to the volume of the probability simplex, as $K$ becomes moderately large. The take home point is that, for a practical choice of $K$, say $K\geq10$, as long as ${\boldsymbol{M}}$ satisfies that each column contains at least $K$ zeros, and the positions of the zeros appear relatively random, it is very likely that it is sufficiently scattered, and thus can be uniquely recovered via solving .

Algorithm
=========

Our identifiability analysis based on the sufficiently scattered condition poses an interesting non-convex optimization problem . First of all, the given co-occurrence probability ${\boldsymbol{\varOmega}}$ may not be exact, therefore it may not be a good idea to put  as a hard constraint. For algorithm design, we propose the following modification to problem . $$\begin{aligned}
\label{prob:alg}
\operatorname*{\textrm{minimize}}_{{\boldsymbol{\varTheta}},{\boldsymbol{M}}}~~~ & 
\sum_{n,\ell=1}^{N}-\varOmega_{n\ell}\log\!\!\sum_{k,j=1}^{K}\!\!M_{nk}\varTheta_{kj}M_{\ell j} 
 + \lambda|\det{\boldsymbol{\varTheta}}| \nonumber\\
\textrm{subject to}~~~ & {\boldsymbol{M}}\geq0, {\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{M}}={\boldsymbol{\mathit{1}}}^{{\!\top\!}}, \\
& {\boldsymbol{\varTheta}}\geq0, {\boldsymbol{\varTheta}}{\boldsymbol{\mathit{1}}}={\boldsymbol{\varTheta}}^{{\!\top\!}}{\boldsymbol{\mathit{1}}}, {\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{\varTheta}}{\boldsymbol{\mathit{1}}}=1. \nonumber\end{aligned}$$ In the loss function of , the first term is the Kullback-Leibler distance between the empirical probability ${\boldsymbol{\varOmega}}$ and the parameterized version ${\boldsymbol{M}}{\boldsymbol{\varTheta}}{\boldsymbol{M}}^{{\!\top\!}}$ (ignoring a constant), and the second term is our identifiability-driven regularization. We need to tune the parameter $\lambda$ to yield good estimation results. However, intuitively we should use a $\lambda$ with a relatively small value. Suppose ${\boldsymbol{\varOmega}}$ is sufficiently accurate, then the priority is to minimize the difference between ${\boldsymbol{\varOmega}}$ and ${\boldsymbol{M}}{\boldsymbol{\varTheta}}{\boldsymbol{M}}^{{\!\top\!}}$; when there exist equally good fits, then the second term comes into play and helps us pick out a solution that is *sufficiently scattered*.

Noticing that the constraints of  are all convex, but not the loss function, we propose to design an iterative algorithm to solve  using successive convex approximation. At iteration $r$ when the updates are ${\boldsymbol{\varTheta}}^r$ and ${\boldsymbol{M}}^r$, we define $$\begin{aligned}
\label{eq:posterior}
\varPi_{n\ell kj}^r = M_{nk}^r\varTheta_{kj}^rM_{\ell j}^r \bigg/ 
\sum_{\kappa,\iota=1}^{K}M_{n\kappa}^r\varTheta_{\kappa\iota}^rM_{\ell\iota}^r.\end{aligned}$$ Obviously, $\varPi_{n\ell kj}^r\geq0$ and $\sum_{k,j=1}^{K}\varPi_{n\ell kj}^r=1$, which defines a probability distribution for fixed $n$ and $\ell$. Using Jensen’s inequality [@jensen1906fonctions], we have that $$\begin{aligned}
\label{eq:upper1}
-\varOmega_{n\ell}\log\sum_{k,j=1}^{K}M_{nk}\varTheta_{kj}M_{\ell j} \leq \sum_{k,j=1}^{K}-\varOmega_{n\ell}\varPi_{n\ell kj}^r
\left(\log M_{nk} + \log\varTheta_{kj} + \log M_{\ell j} - \log\varPi_{n\ell kj}^r\right)\end{aligned}$$ which defines a convex and locally tight upperbound for the first term in the loss function of . Regarding the second term in the loss of , we propose to simply take the linear approximation $$\begin{aligned}
\label{eq:approx2}
|\!\det\!{\boldsymbol{\varTheta}}| \approx |\!\det\!{\boldsymbol{\varTheta}}^r| + 
|\!\det\!{\boldsymbol{\varTheta}}^r|\operatorname{Tr}\!\left( ({\boldsymbol{\varTheta}}^r)^{\!-\!1\!}({\boldsymbol{\varTheta}}\!-\!{\boldsymbol{\varTheta}}^r) \right) \end{aligned}$$

Combining  and , our successive convex approximation algorithm tries to solve the following convex problem at iteration $r$: $$\begin{aligned}
\label{prob:iter}
\operatorname*{\textrm{minimize}}_{{\boldsymbol{\varTheta}},{\boldsymbol{M}}}~~ & \sum_{n,\ell=1}^{N}\sum_{k,j=1}^{K}
-\varOmega_{n\ell}\varPi_{n\ell kj}^r \left(\log M_{nk} + \log M_{\ell j} + \log\varTheta_{kj} \right) 
+ \lambda\sum_{k,j=1}^{K}\varXi_{kj}^r\varTheta_{kj} \\
\textrm{subject to}~~ & {\boldsymbol{M}}\geq0, {\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{M}}={\boldsymbol{\mathit{1}}}^{{\!\top\!}}, \nonumber \\
& {\boldsymbol{\varTheta}}\geq0, {\boldsymbol{\varTheta}}{\boldsymbol{\mathit{1}}}={\boldsymbol{\varTheta}}^{{\!\top\!}}{\boldsymbol{\mathit{1}}}, {\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{\varTheta}}{\boldsymbol{\mathit{1}}}=1,  \nonumber\end{aligned}$$ where we define $\bm{\varXi}^r=|\det{\boldsymbol{\varTheta}}^r|({\boldsymbol{\varTheta}}^r)^{-{{\!\top\!}}}$. Problem  decouples with respect to ${\boldsymbol{M}}$ and ${\boldsymbol{\varTheta}}$, so we can work out their updates individually.

The update of ${\boldsymbol{M}}$ admits a simple closed form solution, which can be derived via checking the KKT conditions. We denote the dual variable corresponding to ${\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{M}}={\boldsymbol{\mathit{1}}}^{{\!\top\!}}$ as ${\boldsymbol{\mu}}\in{\mathbf{R}}^K$. Setting the gradient of the Lagrangian with respect to $M_{nk}$ equal to zero, we have $$M_{nk} = 
\sum_{\ell=1}^N\sum_{j=1}^{K}\left(\varOmega_{n\ell}\varPi_{n\ell kj}^r+\varOmega_{\ell n}\varPi_{\ell njk}^r\right) \bigg/ \mu_k$$ and ${\boldsymbol{\mu}}$ should be chosen so that the constraint ${\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{M}}={\boldsymbol{\mathit{1}}}^{{\!\top\!}}$ is satisfied, which amounts to a simple re-scaling.

The update of ${\boldsymbol{\varTheta}}$ is not as simple as a closed form expression, but it can still be obtained very efficiently. Noticing that the nonnegativity constraint is implicitly implied by the individual $\log$ functions in the loss function, we propose to solve it using Newton’s method with equality constraints [@boyd2004convex §10.2]. Although Newton’s method requires solving a linear system of equations with $K^2$ number of variables in each iteration, there is special structure we can exploit to reduce the per-iteration complexity down to $\mathcal{O}(K^3)$: The Hessian of the loss function of  is diagonal, and the linear equality constraints are highly structured; using block elimination [@boyd2004convex §10.4.2], we ultimately only need to solve a positive definite linear system with $K$ variables. Together with the quadratic convergence rate of Newton’s method, the complexity of updating ${\boldsymbol{\varTheta}}$ is $\mathcal{O}(K^3\log\log\frac{1}{\varepsilon})$, where $\varepsilon$ is the desired accuracy for the ${\boldsymbol{\varTheta}}$ update. Noticing that the complexity of a naive implementation of Newton’s method would be $\mathcal{O}(K^6\log\log\frac{1}{\varepsilon})$, the difference is big for moderately large $K$. The in-line implementation of this tailored Newton’s method <span style="font-variant:small-caps;">ThetaUpdate</span> and the detailed derivation can be found in the supplementary material.

$\lambda>0$ initialize ${\boldsymbol{M}}$ using [@huang2016nips] initialize ${\boldsymbol{\varTheta}}\leftarrow \frac{1}{K(K+1)}({\boldsymbol{I}}+{\boldsymbol{\mathit{1}}}\!{\boldsymbol{\mathit{1}}}^{{\!\top\!}})$ $\widetilde{{\boldsymbol{\varOmega}}}\leftarrow{\boldsymbol{\varOmega}}\big/{\boldsymbol{M}}{\boldsymbol{\varTheta}}{\boldsymbol{M}}^{{{\!\top\!}}}$ $\widetilde{{\boldsymbol{M}}}\leftarrow {\boldsymbol{M}}\ast \left( \widetilde{{\boldsymbol{\varOmega}}}{\boldsymbol{M}}{\boldsymbol{\varTheta}}^{{{\!\top\!}}} + \widetilde{{\boldsymbol{\varOmega}}}^{{\!\top\!}}{\boldsymbol{M}}{\boldsymbol{\varTheta}}\right)$ $\widetilde{{\boldsymbol{\varTheta}}}\leftarrow{\boldsymbol{M}}^{{{\!\top\!}}}\widetilde{{\boldsymbol{\varOmega}}}{\boldsymbol{M}}$ $\widetilde{{\boldsymbol{M}}}\leftarrow \widetilde{{\boldsymbol{M}}} \operatorname{Diag}({\boldsymbol{\mathit{1}}}^{{\!\top\!}}\widetilde{{\boldsymbol{M}}})^{-1}$ $\widetilde{{\boldsymbol{\varTheta}}}\leftarrow$ <span style="font-variant:small-caps;">ThetaUpdate</span> $({\boldsymbol{M}},{\boldsymbol{\varTheta}})\leftarrow$ Amijo line search between $({\boldsymbol{M}},{\boldsymbol{\varTheta}})$ and $(\widetilde{{\boldsymbol{M}}},\widetilde{{\boldsymbol{\varTheta}}})$ ${\boldsymbol{M}}$ and ${\boldsymbol{\varTheta}}$

The entire proposed algorithm to solve Problem  is summarized in Algorithm \[alg:hmm\_id\]. Notice that there is an additional line-search step to ensure decrease of the loss function. The constraint set of is convex, so the line-search step will not incur infeasibility. Computationally, we find that any operation that involves $\varPi_{n\ell kj}^r$ can be carried out succinctly by defining the intermediate matrix $\widetilde{{\boldsymbol{\varOmega}}}={\boldsymbol{\varOmega}}/{\boldsymbol{M}}{\boldsymbol{\varTheta}}{\boldsymbol{M}}^{{{\!\top\!}}}$, where “$/$” denotes element-wise division between two matrices of the same size. The per-iteration complexity of Algorithm \[alg:hmm\_id\] is completely dominated by the operations that involve computing with $\widetilde{{\boldsymbol{\varOmega}}}$, notably comparing with that of <span style="font-variant:small-caps;">Theta-Update</span>. In terms of initialization, which is important since we are optimizing a non-convex problem, we propose to use the method by @huang2016nips to obtain an initialization for ${\boldsymbol{M}}$; for ${\boldsymbol{\varTheta}}$, it is best if we start with a feasible point (so that the Newton’s iterates will remain feasible), and a simple choice is scaling the matrix ${\boldsymbol{I}}+{\boldsymbol{\mathit{1}}}\!{\boldsymbol{\mathit{1}}}^{{\!\top\!}}$ to sum up to one. Finally, we show that this algorithm converges to a stationary point of Problem , with proof relegated to the supplementary material based on [@razaviyayn2013unified].

Assume <span style="font-variant:small-caps;">ThetaUpdate</span> solves Problem  with respect to ${\boldsymbol{\varTheta}}$ exactly, then Algorithm \[alg:hmm\_id\] converges to a stationary point of Problem .

![The total variation difference between the ground truth and estimated transition probability (top) and emission probability (bottom). The total variation difference of the emission probabilities is calculated as $\frac{1}{2K}\|{\boldsymbol{M}}_\natural-{\boldsymbol{M}}_\star\|_1$, since each column of the matrices indicates a (conditional) probability, and the total variation difference is equal to one half of the $L_1$-norm; and similarly for that of the transition probabilities after rescaling the rows of ${\boldsymbol{\varOmega}}_\natural$ and ${\boldsymbol{\varOmega}}_\star$ to sum up to one. The result is averaged over 10 random problem instances.[]{data-label="fig:sim_exp1"}](sim_exp1){width="\linewidth"}

Validation on Synthetic Data
============================

In this section we validate the identifiability performance on synthetic data. In this case, the underlying transition and emission probabilities are generated synthetically, and we compare them with the estimated ones to evaluate performance. The simulations are conducted in MATLAB using the HMM toolbox, which includes functions to generate observation sequences given transition and emission probabilities, as well as an implementation of the Baum-Welch algorithm [@baum1970maximization], i.e., the EM algorithm, to estimate the transition and emission probabilities using the observations. Unfortunately, even for some moderate problem sizes we considered, the streamlined MATLAB implementation of the Baum-Welch algorithm was not able to execute within reasonable amount of time, so its performance is not included here. For the baselines, we compare with the plain NMF approach using multiplicative update [@Vanluyten2008] and the tensor CPD approach [@Sharan2017] using simultaneous diagonalization with Tensorlab [@tensorlab3.0]. Since we work with empirical distributions instead of exact probabilities, the result of the simultaneous diagonalization is not going to be optimal. We therefore use it to initialize the EM algorithm for fitting a nonnegative tensor factorization with KL divergence loss [@shashanka2008probabilistic] for refinement.

We focus on the cases when the number of hidden states $K$ is smaller than the number observed states $N$. As we explained in the introduction, even for this seemingly easier case, it is not known that we can guarantee unique recovery of the HMM parameters *just from the pair-wise co-occurrence probability*. What is known is that the tensor CPD approach is able to guarantee identifiability given exact triple-occurrence probability. We will demonstrate in this section that it is much harder to obtain accurate triple-occurrence probability comparing with the co-occurrence probability. As a result, if the sufficiently scattered assumption holds for the emission probability, the estimated parameters obtained from our method are always more accurate than those obtained from tensor CPD.

Fixing $N=100$ and $K=20$, the transition probabilities are synthetically generated from a random exponential matrix of size $K\times K$ followed by row-normalization; for the emission probabilities, approximately 50% of the entries in the $N\times K$ random exponential matrices are set to zero before normalizing the columns, which is shown to satisfy the sufficiently scattered condition with very high probability [@huang2015principled]. We let the number of HMM realizations go from $10^6$ to $10^8$, and compare the estimation error for the transition matrix and emission matrix by the aforementioned methods. We show the total variation distance between the ground truth probabilities $\Pr[X_{t+1}|X_t]$ and $\Pr[Y_t|X_t]$ and their estimations $\widehat{\Pr}[X_{t+1}|X_t]$ and $\widehat{\Pr}[Y_t|X_t]$ using various methods. The result is shown in Figure \[fig:sim\_exp1\]. As we can see, the proposed method indeed works best, obtaining almost perfect recovery when sample size is above $10^8$. The CPD based method does not work as well since it cannot obtain accurate estimates of the third-order statistics that it needs. Initialized by CPD, EM improves from CPD but the performance is still far away from the proposed method. NMF is not working well since it does not have identifiability in this case.

Application: Hidden Topic Markov Model
======================================

Analyzing text data is one of the core application domains of machine learning. There are two prevailing approaches to model text data. The classical bag-of-words model assumes that each word is *independently* drawn from certain multinomial distributions. These distributions are different across documents, but can be efficiently summarized by a small number of *topics*, again mathematically modeled as distributions over words; this task is widely known as *topic modeling* [@hofmann2001unsupervised; @blei2003latent]. However, it is obvious that the bag-of-words representation is oversimplified. The $n$-gram model, on the other hand, assumes that words are conditionally dependent up to a window-length of $n$. This seems to be a much more realistic model, although the choice of $n$ is totally unclear, and is often dictated by memory and computational limitations in practice—since the size of the joint distribution grows exponentially with $n$. What is more, it is somewhat difficult to extract “topics” from this model, despite some preliminary attempts [@wallach2006topic; @wang2007topical].

We propose to model a document as the realization of a HMM, in which the topics are hidden states emitting words, and the states are evolving according to a Markov chain, hence the name *hidden topic Markov model* (HTMM). For a set of documents, this means we are working with a *collection* of HMMs. Similar to other topic modeling works, we assume that the topic matrix is shared among all documents, meaning all the given HMMs share the same emission probability. For the bag-of-words model, each document has a specific topic distribution ${\boldsymbol{p}}_d$, whereas for our model, each document has its own *topic transition probability* ${\boldsymbol{\varTheta}}_d$; as per our previous discussion, the row-sum and column-sum of ${\boldsymbol{\varTheta}}_d$ are the same, which are also the topic probability for the specific document. The difference is the Markovian assumption on the topics rather than the over-simplifying independence assumption.

We can see some immediate advantages for the HTMM. Since the Markovian assumption is only imposed on the topics, which are not exposed to us, the observations (words) are not independent from each other, which agrees with our intuition. On the other hand, we now understand that although word dependencies exist for a wide neighborhood, we only need to work with pair-wise co-occurrence probabilities, or 2-grams. This releases us from picking a window length $n$ in the $n$-gram model, while maintaining dependencies between words well beyond a neighborhood of $n$ words. It also includes the bag-of-words assumption as a special case: If the topics of the words are indeed independent, this just means that the transition probability has the special form ${\boldsymbol{\mathit{1}}}{\boldsymbol{p}}_d^{{\!\top\!}}$. The closest work to ours is by @gruber2007hidden, which is also termed hidden topic Markov model. However, they make a simplifying assumption that the transition probability takes the form $(1-\epsilon){\boldsymbol{I}}+ \epsilon{\boldsymbol{\mathit{1}}}{\boldsymbol{p}}_d^{{\!\top\!}}$, meaning the topic of the word is either the same as the previous one, or independently drawn from ${\boldsymbol{p}}_d$. Both models are special cases of our general HTMM.

In order to learn the shared topic matrix ${\boldsymbol{M}}$, we can use the co-occurrence statistics for the entire corpus: Denote the co-occurrence statistics for the $d$-th document as ${\boldsymbol{\varOmega}}_d$, then $\operatorname{E}{\boldsymbol{\varOmega}}_d = {\boldsymbol{M}}{\boldsymbol{\varTheta}}_d{\boldsymbol{M}}^{{\!\top\!}}$; consequently $${\boldsymbol{\varOmega}}= \frac{1}{\sum_{d=1}^{D}L_d}\sum_{d=1}^{D}L_d{\boldsymbol{\varOmega}}_d,$$ which is an unbiased estimator for $${\boldsymbol{M}}{\boldsymbol{\varTheta}}{\boldsymbol{M}}^{{\!\top\!}}= \frac{1}{\sum_{d=1}^{D}L_d}\sum_{d=1}^{D}L_d{\boldsymbol{M}}{\boldsymbol{\varTheta}}_d{\boldsymbol{M}}^{{\!\top\!}},$$ where $L_d$ is the length of the $d$-th document and ${\boldsymbol{\varTheta}}$ is conceptually a weighted average of all the topic-transition matrices. Then we may apply Algorithm \[alg:hmm\_id\] to learn the topic matrix. We illustrate the performance of our HTMM by comparing it to three popular bag-of-words topic modeling approaches: pLSA [@hofmann2001unsupervised], LDA [@blei2003latent], and FastAnchor [@arora2013practical], which guarantees identifiability if every topic has a characteristic *anchor word*. Our HTMM model guarantees identifiability if the topic matrix is *sufficiently scattered*, which is a more relaxed condition than the anchor word one. On the Reuters21578 data set obtained at [@reuters21578], we use the raw document to construct the word co-occurrence statistics, as well as bag-of-words representations for each document for the baseline algorithms. We use the version in which the stop-words have been removed, which makes the HTMM model more plausible since any syntactic dependencies have been removed, leaving only semantic dependencies. The vocabulary size of Reuters21578 is around $200,000$, making any method relying on triple-occurrences impossible to implement, and that is why tensor-based methods are not compared here.

![Coherence of the topics.[]{data-label="fig:reuters_coh"}](reuters_coh){width="0.55\linewidth"}

Because of page limitations, we only show the quality of the topics learned by various methods in terms of coherence. Simply put, a higher coherence means more meaningful topics, and the concrete definition can be found in [@arora2013practical] and in the supplementary material. In Figure \[fig:reuters\_coh\], we can see that for different number of topics we tried on the entire dataset, HTMM consistently produces topics with the highest coherence. Additional evaluations can be found in the supplementary material.

Conclusion
==========

We presented an algorithm for learning hidden Markov models in an unsupervised setting, i.e., using only a sequence of observations. Our approach is guaranteed to uniquely recover the ground-truth HMM structure using only pairwise co-occurrence probabilities of the observations, under the assumption that the emission probability is *sufficiently scattered*. Unlike EM, the complexity of the proposed algorithm does not grow with the length of the observation sequence. Compared to tensor-based methods for HMM learning, our approach only requires reliable estimates of pairwise co-occurrence probabilities, which are easier to obtain. We applied our method to topic modeling, assuming each document is a realization of a HMM rather than a simpler bag-of-words model, and obtained improved topic coherence results. We refer the reader to the supplementary material for detailed proofs of the propositions and additional experimental results.

Appendix {#appendix .unnumbered}
========

Proof of Proposition 1
======================

For categorical probabilities $\bm{p}$ and $\bm{q}$, their Kullback-Leiber divergence is defined as $$D_\textup{KL}(\bm{p}\|\bm{q}) = \sum_{n=1}^{N} p_n\log\frac{p_n}{q_n},$$ and their total variation distance is defined as $$D_\textup{TV}(\bm{p}\|\bm{q}) = \frac{1}{2}\sum_{n=1}^{N} |p_n-q_n|.$$

The key to prove Proposition 1 is the fact that the cooccurrence probability ${\boldsymbol{\varOmega}}$ can be obtained by marginalizing $X_{t-1}$ in the triple-occurrence probability $\underline{{\boldsymbol{\varOmega}}_3}$, i.e., $${\boldsymbol{\varOmega}}(i,j) = \sum_{n=1}^{N}\underline{{\boldsymbol{\varOmega}}_3}(n,i,j).$$ Similarly, this holds for the cumulative estimates described in §2 of the main paper as well, $$\widehat{{\boldsymbol{\varOmega}}}(i,j) = \sum_{n=1}^{N}\widehat{\underline{{\boldsymbol{\varOmega}}_3}}(n,i,j).$$

Using the log sum inequality, we have that $${\boldsymbol{\varOmega}}(i,j)\log\frac{{\boldsymbol{\varOmega}}(i,j)}{\widehat{{\boldsymbol{\varOmega}}}(i,j)} \leq
\sum_{n=1}^{N}\underline{{\boldsymbol{\varOmega}}_3}(n,i,j)
\log\frac{\underline{{\boldsymbol{\varOmega}}_3}(n,i,j)}{\widehat{\underline{{\boldsymbol{\varOmega}}_3}}(n,i,j)}.$$ Summing both sides over $i$ and $j$, we result in $$D_\textup{KL}(\widehat{{\boldsymbol{\varOmega}}}\|{\boldsymbol{\varOmega}}) \leq 
D_\textup{KL}(\widehat{\underline{{\boldsymbol{\varOmega}}_3}}\|\underline{{\boldsymbol{\varOmega}}_3})$$

Using Hölder’s inequality with $L_1$-norm and $L_\infty$-norm, we have that $$|{\boldsymbol{\varOmega}}(i,j)-\widehat{{\boldsymbol{\varOmega}}}(i,j)| \leq 
\sum_{n=1}^{N}|\underline{{\boldsymbol{\varOmega}}_3}(n,i,j) - \widehat{\underline{{\boldsymbol{\varOmega}}_3}}(n,i,j)|.$$ Summing both sides over $i$ and $j$ and then dividing by 2, we obtain $$D_\textup{TV}(\widehat{{\boldsymbol{\varOmega}}}\|{\boldsymbol{\varOmega}}) \leq 
D_\textup{TV}(\widehat{\underline{{\boldsymbol{\varOmega}}_3}}\|\underline{{\boldsymbol{\varOmega}}_3})$$ **Q.E.D.**

Proof of Proposition 2
======================

The volume of a hyper-ball in ${\mathbf{R}}^{n}$ with radius $R$ is $$\frac{\pi^{\frac{n}{2}}}{\Gamma(\frac{n}{2}+1)}R^{n}.$$ The elliptical cone ${\mathcal{C}}^*=\{ {\boldsymbol{x}}:\|{\boldsymbol{x}}\|\leq{\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{x}}/\sqrt{K-1} \}$ intersecting with the hyperplane ${\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{x}}=1$ is a hyperball in ${\mathbf{R}}^{K-1}$ with radius $\sqrt{\frac{1}{K(K-1)}}$. Therefore, the volume of the inner-ball is $$V_b = \frac{\pi^{\frac{K-1}{2}}}{\Gamma(\frac{K+1}{2})}(K(K-1))^{-\frac{K-1}{2}}.$$

The nonnegative orthan intersecting with ${\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{x}}=1$ is a regular simplex in ${\mathbf{R}}^{K-1}$ with side length $\sqrt{2}$. Its volume is $$V_s = \frac{\sqrt{K}}{(K-1)!} = \frac{\sqrt{K}}{\Gamma(K)}.$$

Their ratio is $$\begin{aligned}
\frac{V_b}{V_s} 
&= \frac{\frac{\pi^{\frac{K-1}{2}}}{\Gamma(\frac{K+1}{2})}(K(K-1))^{-\frac{K-1}{2}}}
{\frac{\sqrt{K}}{\Gamma(K)}} \\
&= \frac{1}{\sqrt{K}}\left(\frac{\pi}{K(K-1)}\right)^{\frac{K-1}{2}}
\frac{\Gamma(K)}{\Gamma(\frac{K+1}{2})}\\
&= \frac{1}{\sqrt{K}}\left(\frac{\pi}{K(K-1)}\right)^{\frac{K-1}{2}}
\frac{\Gamma(\frac{K}{2})}{2^{1-K}\sqrt{\pi}}\\
&= \frac{1}{\sqrt{\pi K}}\left(\frac{4\pi}{K(K-1)}\right)^{\frac{K-1}{2}}
\Gamma\left(\frac{K}{2}\right)\end{aligned}$$ **Q.E.D.**

This function of volume ratio is plotted in Figure \[fig:vol\_ratio\]. As we can see, as $K$ increases, the volume ratio indeed goes to zero at a super-exponential rate.

![The volume ratio between the hyperball obtained by intersecting $\mathcal{C}$ and the hyperplane ${\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{x}}=1$ and the probability simplex, as $K$ increases.[]{data-label="fig:vol_ratio"}](vol_ratio){width=".5\linewidth"}

Derivation of <span style="font-variant:small-caps;">ThetaUpdate</span>
=======================================================================

It is described in [@boyd2004convex §10.2] that for solving a convex equality constrained problem $$\begin{aligned}
\operatorname*{\textrm{minimize}}_{x}~ & f(x) \\
\text{subject to}~ & Ax=b\end{aligned}$$ using the Newton’s method, we start at a feasible point $x$, and the iterative update takes the form $x\leftarrow x - \alpha\Delta_\text{nt}x$, where the Newton direction is calculated from solving the KKT system $$\begin{bmatrix}
\nabla^2f(x) & A^{{\!\top\!}}~ \\ A & 0
\end{bmatrix}
\begin{bmatrix}
\Delta_\text{nt}x \\ d
\end{bmatrix} = 
\begin{bmatrix}
-\nabla f(x) \\ 0
\end{bmatrix}.$$

Assuming $\nabla^2f(x)\succ0$ and $A$ has full row rank, then the KKT system can be solved via elimination, as described in [@boyd2004convex Algorithm 10.3]. Suppose $A\in{\mathbf{R}}^{m\times n}$, if $\nabla^2f(x)$ is diagonal, the cost of calculating $\Delta_\text{nt}x$ is dominated by forming and inverting the matrix $ADA^{{\!\top\!}}$ with $D$ being diagonal.

Now we follow the steps of [@boyd2004convex Algorithm 10.3] to derive explicit Newton iterates for solving (11). First, we re-write the part of (11) that involve ${\boldsymbol{\varTheta}}$ here: $$\begin{aligned}
\operatorname*{\textrm{minimize}}_{{\boldsymbol{\varTheta}}>0}~~ & \sum_{n,\ell=1}^{N}\sum_{k,j=1}^{K}
-\varOmega_{n\ell}\varPi_{n\ell kj}^r \log\varTheta_{kj} + \lambda\sum_{k,j=1}^{K}\varXi_{kj}^r\varTheta_{kj} \\
\textrm{subject to}~~ & {\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{\varTheta}}{\boldsymbol{\mathit{1}}}=1, {\boldsymbol{\varTheta}}{\boldsymbol{\mathit{1}}}={\boldsymbol{\varTheta}}^{{\!\top\!}}{\boldsymbol{\mathit{1}}}.\end{aligned}$$

Let ${\boldsymbol{\theta}}=\operatorname{vec}({\boldsymbol{\varTheta}})$, then equality constraint has the form $\bm{A}{\boldsymbol{\theta}}=\bm{b}$ where $$\bm{A} = \begin{bmatrix}
{\boldsymbol{\mathit{1}}}^{{\!\top\!}}\otimes{\boldsymbol{\mathit{1}}}^{{\!\top\!}}\\
{\boldsymbol{\mathit{1}}}^{{\!\top\!}}\otimes{\boldsymbol{I}}- {\boldsymbol{I}}\otimes{\boldsymbol{\mathit{1}}}^{{\!\top\!}}\end{bmatrix}.$$ Matrix $\bm{A}$ does not have full row rank, because the last row of $\bm{A}$ is implied by the rest. Therefore, we can discard the last equality constraint. We will keep it when calculating matrix multiplications for simpler expression, and discard the corresponding entry or column/row for other operations.

Obviously $\bm{A\theta}$ has the form $$\bm{A\theta} =
\begin{bmatrix}
{\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{\varTheta}}{\boldsymbol{\mathit{1}}}\\
{\boldsymbol{\varTheta}}{\boldsymbol{\mathit{1}}}- {\boldsymbol{\varTheta}}^{{\!\top\!}}{\boldsymbol{\mathit{1}}}\end{bmatrix},$$ which costs $\mathcal{O}(K^2)$ flops. For a slightly more complicated multiplication $$\bm{A}\operatorname{Diag}({\boldsymbol{\theta}})\bm{A}^{{\!\top\!}}= \begin{bmatrix}
{\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{\varTheta}}{\boldsymbol{\mathit{1}}}& {\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{R}}^{{\!\top\!}}- {\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{\varTheta}}\\
~{\boldsymbol{\varTheta}}{\boldsymbol{\mathit{1}}}- {\boldsymbol{\varTheta}}^{{\!\top\!}}{\boldsymbol{\mathit{1}}}& \operatorname{Diag}({\boldsymbol{\varTheta}}{\boldsymbol{\mathit{1}}}+ {\boldsymbol{\varTheta}}^{{\!\top\!}}{\boldsymbol{\mathit{1}}}) - {\boldsymbol{\varTheta}}- {\boldsymbol{\varTheta}}^{{\!\top\!}}~ 
\end{bmatrix},$$ which also costs $\mathcal{O}(K^2)$ flops to compute. For $[~d_0~\bm{d}^{{\!\top\!}}~]^{{\!\top\!}}\in{\mathbf{R}}^{K+1}$, $$\bm{A}^{{\!\top\!}}[~d_0~\bm{d}^{{\!\top\!}}~]^{{\!\top\!}}= 
\operatorname{vec}\left(d_0{\boldsymbol{\mathit{1}}}{\boldsymbol{\mathit{1}}}^{{\!\top\!}}+ \bm{d}{\boldsymbol{\mathit{1}}}^{{\!\top\!}}- {\boldsymbol{\mathit{1}}}\bm{d}^{{\!\top\!}}\right).$$

At point ${\boldsymbol{\theta}}$, the negative gradient is $-\nabla f({\boldsymbol{\theta}}) = \operatorname{vec}(\bm{G})$ where $$G_{kj} = \frac{\sum_{n,\ell=1}^{N}\varOmega_{n\ell}\varPi_{n\ell kj}^r}{\varTheta_{kj}}
- \lambda\varXi_{kj}^r,$$ and the inverse of the Hessian $\left( \nabla^2f({\boldsymbol{\theta}}) \right)^{-1} = \operatorname{Diag}(\operatorname{vec}(\bm{R}))$ where $$R_{kj} = \frac{\varTheta_{kj}^2}
{\sum_{n,\ell=1}^{N}\varOmega_{n\ell}\varPi_{n\ell kj}^r}.$$

Let $${\boldsymbol{H}}= \begin{bmatrix}
{\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{R}}{\boldsymbol{\mathit{1}}}& {\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{R}}^{{\!\top\!}}- {\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{R}}\\
~{\boldsymbol{R}}{\boldsymbol{\mathit{1}}}- {\boldsymbol{R}}^{{\!\top\!}}{\boldsymbol{\mathit{1}}}& \operatorname{Diag}({\boldsymbol{R}}{\boldsymbol{\mathit{1}}}+ {\boldsymbol{R}}^{{\!\top\!}}{\boldsymbol{\mathit{1}}}) - {\boldsymbol{R}}- {\boldsymbol{R}}^{{\!\top\!}}~ 
\end{bmatrix}$$ and then delete the last column and row of ${\boldsymbol{H}}$, and $$S_{kj} = R_{kj} G_{kj}$$ $${\boldsymbol{g}}= \begin{bmatrix}
{\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{S}}{\boldsymbol{\mathit{1}}}\\ ~{\boldsymbol{S}}{\boldsymbol{\mathit{1}}}- {\boldsymbol{S}}^{{\!\top\!}}{\boldsymbol{\mathit{1}}}~
\end{bmatrix}$$ and then delete the last entry of ${\boldsymbol{g}}$. We can first solve for $\bm{d}$ by $$\bm{d} = {\boldsymbol{H}}^{-1}{\boldsymbol{g}}= [~d_0~\widetilde{\bm{d}}^{~{{\!\top\!}}}~]^{{\!\top\!}}.$$ Then we append a zero at the end of $\bm{d}$ and define $$[~\bm{d}^{{\!\top\!}}~0~]^{{\!\top\!}}\rightarrow \bm{d} = [~d_0~\widetilde{\bm{d}}^{\,{{\!\top\!}}}~]^{{\!\top\!}}.$$

The Newton direction $\Delta_\text{nt}{\boldsymbol{\theta}}$ can then be obtained via $$\Delta_\text{nt}{\boldsymbol{\theta}}= \left(\nabla^2f({\boldsymbol{\theta}})\right)^{-1}
\left( \bm{A}^{{\!\top\!}}\bm{d} + \nabla f({\boldsymbol{\theta}}) \right).$$

In matrix form, it is equivalent to $$\Delta_\text{nt}{\boldsymbol{\varTheta}}= {\boldsymbol{R}}\ast \left(d_0{\boldsymbol{\mathit{1}}}\!{\boldsymbol{\mathit{1}}}^{{\!\top\!}}+ \widetilde{\bm{d}}{\boldsymbol{\mathit{1}}}^{{\!\top\!}}- {\boldsymbol{\mathit{1}}}\widetilde{\bm{d}}^{~{{\!\top\!}}} - \bm{G}\right).$$

The in-line implementation of <span style="font-variant:small-caps;">ThetaUpdate</span> is given here.

${\boldsymbol{\varTheta}},\widetilde{{\boldsymbol{\varTheta}}},\lambda,\rho$ ${\boldsymbol{\varXi}}\leftarrow|\det{\boldsymbol{\varTheta}}|{\boldsymbol{\varTheta}}^{-{{\!\top\!}}}$ ${\boldsymbol{G}}\leftarrow \widetilde{{\boldsymbol{\varTheta}}}\big/{\boldsymbol{\varTheta}}- \lambda{\boldsymbol{\varXi}}$ ${\boldsymbol{R}}\leftarrow {\boldsymbol{\varTheta}}\ast{\boldsymbol{\varTheta}}\big/\widetilde{{\boldsymbol{\varTheta}}}$ $\displaystyle {\boldsymbol{H}}\leftarrow \begin{bmatrix}
{\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{R}}{\boldsymbol{\mathit{1}}}& {\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{R}}^{{\!\top\!}}\!-\! {\boldsymbol{\mathit{1}}}^{{\!\top\!}}{\boldsymbol{R}}\\
~{\boldsymbol{R}}{\boldsymbol{\mathit{1}}}\!-\! {\boldsymbol{R}}^{{\!\top\!}}{\boldsymbol{\mathit{1}}}& \operatorname{Diag}({\boldsymbol{R}}{\boldsymbol{\mathit{1}}}\!+\! {\boldsymbol{R}}^{{\!\top\!}}{\boldsymbol{\mathit{1}}}) \!-\! {\boldsymbol{R}}\!-\! {\boldsymbol{R}}^{{\!\top\!}}~ 
\end{bmatrix}$ delete the last column and row of ${\boldsymbol{H}}$ $\displaystyle {\boldsymbol{g}}\leftarrow\begin{bmatrix}
{\boldsymbol{\mathit{1}}}^{{\!\top\!}}({\boldsymbol{R}}\ast{\boldsymbol{G}}){\boldsymbol{\mathit{1}}}\\ ~({\boldsymbol{R}}\ast{\boldsymbol{G}}){\boldsymbol{\mathit{1}}}- ({\boldsymbol{R}}\ast{\boldsymbol{G}})^{{\!\top\!}}{\boldsymbol{\mathit{1}}}~
\end{bmatrix}$ delete the last entry of ${\boldsymbol{g}}$ $\bm{d} \leftarrow {\boldsymbol{H}}^{-1}{\boldsymbol{g}}$ $[~d_0~\widetilde{\bm{d}}^{~{{\!\top\!}}}~]^{{\!\top\!}}\leftarrow[~\bm{d}^{{\!\top\!}}~0~]^{{\!\top\!}}$ $
\Delta_\text{nt}{\boldsymbol{\varTheta}}= {\boldsymbol{R}}\ast \left(d_0{\boldsymbol{\mathit{1}}}\!{\boldsymbol{\mathit{1}}}^{{\!\top\!}}+ \widetilde{\bm{d}}{\boldsymbol{\mathit{1}}}^{{\!\top\!}}- {\boldsymbol{\mathit{1}}}\widetilde{\bm{d}}^{~{{\!\top\!}}} - \bm{G}\right)
$ ${\boldsymbol{\varTheta}}\leftarrow{\boldsymbol{\varTheta}}-\Delta_\text{nt}{\boldsymbol{\varTheta}}$ ${\boldsymbol{\varTheta}}$

Proof of Proposition 3
======================

The form of Algorithm 1 falls exactly into the framework of block successive convex approximation (BSCA) algorithm proposed by [@razaviyayn2013unified] with only one block of variables. Invoking [@razaviyayn2013unified Theorem 4], we have that every limit point of Algorithm 1 is a stationary point of Problem (7). Additionally, since the constraint set of Problem (7) is compact, *any* sub-sequence has a limit point, which is also a stationary point. This proves that Algorithm 1 converges to a stationary point of Problem (7). **Q.E.D.**

Additional Synthetic Experiments
================================

In this section we conduct a similar synthetic experiment to identify HMM parameters, but with a much smaller problem size, so that we can include the classical Baum-Welch algorithm [@baum1970maximization] as another baseline. Fixing $N=16$ and $K=4$, the transition probabilities are synthetically generated from a random exponential matrix of size $K\times K$ followed by row-normalization; for the emission probabilities, the top $K\times K$ part of the $N\times K$ random exponential matrices are set to be the identity matrix before column normalization, so that it is guaranteed to be sufficiently scattered. We let the number of HMM realizations go from $10^3$ to $10^5$, and compare the estimation error for the transition matrix and emission matrix by the aforementioned methods. We show the total variation distance between the ground truth probabilities $\Pr[X_{t+1}|X_t]$ and $\Pr[Y_t|X_t]$ and their estimations $\widehat{\Pr}[X_{t+1}|X_t]$ and $\widehat{\Pr}[Y_t|X_t]$ using various methods. The result is shown in Figure \[fig:suppl\_small\].

Similar to the experiment shown in the main paper, the proposed method works the best in terms of estimating the HMM parameters, without sacrificing too much computational times. Much to one’s surprise, the Baum-Welch algorithm is not working very well in terms of estimation error. This is possibly because we limit the maximum number of EM iterations to be 500 (default setting of the MATLAB implementation), which may not be enough for convergence. What is expected is that the computational time of Baum-Welch grows linearly with respect to the length of the HMM observations, while other methods are independent from it.

![The total variation difference between the ground truth and estimated transition probability (left) and emission probability (middle), and the elapsed time (right) for $N=16$ and $K=4$. The total variation difference of the emission probabilities is calculated as $\frac{1}{2K}\|{\boldsymbol{M}}_\natural-{\boldsymbol{M}}_\star\|_1$, since each column of the matrices indicates a (conditional) probability, and the total variation difference is equal to one half of the $L_1$-norm; and similarly for that of the transition probabilities after rescaling the rows of ${\boldsymbol{\varOmega}}_\natural$ and ${\boldsymbol{\varOmega}}_\star$ to sum up to one. The result is averaged over 10 random problem instances.[]{data-label="fig:suppl_small"}](suppl_small){width="1.3\linewidth"}

![The elapsed time for the synthetic experiment with $N=100$ and $K=20$ as in the main paper.[]{data-label="fig:suppl_time"}](suppl_time){width=".5\linewidth"}

An interesting remark is that when $T=12,800$, the per-iteration elapsed time of Baum-Welch is approximately 1 second. Recall that each iteration of Baum-Welch calls for the forward-backward algorithm, with complexity $\mathcal{O}(K^2T)$. This means for the problem size considered in the main paper, each iteration of Baum-Welch takes approximately 4 minutes to 7 hours, depending on the realization length. This is clearly not feasible in practice.

We also present the elapsed time of the four algorithms excluding the Baum-Welch algorithm for the case considered in the main paper, i.e., $N=100$ and $K=20$. Similar to the timing result shown in Figure \[fig:suppl\_small\], the proposed method takes the longest time compared to the other three, but not significantly; also recall that the propose method works considerably better in terms of estimation accuracy.

Additional HTMM Evaluations
===========================

In the main body of the paper we showed that HTMM is able to learn topics with higher quality using pairwise word cooccurrences. The quality of topics is evaluated using coherence, which is defined as follows. For each topic, a set of words $\mathcal{V}$ is picked (here we pick the top 20 words with the highest probability of appearing). We calculate the number of documents each word $v_1$ appears $\text{freq}(v_1)$ and the number of documents two words $v_1$ and $v_2$ both appear $\text{freq}(v_1,v_2)$. The coherence of that topic is calculated as $$\sum_{v_1,v_2\in\mathcal{V}} \log\left(\frac{\text{freq}(v_1,v_2)+\epsilon}{\text{freq}(v_1)}\right) .$$ The intuition is that if both $v_1$ and $v_2$ both have high probability of appearing in a topic, then they have high probability of co-occurring in a document as well; hence a higher value of coherence indicates a more indicative topic. The coherence of the individual topics are then averaged to get the coherence for the entire topic matrix.

![The perplexity of different models as number of topics $K$ increases.[]{data-label="fig:perp"}](reuters_perp){width=".5\linewidth"}

Here we show some more evaluation results. Using the learned topic matrix, we can see how it fits the data directly from *perplexity*, defined as [@blei2003latent] $$\exp\left( -\frac{\sum_d\log p(\text{doc}_d)}{\sum_d L_d} \right).$$ A smaller perplexity means the probability model fits the data better. As seen in Figure \[fig:perp\], HTMM gives the smallest perplexity. Notice that since HTMM takes word ordering into account, it is not fair for the other methods to take the bag-of-words representation of the documents. The bag-of-words model is essentially multinomial, whose pdf includes a scaling factor $\frac{n!}{n_1!...n_K!}$ for different combinations of observation orderings. In our case this factor is not included since we *do* know the word ordering in each document. For HTMM the log-likelihood is calculated efficiently using the forward algorithm.

This result is not surprising. Even using the same topic matrix, a bag-of-words model tries to find a $K$-dimensional representation for each document, whereas HTMM looks for a $K^2$-dimension representation. One may wonder if it is causing over-fitting, but we argue that it is not. First of all, we have see that in terms of coherence, HTMM learns a topic matrix with higher quality. For learning feature representations for each document, we showcase the following result. Once we have the topic-word probabilities and topic weights or topic transition probability, we can infer the underlying topic for each word. For bag-of-words models, each word only has one most probable topic in a document, no matter where it appears. For HTMM, once we learn the transition and emission probability, the topic of each word can be optimally estimated using the Viterbi algorithm. For one specific news article from the Reuters21578 data set, the topic inference given by pLSA is:

The word topic inference given by HTMM is:

As we can see, HTMM gets much more consistent and smooth inferred topics, which agrees with human understandings.

Learning HMMs from Triple-occurrences
=====================================

Finally, we show a stronger identifiability result for learning HMMs using triple-occurrence probabilities.

Consider a HMM with $K$ hidden states and $N$ observable states. Suppose the emission probability $\Pr[Y_t|X_t]$ is generic (meaning probabilities not satisfying this condition form a set with Lebesgue measure zero), the transition probabilities $\Pr[X_{t+1}|X_t]$ are linearly independent from each other, and each conditional probability $\Pr[X_{t+1}|X_t]$ contains no more than $N/2$ nonzeros. Then this HMM can be uniquely identified from its triple-occurrence probability $\Pr[Y_{t-1},Y_t,Y_{t+1}]$, up to permutation of the hidden states, for $K\leq \frac{N^2}{16}$.

It is clear that identifiability holds when $K\leq N$, so we focus on the case that $N<K\leq \frac{N^2}{16}$.

As we explained in §1.1, the triple-occurrence probability can be factored into $$\Pr[Y_{t-1},Y_t,Y_{t+1}] = \sum_{k=1}^{K}\Pr[X_t=x_k]\Pr[Y_{t-1}|X_t=x_k] 
 \Pr[Y_{t}|X_t=x_k]\Pr[Y_{t+1}|X_t=x_k].$$ Using tensor notations, this is equivalent to $$\underline{{\boldsymbol{\varOmega}}_3} = {\text{\textlbrackdbl}\bm{p};\bm{L},{\boldsymbol{M}},\bm{N}\text{\textrbrackdbl}},$$ where $$\begin{aligned}
p_k & = \Pr[X_t=x_k], \\
L_{nk} &= \Pr[Y_{t-1}=y_n|X_t=x_k], \\
N_{nk} &= \Pr[Y_{t+1}=y_n|X_t=x_k].\end{aligned}$$ Let $\overline{{\boldsymbol{\varTheta}}}$ denote the row scaled version of ${\boldsymbol{\varTheta}}$ so that each row sums to one, then $\overline{{\boldsymbol{\varTheta}}}$ denotes the transition probability. Then we have $$\label{eq:L=MT}
\bm{L}={\boldsymbol{M}}\overline{{\boldsymbol{\varTheta}}}^{{\!\top\!}}.$$

Since ${\boldsymbol{M}}$ is generic and $\overline{{\boldsymbol{\varTheta}}}$ is full rank, both $\bm{L}$ and $\bm{N}$ are generic as well. The latest tensor identifiability result by @chiantini2012generic [Theorem 1.1] shows that for a $N\times N\times N$ tensor with generic factors, the CPD $\underline{{\boldsymbol{\varOmega}}_3} = {\text{\textlbrackdbl}\bm{p};\bm{L},{\boldsymbol{M}},\bm{N}\text{\textrbrackdbl}}$ is essentially unique if $$K\leq 2^{ 2\lfloor\log_2N\rfloor - 2 },$$ or with a slightly worse bound $$K\leq \frac{N^2}{16}.$$

This does not mean that any non-singular $\overline{{\boldsymbol{\varTheta}}}$ can be uniquely recovered in this case. Equation  is under-determined. A natural assumption to achieve identifiability is that each row of $\overline{{\boldsymbol{\varTheta}}}$, i.e., each conditional transition probability $\Pr[X_{t+1}|X_t]$, can take at most $N/2$ non-zeros. In the context of HMM, this means that at a particular hidden state, there are only a few possible states for the next step, which is very reasonable. For a generic ${\boldsymbol{M}}$, $$\operatorname{spark}({\boldsymbol{M}}) = \operatorname{krank}({\boldsymbol{M}}) + 1 = N+1.$$ @Donoho2003 showed that for such a ${\boldsymbol{M}}$, and a vector $\bm{\theta}$ with at most $N/2$ nonzeros, $\bm{\theta}$ is the unique solution with at most $N/2$ nonzeros that satisfies ${\boldsymbol{M}}\bm{\theta} = \bm{\ell}$. Therefore, if we seek for the sparsest solution to the linear equation , we can uniquely recover ${\boldsymbol{\varTheta}}$ as well.

[^1]: In the supplementary material, we prove that if the emission probability is *generic* and the transition probability is *sparse*, the HMM can be uniquely identified from triple-occurrence probability for using the latest tensor identifiability result [@chiantini2012generic].

[^2]: In some literature $\circ$ is written as the Kronecker product $\otimes$. Strictly speaking, the Kronecker product of three vectors is a very long vector, not a three-way array. For this reason, we chose to use $\circ$ instead of $\otimes$.
