Fast Conical Hull Algorithms for Near-separable Non-negative Matrix Factorization

Introduction

A data matrix $\mathbf{X}$ of size $m\times n$ is said to admit a Non-negative Matrix Factorization (NMF) with inner-dimension $r$, if $\mathbf{X}$ can be expressed as $\mathbf{X}=\mathbf{W}\mathbf{H}$ where $\mathbf{W},\mathbf{H}$ are two non-negative matrices of dimensions $m\times r$ and $r\times n$ respectively. In many applications, a compact (i.e., small $r$) approximate NMF tends to provide a natural and interpretable part-based decomposition of the data (Lee & Seung, 1999), often more appealing than other low-rank factorizations. NMFs arise pervasively in a variety of signal separation problems, such as modeling topics in text and hyperspectral image analysis (Cichocki et al., 2009).

Separability Assumption: The entire dataset, i.e. all columns of $\mathbf{X}$, reside in a cone generated by a small subset of $r$ columns of $\mathbf{X}$.

In algebraic terms, $\mathbf{X}=\mathbf{W}\mathbf{H}=\mathbf{X}_{A}\mathbf{H}$ so that the $r$ columns of $\mathbf{W}$ are hidden among the columns of $\mathbf{X}$ (indexed by an unknown subset of indices $A$). Equivalently, a corresponding subset of $r$ columns of $\mathbf{H}$ happen to constitute the $r\times r$ identity matrix. We refer to these columns as anchors (Arora et al., 2012). Informally, in the context of topic modeling problems where $\mathbf{X}$ is a document-word matrix and $\mathbf{W},\mathbf{H}$ are document-topic and topic-term associations respectively, the separability assumption equivalently posits the existence of special anchor words in the vocabulary, whose occurence uniquely identifies the presence of a topic, and whose usage across the corpus is collectively predictive of the usage of all the other words. The separability assumption was investigated earlier by Donoho & Stodden (2003) who showed that it implied uniqueness of the NMF solution, modulo permutation and scaling. In order to place our contributions in the right context, we first briefly provide a flavor of recently proposed separable NMF algorithms.

Related Work: Assuming that the columns of $\mathbf{X}$ are normalized to have unit $l_{1}$-norm, the separable NMF problem reduces to that of finding the extreme points (that is, points inexpressible as convex combinations of other points) of the convex hull of the columns (Arora et al., 2012). A Linear Program (LP) can be setup to attempt to express a given column as a convex combination of the other columns. If this LP declares infeasibility, an extreme point is identified. This approach (Arora et al., 2012, Section 5) requires solving $n$ feasibility LP’s each involving $n-1$ variables which is not scalable for many problems of interest. A noise-robust version of the procedure further requires knowledge of parameters that are hard to estimate apriori. Bittorf et al. (2012) formulate a single LP whose solution resolves the exactly separable NMF problem. An extension is also developed for noise-robustness. Instead of invoking a general LP solver, a specialized algorithm is derived based on an incremental stochastic gradient descent procedure, and its parallel (multithreaded) implementation is benchmarked on large datasets. On the other hand, this algorithm requires estimates of primal and dual step sizes, converges only asymptotically, and does not explicitly exploit the sparsity of the final solution. Gillis & Vavasis (2012) develop a highly scalable approach closely related to rank-revealing QR factorizations for column subset selection. A perturbation analysis of this algorithm under noise is also presented. In Esser et al. (2012), column subset selection is cast essentially as a form of multivariate regression with row-sparsity inducing norms, e.g., see Bien et al. (2010). Algorithms derived in this framework are asymptotically convergent, and sensitive to near-duplicate columns, making it necessary to perform certain adhoc preprocessing steps.

Contributions: We present a new family of highly scalable and empirically noise-robust algorithms for separable NMFs, with several favorable properties:

The algorithms produce a correct solution for the separable case after exactly $r$ iterations. They require no additional parameters. They are closely related to convex and conical hull finding procedures proposed in the computational geometry literature (Clarkson, 1994; Dula et al., 1998). Computationally, the algorithms bear some resemblence to simultaneous Orthogonal Matching Pursuit (Buhlmann & Geer, 2010; Tropp et al., 2006) for sparse greedy reconstruction of multiple target variables from the same subset of input variables. We also derive a variant based on this connection that performs quite well under noise.

Under controlled noise conditions in synthetic datasets and on real-world topic modeling problems, our algorithms consistently outperform other separable NMF techniques with respect to multiple performance metrics. Our methods are highly competitive with existing non-convex NMF algorithms, but are free of sub-optimal local minima and associated initialization issues.

The solution for $(r-1)$ target anchors is contained in the solution for $r$ target anchors (unlike non-convex NMF methods), which makes it easier to do model selection on real-world datasets by keeping track of performance on a validation set.

The algorithms are highly scalable and have small memory footprint. The sparsity of the data, the intermediate variables and the final solution is carefully exploited in a high-performance parallel and distributed implementation which scales excellently on both shared- and distributed-memory machines. For example, a twitter corpus with 125-thousand tweets can be factorized for $r=100$ in less than 10 seconds on a commodity 8-core machine.

Unlike all existing algorithms, no column normalization is needed. Such normalization interferes with the TFIDF weightings routinely used in text modeling applications, leading to performance loss.

Unlike Esser et al. (2012), the algorithms do not require any special preprocessing to eliminate duplicate or near-duplicate columns.

Fast Conical Hull Algorithms

An informal description: Figure 2 provides some geometric intuition underlying the proposed approach. The algorithm executes $r$ iterations. In each iteration a new anchor column is identified. This corresponds to expanding a cone one extreme ray at a time, until the entire dataset is eventually contained in the cone defined by the full set of anchors. Figure 2 illustrates one step of the algorithm where there is an existing cone defined by three extreme rays (marked 1 to 3). To identify the next extreme ray, the algorithm picks a point outside the current cone (a green point) and projects it to the current cone to compute a residual vector (we call this the projection step). This residual vector separates the current cone from at least one non-selected extreme ray that can be found by maximizing a specific selection criteria (we call this the detection step). Intuitively, the algorithm picks a face of the current cone (spanned by rays 1 and 3 in Figure 2) that “sees” exterior points and rotates this face towards the exterior until it hits the “last” point. In the example shown in Figure 2, ray 4 is identified as a new extreme ray.

These geometric intuitions are inspired by Clarkson (1994); Dula et al. (1998) who present LP-based algorithms for general convex and conical hull problems. Their algorithms are also directly applicable in our NMF setting, provided the data satisfies the separability assumption exactly. In this case, the residual of any single exterior point can be used to correctly expand the cone as described above. However, anchor detection criteria derived from multiple residuals demonstrates radically superior noise robustness, as we report in the experimental section. The emphasis on scalability and noise-robustness thus leads us to a new family of algorithms whose implementation (and associated proof of correctness) is distinct from prior work.

Algorithm 1 details the steps of the proposed family of algorithms which we call Xray . Each iteration consists of two steps: (i) a detection step that finds a column(s) of $\mathbf{X}$ to be added as an anchor, and (ii) a projection step where all data points are projected onto the current cone to get the residuals. Projection is done by solving simultaneous nonnegative least squares problem using Algorithm 2. Every residual vector $\mathbf{R}_{i}$ obtained after the projection step is normal to one of the faces of the current cone. In the selection step, we pick a face of the current cone (identified by its normal $\mathbf{R}_{i}$), normalize all the data points to lie on the hyperplane $\mathbf{p}^{T}\mathbf{x}=1$ $\left(\mathbf{Y}_{j}=\frac{\mathbf{X}_{j}}{\mathbf{p}^{T}\mathbf{X}_{j}}\right)$ for a strictly positive vector $\mathbf{p}$, and expand the current cone by selecting an extreme ray that maximizes the inner product $\mathbf{R}_{i}^{T}\mathbf{Y}_{j}$. The selection step can be implemented in various ways - some options are listed in Algorithm 1.

To show that Xray correctly identifies all the extreme rays, we need the following lemmas.

The residual matrix $\mathbf{R}$, obtained after projection of columns of $\mathbf{X}$ onto the current cone satisfies $\mathbf{R}^{T}\mathbf{X}_{A}\leq 0$, where $\mathbf{X}_{A}$ are the extreme rays of the current cone.

For any point $\mathbf{X}_{i}$ exterior to the current cone, we have $\mathbf{R}_{i}^{T}\mathbf{X}_{i}>0$, where $\mathbf{R}_{i}$ is the residual of $\mathbf{X}_{i}$ obtained by projecting it onto the current cone.

Using the above two lemmas, we prove the following theorem regarding the correctness of Algorithm 1.

The data point $\mathbf{X}_{j^{*}}$ added at each iteration in the Detection step of Algorithm 1, if the maximizer in Eqn. 1 is unique, is an extreme ray of $\mathcal{C}$ that has not been selected in previous iterations.

Let the index set $A$ identify all the extreme rays of $\mathcal{C}$. Under the separability assumption, we have $\mathbf{X}=\mathbf{X}_{A}\mathbf{H}$. Let the index set $A^{t}$ identify the extreme rays of the current cone $\mathcal{C}^{t}$.

Let $\mathbf{Y}_{j}=\frac{\mathbf{X}_{j}}{\mathbf{p}^{T}\mathbf{X}_{j}}$ and $\mathbf{Y}_{A}=\mathbf{X}_{A}[diag(\mathbf{p}^{T}\mathbf{X}_{A})]^{-1}$ (since $\mathbf{p}$ is strictly positive, the inverse exists). Hence $\mathbf{Y}_{j}=$ $\mathbf{Y}_{A}\frac{[diag(\mathbf{p}^{T}\mathbf{X}_{A})]\mathbf{H}_{j}}{\mathbf{p}^{T}\mathbf{X}_{j}}$. Let $\mathbf{C}_{j}=\frac{[diag(\mathbf{p}^{T}\mathbf{X}_{A})]\mathbf{H}_{j}}{\mathbf{p}^{T}\mathbf{X}_{j}}$. We also have $\mathbf{p}^{T}\mathbf{Y}_{j}=1$ and $\mathbf{p}^{T}\mathbf{Y}_{A}=\mathbf{1}^{T}$. Hence, we have $1=\mathbf{p}^{T}\mathbf{Y}_{j}=\mathbf{p}^{T}\mathbf{Y}_{A}\mathbf{C}_{j}=\mathbf{1}^{T}\mathbf{C}_{j}$.

Using Lemma 2.1, Lemma 2.2 and the fact that $\mathbf{p}$ is strictly positive, we have $\max_{1\leq j\leq n}\mathbf{R}_{i}^{T}\mathbf{Y}_{j}=\max_{j\notin A^{t}}\mathbf{R}_{i}^{T}\mathbf{Y}_{j}$. Indeed, for all $j\in A^{t}$ we have $\mathbf{R}_{i}^{t}\mathbf{Y}_{j}\leq 0$ using Lemma 2.1 and there is at least one $j=i\notin A^{t}$ for which $\mathbf{R}_{i}^{t}\mathbf{Y}_{j}>0$ using Lemma 2.2. Hence the maximum lies in the set $\{j:j\notin A^{t}\}$.

Further, we have $\max_{j\notin A^{t}}\mathbf{R}_{i}^{T}\mathbf{Y}_{j}=\max_{j\notin A^{t}}\mathbf{R}_{i}^{T}\mathbf{Y}_{A}\mathbf{C}_{j}\leq\max_{j\in(A\setminus A^{t})}\mathbf{R}_{i}^{\prime}\mathbf{Y}_{j}$. The second inequality is the result of the fact that $\lVert\mathbf{C}_{j}\rVert_{1}=1$ and $\mathbf{C}_{j}\geq 0$. This implies that if there is a unique maximum at a $j^{*}=\mathop{\rm arg\,max}_{j\notin A^{t}}\mathbf{R}_{i}^{T}\mathbf{Y}_{j}$, then $\mathbf{X}_{j^{*}}$ is generator of an extreme ray of the cone $\mathcal{C}$. ∎

Exterior Point Selection: It can be noted that residual of any point exterior to the current cone (i.e., any $\mathbf{R}_{i}\neq 0$) can be used in the selection step of Algorithm 1. This gives us multiple ways of expanding the current cone depending on which $i$ is chosen - all of which solve the separable problem but may behave very differently in the presence of noise. Some natural options are listed in Algorithm 1: choosing a random exterior point (Eqn. 2), one with maximum residual norm (Eqn. 3) or one which defines a normal to a supporting hyperplane of the current cone which “sees” maximum “mass” of points in its positive halfspace, as measured by Eqn. 4. In the experiments, we will refer to these variants as Xray (rand), Xray (max) and Xray (dist) respectively.

A Greedy variation for noisy data: In high dimensional noisy data almost all the points may masquerade as anchors. A natural choice is to expand the current cone greedily by selecting a point that best minimizes the current residual, i.e., $j^{*}=\mathop{\rm arg\,min}_{j}\min_{\mathbf{b}>0}\lVert\mathbf{R}-\mathbf{X}_{j}\mathbf{b}^{T}\rVert_{F}^{2}$. This selection criterion simplifies to Eqn.5 in Algorithm 1 (referred as Xray (greedy) henceforth). One may view this approach as implementing a nonnegative variant of simultaneous orthogonal matching pursuit (Tropp et al., 2006), which is a greedy approach to the problem of sparse regression of multiple response variables on the same subset of explanatory variables, i.e., for solving $\min_{\mathbf{B}\geq 0}\lVert\mathbf{X}-\mathbf{X}\mathbf{B}\rVert_{F}^{2}$ s.t. $\lVert\mathbf{B}\rVert_{0,1}=r$ where $\|\mathbf{B}\|_{0,1}$ pseudo-norm counts the number of non-zero rows in $\mathbf{B}$. In the context of separable NMF, both response variables and explanatory variables are the columns data matrix $\mathbf{X}$. A relaxed version of this problem is solved in Esser et al. (2012) ($\min_{\mathbf{B}\geq 0}\lVert\mathbf{X}-\mathbf{X}\mathbf{B}\rVert_{F}^{2}+\lambda\lVert\mathbf{B}\rVert_{1,\infty}$). It is also possible to have $\lVert\mathbf{B}\rVert_{1,2}$ penalized variant (Tropp, 2006; Bien et al., 2010) which is natural for sparse multivariate regression problems. Note that the greedy approach is not guaranteed to solve the separable NMF problem, but may perform well in the noisy settings as we observe in our experiments. Intuitively, this variant is concerned with greedily optimizing all residuals on average at every iteration, instead of making a decision based on the residual of a single, albeit well-chosen, exterior point.

2 Scalability and Parallelization

Here we describe various implementation details that allow us to gracefully scale to large sparse datasets (e.g., document-term matrices). The detection step can be parallelized by scoring the candidate anchors simultaneously. Likewise, the projection step involves solving Eqn. 6, which is separable in the columns of $\mathbf{B}$ and hence can be optimized in parallel.

Detection Step: We avoid materializing the dense residual matrix $\mathbf{R}$ in the evaluation of the anchor selection criteria. Instead, we score candidate anchors on-the-fly as we compute (but not explicitly materialize) a matrix $\mathbf{Q}=\left(\mathbf{R}^{T}\mathbf{X}\right)_{+}=\left(\mathbf{C}-(\mathbf{C}_{A}\mathbf{H})^{T}\right)_{+}$ where $\mathbf{C}=\mathbf{X}^{T}\mathbf{X}$ denotes a covariance matrix (word-by-word for topic modeling applications). Here, the potential sparsity, symmetry of the covariance matrix $\mathbf{C}$ as well as the non-negativity of $\mathbf{H}$ can be further exploited. For example, if $\mathbf{C}_{ij}=0$, the corresponding entry in the product $(\mathbf{C}_{A}\mathbf{H})$ need not be computed, since the resulting negative value is anyway reset to zero by the $(\cdot)_{+}$ thresholding operator. On a $P$ core machine, the selection criteria may be evaluated in $O(\frac{nnz(\mathbf{C})r}{P})$ time where $nnz(\mathbf{C})$ is the number of non-zeros in $\mathbf{C}$. If $\mathbf{C}$ is dense, we compute $\mathbf{Q}$ using parallel dense BLAS-3 operations. The one time computation of $\mathbf{C}$ is done via a parallel aggregation of rank-one outer-product terms defined by the rows of $\mathbf{X}$.

Projection Step: Algorithm 2 gives the steps of a cyclic block coordinate descent algorithm organized around very light-weight incremental sparsity-exploiting updates for solving Eqn. 6 (derivation omitted for brevity). The algorithm can be invoked in parallel on columns of $\mathbf{X}$ to compute the corresponding columns of $\mathbf{B}$. The previous value of $\mathbf{B}$ is used to warm start the optimization and typically a very small number of iterations is needed for convergence.

Empirical Observations

Here, we report extensive comparisons on synthetic and medium-scale topic modeling problems, and benchmark our parallel implementation on large text datasets on multicore machines and distributed systems. We compare with the methods proposed in Bittorf et al. (2012) (abbrv. as Hottopixx ) and Gillis & Vavasis (2012) (abbrv. as GV), as well as traditional NMFs based on alternating optimization (Cichocki et al., 2009). The source codes for Hottopixx and GV were taken from the respective authors’ websites. In comparisons with Esser et al. (2012), it was observed that it tends to select near-duplicate anchors, as also mentioned in Esser et al. (2012). This characteristic causes it to consistently perform less favorably compared to other methods unless the data is preprocessed in an adhoc fashion to remove similar columns of $\mathbf{X}$; hence we do not include it in our list of baselines. We also do not compare with Arora et al. (2012) since Hottopixx reportedly performs better (Bittorf et al., 2012) and the algorithm requires parameters which are hard to guess apriori.

2 Medium-scale Topic modeling problems

Classification experiments: Figure 4 shows the classification accuracy results obtained with the features (columns of the document-term matrix restricted to anchor words) selected by different methods on the three datasets. Black dotted line is the classification accuracy with full features (all the words). We use $5\%$ of the documents for training and the rest $95\%$ for testing to emulate a semi-supervised learning scenario where we view various methods as inducing a topical representation based on all (unlabeled) data. We use multiclass SVM classifier as implemented in LIBLINEAR (Fan et al., 2008) and use four-fold cross validation to select the parameter $C$. Among separable NMF techniques, the proposed Xray (greedy) and Xray (dist) (with exception on Reuters) outperform Hottopixx and GV on all the three datasets, more so on TDT. On average, traditional NMFs with local optimization perform quite well on these datasets especially when $r$ is small, but can show significant performance variance (shown as error bars) with respect to initialization. As the number of topics increases, the performance gap between the proposed methods and the local optimization method rapidly diminishes. In this regime our techniques are a viable alternative to local optimization methods, and have the advantage of being local-minima-free, i.e., eliminating uncertainty with respect to initialization and therefore not requiring multiple runs.

Clustering experiments: We also evaluated clustering performance by assigning a cluster label to each document based on the maximum element in the corresponding row of $\mathbf{W}$. We refine the solution with a few iterations of alternating optimization. Figure 5 shows the clustering performance in terms of Normalized Mutual Information (NMI) as these iterations proceed. We also show the NMI obtained with local search method after it has converged to a local optimum (averaged NMI from ten runs with different random initializations is shown; error-bar indicates the variation around the average). Again, the proposed Xray methods are among the best performing methods in terms of clustering performance and do not require multiple runs as traditional NMFs do.

Quality of anchor words: Qualitatively, we found that anchor words selected by the proposed Xray methods tend to be more representative of the topics compared to those selected by Hottopixx and GV. Table 2 shows top words and anchors for a few topics (Lewinsky scandal, Iraq nuclear program, National Tobacco Settlement, Indonesia riots of 1998 and Columbia space shuttle) extracted from the TDT dataset.

3 Large-scale Experiments

We implemented a shared- and distributed-memory parallel version of Xray in C++. That is, our implementation can exploit parallelism when running on multi-core machines, or on clusters of multi-core machines. For shared-memory parallelism, we use PFunc (Kambadur et al., 2009), a lightweight and portable library that provides C and C++ APIs to express task parallelism. For distributed-memory parallelism, we use MPI\urlhttp://www.mpi-forum.org/, a popular library specification for message-passing that is used extensively in high-performance computing.

To test the shared-memory performance and scalability of Xray , we ran experiments on daniel, a dual-socket, quad-core Intel® Xeon™ X5570 machine with 64GB of RAM running Linux Kernel 2.6.35-24 (total 8 cores). For compilation, we used GCC v4.4.5 with: “-O3 -fomit-frame- pointer -funroll-loops” in addition to PFunc 1.02, OpenMPI 1.4.5 and untuned ATLAS BLAS. We ran large-scale experiments on three datasets: RCV1 (Lewis et al., 2004), co-occurence matrix of people and places from ClueWeb09 dataset (Lemur, ), and IBM Twitter (IBMT) dataset. The statistics relating to these three large datasets are presented in Table 3.We report scalability results for Xray (greedy) - other variants are computationally very similar.

Figure 6 depicts the multi-threaded performance of our implementation on daniel while detecting 100 topics. Our implementation is able to factorize RCV1 in 409 seconds on 8 cores and achieve $4.2x$ speedup over 8 threads when compared to the sequential implementation. Similarly, for IBMT we achieve $4.5x$ speedup, while completing the factorization in $9.8$ seconds on 8 cores. For the dense $\mathbf{X}^{T}\mathbf{X}$ case, we are able to factorize PPL2 in 1147 seconds with just 8 cores. We believe that further speedup improvements can be demonstrated on these problems by (a) optimizing the data layout of various sparse matrices to alleviate memory contention amongst threads, and (b) in dense problems such as PPL2, by using a version of BLAS tuned to our architecture and by reorganizing our implementation around more BLAS-3 operations that have better memory to compute ratio than BLAS-1 or 2 operations. Our implementation showed good scalability on distributed-memory machines as well (details omitted for brevity).

To compare our performance against the state-of-the-art Hottopixx algorithm (Bittorf et al., 2012), we ran their algorithm on daniel with the options “--dual 0.01 --epochs 10 --splits 8 --hott --normse 1 --primal 1e-6” set in close consultation with the authors. A detailed comparison is shown in Table 4. A head-to-head comparison is difficult because of the different performance characteristics of Hottopixx and Xray . For example, Hottopixx ’s per-epoch runtime is not dependent on $r$, the number of topics, but it’s accuracy is dependent on $E$, the number of epochs, while our methods execute exactly $r$ iterations, where each iteration has a superlinear dependence on $r$. Nonetheless, for all three datasets, we see that Xray performs better than Hottopixx even when Hottopixx is run only for 5 epochs. In particular, for the sparse datasets IBMT and RCV1, Xray runs to completion in significantly shorter amount of time than Hottopixx .

Conclusions and Future Work

Our methods perform favorably in comparison to other recently proposed separable NMF algorithms and offer highly scalable local-minima-free alternatives to existing local optimization techniques. Future work includes a formal noise analysis of the proposed algorithms, investigating the streaming setting where documents or words arrive in an online fashion, and using our models for social media content analysis.

Acknowledgments: We thank Victor Bittorf and Ben Recht for graciously providing their code, associated parameters and technical support. We thank Haim Avron, Christos Boutsidis, Ken Clarkson, Rick Lawrence and Ankur Moitra for insightful and enthusiastic discussions.

Research was sponsored by the U.S. Defense Advanced Research Projects Agency (DARPA) under the Social Media in Strategic Communication (SMISC) program, Agreement Number W911NF-12-C-0028. The views and conclusions contained in this document are those of the author(s) and should not be interpreted as representing the official policies, either expressed or implied, of the U.S. Defense Advanced Research Projects Agency or the U.S. Government. The U.S. Government is authorized to reproduce and distribute reprints for Government purposes notwithstanding any copyright notation hereon.

References