LSRN: A Parallel Iterative Solver for Strongly Over- or Under-Determined Systems

Xiangrui Meng, Michael A. Saunders, Michael W. Mahoney

Introduction

Randomized algorithms have become indispensable in many areas of computer science, with applications ranging from complexity theory to combinatorial optimization, cryptography, and machine learning. Randomization has also been used in numerical linear algebra (for instance, the initial vector in the power iteration is chosen at random so that almost surely it has a nonzero component along the direction of the dominant eigenvector), yet most well-developed matrix algorithms, e.g., matrix factorizations and linear solvers, are deterministic. In recent years, however, motivated by large data problems, very nontrivial randomized algorithms for very large matrix problems have drawn considerable attention from researchers, originally in theoretical computer science and subsequently in numerical linear algebra and scientific computing. By randomized algorithms, we refer in particular to random sampling and random projection algorithms . For a comprehensive overview of these developments, see the review of Mahoney , and for an excellent overview of numerical aspects of coupling randomization with classical low-rank matrix factorization methods, see the review of Halko, Martinsson, and Tropp .

If we let $r=\mathop{\operator@font rank}\nolimits(A)\leq\min(m,n)$ , then recall that if $r<n$ (the LS problem is under-determined or rank-deficient), then (1) has an infinite number of minimizers. In that case, the set of all minimizers is convex and hence has a unique element having minimum length. On the other hand, if $r=n$ so the problem has full rank, there exists only one minimizer to (1) and hence it must have the minimum length. In either case, we denote this unique min-length solution to (1) by $x^{*}$ . That is,

LS problems of this form have a long history, tracing back to Gauss, and they arise in numerous applications. The demand for faster LS solvers will continue to grow in light of new data applications and as problem scales become larger and larger.

In this paper, we describe an LS solver called LSRN for these strongly over- or under-determined, and possibly rank-deficient, systems. LSRN uses random normal projections to compute a preconditioner matrix such that the preconditioned system is provably extremely well-conditioned. Importantly for large-scale applications, the preconditioning process is embarrassingly parallel, and it automatically speeds up with sparse matrices and fast linear operators. LSQR or the Chebyshev semi-iterative (CS) method can be used at the iterative step to compute the min-length solution within just a few iterations. We show that the latter method is preferred on clusters with high communication cost.

Because of its provably-good conditioning properties, LSRN has a fully predictable run-time performance, just like direct solvers, and it scales well in parallel environments. On large dense systems, LSRN is faster than LAPACK’s DGELSD for strongly over-determined problems, and is much faster for strongly under-determined problems, although solvers using fast random projections, like Blendenpik , are still slightly faster in both cases. On sparse systems, LSRN runs significantly faster than competing solvers, for both the strongly over- or under-determined cases.

In section 2 we describe existing deterministic LS solvers and recent randomized algorithms for the LS problem. In section 3 we show how to do preconditioning correctly for rank-deficient LS problems, and in section 4 we introduce LSRN and discuss its properties. Section 5 describes how LSRN can handle Tikhonov regularization for both over- and under-determined systems, and in section 6 we provide a detailed empirical evaluation illustrating the behavior of LSRN.

Least squares solvers

In this section we discuss related work, including deterministic direct and iterative methods as well as recently developed randomized methods, for computing solutions to LS problems, and we discuss how our results fit into this broader context.

A third way to solve (2) is by computing the min-length solution to the normal equation $A^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}Ax=A^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}b$ , namely

It is easy to verify the correctness of (3) by replacing $A$ by its economy-sized SVD $U\Sigma V^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}$ . If $r=\min(m,n)$ , a Cholesky factorization of either $A^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}A$ (if $m\geq n$ ) or $AA^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}$ (if $m\leq n$ ) solves (3) nicely. If $r<\min(m,n)$ , we need the eigensystem of $A^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}A$ or $AA^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}$ to compute $x^{*}$ . The normal equation approach is the least expensive among the three direct approaches we have mentioned, especially when $m\gg n$ or $m\ll n$ , but it is also the least accurate one, especially on ill-conditioned problems. See Chapter 5 of Golub and Van Loan for a detailed analysis.

Instead of these direct methods, we can use iterative methods to solve (1). If all the iterates $\{x^{(k)}\}$ are in $\text{range}(A^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}})$ and if $\{x^{(k)}\}$ converges to a minimizer, it must be the minimizer having minimum length, i.e., the solution to (2). This is the case when we use a Krylov subspace method starting with a zero vector. For example, the conjugate gradient (CG) method on the normal equation leads to the min-length solution (see Paige and Saunders ). In practice, CGLS , LSQR are preferable because they are equivalent to applying CG to the normal equation in exact arithmetic but they are numerically more stable. Other Krylov subspace methods such as the CS method and LSMR can solve (1) as well.

Importantly, however, it is in general hard to predict the number of iterations for CG-like methods. The convergence rate is affected by the condition number of $A^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}A$ . A classical result [16, p.187] states that

2 Randomized methods

In 2007, Drineas, Mahoney, Muthukrishnan, and Sarlós introduced two randomized algorithms for the LS problem, each of which computes a relative-error approximation to the min-length solution in $\mathcal{O}(mn\log n)$ time, when $m\gg n$ . Both of these algorithms apply a randomized Hadamard transform to the columns of $A$ , thereby generating a problem of smaller size, one using uniformly random sampling and the other using a sparse random projection. They proved that, in both cases, the solution to the smaller problem leads to relative-error approximations of the original problem. The accuracy of the approximate solution depends on the sample size; and to have relative precision $\varepsilon$ , one should sample $\mathcal{O}(n/\varepsilon)$ rows after the randomized Hadamard transform. This is suitable when low accuracy is acceptable, but the $\varepsilon$ dependence quickly becomes the bottleneck otherwise. Using those algorithms as preconditioners was also mentioned in . This work laid the ground for later algorithms and implementations.

Later, in 2008, Rokhlin and Tygert described a related randomized algorithm for over-determined systems. They used a randomized transform named SRFT that consists of $m$ random Givens rotations, a random diagonal scaling, a discrete Fourier transform, and a random sampling. They considered using their method as a preconditioning method, and they showed that to get relative precision $\varepsilon$ , only $\mathcal{O}(n\log(1/\varepsilon))$ samples are needed. In addition, they proved that if the sample size is greater than $4n^{2}$ , the condition number of the preconditioned system is bounded above by a constant. Although choosing this many samples would adversely affect the running time of their solver, they also illustrated examples of input matrices for which the $4n^{2}$ sample bound was weak and for which many fewer samples sufficed.

Then, in 2010, Avron, Maymounkov, and Toledo implemented a high-precision LS solver, called Blendenpik, and compared it to LAPACK’s DGELS and to LSQR with no preconditioning. Blendenpik uses a Walsh-Hadamard transform, a discrete cosine transform, or a discrete Hartley transform for blending the rows/columns, followed by a random sampling, to generate a problem of smaller size. The $R$ factor from the QR factorization of the smaller matrix is used as the preconditioner for LSQR. Based on their analysis, the condition number of the preconditioned system depends on the coherence or statistical leverage scores of $A$ , i.e., the maximal row norm of $U$ , where $U$ is an orthonormal basis of $\text{range}(A)$ . We note that a solver for under-determined problems is also included in the Blendenpik package.

In 2011, Coakley, Rokhlin, and Tygert described an algorithm that is also based on random normal projections. It computes the orthogonal projection of any vector $b$ onto the null space of $A$ or onto the row space of $A$ via a preconditioned normal equation. The algorithm solves the over-determined LS problem as an intermediate step. They show that the normal equation is well-conditioned and hence the solution is reliable. For an over-determined problem of size $m\times n$ , the algorithm requires applying $A$ or $A^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}$ $3n+6$ times, while LSRN needs approximately $2n+200$ matrix-vector multiplications under the default setting. Asymptotically, LSRN will become faster as $n$ increases beyond several hundred. See section 4.3 for further complexity analysis of LSRN.

3 Relationship with our contributions

All prior approaches assume that $A$ has full rank, and for those based on iterative solvers, none provides a tight upper bound on the condition number of the preconditioned system (and hence the number of iterations). For LSRN, Theorem 2 ensures that the min-length solution is preserved, independent of the rank, and Theorems 6 and 7 provide bounds on the condition number and number of iterations, independent of the spectrum of $A$ . In addition to handling rank-deficiency well, LSRN can even take advantage of it, resulting in a smaller condition number and fewer iterations.

Some prior work on the LS problem has explored “fast” randomized transforms that run in roughly $\mathcal{O}(mn\log m)$ time on a dense matrix $A$ , while the random normal projection we use in LSRN takes $\mathcal{O}(mn^{2})$ time. Although this could be an issue for some applications, the use of random normal projections comes with several advantages. First, if $A$ is a sparse matrix or a linear operator, which is common in large-scale applications, then the Hadamard-based fast transforms are no longer “fast”. Second, the random normal projection is easy to implement using threads or MPI, and it scales well in parallel environments. Third, the strong symmetry of the standard normal distribution helps give the strong high probability bounds on the condition number in terms of sample size. These bounds depend on nothing but $s/r$ , where $s$ is the sample size. For example, if $s=4r$ , Theorem 6 ensures that, with high probability, the condition number of the preconditioned system is less than $3$ .

This last property about the condition number of the preconditioned system makes the number of iterations and thus the running time of LSRN fully predictable like for a direct method. It also enables use of the CS method, which needs only one level-1 and two level-2 BLAS operations per iteration, and is particularly suitable for clusters with high communication cost because it doesn’t have vector inner products that require synchronization between nodes. Although the CS method has the same theoretical upper bound on the convergence rate as CG-like methods, it requires accurate bounds on the singular values in order to work efficiently. Such bounds are generally hard to come by, limiting the popularity of the CS method in practice, but they are provided for the preconditioned system by our Theorem 6, and we do achieve high efficiency in our experiments.

Preconditioning for linear least squares

In light of (4), much effort has been made to transform a linear system into an equivalent system with reduced condition number. This preconditioning, for a square linear system $Bx=d$ of full rank, usually takes one of the following forms:

Clearly, the preconditioned system is consistent with the original one, i.e., has the same $x^{*}$ as the unique solution, if the preconditioners $M$ and $N$ are nonsingular.

For the general LS problem (2), preconditioning needs better handling in order to produce the same min-length solution as the original problem. For example, if we apply left preconditioning to the LS problem $\min_{x}\|Ax-b\|_{2}$ , the preconditioned system becomes $\min_{x}\|M^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}Ax-M^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}b\|_{2}$ , and its min-length solution is given by

Similarly, the min-length solution to the right preconditioned system is given by

The following lemma states the necessary and sufficient conditions for $A^{\dagger}=N(AN)^{\dagger}$ or $A^{\dagger}=(M^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}A)^{\dagger}M^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}$ to hold. Note that these conditions holding certainly imply that $x_{\text{right}}^{*}=x^{*}$ and $x_{\text{left}}^{*}=x^{*}$ , respectively.

$A^{\dagger}=N(AN)^{\dagger}$ if and only if $\text{range}(NN^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}A^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}})=\text{range}(A^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}})$ ,

$A^{\dagger}=(M^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}A)^{\dagger}M^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}$ if and only if $\text{range}(MM^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}A)=\text{range}(A)$ .

Let $r=\mathop{\operator@font rank}\nolimits(A)$ and $U\Sigma V^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}$ be $A$ ’s economy-sized SVD as in section 2.1, with $A^{\dagger}=V\Sigma^{-1}U^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}$ . Before continuing our proof, we reference the following facts about the pseudoinverse:

$B^{\dagger}=B^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}(BB^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}})^{\dagger}$ for any matrix $B$ ,

For any matrices $B$ and $C$ such that $BC$ is defined, $(BC)^{\dagger}=C^{\dagger}B^{\dagger}$ if (i) $B^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}B=I$ or (ii) $CC^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}=I$ or (iii) $B$ has full column rank and $C$ has full row rank.

Now let’s prove the “if” part of the first statement. If $\text{range}(NN^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}A^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}})=\text{range}(A^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}})=\text{range}(V)$ , we can write $NN^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}A^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}$ as $VZ$ where $Z$ has full row rank. Then,

Conversely, if $N(AN)^{\dagger}=A^{\dagger}$ , we know that $\text{range}(N(AN)^{\dagger})=\text{range}(A^{\dagger})=\text{range}(V)$ and hence $\text{range}(V)\subseteq\text{range}(N)$ . Then we can decompose $N$ as $(V\ V_{c})\hbox{\scriptsize$ \begin{pmatrix}Z\\ Z_{c}\end{pmatrix} $}=VZ+V_{c}Z_{c}$ , where $V_{c}$ is orthonormal, $V^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}V_{c}=0$ , and $\smash[tb]{\begin{pmatrix}Z\\ Z_{c}\end{pmatrix}}$ has full row rank. Then,

Multiplying by $V_{c}^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}$ on the left and $U\Sigma$ on the right, we get $Z_{c}Z^{\dagger}=0$ , which is equivalent to $Z_{c}Z^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}=0$ . Therefore,

where we used the facts that $Z$ has full row rank and hence $ZZ^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}$ is nonsingular, $\Sigma$ is nonsingular, and $U$ has full column rank.

To prove the second statement, let us take $B=A^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}$ . By the first statement, we know $B^{\dagger}=M(BM)^{\dagger}$ if and only if $\text{range}(MM^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}B^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}})=\text{range}(B^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}})$ , which is equivalent to saying $A^{\dagger}=(M^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}A)^{\dagger}M^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}$ if and only if $\text{range}(MM^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}A)=\text{range}(A)$ . ∎

Although Lemma 1 gives the necessary and sufficient condition, it does not serve as a practical guide for preconditioning LS problems. In this work, we are more interested in a sufficient condition that can help us build preconditioners. To that end, we provide the following theorem.

$x_{\text{right}}^{*}=x^{*}$ if $\text{range}(N)=\text{range}(A^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}})$ ,

$x_{\text{left}}^{*}=x^{*}$ if $\text{range}(M)=\text{range}(A)$ .

Let $r=\mathop{\operator@font rank}\nolimits(A)$ and $U\Sigma V^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}$ be $A$ ’s economy-sized SVD. If $\text{range}(N)=\text{range}(A^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}})=\text{range}(V)$ , we can write $N$ as $VZ$ , where $Z$ has full row rank. Therefore,

By Lemma 1, $A^{\dagger}=N(AN)^{\dagger}$ and hence $x_{\text{left}}^{*}=x^{*}$ . The second statement can be proved by similar arguments. ∎

Algorithm LSRN

In this section we present LSRN, an iterative solver for solving strongly over- or under-determined systems, based on “random normal projection”. To construct a preconditioner we apply a transformation matrix whose entries are independent random variables drawn from the standard normal distribution. We prove that the preconditioned system is almost surely consistent with the original system, i.e., both have the same min-length solution. At least as importantly, we prove that the spectrum of the preconditioned system is independent of the spectrum of the original system; and we provide a strong concentration result on the extreme singular values of the preconditioned system. This concentration result enables us to predict the number of iterations for CG-like methods, and it also enables use of the CS method, which requires an accurate bound on the singular values to work efficiently.

Algorithm 1 shows the detailed procedure of LSRN to compute the min-length solution to a strongly over-determined problem, and Algorithm 2 shows the detailed procedure for a strongly under-determined problem. We refer to these two algorithms together as LSRN. Note that they only use the input matrix $A$ for matrix-vector and matrix-matrix multiplications, and thus $A$ can be a dense matrix, a sparse matrix, or a linear operator. In the remainder of this section we focus on analysis of the over-determined case. We emphasize that analysis of the under-determined case is quite analogous.

2 Theoretical properties

The use of random normal projection offers LSRN some nice theoretical properties. We start with consistency.

In Algorithm 1, we have $\hat{x}=A^{\dagger}b$ almost surely.

Let $r=\mathop{\operator@font rank}\nolimits(A)$ and $U\Sigma V^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}$ be $A$ ’s economy-sized SVD. We have

and hence by Theorem 2 we have $\hat{x}=A^{\dagger}b$ almost surely. ∎

A more interesting property of LSRN is that the spectrum (the set of singular values) of the preconditioned system is solely associated with a random matrix of size $s\times r$ , independent of the spectrum of the original system.

In Algorithm 1, the spectrum of $AN$ is the same as the spectrum of $G_{1}^{\dagger}=(GU)^{\dagger}$ , independent of $A$ ’s spectrum.

which gives $AN$ ’s SVD. Therefore, $AN$ ’s singular values are $\text{diag}(\Sigma_{1}^{-1})$ , the same as $G_{1}^{\dagger}$ ’s spectrum, but independent of $A$ ’s. ∎

We know that $G_{1}=GU$ is a random matrix whose entries are independent random variables following the standard normal distribution. The spectrum of $G_{1}$ is a well-studied problem in Random Matrix Theory, and in particular the properties of extreme singular values have been studied. Thus, the following lemma is important for us. We use $\mathcal{P}(\cdot)$ to refer to the probability that a given event occurs.

(Davidson and Szarek ) Consider an $s\times r$ random matrix $G_{1}$ with $s\geq r$ , whose entries are independent random variables following the standard normal distribution. Let the singular values be $\sigma_{1}\geq\cdots\geq\sigma_{r}$ . Then for any $t>0$ ,

With the aid of Lemma 5, it is straightforward to obtain the concentration result of $\sigma_{1}(AN)$ , $\sigma_{r}(AN)$ , and $\kappa(AN)$ as follows.

In Algorithm 1, for any $\alpha\in(0,1-\sqrt{r/s})$ , we have

In order to estimate the number of iterations for CG-like methods, we can now combine (4) and (7).

In exact arithmetic, given a tolerance $\varepsilon>0$ , a CG-like method applied to the preconditioned system $\min_{y}\|ANy-b\|_{2}$ with $y^{(0)}=0$ converges within $(\log\varepsilon-\log 2)/\log(\alpha+\sqrt{r/s})$ iterations in the sense that

holds with probability at least $1-2e^{-\alpha^{2}s/2}$ for any $\alpha\in(0,1-\sqrt{s/r})$ , where $\hat{y}_{\textrm{\tiny CG}}$ is the approximate solution returned by the CG-like solver and $y^{*}=(AN)^{\dagger}b$ . Let $\hat{x}_{\textrm{\tiny CG}}=N\hat{y}_{\textrm{\tiny CG}}$ be the approximate solution to the original problem. Since $x^{*}=Ny^{*}$ , (8) is equivalent to

where $\hat{r}_{\textrm{\tiny CG}}=b-A\hat{x}_{\textrm{\tiny CG}}$ and $r^{*}=b-Ax^{*}$ .

In addition to allowing us to bound the number of iterations for CG-like methods, the result given by (6) also allows us to use the CS method. This method needs only one level-1 and two level-2 BLAS operations per iteration; and, importantly, because it doesn’t have vector inner products that require synchronization between nodes, this method is suitable for clusters with high communication cost. It does need an explicit bound on the singular values, but once that bound is tight, the CS method has the same theoretical upper bound on the convergence rate as other CG-like methods. Unfortunately, in many cases, it is hard to obtain such an accurate bound, which prevents the CS method becoming popular in practice. In our case, however, (6) provides a probabilistic bound with very high confidence. Hence, we can employ the CS method without difficulty. For completeness, Algorithm 3 describes the CS method we implemented for solving LS problems. For discussion of its variations, see Gutknecht and Rollin .

3 Running time complexity

where lower-order terms are ignored. Here, $\text{flops}(\text{randn})$ is the average flop count to generate a sample from the standard normal distribution, while $\text{flops}(Av)$ and $\text{flops}(A^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}u)$ are the flop counts for the respective matrix-vector products. If $A$ is a dense matrix, then we have $\text{flops}(Av)=\text{flops}(A^{\raisebox{0.62997pt}{\hbox{\tiny$ T $}}}u)=2mn$ . Hence, the total cost becomes

where $T_{\text{randn}}$ , $T_{\text{mult}}$ , and $T_{\text{iter}}$ are the running times for the respective stages if LSRN runs on a single core, $T_{\text{svd}}^{\text{mt},p}$ is the running time of SVD using $p$ cores, and communication cost among threads is ignored. Hence, multi-threaded LSRN has very good scalability with near-linear speedup.

Tikhonov regularization

We point out that it is easy to extend LSRN to handle certain types of Tikhonov regularization, also known as ridge regression. Recall that Tikhonov regularization involves solving the problem

This is an over-determined problem of size $(m+n)\times n$ . If $m\gg n$ , then we certainly have $m+n\gg n$ . Therefore, if $m\gg n$ , we can directly apply LSRN to (14) in order to solve (13). On the other hand, if $m\ll n$ , then although (14) is still over-determined, it is “nearly square,” in the sense that $m+n$ is only slightly larger than $n$ . In this regime, random sampling methods and random projection methods like LSRN do not perform well. In order to deal with this regime, note that (13) is equivalent to

where $r=b-Ax$ is the residual vector. (Note that we use $r$ to denote the matrix rank in a scalar context and the residual vector in a vector context.) By introducing $z=Wx$ and assuming that $W$ is non-singular, we can re-write the above problem as

i.e., as computing the min-length solution to

Note that (15) is an under-determined problem of size $m\times(m+n)$ . Hence, if $m\ll n$ , we have $m\ll m+n$ and we can use LSRN to compute the min-length solution to (15), denoted by $\begin{pmatrix}z^{*}\\ r^{*}\end{pmatrix}$ . The solution to the original problem (13) is then given by $x^{*}=W^{-1}z^{*}$ . Here, we assume that $W^{-1}$ is easy to apply, as is the case when $W=\lambda I_{n}$ , so that $AW^{-1}$ can be treated as an operator. The equivalence between (13) and (15) was first established by Herman, Lent, and Hurwitz .

In most applications of regression analysis, the amount of regularization, e.g., the optimal regularization parameter, is unknown and thus determined by cross-validation. This requires solving a sequence of LS problems where only $W$ differs. For over-determined problems, we only need to perform a random normal projection on $A$ once. The marginal cost to solve for each $W$ is the following: a random normal projection on $W$ , an SVD of size $\lceil\gamma n\rceil\times n$ , and a predictable number of iterations. Similar results hold for under-determined problems when each $W$ is a multiple of the identity matrix.

Numerical experiments

We implemented our LS solver LSRN and compared it with competing solvers: LAPACK’s DGELSD, MATLAB’s backslash, and Blendenpik by Avron, Maymounkov, and Toledo . MATLAB’s backslash uses different algorithms for different problem types. For sparse rectangular systems, as stated by Tim Davishttp://www.cise.ufl.edu/research/sparse/SPQR/, “SuiteSparseQR is now QR in MATLAB 7.9 and $x=A\backslash b$ when $A$ is sparse and rectangular.” Table 1 summarizes the properties of those solvers. We report our empirical results in this section.

The experiments were performed on either a local shared-memory machine or a virtual cluster hosted on Amazon’s Elastic Compute Cloud (EC2). The shared-memory machine has $12$ Intel Xeon CPU cores at clock rate 2GHz with 128GB RAM. The virtual cluster consists of $20$ m1.large instances configured by a third-party tool called StarClusterhttp://web.mit.edu/stardev/cluster/. An m1.large instance has $2$ virtual cores with $2$ EC2 Compute Units“One EC2 Compute Unit provides the equivalent CPU capacity of a 1.0-1.2 GHz 2007 Opteron or 2007 Xeon processor.” from http://aws.amazon.com/ec2/faqs/ each. To attain top performance on the shared-memory machine, we implemented a multi-threaded version of LSRN in C, and to make our solver general enough to handle large problems on clusters, we also implemented an MPI version of LSRN in Python with NumPy, SciPy, and mpi4py. Both packages are available for downloadhttp://www.stanford.edu/group/SOL/software/lsrn.html. We use the multi-threaded implementation to compare LSRN with other LS solvers and use the MPI implementation to explore scalability and to compare iterative solvers under a cluster environment. To generate values from the standard normal distribution, we adopted the code from Marsaglia and Tsang and modified it to use threads; this can generate a billion samples in less than two seconds on the shared-memory machine. We also modified Blendenpik to call multi-threaded FFTW routines. Blendenpik’s default settings were used, i.e., using randomized discrete Fourier transform and sampling $4\min(m,n)$ rows/columns. All LAPACK’s LS solvers, Blendenpik, and LSRN are linked against MATLAB’s own multi-threaded BLAS and LAPACK libraries. So, in general, this is a fair setup because all the solvers can use multi-threading automatically and are linked against the same BLAS and LAPACK libraries. The running times were measured in wall-clock times.

2 κ(AN)𝜅𝐴𝑁\kappa(AN) and number of iterations

Recall that Theorem 6 states that $\kappa(AN)$ , the condition number of the preconditioned system, is roughly bounded by $(1+\sqrt{r/s})/(1-\sqrt{r/s})$ when $s$ is large enough such that we can ignore $\alpha$ in practice. To verify this statement, we generate random matrices of size $10^{4}\times 10^{3}$ with condition numbers ranged from $10^{2}$ to $10^{8}$ . The left figure in Figure 1 compares $\kappa(AN)$ with $\kappa_{+}(A)$ , the effective condition number of $A$ , under different choices of $s$ and $r$ . We take the largest value of $\kappa(AN)$ in $10$ independent runs as the $\kappa(AN)$ in the plot. For each pair of $s$ and $r$ , the corresponding estimate $(1+\sqrt{r/s})/(1-\sqrt{r/s})$ is drawn in a dotted line of the same color, if not overlapped with the solid line of $\kappa(AN)$ . We see that $(1+\sqrt{r/s})/(1-\sqrt{r/s})$ is indeed an accurate estimate of the upper bound on $\kappa(AN)$ . Moreover, $\kappa(AN)$ is not only independent of $\kappa_{+}(A)$ , but it is also quite small. For example, we have $(1+\sqrt{r/s})/(1-\sqrt{r/s})<6$ if $s>2r$ , and hence we can expect super fast convergence of CG-like methods.

Based on Theorem 7, the number of iterations should be less than $(\log\varepsilon-\log 2)/\log\sqrt{r/s}$ , where $\varepsilon$ is a given tolerance. In order to match the accuracy of direct solvers, we set $\varepsilon=10^{-14}$ . The right figure in Figure 1 shows the number of LSQR iterations for different combinations of $r/s$ and $\kappa_{+}(A)$ . Again, we take the largest iteration number in $10$ independent runs for each pair of $r/s$ and $\kappa_{+}(A)$ . We also draw the theoretical upper bound $(\log\varepsilon-\log 2)/\log\sqrt{r/s}$ in a dotted line. We see that the number of iterations is basically a function of $r/s$ , independent of $\kappa_{+}(A)$ , and the theoretical upper bound is very good in practice. This confirms that the number of iterations is fully predictable given $\gamma$ .

3 Tuning the oversampling factor γ𝛾\gamma

We experimented with various LS problems. The best choice of $\gamma$ ranges from $1.6$ to $2.5$ , depending on the type and the size of the problem. We also note that, when $\gamma$ is given, the running time of the iteration stage is fully predictable. Thus we can initialize LSRN by measuring randn/sec and flops/sec for matrix-vector multiplication, matrix-matrix multiplication, and SVD, and then determine the best value of $\gamma$ by minimizing the total running time (12). For simplicity, we set $\gamma=2.0$ in all later experiments; although this is not the optimal setting for all cases, it is always a reasonable choice.

4 Dense least squares

As the state-of-the-art dense linear algebra library, LAPACK provides several routines for solving LS problems, e.g., DGELS, DGELSY, and DGELSD. DGELS uses QR factorization without pivoting, which cannot handle rank-deficient problems. DGELSY uses QR factorization with pivoting, which is more reliable than DGELS on rank-deficient problems. DGELSD uses SVD. It is the most reliable routine, and should be the most expensive as well. However, we find that DGELSD actually runs much faster than DGELSY on strongly over- or under-determined systems on the shared-memory machine. It may be because of better use of multi-threaded BLAS, but we don’t have a definitive explanation.

Figure 3 compares the running times of LSRN and competing solvers on randomly generated full-rank dense strongly over- or under-determined problems. We set the condition numbers to $10^{6}$ for all problems. Note that DGELS and DGELSD almost overlapped. The results show that Blendenpik is the winner. For small-sized problems ( $m\leq 3e4$ ), the follow-ups are DGELS and DGELSD. When the problem size goes larger, LSRN becomes faster than DGELS/DGELSD. DGELSY is always slower than DGELS/DGELSD, but still faster than MATLAB’s backslash. The performance of LAPACK’s solvers decreases significantly for under-determined problems. We monitored CPU usage and found that they couldn’t fully use all the CPU cores, i.e., they couldn’t effectively call multi-threaded BLAS. Though still the best, the performance of Blendenpik also decreases. LSRN’s performance does not change much.

LSRN is also capable of solving rank-deficient problems, and in fact it takes advantage of any rank-deficiency (in that it finds a solution in fewer iterations). Figure 4 shows the results on over- and under-determined rank-deficient problems generated the same way as in previous experiments, except that we set $r=800$ . DGELSY and DGELSD remain the same speed on over-determined problems as in full-rank cases, respectively, and run slightly faster on under-determined problems. LSRN’s running times reduce to $93$ seconds on the problem of size $10^{6}\times 10^{3}$ , from $100$ seconds on its full-rank counterpart.

We see that, for strongly over- or under-determined problems, DGELSD is the fastest and most reliable routine among the LS solvers provided by LAPACK. However, it (or any other LAPACK solver) runs much slower on under-determined problems than on over-determined problems, while LSRN works symmetrically on both cases. Blendenpik is the fastest dense least squares solver in our tests. Though it is not designed for solving rank-deficient problems, Blendenpik should be modifiable to handle such problems following Theorem 2. We also note that Blendenpik’s performance depends on the distribution of the row norms of $U$ . We generate test problems randomly so that the row norms of $U$ are homogeneous, which is ideal for Blendenpik. When the row norms of $U$ are heterogeneous, Blendenpik’s performance may drop. See Avron, Maymounkov, and Toledo for a more detailed analysis.

5 Sparse least squares

In LSRN, $A$ is only involved in the computation of matrix-vector and matrix-matrix multiplications. Therefore LSRN accelerates automatically when $A$ is sparse, without exploring $A$ ’s sparsity pattern. LAPACK does not have any direct sparse LS solver. MATLAB’s backslash uses SuiteSparseQR by Tim Davis when $A$ is sparse and rectangular; this requires explicit knowledge of $A$ ’s sparsity pattern to obtain a sparse QR factorization.

We generated sparse LS problems using MATLAB’s “sprandn” function with density $0.01$ and condition number $10^{6}$ . All problems have full rank. Figure 5 shows the results on over-determined problems. LAPACK’s solvers and Blendenpik basically perform the same as in the dense case. DGELSY is the slowest among the three. DGELS and DGELSD still overlap with each other, faster than DGELSY but slower than Blendenpik. We see that MATLAB’s backslash handles sparse problems very well. On the $10^{6}\times 10^{3}$ problem, backslash’s running time reduces to $55$ seconds, from $273$ seconds on the dense counterpart. The overall performance of MATLAB’s backslash is better than Blendenpik’s. LSRN’s curve is very flat. For small problems ( $m\leq 10^{5}$ ), LSRN is slow. When $m>10^{5}$ , LSRN becomes the fastest solver among the six. LSRN takes only $23$ seconds on the over-determined problem of size $10^{6}\times 10^{3}$ . On large under-determined problems, LSRN still leads by a huge margin.

LSRN makes no distinction between dense and sparse problems. The speedup on sparse problems is due to faster matrix-vector and matrix-matrix multiplications. Hence, although no test was performed, we expect a similar speedup on fast linear operators as well. Also note that, in the multi-threaded implementation of LSRN, we use a naive multi-threaded routine for sparse matrix-vector and matrix-matrix multiplications, which is far from optimized and thus leaves room for improvement.

6 Real-world problems

In this section, we report results on some real-world large data problems. The problems are summarized in Table 2, along with running times.

landmark and rail4284 are from the University of Florida Sparse Matrix Collection . landmark originated from a rank-deficient LS problem. rail4284 has full rank and originated from a linear programming problem on Italian railways. Both matrices are very sparse and have structured patterns. MATLAB’s backslash (SuiteSparseQR) runs extremely fast on these two problems, though it doesn’t guarantee to return the min-length solution. Blendenpik is not designed to handle the rank-deficient landmark, and it unfortunately runs out of memory (OOM) on rail4284. LSRN takes 17.55 seconds on landmark and 136.0 seconds on rail4284. DGELSD is slightly slower than LSRN on landmark and much slower on rail4284.

We see that, though both methods taking advantage of sparsity, MATLAB’s backslash relies heavily on the sparsity pattern, and its performance is unpredictable until the sparsity pattern is analyzed, while LSRN doesn’t rely on the sparsity pattern and always delivers predictable performance and, moreover, the min-length solution.

7 Scalability and choice of iterative solvers on clusters

Ideally, from the complexity analysis (12), when we double $n$ and double the number of cores, the increase in running time should be a constant if the cluster is homogeneous and has perfect load balancing (which we have observed is not true on Amazon EC2). For LSRN with CS, from tnimg_10 to tnimg_20 the running time increases $27.6$ seconds, and from tnimg_20 to tnimg_40 the running time increases $34.7$ seconds. We believe the difference between the time increases is caused by the heterogeneity of the cluster, because Amazon EC2 doesn’t guarantee the connection speed among nodes. From tnimg_4 to tnimg_40, the problem scale is enlarged by a factor of $10$ while the running time only increases by a factor of $50\%$ . The result still demonstrates LSRN’s good scalability. We also compare the performance of LSQR and CS as the iterative solvers in LSRN. For all problems LSQR converges in $84$ iterations and CS converges in $106$ iterations. However, LSQR is slower than CS. The communication cost saved by CS is significant on those tests. As a result, we recommend CS as the default LSRN iterative solver for cluster environments. Note that to reduce the communication cost on a cluster, we could also consider increasing $\gamma$ to reduce the number of iterations.

Conclusion

We developed LSRN, a parallel solver for strongly over- or under-determined, and possibly rank-deficient, systems. LSRN uses random normal projection to compute a preconditioner matrix for an iterative solver such as LSQR and the Chebyshev semi-iterative (CS) method. The preconditioning process is embarrassingly parallel and automatically speeds up on sparse matrices and fast linear operators, and on rank-deficient data. We proved that the preconditioned system is consistent and extremely well-conditioned, and derived strong bounds on the number of iterations of LSQR or the CS method, and hence on the total running time. On large dense systems, LSRN is competitive with the best existing solvers, and it runs significantly faster than competing solvers on strongly over- or under-determined sparse systems. LSRN is easy to implement using threads or MPI, and it scales well in parallel environments.

Acknowledgements

After completing the initial version of this manuscript, we learned of the LS algorithm of Coakley et al. . We thank Mark Tygert for pointing us to this reference. We are also grateful to Lisandro Dalcin, the author of mpi4py, for his own version of the MPI_Barrier function to prevent idle processes from interrupting the multi-threaded SVD process too frequently.