A Practical Algorithm for Topic Modeling with Provable Guarantees

Introduction

Topic modeling is a popular method that learns thematic structure from large document collections without human supervision. The model is simple: documents are mixtures of topics, which are modeled as distributions over a vocabulary [Blei, 2012]. Each word token is generated by selecting a topic from a document-specific distribution, and then selecting a specific word from that topic-specific distribution. Posterior inference over document-topic and topic-word distributions is intractable — in the worst case it is NP-hard even for just two topics [Arora et al., 2012b]. As a result, researchers have used approximate inference techniques such as singular value decomposition [Deerwester et al., 1990], variational inference [Blei et al., 2003], and MCMC [Griffiths & Steyvers, 2004].

Recent work in theoretical computer science focuses on designing provably efficient algorithms for topic modeling. These treat the topic modeling problem as one of statistical recovery: assuming the data was generated perfectly from the hypothesized model using an unknown set of parameter values, the goal is to recover the model parameters in polynomial time given a reasonable number of samples.

Arora et al. [2012b] present an algorithm that provably recovers the parameters of topic models provided that the topics meet a certain separability assumption [Donoho & Stodden, 2003]. Separability requires that every topic contains at least one anchor word that has non-zero probability only in that topic. If a document contains this anchor word, then it is guaranteed that the corresponding topic is among the set of topics used to generate the document. The algorithm proceeds in two steps: first it selects anchor words for each topic; and second, in the recovery step, it reconstructs topic distributions given those anchor words. The input for the algorithm is the second-order moment matrix of word-word co-occurrences.

Anandkumar et al. present a provable algorithm based on third-order moments that does not require separability, but, unlike the algorithm of Arora et al., assumes that topics are not correlated. Although standard topic models like LDA [Blei et al., 2003] assume that the choice of topics used to generate the document are uncorrelated, there is strong evidence that topics are dependent [Blei & Lafferty, 2007, Li & McCallum, 2007]: economics and politics are more likely to co-occur than economics and cooking.

Both algorithms run in polynomial time, but the bounds that have been proven on their sample complexity are weak and their empirical runtime performance is slow. The algorithm presented by Arora et al. [2012b] solves numerous linear programs to find anchor words. Bittorf et al. and Gillis reduce the number of linear programs needed. All of these algorithms infer topics given anchor words using matrix inversion, which is notoriously unstable and noisy: matrix inversion frequently generates negative values for topic-word probabilities.

In this paper we present three contributions. First, we replace linear programming with a combinatorial anchor selection algorithm. So long as the separability assumption holds, we prove that this algorithm is stable in the presence of noise and thus has polynomial sample complexity for learning topic models. Second, we present a simple probabilistic interpretation of topic recovery given anchor words that replaces matrix inversion with a new gradient-based inference method. Third, we present an empirical comparison between recovery-based algorithms and existing likelihood-based topic inference. We study both the empirical sample complexity of the algorithms on synthetic distributions and the performance of the algorithms on real-world document corpora. We find that our algorithm performs as well as collapsed Gibbs sampling on a variety of metrics, and runs at least an order of magnitude faster.

Our algorithm both inherits the provable guarantees from Arora et al. [2012a, b] and also results in simple, practical implementations. We view our work as a step toward bridging the gap between statistical recovery approaches to machine learning and maximum likelihood estimation, allowing us to circumvent the computational intractability of maximum likelihood estimation yet still be robust to model error.

Background

We consider the learning problem for a class of admixture distributions that are frequently used for probabilistic topic models. Examples of such distributions include latent Dirichlet allocation [Blei et al., 2003], correlated topic models [Blei & Lafferty, 2007], and Pachinko allocation [Li & McCallum, 2007]. We denote the number of words in the vocabulary by $V$ and the number of topics by $K$. Associated with each topic $k$ is a multinomial distribution over the words in the vocabulary, which we will denote as the column vector $A_{k}$ of length $V$. Each of these topic models postulates a particular prior distribution $\tau$ over the topic distribution of a document. For example, in latent Dirichlet allocation (LDA) $\tau$ is a Dirichlet distribution, and for the correlated topic model $\tau$ is a logistic Normal distribution. The generative process for a document $d$ begins by drawing the document’s topic distribution $W_{d}\sim\tau$. Then, for each position $i$ we sample a topic assignment $z_{i}\sim W_{d}$, and finally a word $w_{i}\sim A_{z_{i}}$.

We can combine the column vectors $A_{k}$ for each of the $K$ topics to obtain the word-topic matrix $A$ of dimension $V\times K$. We can similarly combine the column vectors $W_{d}$ for $M$ documents to obtain the topic-document matrix $W$ of dimension $K\times M$. We emphasize that $W$ is unknown and stochastically generated: we can never expect to be able to recover it. The learning task that we consider is to find the word-topic matrix $A$. For the case when $\tau$ is Dirichlet (LDA), we also show how to learn hyperparameters of $\tau$.

Maximum likelihood estimation of the word-topic distributions is NP-hard even for two topics [Arora et al., 2012b], and as a result researchers typically use approximate inference. The most popular approaches are variational inference [Blei et al., 2003], which optimizes an approximate objective, and Markov chain Monte Carlo [McCallum, 2002], which asymptotically samples from the posterior distribution but has no guarantees of convergence.

Arora et al. [2012b] present an algorithm that provably learns the parameters of a topic model given samples from the model, provided that the word-topic distributions satisfy a condition called separability:

The word-topic matrix $A$ is $p$-separable for $p>0$ if for each topic $k$, there is some word $i$ such that $A_{i,k}\geq p$ and $A_{i,k^{\prime}}=0$ for $k^{\prime}\neq k$.

There is a polynomial time algorithm that learns the parameters of a topic model if the number of documents is at least

where $p$ is defined above, $\gamma$ is the condition number of $R$, and $a=\max_{k,k^{\prime}}\alpha_{k}/\alpha_{k^{\prime}}$. The algorithm learns the word-topic matrix $A$ and the topic-topic covariance matrix $R$ up to additive error $\epsilon$.

Unfortunately, this algorithm is not practical. Its running time is prohibitively large because it solves $V$ linear programs, and its use of matrix inversion makes it unstable and sensitive to noise. In this paper, we will give various reformulations and modifications of this algorithm that alleviate these problems altogether.

A Probabilistic Approach to Exploiting Separability

The Arora et al. [2012b] algorithm has two steps: anchor selection, which identifies anchor words, and recovery, which recovers the parameters of $A$ and of $\tau$. Both anchor selection and recovery take as input the matrix $Q$ (of size $V\times V$) of word-word co-occurrence counts, whose construction is described in the supplementary material. $Q$ is normalized so that the sum of all entries is $1$. The high-level flow of our complete learning algorithm is described in Algorithm 1, and follows the same two steps. In this section we will introduce a new recovery method based on a probabilistic framework. We defer the discussion of anchor selection to the next section, where we provide a purely combinatorial algorithm for finding the anchor words.

where $D$ is a diagonal matrix of size $K\times K$. Next, it solves for $A$ and $R$ using the algebraic manipulations outlined in Algorithm 2.

The use of matrix inversion in Algorithm 2 results in substantial imprecision in the estimates when we have small sample sizes. The returned $A$ and $R$ matrices can even contain small negative values, requiring a subsequent projection onto the simplex. As we will show in Section 5, the original recovery method performs poorly relative to a likelihood-based algorithm. Part of the problem is that the original recover algorithm uses only $K$ rows of the matrix $Q$ (the rows for the anchor words), whereas $Q$ is of dimension $V\times V$. Besides ignoring most of the data, this has the additional complication that it relies completely on co-occurrences between a word and the anchors, and this estimate may be inaccurate if both words occur infrequently.

Here we adopt a new probabilistic approach, which we describe below after introducing some notation. Consider any two words in a document and call them $w_{1}$ and $w_{2}$, and let $z_{1}$ and $z_{2}$ refer to their topic assignments. We will use $A_{i,k}$ to index the matrix of word-topic distributions, i.e. $A_{i,k}=p(w_{1}=i|z_{1}=k)=p(w_{2}=i|z_{2}=k)$. Given infinite data, the elements of the $Q$ matrix can be interpreted as $Q_{i,j}=p(w_{1}=i,w_{2}=j)$. The row-normalized $Q$ matrix, denoted $\bar{Q}$, which plays a role in both finding the anchor words and the recovery step, can be interpreted as a conditional probability $\bar{Q}_{i,j}=p(w_{2}=j|w_{1}=i)$.

Denoting the indices of the anchor words as $\mathbf{S}=\{s_{1},s_{2},...,s_{K}\}$, the rows indexed by elements of $\mathbf{S}$ are special in that every other row of $\bar{Q}$ lies in the convex hull of the rows indexed by the anchor words. To see this, first note that for an anchor word $s_{k}$,

where (1) uses the fact that in an admixture model $w_{2}\bot w_{1}\mid z_{1}$, and (2) is because $p(z_{1}=k|w_{1}=s_{k})=1$. For any other word $i$, we have

Denoting the probability $p(z_{1}=k|w_{1}=i)$ as $C_{i,k}$, we have $\bar{Q}_{i,j}=\sum_{k}C_{i,k}\bar{Q}_{s_{k},j}$. Since $C$ is non-negative and $\sum_{k}C_{i,k}=1$, we have that any row of $\bar{Q}$ lies in the convex hull of the rows corresponding to the anchor words. The mixing weights give us $p(z_{1}|w_{1}=i)$! Using this together with $p(w_{1}=i)$, we can recover the $A$ matrix simply by using Bayes’ rule:

Finally, we observe that $p(w_{1}=i)$ is easy to solve for since $\sum_{j}Q_{i,j}=\sum_{j}p(w_{1}=i,w_{2}=j)=p(w_{1}=i)$.

Our new algorithm finds, for each row of the empirical row normalized co-occurrence matrix, $\hat{Q}_{i}$, the coefficients $p(z_{1}|w_{1}=i)$ that best reconstruct it as a convex combination of the rows that correspond to anchor words. This step can be solved quickly and in parallel (independently) for each word using the exponentiated gradient algorithm. Once we have $p(z_{1}|w_{1})$, we recover the $A$ matrix using Bayes’ rule. The full algorithm using KL divergence as an objective is found in Algorithm 3. Further details of the exponentiated gradient algorithm are given in the supplementary material.

One reason to use KL divergence as the measure of reconstruction error is that the recovery procedure can then be understood as maximum likelihood estimation. In particular, we seek the parameters $p(w_{1})$, $p(z_{1}|w_{1})$, $p(w_{2}|z_{1})$ that maximize the likelihood of observing the word co-occurence counts, $\hat{Q}$. However, the optimization problem does not explicitly constrain the parameters to correspond an admixture model.

We can also define a similar algorithm using quadratic loss, which we call RecoverL2. This formulation has the extremely useful property that both the objective and gradient can be kernelized so that the optimization problem is independent of the vocabulary size. To see this, notice that the objective can be re-written as

where $\overline{Q}_{\mathbf{S}}\overline{Q}^{T}_{\mathbf{S}}$ is $K\times K$ and can be computed once and used for all words, and $\overline{Q}_{\mathbf{S}}\overline{Q}_{i}^{T}$ is $K\times 1$ and can be computed once prior to running the exponentiated gradient algorithm for word $i$.

To recover the $R$ matrix for an admixture model, recall that $Q=ARA^{T}$. This may be an over-constrained system of equations with no solution for $R$, but we can find a least-squares approximation to $R$ by pre- and post-multiplying $Q$ by the pseudo-inverse $A^{\dagger}$. For the special case of LDA we can learn the Dirichlet hyperparameters. Recall that in applying Bayes’ rule we calculated $p(z_{1})=\sum_{i^{\prime}}p(z_{1}|w_{1}=i^{\prime})p(w_{1}=i^{\prime})$. These values for $p(z)$ specify the Dirichlet hyperparameters up to a constant scaling. This constant could be recovered from the $R$ matrix [Arora et al., 2012b], but in practice we find it is better to choose it using a grid search to maximize the likelihood of the training data.

We will see in Section 5 that our nonnegative recovery algorithm performs much better on a wide range of performance metrics than the recovery algorithm in Arora et al. [2012b]. In the supplementary material we show that it also inherits the theoretical guarantees of Arora et al. [2012b]: given polynomially many documents, our algorithm returns an estimate $\hat{A}$ at most $\epsilon$ from the true word-topic matrix $A$.

A Combinatorial Algorithm for Finding Anchor Words

Here we consider the anchor selection step of the algorithm where our goal is to find the anchor words. In the infinite data case where we have infinitely many documents, the convex hull of the rows in $\overline{Q}$ will be a simplex where the vertices of this simplex correspond to the anchor words. Since we only have a finite number of documents, the rows of $\overline{Q}$ are only an approximation to their expectation. We are therefore given a set of $V$ points $d_{1},d_{2},...d_{V}$ that are each a perturbation of $a_{1},a_{2},...a_{V}$ whose convex hull $P$ defines a simplex. We would like to find an approximation to the vertices of $P$. See Arora et al. [2012a] and Arora et al. [2012b] for more details about this problem.

Arora et al. [2012a] give a polynomial time algorithm that finds the anchor words. However, their algorithm is based on solving $V$ linear programs, one for each word, to test whether or not a point is a vertex of the convex hull. In this section we describe a purely combinatorial algorithm for this task that avoids linear programming altogether. The new “FastAnchorWords” algorithm is given in Algorithm 4. To find all of the anchor words, our algorithm iteratively finds the furthest point from the subspace spanned by the anchor words found so far.

Since the points we are given are perturbations of the true points, we cannot hope to find the anchor words exactly. Nevertheless, the intuition is that even if one has only found $r$ points $S$ that are close to $r$ (distinct) anchor words, the point that is furthest from $\operatorname{span}(S)$ will itself be close to a (new) anchor word. The additional advantage of this procedure is that when faced with many choices for a next anchor word to find, our algorithm tends to find the one that is most different than the ones we have found so far.

The main contribution of this section is a proof that the FastAnchorWords algorithm succeeds in finding $K$ points that are close to anchor words. To precisely state the guarantees, we recall the following definition from [Arora et al., 2012a]:

In most reasonable settings the parameters of the topic model define lower bounds on the robustness of the polytope $P$. For example, in LDA, this lower bound is based on the largest ratio of any pair of hyper-parameters in the model [Arora et al., 2012b]. Our goal is to find a set of points that are close to the vertices of the simplex, and to make this precise we introduce the following definition:

Let $a_{1},a_{2},...a_{V}$ be a set of points whose convex hull $P$ is a simplex with vertices $v_{1},v_{2},...v_{K}$. Then we say $a_{i}$ $\epsilon$-covers $v_{j}$ if when $a_{j}$ is written as a convex combination of the vertices as $a_{i}=\sum_{j}c_{j}v_{j}$, then $c_{j}\geq 1-\epsilon$. Furthermore we will say that a set of $K$ points $\epsilon$-covers the vertices if each vertex is $\epsilon$ covered by some point in the set.

We will prove the following theorem: suppose there is a set of points $\mathcal{A}=a_{1},a_{2},...a_{V}$ whose convex hull $P$ is $\gamma$-robust and has vertices $v_{1},v_{2},...v_{K}$ (which appear in $\mathcal{A}$) and that we are given a perturbation $d_{1},d_{2},...d_{V}$ of the points so that for each $i$, $\|a_{i}-d_{i}\|\leq\epsilon$, then:

This new algorithm not only helps us avoid linear programming altogether in inferring the parameters of a topic model, but also can be used to solve the nonnegative matrix factorization problem under the separability assumption, again without resorting to linear programming. Our analysis rests on the following lemmas, whose proof we defer to the supplementary material. Suppose the algorithm has found a set $S$ of $k$ points that are each $\delta$-close to distinct vertices in $\{v_{1},v_{2},...,v_{K}\}$ and that $\delta<\gamma/20K$.

There is a vertex $v_{i}$ whose distance from $\operatorname{span}(S)$ is at least $\gamma/2$.

The proof of this lemma is based on a volume argument, and the connection between the volume of a simplex and the determinant of the matrix of distances between its vertices.

The point $d_{j}$ found by the algorithm must be $\delta=O(\epsilon/\gamma^{2})$ close to some vertex $v_{i}$.

This lemma is used to show that the error does not accumulate too badly in our algorithm, since $\delta$ only depends on $\epsilon$, $\gamma$ (not on the $\delta$ used in the previous step of the algorithm). This prevents the error from accumulating exponentially in the dimension of the problem, which would be catastrophic for our proof.

After running the first phase of our algorithm, we run a cleanup phase (the second loop in Alg. 4) that can reduce the error in our algorithm. When we have $K-1$ points close to $K-1$ vertices, only one of the vertices can be far from their span. The farthest point must be close to this missing vertex. The following lemma shows that this cleanup phase can improve the guarantees of Lemma A.2:

Suppose $|S|=K-1$ and each point in $S$ is $\delta=O(\epsilon/\gamma^{2})<\gamma/20K$ close to some vertex $v_{i}$, then the farthest point $v_{j}^{\prime}$ found by the algorithm is $1-O(\epsilon/\gamma)$ close to the remaining vertex.

This algorithm is a greedy approach to maximizing the volume of the simplex. The larger the volume is, the more words per document the resulting model can explain. Better anchor word selection is an open question for future work. We have experimented with a variety of other heuristics for maximizing simplex volume, with varying degrees of success.

Related work. The separability assumption has also been studied under the name “pure pixel assumption” in the context of hyperspectral unmixing. A number of algorithms have been proposed that overlap with ours – such as the VCA [Nascimento & Dias, 2004] algorithm (which differs in that there is no clean-up phase) and the N-FINDR [Gomez et al., 2007] algorithm which attempts to greedily maximize the volume of a simplex whose vertices are data points. However these algorithms have only been proven to work in the infinite data case, and for our algorithm we are able to give provable guarantees even when the data points are perturbed (e.g., as the result of sampling noise). Recent work of Thurau et al. and Kumar et al. follow the same pattern as our paper, but use non-negative matrix factorization under the separability assumption. While both give applications to topic modeling, in realistic applications the term-by-document matrix is too sparse to be considered a good approximation to its expectation (because documents are short). In contrast, our algorithm works with the Gram matrix $Q$ so that we can give provable guarantees even when each document is short.

Experimental Results

We compare three parameter recovery methods, Recover [Arora et al., 2012b], RecoverKL and RecoverL2 to a fast implementation of Gibbs sampling [McCallum, 2002].We were not able to obtain Anandkumar et al. ’s implementation of their algorithm, and our own implementation is too slow to be practical. Linear programming-based anchor word finding is too slow to be comparable, so we use FastAnchorWords for all three recovery algorithms. Using Gibbs sampling we obtain the word-topic distributions by averaging over 10 saved states, each separated by 100 iterations, after 1000 burn-in iterations.

We train models on two synthetic data sets to evaluate performance when model assumptions are correct, and real documents to evaluate real-world performance. To ensure that synthetic documents resemble the dimensionality and sparsity characteristics of real data, we generate semi-synthetic corpora. For each real corpus, we train a model using MCMC and then generate new documents using the parameters of that model (these parameters are not guaranteed to be separable).

We use two real-world data sets, a large corpus of New York Times articles (295k documents, vocabulary size 15k, mean document length 298) and a small corpus of NIPS abstracts (1100 documents, vocabulary size 2500, mean length 68). Vocabularies were pruned with document frequency cutoffs. We generate semi-synthetic corpora of various sizes from models trained with $K=100$ from NY Times and NIPS, with document lengths set to 300 and 70, respectively, and with document-topic distributions drawn from a Dirichlet with symmetric hyperparameters $0.03$.

The parameter $\epsilon=0.01$ is used to avoid taking the $\log$ of zero for words that never co-occur [Stevens et al., 2012]. This metric has been shown to correlate well with human judgments of topic quality. If we perfectly reconstruct topics, all the high-probability words in a topic should co-occur frequently, otherwise, the model may be mixing unrelated concepts. Coherence measures the quality of individual topics, but does not measure redundancy, so we measure inter-topic similarity. For each topic, we gather the set of the $N$ most probable words. We then count how many of those words do not appear in any other topic’s set of $N$ most probable words. Some overlap is expected due to semantic ambiguity, but lower numbers of unique words indicate less useful models.

2 Efficiency

The Recover algorithms, in Python, are faster than a heavily optimized Java Gibbs sampling implementation [Yao et al., 2009].

Fig. 1 shows the time to train models on synthetic corpora on a single machine. Gibbs sampling is linear in the corpus size. RecoverL2 is also linear ($\rho=0.79$), but only varies from 33 to 50 seconds. Estimating $Q$ is linear, but takes only 7 seconds for the largest corpus. FastAnchorWords takes less than 6 seconds for all corpora.

3 Semi-synthetic documents

4 Effect of separability

5 Effect of correlation

The theoretical guarantees of the new algorithms apply even if topics are correlated. To test how algorithms respond to correlation, we generated new synthetic corpora from the same $K=100$ model trained on NY Times articles. Instead of a symmetric Dirichlet distribution, we use a logistic normal distribution with a block-structured covariance matrix. We partition topics into 10 groups. For each pair of topics in a group, we add a non-zero off-diagonal element to the covariance matrix. This block structure is not necessarily realistic, but shows the effect of correlation. Results for two levels of covariance ($\rho=0.05,\rho=0.1$) are shown in Fig. 5.

6 Real documents

The new algorithms produce comparable quantitative and qualitative results on real data. Fig. 6 shows three metrics for both corpora. Error bars show the distribution of log probabilities across held-out documents (top panel) and coherence and unique words across topics (center and bottom panels). Held-out sets are 230 documents for NIPS and 59k for NY Times. For the small NIPS corpus we average over 5 non-overlapping train/test splits. The matrix-inversion in Recover failed for the smaller corpus (NIPS). In the larger corpus (NY Times), Recover produces noticeably worse held-out log probability per token than the other algorithms. Gibbs sampling produces the best average held-out probability ($p<0.0001$ under a paired $t$-test), but the difference is within the range of variability between documents. We tried several methods for estimating hyperparameters, but the observed differences did not change the relative performance of algorithms.

Gibbs sampling has worse coherence than the Recover algorithms, but produces more unique words per topic. These patterns are consistent with semi-synthetic results for similarly sized corpora (details are in supplementary material).

Conclusions

References

Appendix A Proof for Anchor-Words Finding Algorithm

Recall that the correctness of the algorithm depends on the following Lemmas:

There is a vertex $v_{i}$ whose distance from $\operatorname{span}(S)$ is at least $\gamma/2$.

The point $\Delta_{j}$ found by the algorithm must be $\delta=O(\epsilon/\gamma^{2})$ close to some vertex $v_{i}$.

In order to prove Lemma A.1, we use a volume argument. First we show that the volume of a robust simplex cannot change by too much when the vertices are perturbed.

Suppose $\{v_{1},v_{2},...,v_{K}\}$ are the vertices of a $\gamma$-robust simplex $S$. Let $S^{\prime}$ be a simplex with vertices $\{v_{1}^{\prime},v_{2}^{\prime},...,v_{K}^{\prime}\}$, each of the vertices $v_{i}^{\prime}$ is a perturbation of $v_{i}$ and $\left\lVert v_{i}^{\prime}-v_{i}\right\rVert_{2}\leq\delta$. When $10\sqrt{K}\delta<\gamma$ the volume of the two simplices satisfy

Proof: As the volume of a simplex is proportional to the determinant of a matrix whose columns are the edges of the simplex, we first show the following perturbation bound for determinant.

Let $A$, $E$ be $K\times K$ matrices, the smallest eigenvalue of $A$ is at least $\gamma$, the Frobenius norm $\left\lVert E\right\rVert_{F}\leq\sqrt{K}\delta$, when $\gamma>5\sqrt{K}\delta$ we have

Proof: Since $\det(AB)=\det(A)\det(B)$, we can multiply both $A$ and $A+E$ by $A^{-1}$. Hence $\det(A+E)/\det(A)=\det(I+A^{-1}E)$.

The Frobenius norm of $A^{-1}E$ is bounded by

Let the eigenvalues of $A^{-1}E$ be $\lambda_{1},\lambda_{2},...,\lambda_{K}$, then by definition of Frobenius Norm $\sum_{i=1}^{K}\lambda_{i}^{2}\leq\left\lVert A^{-1}E\right\rVert_{F}^{2}\leq K\delta^{2}/\gamma^{2}$. The eigenvalues of $I+A^{-1}E$ are just $1+\lambda_{1},1+\lambda_{2},...,1+\lambda_{K}$, and the determinant $\det(I+A^{-1}E)=\prod_{i=1}^{K}(1+\lambda_{i})$. Hence it suffices to show

To do this we apply Lagrangian method and show the minimum is only obtained when all $\lambda_{i}$’s are equal. The optimal value must be obtained at a local optimum of

Taking partial derivatives with respect to $\lambda_{i}$’s, we get the equations $-\lambda_{i}(1+\lambda_{i})=-\prod_{i=1}^{K}(1+\lambda_{i})/2C$ (here using $\sqrt{K}\delta/\gamma$ is small so $1+\lambda_{i}>1/2>0$). The right hand side is a constant, so each $\lambda_{i}$ must be one of the two solutions of this equation. However, only one of the solution is larger than $1/2$, therefore all the $\lambda_{i}$’s are equal. $\blacksquare$

For the lower bound, we can project the perturbed subspace to the $K-1$ dimensional space. Such a projection cannot increase the volume and the perturbation distances only get smaller. Therefore we can apply the claim directly, the columns of $A$ are just $v_{i+1}-v_{1}$ for $i=1,2,...,K-1$; columns of $E$ are just $v_{i+1}^{\prime}-v_{i+1}-(v_{1}^{\prime}-v_{1})$. The smallest eigenvalue of $A$ is at least $\gamma$ because the polytope is $\gamma$ robust, which is equivalent to saying after orthogonalization each column still has length at least $\gamma$. The Frobenius norm of $E$ is at most $2\sqrt{K-1}\delta$. We get the lower bound directly by applying the claim.

For the upper bound, swap the two sets $S$ and $S^{\prime}$ and use the argument for the lower bound. The only thing we need to show is that the smallest eigenvalue of the matrix generated by points in $S^{\prime}$ is still at least $\gamma/2$. This follows from Wedin’s Theorem [Wedin, 1972] and the fact that $\left\lVert E\right\rVert\leq\left\lVert E\right\rVert_{F}\leq\sqrt{K}\delta\leq\gamma/2$. $\blacksquare$

Proof: The first case is for the first step of the algorithm, when we try to find the farthest point to the origin. Here essentially $S=\{\vec{0}\}$. For any two vertices $v_{1},v_{2}$, since the simplex is $\gamma$ robust, the distance between $v_{1}$ and $v_{2}$ is at least $\gamma$. Which means $\mbox{dis}(\vec{0},v_{1})+\mbox{dis}(\vec{0},v_{2})\geq\gamma$, one of them must be at least $\gamma/2$.

For the later steps, recall that $S$ contains vertices of a perturbed simplex. Let $S^{\prime}$ be the set of original vertices corresponding to the perturbed vertices in $S$. Let $v$ be any vertex in $\{v_{1},v_{2},...,v_{K}\}$ which is not in $S$. Now we know the distance between $v$ and $S$ is equal to $\mbox{vol}(S\cup\{v\})/(|S|-1)\mbox{vol}(S)$. On the other hand, we know $\mbox{vol}(S^{\prime}\cup\{v\})/(|S^{\prime}|-1)\mbox{vol}(S^{\prime})\geq\gamma$. Using Lemma A.3 to bound the ratio between the two pairs $\mbox{vol}(S)/\mbox{vol}(S^{\prime})$ and $\mbox{vol}(S\cup\{v\})/\mbox{vol}(S^{\prime}\cup\{v\})$, we get:

when $\gamma>20K\epsilon^{\prime}$. $\blacksquare$

Since $d_{j}$ is the point found by the algorithm, let us consider the point $a_{j}$ before perturbation. The point $a_{j}$ is inside the simplex, therefore we can write $a_{j}$ as a convex combination of the vertices:

Let $v_{t}$ be the vertex with largest coefficient $c_{t}$. Let $\Delta$ be the largest distance from some vertex to the space spanned by points in $S$ ($\Delta=\max_{l}\mbox{dis}(v_{l},\operatorname{span}(S))$. By Lemma A.1 we know $\Delta>\gamma/2$. Also notice that we are not assuming $\mbox{dis}(v_{t},\operatorname{span}(S))=\Delta$.

Now we rewrite $a_{j}$ as $c_{t}v_{t}+(1-c_{t})w$, where $w$ is a vector in the convex hull of vertices other than $v_{t}$. Observe that $a_{j}$ must be far from $\operatorname{span}(S)$, because $d_{j}$ is the farthest point found by the algorithm. Indeed:

The second inequality is because there must be some point $d_{l}$ that correspond to the farthest vertex $v_{l}$ and have $\mbox{dis}(d_{l},\operatorname{span}(S))\geq\Delta-\epsilon$. Thus as $d_{j}$ is the farthest point $\mbox{dis}(d_{j},\operatorname{span}(S))\geq\mbox{dis}(d_{l},\operatorname{span}(S))\geq\Delta-\epsilon$.

The point $v_{t}$ must be far from $w$ by applying Lemma A.1: consider the set of vertices V^{\prime}=\{v_{i}:v_{i}\mbox{ does not correspond to any point inSand }i\neq t\}. The set $V^{\prime}\cup S$ satisfy the assumptions in Lemma A.1 so there must be one vertex that is far from $\operatorname{span}(V^{\prime}\cup S)$, and it can only be $v_{t}$. Therefore even after projecting to orthogonal subspace of $\operatorname{span}(S)$, $v_{t}$ is still far from any convex combination of $V^{\prime}$. The vertices that are not in $V^{\prime}$ all have very small norm after projecting to orthogonal subspace (at most $\delta_{0}$) so we know the distance of $v_{t}$ and $w$ is at least $\gamma/2-\delta_{0}>\gamma/4$.

Now the problem becomes a two dimensional calculation. When $c_{t}$ is fixed the length of $a_{j}$ is strictly increasing when the distance of $v_{t}$ and $w$ decrease, so we assume the distance is $\gamma/4$. Simple calculation (using essentially just pythagorean theorem) shows

The right hand side is largest when $\Delta=2$ (since the vectors are in unit ball) and the maximum value is $O(\epsilon/\gamma^{2})$. When this value is smaller than $1/K$, we must have $1-c_{t}\leq O(\epsilon/\gamma^{2})$. Thus $c_{t}\geq 1-O(\epsilon/\gamma^{2})$ and $\delta\leq(1-c_{t})+\epsilon\leq O(\epsilon/\gamma^{2})$. $\blacksquare$

The cleanup phase tries to find the farthest point to a subset of $K-1$ vertices, and use that point as the $K$-th vertex. This will improve the result because when we have $K-1$ points close to $K-1$ vertices, only one of the vertices can be far from their span. Therefore the farthest point must be close to the only remaining vertex. Another way of viewing this is that the algorithm is trying to greedily maximize the volume of the simplex, which makes sense because the larger the volume is, the more words/documents the final LDA model can explain.

The following lemma makes the intuitions rigorous and shows how cleanup improves the guarantee of Lemma A.2.

Suppose $|S|=K-1$ and each point in $S$ is $\delta=O(\epsilon/\gamma^{2})<\gamma/20K$ close to some vertex $v_{i}$, then the farthest point $v_{j}^{\prime}$ found by the algorithm is $1-O(\epsilon/\gamma)$ close to the remaining vertex.

Proof: We still look at the original point $a_{j}$ and express it as $\sum_{t=1}^{K}c_{t}v_{t}$. Without loss of generality let $v_{1}$ be the vertex that does not correspond to anything in $S$. By Lemma A.1 $v_{1}$ is $\gamma/2$ far from $\operatorname{span}(S)$. On the other hand all other vertices are at least $\gamma/20r$ close to $\operatorname{span}(S)$. We know the distance $\mbox{dis}(a_{j},\operatorname{span}(S))\geq\mbox{dis}(v_{1},\operatorname{span}(S))-2\epsilon$, this cannot be true unless $c_{t}\geq 1-O(\epsilon/\gamma)$. $\blacksquare$

These lemmas directly lead to the following theorem:

Proof: In the first phase of the algorithm, do induction using Lemma A.2. When $20K\epsilon/\gamma^{2}<\gamma$ Lemma A.2 shows that we find a set of points that $O(\epsilon/\gamma^{2})$-covers the vertices. Now Lemma A.5 shows after cleanup phase the points are refined to $O(\epsilon/\gamma)$-cover the vertices. $\blacksquare$

Appendix B Proof for Nonnegative Recover Procedure

In order to show RecoverL2 learns the parameters even when the rows of $\bar{Q}$ are perturbed, we need the following lemma that shows when columns of $\bar{Q}$ are close to the expectation, the posteriors $c$ computed by the algorithm is also close to the true value.

For a $\gamma$ robust simplex $S$ with vertices $\{v_{1},v_{2},...,v_{K}\}$, let $v$ be a point in the simplex that can be represented as a convex combination $v=\sum_{i=1}^{K}c_{i}v_{i}$. If the vertices of $S$ are perturbed to $S^{\prime}=\{...,v_{i}^{\prime},...\}$ where $\left\lVert v_{i}^{\prime}-v_{i}\right\rVert\leq\delta_{1}$ and $v$ is perturbed to $v^{\prime}$ where $\left\lVert v-v^{\prime}\right\rVert\leq\delta_{2}$. Let $v^{*}$ be the point in $S^{\prime}$ that is closest to $v^{\prime}$, and $v^{*}=\sum_{i=1}^{K}c_{i}^{\prime}v_{i}$, when $10\sqrt{K}\delta_{1}\leq\gamma$ for all $i\in[K]$ $|c_{i}-c_{i}^{\prime}|\leq 4(\delta_{1}+\delta_{2})/\gamma$.

Proof: Consider the point $u=\sum_{i=1}^{K}c_{i}v_{i}^{\prime}$, by triangle inequality: $\left\lVert u-v\right\rVert\leq\sum_{i=1}^{K}c_{i}\left\lVert v_{i}-v_{i}^{\prime}\right\rVert\leq\delta_{1}$. Hence $\left\lVert u-v^{\prime}\right\rVert\leq\left\lVert u-v\right\rVert+\left\lVert v-v^{\prime}\right\rVert\leq\delta_{1}+\delta_{2}$, and $u$ is in $S^{\prime}$. The point $v^{*}$ is the point in $S^{\prime}$ that is closest to $v^{\prime}$, so $\left\lVert v^{*}-v^{\prime}\right\rVert\leq\delta_{1}+\delta_{2}$ and $\left\lVert v^{*}-u\right\rVert\leq 2(\delta_{1}+\delta_{2})$.

Then we need to show when a point ($u$) moves a small distance, its representation also changes by a small amount. Intuitively this is true because $S$ is $\gamma$ robust. By Lemma A.1 when $10\sqrt{K}\delta_{1}<\gamma$, the simplex $S^{\prime}$ is also $\gamma/2$ robust. For any $i$, let $Proj_{i}(v^{*})$ and $Proj_{i}(u)$ be the projections of $v^{*}$ and $u$ in the orthogonal subspace of $\operatorname{span}(S^{\prime}\backslash v_{i}^{\prime})$, then

and this completes the proof. $\blacksquare$

With this lemma it is not hard to show that RecoverL2 has polynomial sample complexity.

When the number of documents $M$ is at least

our algorithm using the conjunction of FastAnchorWords and RecoverL2 learns the $A$ matrix with entry-wise error at most $\epsilon$.

The error in each column of $\bar{Q}$ can be at most $\delta_{2}=\epsilon_{Q}\sqrt{4aK/\epsilon}$. By Theorem A.6 when $20K\delta_{2}/(\gamma p)^{2}<\gamma p$ (which is satisfied when $M=O(aK^{3}\log V/D(\gamma p)^{6}\epsilon)$) , the anchor words found are $\delta_{1}=O(\delta_{2}/(\gamma p))$ close to the true anchor words. Hence by Lemma B.1 every entry of $C$ has error at most $O(\delta_{2}/(\gamma p)^{2})$.

With such number of documents, all the word probabilities $p(w=i)$ are estimated more accurately than the entries of $C_{i,j}$, so we omit their perturbations here for simplicity. When we apply the Bayes rule, we know $A_{i,k}=C_{i,k}p(w=i)/p(z=k)$, where $p(z=k)$ is $\alpha_{k}$ which is lower bounded by $1/aK$. The numerator and denominator are all related to entries of $C$ with positive coefficients sum up to at most 1. Therefore the errors $\delta_{num}$ and $\delta_{denom}$ are at most the error of a single entry of $C$, which is bounded by $O(\delta_{2}/(\gamma p)^{2})$. Applying Taylor’s Expansion to $(p(z=k,w=i)+\delta_{num})/(\alpha_{k}+\delta_{denom})$, the error on entries of $A$ is at most $O(aK\delta_{2}/(\gamma p)^{2})$. When $\epsilon_{Q}\leq O((\gamma p)^{2}\epsilon^{1.5}/(aK)^{1.5})$, we have $O(aK\delta_{2}/(\gamma p)^{2})\leq\epsilon$, and get the desired accuracy of $A$. The number of document required is $M=O((aK)^{3}\log V/D\epsilon^{3}(\gamma p)^{4})$.

The sample complexity for $R$ can then be bounded using matrix perturbation theory. $\blacksquare$

Appendix C Empirical Results

Appendix D Algorithmic Details

D.2 Exponentiated gradient algorithm

The optimization problem that arises in RecoverKL and RecoverL2 has the following form,

where $d(\cdot,\cdot)$ is a Bregman divergence, $x$ is a vector of length $K$, and $T$ is a matrix of size $V\times K$. We solve this optimization problem using the Exponentiated Gradient algorithm [Kivinen & Warmuth, 1995], described in Algorithm 5. In our experiments we show results using both squared Euclidean distance and KL divergence for the divergence measure. Stepsizes are chosen with a line search to find an $\eta$ that satisfies the Wolfe and Armijo conditions (For details, see Nocedal & Wright ). We test for convergence using the KKT conditions. Writing the KKT conditions for our constrained minimization problem:

Stationarity: $\nabla_{x}d(b,Tx^{*})-\vec{\lambda}+\mu{\bf 1}$ = 0

Primal Feasibility: $x^{*}\geq 0$, $|x^{*}|_{1}=1$

Complementary Slackness: $\lambda_{i}x_{i}^{*}=0$

For every iterate of $x$ generated by Exponentiated Gradient, we set $\lambda,\mu$ to satisfy conditions 1-3. This gives the following equations:

By construction conditions 1-3 are satisfied (note that the multiplicative update and the projection step ensure that $x$ is always primal feasible). Convergence is tested by checking whether the final KKT condition holds within some tolerance. Since $\lambda$ and $x$ are nonnegative, we check complimentary slackness by testing whether $\lambda^{T}x<\epsilon$. This convergence test can also be thought of as testing the value of the primal-dual gap, since the Lagrangian function has the form: $L(x,\lambda,\mu)=d(b,Tx)-\lambda^{T}x+\mu(x^{T}{\bf 1}-1)$, and $(x^{T}{\bf 1}-1)$ is zero at every iteration.

The running time of RecoverL2 is the time of solving $V$ small ($K\times K$) quadratic programs. Especially when using Exponentiated Gradient to solve the quadratic program, each word requires $O(KV)$ time for preprocessing and $O(K^{2})$ per iteration. The total running time is $O(KV^{2}+K^{2}VT)$ where $T$ is the average number of iterations. The value of $T$ is about $100-1000$ depending on data sets.