Factoring nonnegative matrices with linear programs

Introduction

Nonnegative matrix factorization (NMF) is a popular approach for selecting features in data . Many machine-learning and data-mining software packages (including Matlab , R , and Oracle Data Mining ) now include heuristic computational methods for NMF. Nevertheless, we still have limited theoretical understanding of when these heuristics are correct.

The difficulty in developing rigorous methods for NMF stems from the fact that the problem is computationally challenging. Indeed, Vavasis has shown that NMF is NP-Hard ; see for further worst-case hardness results. As a consequence, we must instate additional assumptions on the data if we hope to compute nonnegative matrix factorizations in practice.

In this spirit, Arora, Ge, Kannan, and Moitra (AGKM) have exhibited a polynomial-time algorithm for NMF that is provably correct—provided that the data is drawn from an appropriate model, based on ideas from . The AGKM result describes one circumstance where we can be sure that NMF algorithms are capable of producing meaningful answers. This work has the potential to make an impact in machine learning because proper feature selection is an important preprocessing step for many other techniques. Even so, the actual impact is damped by the fact that the AGKM algorithm is too computationally expensive for large-scale problems and is not tolerant to departures from the modeling assumptions. Thus, for NMF, there remains a gap between the theoretical exercise and the actual practice of machine learning.

The present work presents a scalable, robust algorithm that can successfully solve the NMF problem under appropriate hypotheses. Our first contribution is a new formulation of the nonnegative feature selection problem that only requires the solution of a single linear program. Second, we provide a theoretical analysis of this algorithm. This argument shows that our method succeeds under the same modeling assumptions as the AGKM algorithm with an additional margin constraint that is common in machine learning. We prove that if there exists a unique, well-defined model, then we can recover this model accurately; our error bound improves substantially on the error bound for the AGKM algorithm in the high SNR regime. One may argue that NMF only “makes sense” (i.e., is well posed) when a unique solution exists, and so we believe our result has independent interest. Furthermore, our algorithm can be adapted for a wide class of noise models.

In addition to these theoretical contributions, our work also includes a major algorithmic and experimental component. Our formulation of NMF allows us to exploit methods from operations research and database systems to design solvers that scale to extremely large datasets. We develop an efficient stochastic gradient descent (SGD) algorithm that is (at least) two orders of magnitude faster than the approach of AGKM when both are implemented in Matlab. We describe a parallel implementation of our SGD algorithm that can robustly factor matrices with $10^{5}$ features and $10^{6}$ examples in a few minutes on a multicore workstation.

Our formulation of NMF uses a data-driven modeling approach to simplify the factorization problem. More precisely, we search for a small collection of rows from the data matrix that can be used to express the other rows. This type of approach appears in a number of other factorization problems, including rank-revealing QR , interpolative decomposition , subspace clustering , dictionary learning , and others. Our computational techniques can be adapted to address large-scale instances of these problems as well.

Separable Nonnegative Matrix Factorizations and Hott Topics

Notation. For a matrix $\bm{M}$ and indices $i$ and $j$, we write $\bm{M}_{i\cdot}$ for the $i$th row of $\bm{M}$ and $\bm{M}_{\cdot j}$ for the $j$th column of $\bm{M}$. We write $M_{ij}$ for the $(i,j)$ entry.

Let $\bm{Y}$ be a nonnegative $f\times n$ data matrix with columns indexing examples and rows indexing features. Exact NMF seeks a factorization $\bm{Y}=\bm{F}\bm{W}$ where the feature matrix $\bm{F}$ is $f\times r$, where the weight matrix $\bm{W}$ is $r\times n$, and both factors are nonnegative. Typically, $r\ll\min\{f,n\}$.

Unless stated otherwise, we assume that each row of the data matrix $\bm{Y}$ is normalized so it sums to one. Under this hypothesis, we may also assume that each row of $\bm{F}$ and of $\bm{W}$ also sums to one .

It is notoriously difficult to solve the NMF problem. Vavasis showed that it is NP-complete to decide whether a matrix admits a rank-$r$ nonnegative factorization . AGKM proved that an exact NMF algorithm can be used to solve 3-SAT in subexponential time .

The literature contains some mathematical analysis of NMF that can be used to motivate algorithmic development. Thomas developed a necessary and sufficient condition for the existence of a rank-$r$ NMF. More recently, Donoho and Stodden obtained a related sufficient condition for uniqueness. AGKM exhibited an algorithm that can produce a nonnegative matrix factorization under a weaker sufficient condition. To state their results, we need a definition.

These ideas support the uniqueness results of Donoho and Stodden and the AGKM algorithm. Indeed, we can find an NMF of $\bm{Y}$ efficiently if $\bm{Y}$ contains a set of $r$ rows that is simplicial and whose convex hull contains the remaining rows.

An NMF $\bm{Y}=\bm{F}\bm{W}$ is called separable if the rows of $\bm{W}$ are simplicial and there is a permutation matrix $\bm{\Pi}$ such that

To compute a separable factorization of $\bm{Y}$, we must first identify a simplicial set of rows from $\bm{Y}$. Afterward, we compute weights that express the remaining rows as convex combinations of this distinguished set. We call the simplicial rows hott and the corresponding features hott topics.

This model allows us to express all the features for a particular instance if we know the values of the instance at the simplicial rows. This assumption can be justified in a variety of applications. For example, in text, knowledge of a few keywords may be sufficient to reconstruct counts of the other words in a document. In vision, localized features can be used to predict gestures. In audio data, a few bins of the spectrogram may allow us to reconstruct the remaining bins.

While a nonnegative matrix one encounters in practice might not admit a separable factorization, it may be well-approximated by a nonnnegative matrix with separable factorization. AGKM derived an algorithm for nonnegative matrix factorization of a matrix that is well-approximated by a separable factorization. To state their result, we introduce a norm on $f\times n$ matrices:

Main Theoretical Results: NMF by Linear Programming

This paper shows that we can factor an approximately separable nonnegative matrix by solving a linear program. A major advantage of this formulation is that it scales to very large data sets.

Here is the key observation: Suppose that $\bm{Y}$ is any $f\times n$ nonnegative matrix that admits a rank-$r$ separable factorization $\bm{Y}=\bm{FW}$. If we pad $\bm{F}$ with zeros to form an $f\times f$ matrix, we have

We call the matrix $\bm{C}$ factorization localizing. Note that any factorization localizing matrix $\bm{C}$ is an element of the polyhedral set

Thus, to find an exact NMF of $\bm{Y}$, it suffices to find a feasible element of $\bm{C}\in\Phi(\bm{Y})$ whose diagonal is integral. This task can be accomplished by linear programming. Once we have such a $\bm{C}$, we construct $\bm{W}$ by extracting the rows of $\bm{X}$ that correspond to the indices $i$ where $C_{ii}=1$. We construct the feature matrix $\bm{F}$ by extracting the nonzero columns of $\bm{C}$. This approach is summarized in Algorithm 2. In turn, we can prove the following result.

Suppose $\bm{Y}$ is a nonnegative matrix with a rank-$r$ separable factorization $\bm{Y}=\bm{F}\bm{W}$. Then Algorithm 2 constructs a rank-$r$ nonnegative matrix factorization of $\bm{Y}$.

As the theorem suggests, we can isolate the rows of $\bm{Y}$ that yield a simplicial factorization by solving a single linear program. The factor $\bm{F}$ can be found by extracting columns of $\bm{C}$.

Suppose we observe a nonnegative matrix $\bm{X}$ whose rows sum to one. Assume that $\bm{X}=\bm{Y}+\bm{\Delta}$ where $\bm{Y}$ is a nonnegative matrix whose rows sum to one, which has a rank-$r$ separable factorization $\bm{Y}=\bm{FW}$ such that the rows of $\bm{W}$ are $\alpha$-robust simplicial, and where \big{\lVert}{\bm{\Delta}}\big{\rVert}_{\infty,1}\leq\epsilon. Define the polyhedral set

The set $\Phi(\bm{X})$ consists of matrices $\bm{C}$ that approximately locate a factorization of $\bm{X}$. We can prove the following result.

Suppose that $\bm{X}$ satisfies the assumptions stated in the previous paragraph. Furthermore, assume that for every row $\bm{Y}_{j,\cdot}$ that is not hott, we have the margin constraint $\|\bm{Y}_{j,\cdot}-\bm{Y}_{i,\cdot}\|\geq d_{0}$ for all hott rows $i$. Then we can find a nonnegative factorization satisfying \big{\lVert}{\bm{X}-\hat{\bm{F}}\hat{\bm{W}}}\big{\rVert}_{\infty,1}\leq 2\epsilon provided that $\epsilon<\tfrac{\min\{\alpha d_{0},\alpha^{2}\}}{9(r+1)}$. Furthermore, this factorization correctly identifies the hott topics appearing in the separable factorization of $\bm{Y}$.

Algorithm 3 requires the solution of two linear programs. The first minimizes a cost vector over $\Phi_{2\epsilon}(\bm{X})$. This lets us find $\hat{\bm{W}}$. Afterward, the matrix $\hat{\bm{F}}$ can be found by setting

Our robustness result requires a margin-type constraint assuming that the original configuration consists either of duplicate hott topics, or topics that are reasonably far away from the hott topics. On the other hand, under such a margin constraint, we can construct a considerably better approximation than that guaranteed by the AGKM algorithm. Moreover, unlike AGKM, our algorithm does not need to know the parameter $\alpha$.

The proofs of Theorems 3.1 and 3.2 can be found in the appendix. The main idea is to show that we can only represent a hott topic efficiently using the hott topic itself. Some earlier versions of this paper contained incomplete arguments, which we have remedied. For a signifcantly stronger robustness analysis of Algorithm 3, see the recent paper .

Having established these theoretical guarantees, it now remains to develop an algorithm to solve the LP. Off-the-shelf LP solvers may suffice for moderate-size problems, but for large-scale matrix factorization problems, their running time is prohibitive, as we show in Section 5. In Section 4, we turn to describe how to solve Algorithm 3 efficiently for large data sets.

2 Related Work

Localizing factorizations via column or row subset selection is a popular alternative to direct factorization methods such as the SVD. Interpolative decomposition such as Rank-Revealing QR and CUR have favorable efficiency properties as compared to factorizations (such as SVD) that are not based on exemplars. Factorization localization has been used in subspace clustering and has been shown to be robust to outliers .

In recent work on dictionary learning, Esser et al. and Elhamifar et al. have proposed a factorization localization solution to nonnegative matrix factorization using group sparsity techniques . Esser et al. prove asymptotic exact recovery in a restricted noise model, but this result requires preprocessing to remove duplicate or near-duplicate rows. Elhamifar shows exact representative recovery in the noiseless setting assuming no hott topics are duplicated. Our work here improves upon this work in several aspects, enabling finite sample error bounds, the elimination of any need to preprocess the data, and algorithmic implementations that scale to very large data sets.

Incremental Gradient Algorithms for NMF

The rudiments of our fast implementation rely on two standard optimization techniques: dual decomposition and incremental gradient descent. Both techniques are described in depth in Chapters 3.4 and 7.8 of Bertsekas and Tstisklis .

We aim to minimize $\bm{p}^{T}\operatorname{diag}(\bm{C})$ subject to $\bm{C}\in\Phi_{\tau}(\bm{X})$. To proceed, form the Lagrangian

with multipliers $\beta$ and $\bm{w}\geq\bm{0}$. Note that we do not dualize out all of the constraints. The remaining ones appear in the constraint set $\Phi_{0}=\{\bm{C}~{}:~{}\bm{C}\geq\bm{0},\ \operatorname{diag}(\bm{C})\leq 1,\mbox{ and }C_{ij}\leq C_{jj}~{}\mbox{for all}~{}i,j\}$.

Dual subgradient ascent solves this problem by alternating between minimizing the Lagrangian over the constraint set $\Phi_{0}$, and then taking a subgradient step with respect to the dual variables

where $\bm{C}^{\star}$ is the minimizer of the Lagrangian over $\Phi_{0}$. The update of $w_{i}$ makes very little difference in the solution quality, so we typically only update $\beta$.

We minimize the Lagrangian using projected incremental gradient descent. Note that we can rewrite the Lagrangian as

Here, $\operatorname{supp}(\bm{x})$ is the set indexing the entries where $\bm{x}$ is nonzero, and $\mu_{j}$ is the number of nonzeros in row $j$ divided by $n$. The incremental gradient method chooses one of the $n$ summands at random and follows its subgradient. We then project the iterate onto the constraint set $\Phi_{0}$. The projection onto $\Phi_{0}$ can be performed in the time required to sort the individual columns of $\bm{C}$ plus a linear-time operation. The full procedure is described in Appendix B. In the case where we expect a unique solution, we can drop the constraint $C_{ij}\leq C_{jj}$, resulting in a simple clipping procedure: set all negative items to zero and set any diagonal entry exceeding one to one. In practice, we perform a tradeoff. Since the constraint $C_{ij}\leq C_{jj}$ is used solely for symmetry breaking, we have found empirically that we only need to project onto $\Phi_{0}$ every $n$ iterations or so.

This incremental iteration is repeated $n$ times in a phase called an epoch. After each epoch, we update the dual variables and quit after we believe we have identified the large elements of the diagonal of $\bm{C}$. Just as before, once we have identified the hott rows, we can form $\bm{W}$ by selecting these rows of $\bm{X}$. We can find $\bm{F}$ just as before, by solving (2). Note that this minimization can also be computed by incremental subgradient descent. The full procedure, called Hottopixx, is described in Algorithm 4.

For small-scale problems, Hottopixx can be implemented in a few lines of Matlab code. But for the very large data sets studied in Section 5, we take advantage of natural parallelism and a host of low-level optimizations that are also enabled by our formulation. As in any numerical program, memory layout and cache behavior can be critical factors for performance. We use standard techniques: in-memory clustering to increase prefetching opportunities, padded data structures for better cache alignment, and compiler directives to allow the Intel compiler to apply vectorization.

Note that the incremental gradient step (step 6 in Algorithm 4) only modifies the entries of $\bm{C}$ where $\bm{X}_{\cdot k}$ is nonzero. Thus, we can parallelize the algorithm with respect to updating either the rows or the columns of $\bm{C}$. We store $\bm{X}$ in large contiguous blocks of memory to encourage hardware prefetching. In contrast, we choose a dense representation of our localizing matrix $\bm{C}$; this choice trades space for runtime performance.

Each worker thread is assigned a number of rows of $\bm{C}$ so that all rows fit in the shared L3 cache. Then, each worker thread repeatedly scans $\bm{X}$ while marking updates to multiple rows of $\bm{C}$. We repeat this process until all rows of $\bm{C}$ are scanned, similar to the classical block-nested loop join in relational databases .

Experiments

Except for the speedup curves, all of the experiments were run on an identical configuration: a dual Xeon X650 (6 cores each) machine with 128GB of RAM. The kernel is Linux 2.6.32-131.

In small-scale, synthetic experiments, we compared Hottopixx to the AGKM algorithm and the linear programming formulation of Algorithm 3 implemented in Matlab. Both AGKM and Algorithm 3 were run using CVX coupled to the SDPT3 solver . We ran Hottopixx for $50$ epochs with primal stepsize 1e-1 and dual stepsize 1e-2. Once the hott topics were identified, we fit $\bm{F}$ using two cleaning epochs of incremental gradient descent for all three algorithms.

Because we ran over $2000$ experiments with $405$ different parameter settings, it is convenient to use the performance profiles to compare the performance of the different algorithms . Let $\mathcal{P}$ be the set of experiments and $\mathcal{A}$ denote the set of different algorithms we are comparing. Let $Q_{a}(p)$ be the value of some performance metric of the experiment $p\in\mathcal{P}$ for algorithm $a\in\mathcal{A}$. Then the performance profile at $\tau$ for a particular algorithm is the fraction of the experiments where the value of $Q_{a}(p)$ lies within a factor of $\tau$ of the minimal value of $\min_{b\in\mathcal{A}}Q_{b}(p)$. That is,

In a performance profile, the higher a curve corresponding to an algorithm, the more often it outperforms the other algorithms. This gives a convenient way to contrast algorithms visually.

Our performance profiles are shown in Figure 2. The first two figures correspond to experiments with $f=40$ and $n=400$. The third figure is for the synthetic experiments with all other values of $f$ and $n$. In terms of $(\infty,1)$-norm error, the linear programming solver typically achieves the lowest error. However, using SDPT3, it is prohibitively slow to factor larger matrices. On the other hand, Hottopixx achieves better noise performance than the AGKM algorithm in much less time. Moreover, the AGKM algorithm must be fed the values of $\epsilon$ and $\alpha$ in order to run. Hottopixx does not require this information and still achieves about the same error performance.

We also display a graph for running only four epochs (hott (fast)). This algorithm is by far the fastest algorithm, but does not achieve as optimal a noise performance. For very high levels of noise, however, it achieves a lower reconstruction error than the AGKM algorithm, whose performance degrades once $\eta$ approaches or exceeds $1$ (Figure 2(f)). We also provide performance profiles for the root-mean-square error of the nonnegative matrix factorizations (Figure 2 (d) and (e)). The performance is qualitatively similar to that for the $(\infty,1)$-norm.

We also coded Hottopixx in C++, using the design principles described in Section 4.1, and ran on three large data sets. We generated a large synthetic example (jumbo) as above with $r=100$. We generated a co-occurrence matrix of people and places from the ClueWeb09 Dataset , normalized by TFIDF. We also used Hottopixx to select features from the RCV1 data set to recognize the class CCAT . The statistics for these data sets can be found in Table 1.

In Figure 3 (left), we plot the speed-up over a serial implementation. In contrast to other parallel methods that exhibit memory contention , we see superlinear speed-ups for up to 20 threads due to hardware prefetching and cache effects. All three of our large data sets can be trained in minutes, showing that we can scale Hottopixx on both synthetic and real data. Our algorithm is able to correctly identify the hott topics on the jumbo set. For clueweb, we plot the RMSE Figure 3 (middle). This curve rolls off quickly for the first few hundred topics, demonstrating that our algorithm may be useful for dimensionality reduction in Natural Language Processing applications. For RCV1, we trained an SVM on the set of features extracted by Hottopixx and plot the misclassification error versus the number of topics in Figure 3 (right). With $1500$ hott topics, we achieve $7\%$ misclassification error as compared to $5.5\%$ with the entire set of features.

Discussion

This paper provides an algorithmic and theoretical framework for analyzing and deploying any factorization problem that can be posed as a linear (or convex) factorization localizing program. Future work should investigate the applicability of Hottopixx to other factorization localizing algorithms, such as subspace clustering, and should revisit earlier theoretical bounds on such prior art.

The authors would like to thank Sanjeev Arora, Michael Ferris, Rong Ge, Nicolas Gillis, Ankur Moitra, and Stephen Wright for helpful suggestions. BR is generously supported by ONR award N00014-11-1-0723, NSF award CCF-1139953, and a Sloan Research Fellowship. CR is generously supported by NSF CAREER award under IIS-1054009, ONR award N000141210041, and gifts or research awards from American Family Insurance, Google, Greenplum, and Oracle. JAT is generously supported by ONR award N00014-11-1002, AFOSR award FA9550-09-1-0643, and a Sloan Research Fellowship.

References

Appendix A Proofs

Let $\bm{Y}$ be a nonnegative matrix whose rows sum to one. Assume that $\bm{Y}$ admits an exact separable factorization of rank $r$. In other words, we can write $\bm{Y}=\bm{FW}$ where the rows of $\bm{W}$ are $\alpha$-robust simplicial and

for some permutation $\bm{\Pi}$. Let $I$ denote the indices of the rows in $\bm{Y}$ that correspond with the identity matrix in the factorization, which we have called the hott rows. Then we can write each row $j$ that is not hott as a convex combination of the hott rows:

As we have discussed, we may assume that $\sum_{k}M_{jk}=1$ for each $j\notin I$ because each row of $\bm{Y}$ sums to one.

The first lemma offers a stronger bound on the coefficients $M_{jk}$ in terms of the distance between row $j$ and the hott rows.

Proof Let us introduce notation for the quantity of interest: $w_{i}=w_{i}(\bm{c})=\sum_{j\in\mathcal{B}_{\delta}(\bm{X}_{i\cdot})}c_{j}$. We may assume that $w_{i}<1$, or else the result holds trivially. Since the entries of $\bm{c}$ sum to one, we have

To establish the result, we may as well assume that $w_{i}(\bm{c})$ achieves its minimum possible value subject to the constraints that the value of $\bm{c}^{T}\bm{Y}$ is fixed and that $\bm{c}$ is a nonnegative vector that sums to one. We claim that this minimum such $w_{i}$ occurs if and only if $c_{j}=0$ for all $j\in\mathcal{B}_{\delta}(i)\setminus\{i\}$. We complete the proof under this additional surmise.

Let us continue. Owing to the assumption that $\bm{Y}_{i\cdot}$ is no farther than $\tau$ from $\bm{c}^{T}\bm{Y}$, we have

The first line follows when we split the sum over $j$ based on whether or not the components fall in $\mathcal{B}_{\delta}(i)$. Then we apply the property that $c_{j}=0$ for $j\in\mathcal{B}_{\delta}(i)\setminus\{i\}$, and we identify the quantity $w_{i}$. In the last line, we factored out $1-w_{i}$, and we introduced the separable factorization of $\bm{Y}$.

and note that $\pi_{k}\geq 0$. Furthermore,

because the rows of $\bm{M}$ sum to one and because of the definition of $w_{i}$. Lemma A.1 implies that $\pi_{i}$ satisfies the bound

Indeed, the lemma is valid because $\bm{Y}_{j\cdot}$ is at least a distance of $\delta$ away from $\bm{Y}_{i\cdot}$ for every $j\notin\mathcal{B}_{\delta}(i)$.

With these observations, we can continue our calculation from (4):

This result is almost obvious when there are no duplicated rows. Indeed, since the hott topics form a simplicial set and the matrix $\bm{Y}$ admits a separable factorization, the only way we can represent all $r$ hott topics exactly is to have $C_{ii}=1$ for every hott row $i$. This exhausts the trace constraint, and we see that every other diagonal entry $C_{kk}=0$ for every not hott row $k$. The only matrices that are feasible identify the hott rows on the diagonal. They must represent the remaining rows using linear combinations of the hott topics because of the constraints $\bm{CY}=\bm{Y}$ and $C_{ij}\leq C_{jj}$. It follows that the only feasible matrices are factorization localizing matrices.

When there are duplicated rows, the analysis is slightly more delicate. By the same argument as above, all the weight on the diagonal must be concentrated on hott rows. But the objective $\bm{p}^{T}\operatorname{diag}(\bm{C})$ ensures that, out of any set of duplicates of a given topic, we always pick the duplicate row $j$ where $p_{j}$ is smallest; otherwise, we could reduce the objective further. Therefore, the diagonal of $\bm{C}$ identifies all $r$ distinct hott topics, and we select each one duplicate of each topic. As before, the other constraints ensure that the remaining rows are represented with this distinguished choice of hott topic exemplars. Therefore, the only minimizers are factorization localizing matrices that identify each hott topic exactly once.

A.2 Proof of Theorem 3.2

Let $\bm{X}=\bm{Y}+\bm{\Delta}$. The matrix $\bm{X}$ is the observed data, with rows scaled to have unit sum, and the perturbation matrix $\bm{\Delta}$ satisfies \big{\lVert}{\bm{\Delta}}\big{\rVert}_{\infty,1}\leq\epsilon. We assume that $\bm{Y}$ is a nonnegative matrix whose rows sum to one, and we posit that it admits a rank-$r$ separable NMF $\bm{Y}=\bm{FW}$ where $\bm{W}$ is $\alpha$-robust simplicial. We write $I$ for the set of rows corresponding to hott topics in $\bm{Y}$.

Suppose that $\bm{C}_{0}$ is a factorization localizing matrix for the underlying matrix $\bm{Y}$. That is, $\bm{C}_{0}\bm{Y}=\bm{Y}$ and each row of $\bm{C}_{0}$ sums to one. It follows that

Using our decomposition $\bm{X}=\bm{Y}+\bm{\Delta}$, we quickly verify that

The point here is that a factorization localizing matrix for $\bm{Y}$ serves as an approximate factorization localizing matrix for $\bm{X}$.

Our approach for constructing an approximate factorization of $\bm{X}$ requires us to minimize a cost function $\bm{t}^{T}\operatorname{diag}(\bm{C})$ over the constraint set

Note that the factorization localizing matrix $\bm{C}_{0}$ for $\bm{Y}$ is a member of this set, so the optimization problem we solve in Theorem 3.2 is feasible.

Suppose that $\bm{C}\in\Phi_{2\epsilon}(\bm{X})$ is arbitrary. Let us check that the row sums of $\bm{C}$ are not much larger than one. To that end, note that

We have twice used the fact that every row of $\bm{X}$ sums to one. For any row $\bm{c}$ of the matrix $\bm{C}$, this formula yields $\bm{c}^{T}\bm{1}\leq 1+2\epsilon$ since \big{\lVert}{\bm{CX}-\bm{X}}\big{\rVert}_{\infty,1}\leq 2\epsilon. As a consequence,

In particular, every matrix $\bm{C}$ in the set $\Phi_{2\epsilon}(\bm{X})$ has $C_{ii}\geq 1-(8\epsilon+4\epsilon^{2})/\min\{\alpha d_{0},\alpha^{2}\}$ for each hott topic $i$. To ensure that hott topic $i$ has weight $C_{ii}$ greater than $1-1/(r+1)$ for each $i$, we need

Since there are $r$ hott rows, they carry total weight greater than $r(1-1/(r+1))$. Given the trace constraint, that leaves less than $1-1/(r+1)$ for the remaining rows. We see that each of the $r$ hott rows must carry more weight than every row that is not hott, so we can easily identify them.

Once we have identified the set $I$ of hott topics, we simply solve the second linear program

to find a $2\epsilon$-accurate factorization.

To project onto the set $\Phi_{0}$, note that we can compute the projection one column at a time. Moreover, the projection for each individual column amounts to (after permuting the entries of the column),

Assume, again without loss of generality, that we want to project a vector $\bm{z}$ with $z_{2}\geq z_{3}\geq\ldots\geq z_{n}$. Then we need to solve the quadratic program

The optimal solution can be found as follows. Let $k_{c}$ be the largest $k\in\{2,\ldots,f\}$ such that

where $\Pi_{}$ denotes the projection onto the interval $$. Set

Then $\hat{\bm{x}}$ is the optimal solution. A linear time algorithm for computing $\hat{x}$ is given by Algorithm 5

To prove that $\hat{\bm{x}}$ is optimal, define

$y_{i}$ is the gradient of $\tfrac{1}{2}\|\bm{x}-\bm{z}\|^{2}$ at $\hat{x}$. Consider the LP

$\hat{\bm{x}}$ is an optimal solution for this LP because the cost is negative on the negative entries, on the nonnegative entries that are larger than $k_{c}$, positive for $2\leq k\leq k_{c}$, and nonpositive for $k=1$. Hence, by the minimum principle, $\hat{\bm{x}}$ is also a solution of (8).