Random Projections for the Nonnegative Least-Squares Problem

Introduction

The Nonnegative Least Squares (NNLS) problem is a constrained least-squares regression problem where the variables are allowed to take only nonnegative values. More specifically, the NNLS problem is defined as follows:

NNLS is a quadratic optimization problem with linear inequality constraints. As such, it is a convex optimization problem and thus it is solvable (up to arbitrary accuracy) in polynomial time . In words, NNLS seeks to find the best nonnegative vector $x_{opt}$ in order to approximately express $b$ as a strictly nonnegative linear combination of the columns of $A$, i.e., $b\approx Ax_{opt}$.

The motivation for NNLS problems in data mining and machine learning stems from the fact that given least-squares regression problems on nonnegative data such as images, text, etc., it is natural to seek nonnegative solution vectors. (Examples of data applications are described in .) NNLS is also useful in the computation of the Nonnegative Matrix Factorization , which has received considerable attention in the past few years. Finally, NNLS is the core optimization problem and the computational bottleneck in designing a class of Support Vector Machines . Since modern datasets are often massive, there is continuous need for faster, more efficient algorithms for NNLS.

The following theorem is the main quality-of-approximation result for our randomized NNLS algorithm.

holds with probability at least $0.5$Note that a small number of repetitions of the algorithm suffices to boost its success probability.. The running time of the RandomizedNNLS algorithm is

The latter term corresponds to the time required to exactly solve an NNLS problem on an input matrix of dimensions $r\times d$.

One should compare the running time of our method to $T_{NNLS}\left(n,d\right)$, which corresponds to the time required to solve the NNLS problem exactly. We experimentally evaluate our approach on 3,000 NNLS problems constructed from a large and sparse term-document data collection. On average (see section 4.1), the proposed algorithm achieves a three-fold speedup when compared to a state-of-the-art NNLS solver with a small (approx. $10\%$) loss in accuracy; a two-fold speedup is achieved with a $4\%$ loss in accuracy. Computational savings are more pronounced for NNLS problems with denser input matrices $A$ and vectors $b$ (see section 4.2).

The remainder of the paper is organized as follows. Section 2 reviews basic linear algebraic definitions and discusses related work. In Section 3 we present our randomized algorithm for approximating the NNLS problem, discuss its running time, and give the proof of Theorem 1. Finally, in section 4 we provide an experimental evaluation of our method.

Background and related work

The (non-normalized) $n\times n$ matrix of the Hadamard-Walsh transform $H_{n}$ is defined recursively as follows:

The $n\times n$ normalized matrix of the Hadamard-Walsh transform is equal to $\frac{1}{\sqrt{n}}H_{n}$; hereafter, we will denote this normalized matrix by $H_{n}$ ($n$ is a power of $2$). For simplicity, throughout this paper we will assume that $n$ is a power of two; padding $A$ and $b$ with all-zero rows suffices to remove the assumption. Finally, all logarithms are base two.

2 Algorithms for the NNLS problem

A Random Projection Type Algorithm for the NNLS problem

This section describes our main algorithm for the NNLS problem. Our algorithm employs a randomized Hadamard transform to construct a much smaller NNLS problem and solves this smaller problem exactly using a standard NNLS solver. The approximation accuracy of our algorithm is a function of the size of the small NNLS problem.

In section 3.3 we will argue that, for any $\epsilon\in(0,1/3]$, if we set

2 The proof of Theorem 1

for some $\beta\in(0,1]$, and the parameter $r$ satisfies

for a sufficiently large constant $c_{o}$, then with probability at least $2/3$ all $d$-dimensional vectors $y$ satisfy,

Using Theorem 7 of (originally proven by Rudelson and Virshynin in ), we see that

for a sufficiently large constant $c_{o}$ ($c_{o}$ is not specified in ). Markov’s inequality implies that

with probability at least 2/3. Finally, using equation (9) concludes the proof of the lemma. $\diamond$

2.2 Another useful result

Let $U$ be an $n\times d$ orthogonal matrix ($n\geq 20$ and $n\geq d$). Then, for all $i\in[n]$,

2.3 The proof of Theorem 1

We are now ready to prove Theorem 1. We apply lemma 1 for \Phi=\left[\begin{array}[]{cc}H_{n}DA&-H_{n}Db\\ \end{array}\right]\in R^{n\times(d+1)}, the parameter $r$ of Theorem 1, sampling probabilities $p_{i}=1/n$, for all $i\in[n]$, $\beta=1/\left(4.2\log n\right)$, and $\epsilon^{\prime}=\epsilon/3\in(0,1/3]$, where $\epsilon\in(0,1]$ is the parameter of Theorem 1. Let $U_{\Phi}$ be the $n\times(d+1)$ matrix of the left singular vectors of $\Phi$. Note that $U_{\Phi}$ is exactly equal to $U_{\Phi}=H_{n}DU_{[A\hskip 3.61371pt-b]}$, where $U_{[A\hskip 3.61371pt-b]}$ is the $n\times(d+1)$ matrix of the left singular vectors of $[A\hskip 3.61371pt-b]$. Lemma 2 for $U_{[A\hskip 3.61371pt-b]}$ and our choice of $\beta$, guarantee that for all $i\in[n]$, with probability at least $0.9$

Manipulating equations (14) and (15) we get

3 What is the minimal value of r𝑟r ?

To derive values of $r$ for which the RandomizedNNLS algorithm satisfies the relative error guarantees of Theorem 1, we need to solve equation (2); this is hard since the solution depends on the Lambart $W$ function. Thus, we identify a range of values of $r$ that are sufficient for our purposes. Using the fact that for any $\alpha\geq 4$, and for any $\gamma\geq 2\alpha\log(\alpha)$,

and by setting $\alpha=342c_{o}^{2}(d+1)\log(n)/\epsilon^{2}$ in equation (2) (note that $342c_{o}^{2}(d+1)\log(n)/\epsilon^{2}\geq 4$), it can be proved that every $r$ such that

satisfies the inequality 2 ($c_{o}$ is the constant of Theorem 1).

4 Running time analysis

After this preconditioning step, we employ an NNLS solver on the smaller problem. The computational cost of the NNLS solver on the small problem was denoted as $T_{NNLS}(r,d)$ in Theorem 1. $T_{NNLS}(r,d)$ cannot be specified exactly since theoretical running times for exact NNLS solvers are unknown. In the sequel we comment on the computational costs of some well defined segments of some NNLS solvers.

The NNLS formulation of Definition 1 is a convex quadratic program, and is equivalent to

On the other hand, many standard implementations of NNLS solvers (and in particular those that are based on active set methods) work directly on the formulation of Definition 1. A typical cost of these implementations is of the order $O(nd^{2})$ per iteration. Other approaches, for example the NNLS method of , proceed by computing matrix-vector products of the form $Au$, for an appropriate $d$-dimensional vector $u$, thus cost typically $O(nd)$ time per iteration. In these cases our algorithm needs again $T_{precond}$ preprocessing time, but costs only $O(rd^{2})$ or $O(rd)$ time per iteration, respectively. Again, if given our choice of $r$, the computational savings per iteration are comparable with the $O(d/\log(d))$ speedup described above.

Experimental Evaluation

In this section, we experimentally evaluate our RandomizedNNLS algorithm on (i) large, sparse matrices from a text-mining application, and (ii) random matrices with varying sparsity. We mainly focus on employing the state-of-the-art solver of to solve the small NNLS problem.

Our data come from the Open Directory Project (ODP) , a multilingual open content directory of WWW links that is constructed and maintained by a community of volunteer editors. ODP uses a hierarchical ontology scheme for organizing site listings. Listings on similar topics are grouped into categories, which can then include smaller subcategories. Gabrilovich and Markovitch constructed a benchmark set of 300 term-document matrices from ODP, called TechTC300 (Technion Repository of Text Categorization Datasets ), which they made publicly available. Each term-document matrix of the TechTC300 dataset consists of a total of 150 to 400 documents from two different ODP categories, and a total of 15,000 to 35,000 terms. We chose this dataset because we believe that it does represent an important application area, namely text mining, and we do believe that the results from our experiments will be representative of the potential usefulness of our randomized NNLS algorithm in large, sparse, term-document NNLS problems.

We present average results from 3,000 NNLS problems. More specifically, for each of the 300 matrices of the TechTC300 dataset, we randomly choose a column from the term-document matrix as the vector $b$, we assign the remaining columns of the same term-document matrix to the matrix $A$, and solve the resulting NNLS problem with inputs $A$ and $b$. We repeat this process ten times for each term-document matrix of the TechTC300 dataset, and thus solve a total of 3,000 problems. Whenever an NNLS routine is called, it is initialized with the all-zeros vector. We evaluate the accuracy and the running time of our algorithm when compared to two standard NNLS algorithms. The first one is described in We would like to thank the authors of for providing us with a Matlab implementation of their algorithm., and the second one is the active set method of , implemented as the built-in function lsqnonneg in Matlab. We would also like to emphasize that in the authors compare their approach to other NNLS approaches and conclude that their algorithm is significantly faster. Note that the method of operates on the quadratic programming formulation discussed in Section 3.1The actual implementation involves computations of the form $t=Au$ and $s=A^{T}t$, avoiding the computation and storage of the matrix $A^{T}A$., while lsqnonneg operates on the formulation of Definition 1. Finally, we implemented our RandomizedNNLS algorithm in Matlab. The platform used for the experiments was a 2.0 GHz Pentium IV with 1GB RAM.

Our (average) results are shown in Figure 1. We only focus on the algorithm of , which was significantly faster, running (on average) in five seconds, compared to more than one minute for the lsqnonneg function. We experimented with eight different values of the parameter $r$, which dictates the size of the small subproblem (see the RandomizedNNLS algorithm). More specifically, we set $r$ to $d+i\cdot 50$, for $i=1\ldots 8$, where $d$ is the number of columns in the matrix $A$. Our results verify that ($i$) the RandomizedNNLS algorithm is very accurate, ($ii$) that it reduces the running time of the NNLS method of , and $(iii)$ that there exists a natural tradeoff between the approximation accuracy and the number of sampled rows. Notice, for example, that the running time of the state-of-the-art NNLS solver of can be reduced from two to three times, while the residual error is from $4\%$ up to $10\%$ worse than the optimal residual error.

We briefly comment on the performance of RandomizedNNLS when compared to the lsqnonneg algorithm. As expected, the accuracy results are essentially identical with the method of , since both methods solve the NNLS problem exactly. Our speedup, however, was much more significant, ranging from 14-fold to 10-fold for $r=d+50$ and $r=d+400$ respectively (data not shown).

2 Sparse vs dense NNLS problems

The astute reader might notice that we evaluated the performance of our algorithm in a rather adversarial setting. The TechTC300 data are quite sparse, hence existing NNLS methods would operate on sparse matrices. However, our preprocessing step in the RandomizedNNLS algorithm destroys the sparsity, and the induced subproblem becomes dense. Thus, we are essentially comparing the time required to solve a sparse, large NNLS problem to the time required to solve a dense, small NNLS problem. If the original problem were dense as well, we would expect more pronounced computational savings. In this section we experiment with random matrices of varying density in order to confirm this hypothesis.

First, it is worth noting that the sparsity of the input matrix $A$ and/or the target vector $b$ do not seem to affect the approximation accuracy of the RandomizedNNLS algorithm. This should not come as a surprise since our results in Theorem 1 do not make any assumptions on the inputs $A$ and $b$. Indeed, our experiments in Figure 2 confirm our expectations.

Prior to discussing our experiments on random matrices of varying density, it is worth noting that the NNLS solver of has a running time that is a function of the number of non-zero entries of $A$. Indeed, the method of is an iterative method where the computational bottleneck in the $j$-th iteration involves computations of the form $A^{T}Au$, for a $d$-dimensional vector $u$. implemented their algorithm by computing the two matrix-vector products $Au$ and $A^{T}\left(Au\right)$ separately, thus never forming the matrix $A^{T}A$ and thus taking advantage of the sparsity of $A$. Indeed, NNLS problems with sparse coefficient matrices $A$ are solved faster than NNLS problems with similar-size dense coefficient matrices $A$ by using the method of The authors of performed extensive numerical experiments to verify that observation; for example see the last row of Table 4 on page 14 in and notice that the running time of their method increases as the density of $A$ increases..

In order to measure how the speedup of our RandomizedNNLS algorithm improves as the matrix $A$ and vector $b$ become denser, we designed the following experiment. First, let the density of an NNLS problem denote the percentage of non-zero entries in $A$ and $b$; for example, $density(A,b)=10\%$ means that approximately $0.9(nd+n)$ entries in the $n\times d$ matrix $A$ and the $n\times 1$ vector $b$ are zero. We chose six density parameters ($2\%$, $4\%$, $8\%$, $16\%$, $32\%$, and $64\%$) and generated $100$ NNLS problems for each density parameter. More specifically, we first constructed ten $n\times(d+1)$ random matrices with the target density (the non-zero entries are normally distributed in $$). Then, for each matrix, we randomly selected one column to form the vector$b$and assigned the remaining$d$columns to the matrix$A$. We repeated this selection process ten times for each of the ten$n\times(d+1)$matrices, thus forming a set of$100$NNLS problems with inputs$A$and$b$. We fixed the dimensions to$n=10,000$and$d=300$and we experimented with four values of$r=(d,d+50,d+100,d+150)$. In Figure 2 we present average results over the$100$NNLS problems for each choice of the density parameter. Notice that increasing the density of the inputs$A$and$b$, the computational gains increase as well. On top of that, our method becomes more accurate while the number of the zero entries in$A$and$b$become fewer. Notice for example, on the right plot of Figure 2, when$r=300$, the two extreme cases ($density=2\%$and$density=64\%$) correspond to an$18\%$and a$4\%$loss in accuracy, respectively. Given these two observations as well as the actual times of$T_{x_{opt}}$ (see the caption of Figure 2), we conclude that the random projection ideas empirically seem more promising for dense rather than sparse NNLS problems.

Conclusions

We presented a random projection algorithm for the Nonnegative Least Squares Problem. We experimentally evaluated our algorithm on a large, text-mining dataset, and verified that, as promised in our theoretical findings, practically it does give very accurate approximate solutions, while outperforming two standard NNLS methods in terms of computational efficiency. Future work includes the extension of our theoretical findings of Theorem 1 to NNLS problems with multiple right hand side vectors. An immediate application of this would be the computation of Nonnegative Matrix Factorizations based on Alternating Least Squares type approaches . Finally, notice that, since our analysis is independent of the type of constraints on the vector $x$, our main algorithm can be employed to approximate a least-squares problem with any type of constraints on $x$.

We would like to thank Kristin P. Bennett and Michael W. Mahoney for useful discussions. The first author would like also thank the Institute of Pure and Applied Mathematics of the University of California at Los Angeles for its generous hospitality during the period Sept. 2008 - Dec. 2008, when part of this work was done.

References