Randomized Extended Kaczmarz for Solving Least-Squares

Introduction

The Kaczmarz method is an iterative projection algorithm for solving linear systems of equations [Kac37]. Due to its simplicity, the Kaczmarz method has found numerous applications including image reconstruction, distributed computation and signal processing to name a few [FCM+92, Her80, Nat01, FZ12], see [Cen81] for more applications. The Kaczmarz method has also been rediscovered in the field of image reconstruction and called ART (Algebraic Reconstruction Technique) [GBH70], see also [CZ97, Her80] for additional references. It has been also applied to more general settings, see [Cen81, Table 1] and [Tom55, McC75] for non-linear versions of the Kaczmarz method.

Kaczmarz proved that this process converges to the unique solution for square non-singular matrices [Kac37], but without any attempt to bound the rate of convergence. Bounds on the rate of convergence of the Kaczmarz method are given in [McC75], [Ans84] and [Gal03, Theorem 4.4, p.120]. In addition, an error analysis of the Kaczmarz method under the finite precision model of computation is given in [Kni93, Kni96].

Nevertheless, the Kaczmarz method converges even if the linear system ${\mathsf{A}}\mathbf{x}={\mathbf{b}}$ is overdetermined ($m>n$) and has no solution. In this case and provided that ${\mathsf{A}}$ has full column rank, the Kaczmarz method converges to the least squares estimate. This was first observed by Whitney and Meany [WM67] who proved that the relaxed Kaczmarz method converges provided that the relaxation parameters are within $$and$\lambda_{k}\to 0$, see also [CEG83, Theorem 1], [Tan71] and [HN90] for additional references.

In the literature there was empirical evidence that selecting the rows non-uniformly at random may be more effective than selecting the rows via Kaczmarz’s cyclic manner [HM93, FCM+92]. Towards explaining such an empirical evidence, Strohmer and Vershynin proposed a simple randomized variant of the Kaczmarz method that has exponential convergence in expectation [SV09] assuming that the linear system is solvable; see also [LL10] for extensions to linear constraints. A randomized iterative algorithm that computes a sequence of random vectors $\mathbf{x}^{(0)},\mathbf{x}^{(1)},\ldots$ is said to converge in expectation to a vector $\mathbf{x}^{*}$ if and only if $\EE\left\|\mathbf{x}^{(k)}-\mathbf{x}^{*}\right\|_{2}^{2}\to 0$ as $k\to\infty$, where the expectation is taken over the random choices of the algorithm. Soon after [SV09], Needell analyzed the behavior of the randomized Kaczmarz method for the case of full column rank linear systems that do not have any solution [Nee10]. Namely, Needell proved that the randomized Kaczmarz estimate vector is (in the limit) within a fixed distance from the least squares solution and also that this distance is proportional to the distance of ${\mathbf{b}}$ from the column space of ${\mathsf{A}}$. In other words, Needell proved that the randomized Kaczmarz method is effective for least squares problems whose least squares error is negligible.

In this paper we present a randomized iterative least squares solver (Algorithm 3) that converges in expectation to the minimum Euclidean norm solution of

The proposed algorithm is based on [SV09, Nee10] and inspired by [Pop99]. More precisely the proposed algorithm can be thought of as a randomized variant of Popa’s extended Kaczmarz method [Pop99], therefore we named it as randomized extended Kaczmarz.

In Section 2, we briefly discuss related work on the design of deterministic and randomized algorithms for solving least squares problems. In Section 3, we present a randomized iterative algorithm for projecting a vector onto a subspace (represented as the column space of a given matrix) which may be of independent interest. In addition, we discuss the convergence properties of the randomized Kaczmarz algorithm for solvable systems (Section 3.2) and recall its analysis for non-solvable systems (Section 3.3). In Section 4, we present and analyze the randomized extended Kaczmarz algorithm. Finally, in Section 5 we provide a numerical evaluation of the proposed algorithm.

Least squares solvers

In this section we give a brief discussion on least squares solvers including deterministic direct and iterative algorithms together with recently proposed randomized algorithms. For a detailed discussion on deterministic methods, the reader is referred to [Bj96]. In addition, we place our contribution in context with prior work.

In the literature, several methods have been proposed for solving least squares problems of the form (1). Here we briefly describe a representative sample of such methods including the use of QR factorization with pivoting, the use of the singular value decomposition (SVD) and iterative methods such as Krylov subspace methods applied on the normal equations [Saa03]. LAPACK provides robust implementations of the first two methods; DGELSY uses QR factorization with pivoting and DGELSD uses the singular value decomposition [ABD+90]. For the iterative methods, LSQR is equivalent to applying the conjugate gradient method on the normal equations [PS82] and it is a robust and numerically stable method.

Randomized algorithms

To the best of our knowledge, most randomized algorithms proposed in the theoretical computer science literature for approximately solving least squares are mainly based on the following generic two step procedure: first randomly (and efficiently) project the linear system into sufficiently many dimensions, and second return the solution of the down-sampled linear system as an approximation to the original optimal solution [DMM06, Sar06, CW09, NDT09, MZ11, DMMS11], see also [CW12]. Concentration of measure arguments imply that the optimal solution of the down-sampled system is close to the optimal solution of the original system. The accuracy of the approximate solution using this approach depends on the sample size and to achieve relative accuracy $\varepsilon$, the sample size should depend inverse polynomially on $\varepsilon$. This makes these approaches unsuitable for the high-precision regime of error that is considered here.

A different approach is the so called randomized preconditioning method, see [RT08, AMT10]. The authors of [AMT10] implemented Blendenpik, a high-precision least squares solver. Blendenpik consists of two steps. In the first step, the input matrix is randomly projected and an effective preconditioning matrix is extracted from the projected matrix. In the second step, an iterative least squares solver such as the LSQR algorithm of Paige and Saunders [PS82] is applied on the preconditioned system. Blendenpik is effective for overdetermined and underdetermined problems.

A parallel iterative least squares solver based on normal random projections called LSRN was recently implemented by Meng, Saunders and Mahoney [MSM11]. LSRN consists of two phases. In the first preconditioning phase, the original system is projected using random normal projection from which a preconditioner is extracted. In the second step, an iterative method such as LSQR or the Chebyshev semi-iterative method [GV61] is applied on the preconditioned system. This approach is also effective for over-determined and under-determined least squares problems assuming the existence of a parallel computational environment.

1 Relation with our contribution

In Section 5, we compare the randomized extended Kaczmarz algorithm against DGELSY, DGELSD, Blendenpik. LSRN [MSM11] did not perform well under a setup in which no parallelization is allowed, so we do not include LSRN’s performance. The numerical evaluation of Section 5 indicates that the randomized extended Kaczmarz is effective on the case of sparse, well-conditioned and strongly rectangular (both overdetermined and underdetermined) least squares problems, see Figure 1. Moreover, the randomized extended Kaczmarz algorithm has also comparable performance with LAPACK’s routine for the dense random input matrices, see Figure 2 (notice that the proposed algorithm almost matches Blendenpik’s performance for the underdetermined case, see Figure 2(b)). On the other hand, a preconditioned version of the proposed algorithm does not perform well under the case of ill-conditioned matrices, see Figure 3.

Background

where $p_{j}:=\left\|{\mathsf{A}}_{(j)}\right\|_{2}^{2}/\left\|{\mathsf{A}}\right\|_{\text{\rm F}}^{2}$ for every $i\in{[n]}$ and $q_{i}:=\left\|{\mathsf{A}}^{(i)}\right\|_{2}^{2}/\left\|{\mathsf{A}}\right\|_{\text{\rm F}}^{2}$ for every $i\in{[m]}$. The following fact will be used extensively in the paper.

1 Randomized Approximate Orthogonal Projection

Algorithm 1 is iterative. Initially, it starts with $\mathbf{z}^{(0)}={\mathbf{b}}$. At the $k$-th iteration, the algorithm randomly selects a column ${\mathsf{A}}_{(j)}$ of ${\mathsf{A}}$ for some $j$, and updates $\mathbf{z}^{(k)}$ by projecting it onto the orthogonal complement of the space of ${\mathsf{A}}_{(j)}$. The claim is that randomly selecting the columns of ${\mathsf{A}}$ with probability proportional to their square norms implies that the algorithm converges to ${{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})^{\bot}}}$ in expectation. After $T$ iterations, the algorithm outputs $\mathbf{z}^{(T)}$ and by orthogonality ${\mathbf{b}}-\mathbf{z}^{(T)}$ serves as an approximation for ${{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})}}$. The next theorem bounds the expected rate of convergence for Algorithm 1.

Moreover, each iteration of Algorithm 1 requires in expectation (over the random choices of the algorithm) at most $5\text{C}_{\text{avg}}$ arithmetic operations.

A suggestion for a stopping criterion for Algorithm 1 is to regularly check: $\frac{\left\|{\mathsf{A}}^{\top}\mathbf{z}^{(k)}\right\|_{2}}{\left\|{\mathsf{A}}\right\|_{\text{\rm F}}\left\|\mathbf{z}^{(k)}\right\|_{2}}\leq\varepsilon$ for some given accuracy $\varepsilon>0$. It is easy to see that whenever this criterion is satisfied, it holds that $\left\|{{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})^{\bot}}}-\mathbf{z}^{(k)}\right\|_{2}/\left\|\mathbf{z}^{(k)}\right\|_{2}\leq\varepsilon\kappa_{\textrm{\tiny F}}({\mathsf{A}})$, i.e., ${\mathbf{b}}-\mathbf{z}^{(k)}\approx{{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})}}$.

We devote the rest of this subsection to prove Theorem 2. Define ${\mathsf{P}}(j):=\mathbf{I}_{m}-\frac{{\mathsf{A}}_{(j)}{\mathsf{A}}_{(j)}^{\top}}{\left\|{\mathsf{A}}_{(j)}\right\|_{2}^{2}}$ for every $j\in[n]$. Observe that ${\mathsf{P}}(j){\mathsf{P}}(j)={\mathsf{P}}(j)$, i.e., ${\mathsf{P}}(j)$ is a projector matrix. Let $X$ be a random variable over $\{1,2,\ldots,n\}$ that picks index $j$ with probability $\left\|{\mathsf{A}}_{(j)}\right\|_{2}^{2}/\left\|{\mathsf{A}}\right\|_{\text{\rm F}}^{2}$. It is clear that $\EE[{\mathsf{P}}(X)]=\mathbf{I}_{m}-{\mathsf{A}}{\mathsf{A}}^{\top}/\left\|{\mathsf{A}}\right\|_{\text{\rm F}}^{2}$. Later we will make use of the following fact.

For every vector $\mathbf{u}$ in the column space of ${\mathsf{A}}$, it holds $\left\|\left(\mathbf{I}_{m}-\frac{{\mathsf{A}}{\mathsf{A}}^{\top}}{\left\|{\mathsf{A}}\right\|_{\text{\rm F}}^{2}}\right)\mathbf{u}\right\|_{2}\leq\left(1-\frac{\sigma^{2}_{\min}}{\left\|{\mathsf{A}}\right\|_{\text{\rm F}}^{2}}\right)\left\|\mathbf{u}\right\|_{2}$.

Define ${\mathbf{e}}^{(k)}:=\mathbf{z}^{(k)}-{{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})^{\bot}}}$ for every $k\geq 0$. A direct calculation implies that

Indeed, ${\mathbf{e}}^{(k)}=\mathbf{z}^{(k)}-{{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})^{\bot}}}={\mathsf{P}}(j_{k})\mathbf{z}^{(k-1)}-{{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})^{\bot}}}={\mathsf{P}}(j_{k})({\mathbf{e}}^{(k-1)}+{{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})^{\bot}}})-{{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})^{\bot}}}={\mathsf{P}}(j_{k}){\mathbf{e}}^{(k-1)}$ using the definitions of ${\mathbf{e}}^{(k)}$, $\mathbf{z}^{(k)}$, ${\mathbf{e}}^{(k-1)}$ and the fact that ${\mathsf{P}}(j_{k}){{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})^{\bot}}}={{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})^{\bot}}}$ for any $j_{k}\in{[n]}$. Moreover, it is easy to see that for every $k\geq 0$ ${\mathbf{e}}^{(k)}$ is in the column space of ${\mathsf{A}}$, since ${\mathbf{e}}^{(0)}={\mathbf{b}}-{{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})^{\bot}}}={{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})}}\in\mathcal{R}({\mathsf{A}})$, ${\mathbf{e}}^{(k)}={\mathsf{P}}(j_{k}){\mathbf{e}}^{(k-1)}$ and in addition ${\mathsf{P}}(j_{k})$ is a projector matrix for every $j_{k}\in[n]$.

Let $X_{1},X_{2},\ldots$ be a sequence of independent and identically distributed random variables distributed as $X$. For ease of notation, we denote by $\EE_{k-1}[\cdot]=\EE_{X_{k}}[\cdot\ |\ X_{1},X_{2},\ldots,X_{k-1}]$, i.e., the conditional expectation conditioned on the first $(k-1)$ iteration of the algorithm. It follows that

where we used linearity of expectation, the fact that ${\mathsf{P}}(\cdot)$ is a projector matrix, Cauchy-Schwarz inequality and Fact 3. Repeating the same argument $k-1$ times we get that

Note that ${\mathbf{e}}^{(0)}={\mathbf{b}}-{{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})^{\bot}}}={{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})}}$ to conclude.

2 Randomized Kaczmarz

Strohmer and Vershynin proposed the following randomized variant of Kaczmarz algorithm (Algorithm 2), see [SV09] for more details. The following theorem is a restatement of the main result of [SV09] without imposing the full column rank assumption.

We devote the rest of this subsection to prove Theorem 4 following [SV09]. The proof is based on the following two elementary lemmas which both appeared in [SV09]. However, in our setting, the second lemma is not identical to that in [SV09]. We deferred their proofs to the Appendix.

Assume that ${\mathsf{A}}\mathbf{x}={\mathbf{b}}$ has a solution and use the notation of Algorithm 2, then $\mathbf{x}^{(k+1)}-\mathbf{x}_{\text{\tiny LS}}$ is perpendicular to $\mathbf{x}^{(k+1)}-\mathbf{x}^{(k)}$ for any $k\geq 0$. In particular, in exact arithmetic it holds that $\left\|\mathbf{x}^{(k+1)}-\mathbf{x}_{\text{\tiny LS}}\right\|_{2}^{2}=\left\|\mathbf{x}^{(k)}-\mathbf{x}_{\text{\tiny LS}}\right\|_{2}^{2}-\left\|\mathbf{x}^{(k+1)}-\mathbf{x}^{(k)}\right\|_{2}^{2}$.

The above lemma provides a formula for the error at each iteration. Ideally, we seek to minimize the error at each iteration which is equivalent to maximizing $\left\|\mathbf{x}^{(k+1)}-\mathbf{x}^{(k)}\right\|_{2}$ over the choice of the row projections of the algorithm. The next lemma suggests that by randomly picking the rows of ${\mathsf{A}}$ reduces the error in expectation.

Theorem 4 follows by iterating Lemma 6, we get that

3 Randomized Kaczmarz Applied to Noisy Linear Systems

The analysis of Strohmer and Vershynin is based on the restrictive assumption that the linear system has a solution. Needell made a step further and analyzed the more general setting in which the linear system does not have any solution and ${\mathsf{A}}$ has full column rank [Nee10]. In this setting, it turns out that the randomized Kaczmarz algorithm computes an estimate vector that is within a fixed distance from the solution; the distance is proportional to the norm of the “noise vector” multiplied by $\kappa^{2}_{\textrm{\tiny F}}({\mathsf{A}})$ [Nee10]. The following theorem is a restatement of the main result in [Nee10] with two modifications: the full column rank assumption on the input matrix is dropped and the additive term $\gamma$ of Theorem $2.1$ in [Nee10] is improved to $\left\|{\mathbf{w}}\right\|_{2}^{2}/\left\|{\mathsf{A}}\right\|_{\text{\rm F}}^{2}$. The only technical difference here from [Nee10] is that the full column rank assumption is not necessary, so we defer the proof to the Appendix for completeness.

Randomized Extended Kaczmarz

In the present paper, the main observation is that it is possible to efficiently reduce the norm of the “noisy” part of ${\mathbf{b}}$, ${{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})^{\bot}}}$ (using Algorithm 1) and then apply the randomized Kaczmarz algorithm on a new linear system whose right hand side vector is now arbitrarily close to the column space of ${\mathsf{A}}$, i.e., ${\mathsf{A}}\mathbf{x}\approx{{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})}}$. This idea together with the observation that the least squares solution of the latter linear system is equal (in the limit) to the least squares solution of the original system (see Fact 1) implies a randomized algorithm for solving least squares.

Next we present the randomized extended Kaczmarz algorithm which is a specific combination of the randomized orthogonal projection algorithm together with the randomized Kaczmarz algorithm.

The stopping criterion of Step 8 was decided based on the following analysis. Assume that the termination criteria are met for some $k>0$. Let $\mathbf{z}^{(k)}={{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})^{\bot}}}+{\mathbf{w}}$ for some ${\mathbf{w}}\in\mathcal{R}({\mathsf{A}})$ (which holds by the definition of $\mathbf{z}^{(k)}$). Then,

By re-arranging terms and using the second part of the termination criterion, it follows that $\left\|\mathbf{z}^{(k)}-{{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})^{\bot}}}\right\|_{2}\leq\varepsilon\frac{\left\|{\mathsf{A}}\right\|_{\text{\rm F}}^{2}}{\sigma_{\min}}\left\|\mathbf{x}^{(k)}\right\|_{2}$. Now,

where we used the triangle inequality, the first part of the termination rule together with ${{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})}}={\mathsf{A}}\mathbf{x}_{\text{\tiny LS}}$ and the above discussion. Now, since $\mathbf{x}^{(k)},\mathbf{x}_{\text{\tiny LS}}\in\mathcal{R}({\mathsf{A}}^{\top})$, it follows that

Equation (6) demonstrates that the forward error of REK after termination is bounded.

2 Rate of convergence

The following theorem bounds the expected rate of convergence of Algorithm 3.

After $T>1$ iterations, in exact arithmetic, Algorithm 3 with input ${\mathsf{A}}$ (possibly rank-deficient) and ${\mathbf{b}}$ computes a vector $\mathbf{x}^{(T)}$ such that

For the sake of notation, set $\alpha=1-1/\kappa^{2}_{\textrm{\tiny F}}({\mathsf{A}})$ and denote by $\EE_{k}[\cdot]:=\EE[\cdot\ |\ i_{0},j_{0},i_{1},j_{1},\ldots,i_{k},j_{k}]$, i.e., the conditional expectation with respect to the first $k$ iterations of Algorithm 3. Observe that Steps $5$ and $6$ are independent from Steps $4$ and $7$ of Algorithm 3, so Theorem 2 implies that for every $l\geq 0$

Fix a parameter $k^{*}:=\lfloor T/2\rfloor$. After the $k^{*}$-th iteration of Algorithm 3, it follows from Theorem 7 (Inequality (5)) that

Indeed, the randomized Kaczmarz algorithm is executed with input $({\mathsf{A}},{\mathbf{b}}-\mathbf{z}^{(k^{*}-1)})$ and current estimate vector $\mathbf{x}^{(k^{*}-1)}$. Set $\mathbf{y}={{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})}}$ and ${\mathbf{w}}={{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})^{\bot}}}-\mathbf{z}^{(k^{*}-1)}$ in Theorem 7 and recall that $\mathbf{x}_{\text{\tiny LS}}={{\mathsf{A}}}^{\dagger}{\mathbf{b}}={{\mathsf{A}}}^{\dagger}{{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})}}={{\mathsf{A}}}^{\dagger}\mathbf{y}$.

Now, averaging the above inequality over the random variables $i_{1},j_{1},i_{2},j_{2},\ldots,i_{k^{*}-1},j_{k^{*}-1}$ and using linearity of expectation, it holds that

Simplifying the right hand side using the fact that $\sum_{l=0}^{\infty}\alpha^{l}=\frac{1}{1-\alpha}=\kappa^{2}_{\textrm{\tiny F}}({\mathsf{A}})$, it follows

Moreover, observe that for every $l\geq 0$

Now for any $k>0$, similar considerations as Ineq. (8) implies that

To derive the last inequality, consider two cases. If $T$ is even, set $k=k^{*}$, otherwise set $k=k^{*}+1$. In both cases, $(\alpha^{k}+\alpha^{k^{*}})\leq 2\alpha^{k^{*}}$. ∎

3 Theoretical bounds on time complexity

In this section, we discuss the running time complexity of the randomized extended Kaczmarz (Algorithm 3). Recall that REK is a Las-Vegas randomized algorithm, i.e., the algorithm always outputs an “approximately correct” least squares estimate (satisfying (6)) but its runnning time is a random variable. Given any fixed accuracy parameter $\varepsilon>0$ and any fixed failure probability $0<\delta<1$ we bound the number of iterations required by the algorithm to terminate with probability at least $1-\delta$.

Fix an accuracy parameter $0<\varepsilon<2$ and failure probability $0<\delta<1$. In exact arithmetic, Algorithm 3 terminates after at most

iterations with probability at least $1-\delta$.

Denote $\alpha:=1-1/\kappa^{2}_{\textrm{\tiny F}}({\mathsf{A}})$ for notational convenience. It suffices to prove that with probability at least $1-\delta$ the conditions of Step 8 of Algorithm 3 are met. Instead of proving this, we will show that:

With probability at least $1-\delta/2$: $\left\|({\mathbf{b}}-\mathbf{z}^{(T^{*})})-{{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})}}\right\|_{2}\leq\varepsilon\left\|{{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})}}\right\|_{2}/4$.

With probability at least $1-\delta/2$: $\left\|\mathbf{x}^{(T^{*})}-\mathbf{x}_{\text{\tiny LS}}\right\|_{2}\leq\varepsilon\left\|\mathbf{x}_{\text{\tiny LS}}\right\|_{2}/4$.

Later we prove that Items (1) and (2) imply the Lemma. First we prove Item (1). By the definition of the algorithm,

the first equality follows since ${\mathbf{b}}-{{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})}}={{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})^{\bot}}}$, the second inequality is Markov’s inequality, the third inequality follows by Theorem 2, and the last inequality since $T^{*}\geq\kappa^{2}_{\textrm{\tiny F}}({\mathsf{A}})\ln(\frac{32}{\delta\varepsilon^{2}})$.

the first inequality is Markov’s inequality, the second inequality follows by Theorem 8, and the last inequality follows provided that $T^{*}\geq 2\kappa^{2}_{\textrm{\tiny F}}({\mathsf{A}})\ln\left(\frac{32(1+2\kappa^{2}\left({\mathsf{A}}\right))}{\delta\varepsilon^{2}}\right)$

A union bound on the complement of the above two events (Item (1) and (2)) implies that both events happen with probability at least $1-\delta$. Now we show that conditioning on Items (1) and (2), it follows that REK terminates after $T^{*}$ iterations, i.e.,

We start with the first condition. First, using triangle inequality and Item 2, it follows that

where the first inequality is triangle inequality, the second inequality follows by Item $1$ and ${{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})}}={\mathsf{A}}\mathbf{x}_{\text{\tiny LS}}$, the third and forth inequality follows by Item $2$ and the fifth inequality holds by Inequality (11) and the last inequality follows since $\varepsilon<2$. The second condition follows since

the first equation follows by orthogonality, the second inequality assuming Item (2), the third inequality follows since ${{\mathbf{b}}_{\mathcal{R}({\mathsf{A}})}}={\mathsf{A}}\mathbf{x}_{\text{\tiny LS}}$, the forth inequality follows by (11) and the final inequality since $\varepsilon<2$. ∎

Lemma 9 bounds the number of iterations with probability at least $1-\delta$, next we bound the total number of arithmetic operations in worst case (Eqn. (12)) and in expectation (Eqn. (13)). Let’s calculate the computational cost of REK in terms of floating-point operations (flops) per iteration. For the sake of simplicity, we ignore the additional (negligible) computational overhead required to perform the sampling operations (see Section 5 for more details) and checking for convergence.

Each iteration of Algorithm 3 requires four level-1 BLAS operations (two DDOT operations of size $m$ and $n$, respectively, and two DAXPY operations of size $n$ and $m$, respectively) and additional four flops. In total, $4(m+n)+2$ flops per iteration.

Therefore by Lemma 9, with probability at least $1-\delta$, REK requires at most

Exploiting the (possible) sparsity of ${\mathsf{A}}$, we first show that each iteration of Algorithm 3 requires at most $5(\text{C}_{\text{avg}}+\text{R}_{\text{avg}})$ operations in expectation. For simplicity of presentation, we assume that we have stored ${\mathsf{A}}$ in compressed column sparse format and compressed row sparse format [BBC+87].

By the linearity of expectation and the definitions of $\text{C}_{\text{avg}}$ and $\text{R}_{\text{avg}}$, the expected running time after $T^{*}$ iterations is at most $5T^{*}(\text{C}_{\text{avg}}+\text{R}_{\text{avg}})$. It holds that (recall that $p_{j}=\left\|{\mathsf{A}}_{(j)}\right\|_{2}^{2}/\left\|{\mathsf{A}}\right\|_{\text{\rm F}}^{2}$)

Hence by Lemma 9, with probability at least $1-\delta$, the expected number of arithmetic operations of REK is at most

Implementation and Experimental Results

The proposed algorithm has been entirely implemented in C. We provide three implementation of Algorithm 3: REK-C, REK-BLAS and REK-BLAS-PRECOND. REK-C corresponds to a direct translation of Algorithm 3 to C code. REK-BLAS is an implementation of REK with two additional technical features. First, REK-BLAS uses level-1 BLAS routines for all operations of Algorithm 3 and secondly REK-BLAS additionally stores explicitly the transpose of ${\mathsf{A}}$ for more efficiently memory access of both the rows and columns of ${\mathsf{A}}$ using BLAS. REK-BLAS-PRECOND is an implementation of REK-BLAS that additionally supports upper triangular preconditioning; we used Blendenpik’s preconditioning code to ensure a fair comparison with Blendenpik (see Section 5.2). In the implementations of REK-C and REK-BLAS we check for convergence every $8\min(m,n)$ iterations.

Moreover, all implementations include efficient code that handles sparse input matrices using the compressed column (and row) sparse matrix format [BBC+87].

The sampling operations of Algorithm 3 (Steps 4 and 5) are implemented using the so-called “alias method” for generating samples from any given discrete distribution [Wal77, Vos91]. The alias method, assuming access to a uniform random variable on $$in constant time and linear time preprocessing, generates one sample of the given distribution in constant time [Vos91]. We use an implementation of W. D. Smith that is described in [Smi02] and C’s drand48() to get uniform samples from$$.

2 Experimental Results

We report our experimental results in this section. We compared the randomized extended Kaczmarz (REK-C, REK-BLAS,REK-BLAS-PRECOND) algorithm to LAPACK’s DGELSY and DGELSD least squares solvers, Blendenpick Available at http://www.mathworks.com/matlabcentral/fileexchange/25241-blendenpik. Blendenpik’s default settings were used. (version 1.3, [AMT10]) and MATLAB’s backslash operator. LSRN [MSM11] did not perform well under a setup in which no parallelization is allowed as the one used here, so we do not include LSRN’s performance. DGELSY uses QR factorization with pivoting and DGELSD uses the singular value decomposition. We use MATLAB with version 7.9.0.529 (R2009b). In addition, we use MATLAB’s included BLAS and LAPACK packages and we call LAPACK’s functions from MATLAB using MATLAB’s CMEX technology which allows us to measure only LAPACK’s elapsed time. We should highlight that MATLAB is used as a scripting language and no MATLAB-related overheads have been taken under consideration. Blendenpik requires the FFTW library http://www.fftw.org/; we used FFTW-3.3.3. To match the accuracy of LAPACK’s direct solvers, we fixed $\varepsilon$ in Algorithm 3 to be 10e-14. Moreover, during our experiments we ensured that the residual error of all the competing algorithms were about of the same order of magnitude.

We used a Pentium(R) Dual-Core E5300 (2.60GHz) equipped with 5GB of RAM and compiled our source code using GCC-4.7.2 under Linux operating system. All running times displayed below are measured using the ftime Linux system call by taking the average of the running time of 10 independent executions.

We experimented our algorithm under three different distributions of random input matrices (sparse, dense and ill-conditioned) under the setting of strongly rectangular settings of least squares instances. In all cases we normalized the column norms of the input matrices to unity and generate the right hand side vector ${\mathbf{b}}$ having Gaussian entries of variance one.

We tested our algorithm in the overdetermined setting of random sparse $m\times n$ matrices with $n=800$ and $m=2000,3000,\ldots,20000$ and density $0.25$. We also tested REK-BLAS on the underdetermined case where $m=800$ and $n=2000,3000,\ldots,20000$. In both cases, the density of the sparse matrices was set to $0.25$ (for even sparser matrices REK-BLAS performed even better compared to all other mentioned methods). To generate these sparse matrix ensembles, we used MATLAB’s sprandn function with variance one. The results are depicted in Figure 1. Both plots demonstrate that REK-BLAS is superior on both the underdetermined (Figure 1(b)) and overdetermined case (Figure 1(a)). It is interesting that REK-BLAS performs well in the underdetermined case.

Dense and well-conditioned least squares

In this scenario, we used random overdetermined dense $m\times n$ matrices with much more rows than columns, i.e., we set $n=500$ and $m=1000,2000,\ldots,20000$. We also tested REK-BLAS on the underdetermined case where $m=500$ and $n=1000,2000,\ldots,20000$. We generated this set of matrices using MATLAB’s randn function with variance ten. We depicted the results in Figure 2(a). In the overdetermined case (Figure 2(a)), REK-BLAS is marginally superior compared to LAPACK’s routines whereas REK-C (as a naive implementation of Algorithm 3) is inferior. Blendepik is the winner in this case. Interestingly, REK-BLAS almost matches the performance of Blendenpik in the underdetermined case, see Figure 2(b).

Dense and ill-conditioned least squares

Finally, we tested all algorithms under a particular case of random ill-conditioned dense matrices with $n=500$ and $m=1000,2000,\ldots,20000$. Namely, we used Higham’s randSVD function for generating these matrices [Hig89, Hig96]. More precisely, we set the condition number of these matrices to be $10e6$; set the top singular value to one and the rest to 10e-6. The results are displayed in Figure 3. Unfortunately, in the ill-conditioned setting REK-BLAS-PRECOND is inferior compared to LAPACK’s routines and Blendepik. We also verified the results of [AMT10] that Blendepik is superior compared to LAPACK’s least squares solvers in this setting.

Acknowledgements

We would like to thank the anonymous reviewers for their invaluable comments on an earlier draft of the present manuscript. The first author would like to thank Haim Avron for his technical support on several issues regarding Blendenpik and Philip A. Knight for sharing his unpublished manuscript [Kni96].

References

Appendix

We present the proof of known facts from previous works for completeness.

(of Lemma 5) It suffices to show that $\left\langle{\mathbf{x}^{(k+1)}-\mathbf{x}_{\text{\tiny LS}}},\ {\mathbf{x}^{(k+1)}-\mathbf{x}^{(k)}}\right\rangle=0$. For notational convenience, let $\alpha_{i}:=\frac{b_{i}-\left\langle{\mathbf{x}^{(k)}},\ {{\mathsf{A}}^{(i)}}\right\rangle}{\left\|{\mathsf{A}}^{(i)}\right\|_{2}^{2}}$ for every $i\in{[m]}$. Assume that $\mathbf{x}^{(k+1)}=\mathbf{x}^{(k)}+\alpha_{i_{k}}{\mathsf{A}}^{(i_{k})}$ for some arbitrary $i_{k}\in[m]$. Then,

using the definition of $\mathbf{x}^{(k+1)}$, and the fact that $\left\langle{\mathbf{x}_{\text{\tiny LS}}},\ {{\mathsf{A}}^{(i_{k})}}\right\rangle=b_{i_{k}}$ since $\mathbf{x}_{\text{\tiny LS}}$ is a solution to ${\mathsf{A}}\mathbf{x}={\mathbf{b}}$. Now, by the definition of $\alpha_{i_{k}}$, $\left\langle{\mathbf{x}^{(k+1)}},\ {{\mathsf{A}}^{(i_{k})}}\right\rangle=\left\langle{\mathbf{x}^{(k)}},\ {{\mathsf{A}}^{(i_{k})}}\right\rangle+\alpha_{i_{k}}\left\|{\mathsf{A}}^{(i_{k})}\right\|_{2}^{2}=\left\langle{\mathbf{x}^{(k)}},\ {{\mathsf{A}}^{(i_{k})}}\right\rangle+b_{i_{k}}-\left\langle{\mathbf{x}^{(k)}},\ {{\mathsf{A}}^{(i_{k})}}\right\rangle=b_{i_{k}}$. ∎

(of Lemma 6) In light of Lemma 5, it suffices to show that $\EE_{Z}\left\|\mathbf{x}^{(k+1)}-\mathbf{x}^{(k)}\right\|_{2}^{2}\geq\frac{1}{\kappa^{2}_{\textrm{\tiny F}}({\mathsf{A}})}\left\|\mathbf{x}^{(k)}-\mathbf{x}_{\text{\tiny LS}}\right\|_{2}^{2}$. By the definition of $\mathbf{x}^{(k+1)}$, it follows

By hypothesis, $\mathbf{x}^{(k)}$ is in the row space of ${\mathsf{A}}$ for any $k$ when $\mathbf{x}^{(0)}$ is; in addition, the same is true for $\mathbf{x}_{\text{\tiny LS}}$ by the definition of pseudo-inverse [GL96]. Therefore, $\left\|{\mathsf{A}}(\mathbf{x}_{\text{\tiny LS}}-\mathbf{x}^{(k)})\right\|_{2}\geq\sigma_{\min}\left\|\mathbf{x}_{\text{\tiny LS}}-\mathbf{x}^{(k)}\right\|_{2}$. ∎

(of Theorem 7) As in [Nee10], for any $i\in{[m]}$ define the affine hyper-planes:

Assume for now that at the $k$-th iteration of the randomized Kaczmarz algorithm applied on $({\mathsf{A}},{\mathbf{b}})$, the $i$-th row is selected. Note that $\hat{\mathbf{x}}^{(k)}$ is the projection of $\hat{\mathbf{x}}^{(k-1)}$ on $\mathcal{H}_{i}^{w_{i}}$ by the definition of the randomized Kaczmarz algorithm on input $({\mathsf{A}},{\mathbf{b}})$. Let us denote the projection of $\hat{\mathbf{x}}^{(k-1)}$ on $\mathcal{H}_{i}$ by $\mathbf{x}^{(k)}$. The two affine hyper-planes $\mathcal{H}_{i},\mathcal{H}_{i}^{w_{i}}$ are parallel with common normal ${\mathsf{A}}^{(i)}$, so $\mathbf{x}^{(k)}$ is the projection of $\hat{\mathbf{x}}^{(k)}$ on $\mathcal{H}_{i}$ and the minimum distance between $\mathcal{H}_{i}$ and $\mathcal{H}_{i}^{w_{i}}$ equals $|w_{i}|/\left\|{\mathsf{A}}^{(i)}\right\|_{2}$. In addition, $\mathbf{x}^{*}\in\mathcal{H}_{i}$ since $\left\langle{\mathbf{x}^{*}},\ {{\mathsf{A}}^{(i)}}\right\rangle=y_{i}$, therefore by orthogonality we get that

Since $\mathbf{x}^{(k)}$ is the projection of $\hat{\mathbf{x}}^{(k-1)}$ onto $\mathcal{H}_{i}$ (that is to say, $\mathbf{x}^{(k)}$ is a randomized Kaczmarz step applied on input $({\mathsf{A}},\mathbf{y})$ where the $i$-th row is selected on the $k$-th iteration) and $\hat{\mathbf{x}}^{(k-1)}$ is in the row space of ${\mathsf{A}}$, Lemma 6 tells us that

Note that for given selected row $i$ we have $\left\|\hat{\mathbf{x}}^{(k)}-\mathbf{x}^{(k)}\right\|_{2}^{2}=\frac{w_{i}^{2}}{\left\|{\mathsf{A}}^{(i)}\right\|_{2}^{2}}$; by the distribution of selecting the rows of ${\mathsf{A}}$ we have that

Inequality (5) follows by taking expectation on both sides of Equation (14) and bounding its resulting right hand side using Equations (15) and (16). Applying Inequality (5) inductively, it follows that

where we used that $\mathbf{x}^{(0)}$ is in the row space of ${\mathsf{A}}$. The latter sum is bounded above by $\sum_{i=0}^{\infty}\left(1-\frac{1}{\kappa^{2}_{\textrm{\tiny F}}({\mathsf{A}})}\right)^{i}=\left\|{\mathsf{A}}\right\|_{\text{\rm F}}^{2}/\sigma^{2}_{\min}$. ∎