Source-linked AI summary
Randomized Extended Kaczmarz for Solving Least-Squares
Anastasios Zouzias, Nikolaos Freris
TL;DR
Least-squares systems may be inconsistent, limiting randomized Kaczmarz when the residual is not negligible. The paper removes the noisy component of the right-hand side before applying randomized Kaczmarz, yielding an expected-convergent solver for the minimum-norm least-squares solution. The method is evaluated across sparse, rectangular, and ill-conditioned settings, where performance varies by regime.
Problem
Randomized Kaczmarz analyses addressed solvable systems or were effective mainly when the least-squares error was negligible.
Method
The algorithm reduces the component of b orthogonal to the column space of A, then applies randomized Kaczmarz to the resulting system.
Results
The algorithm converges in expectation to the minimum Euclidean norm least-squares solution, with a convergence guarantee for possibly rank-deficient A.
Takeaways & Limitations
Randomized extended Kaczmarz provides a randomized iterative approach for least-squares problems beyond the negligible-residual setting targeted by randomized Kaczmarz.
Takeaways & Limitations
The experimental results show that performance depends on the problem regime: the method is inferior to LAPACK and Blendenpik for the tested ill-conditioned dense matrices.
Abstract
from arXiv · showhide
We present a randomized iterative algorithm that exponentially converges in expectation to the minimum Euclidean norm least squares solution of a given linear system of equations. The expected number of arithmetic operations required to obtain an estimate of given accuracy is proportional to the square condition number of the system multiplied by the number of non-zeros entries of the input matrix. The proposed algorithm is an extension of the randomized Kaczmarz method that was analyzed by Strohmer and Vershynin.
1 Introduction
The paper develops randomized extensions of Kaczmarz methods for least-squares systems, building on exponential expected convergence results that assumed solvability. Its proposed randomized extended method targets the minimum Euclidean norm least-squares solution.
- Kaczmarz method: Kaczmarz iteratively projects an estimate onto affine hyperplanes defined by selected rows of Ax = b.The standard method cycles through rows, while relaxation parameters control the projection updates.
- Kaczmarz method: For overdetermined inconsistent systems, Kaczmarz can converge to the least-squares estimate under full column rank and suitable diminishing relaxation parameters.The cited condition is relaxation parameters in [0, 2] with λk → 0.
- Randomized Kaczmarz: Randomized Kaczmarz selects rows non-uniformly and has exponential convergence in expectation when the linear system is solvable.Needell extended the analysis to inconsistent full-column-rank systems, obtaining convergence within a distance proportional to b’s distance from the column space of A.
- Randomized extended Kaczmarz: The proposed randomized extended Kaczmarz algorithm is a randomized variant of Popa’s extended Kaczmarz method, based on prior randomized Kaczmarz analyses.The paper presents it as a randomized iterative least-squares solver converging to the minimum Euclidean norm solution.
- Paper organization: The paper organizes related work, randomized subspace projection, convergence analysis for solvable systems, and the proposed least-squares algorithm across its sections.The stated organization identifies these as the main components of the development.
2 Least squares solvers
The paper situates randomized extended Kaczmarz among deterministic direct and iterative solvers, random projection methods, and randomized preconditioning approaches. It compares these methods experimentally on several least-squares settings.
- Deterministic algorithms: Deterministic least-squares solvers include pivoted QR, SVD, and Krylov-subspace methods such as LSQR applied to normal equations.LAPACK provides implementations of the direct QR and SVD approaches, while LSQR is described as numerically stable.
- Randomized algorithms: Many randomized algorithms first project the system into fewer dimensions and then solve the resulting down-sampled system.This two-step procedure is presented as a generic approach in the theoretical computer science literature.
- Randomized algorithms: Randomized preconditioning projects the input matrix, extracts a preconditioner, and applies an iterative solver such as LSQR to the preconditioned system.Blendenpik is given as an implementation of this approach.
- Randomized algorithms: LSRN uses normal random projections for preconditioning and then applies LSQR or a Chebyshev semi-iterative method.Its design is described as a two-phase preconditioning and iterative-solving approach.
- Experimental comparison: The paper compares randomized extended Kaczmarz with DGELSY, DGELSD, and Blendenpik, excluding LSRN because it performed poorly without parallelization.The numerical evaluation focuses on sparse, well-conditioned, and strongly rectangular least-squares problems.
3 Background
This section develops randomized projection methods for least-squares problems, including approximate projection onto a matrix column space and randomized row projections. It establishes expected convergence results and discusses limitations of prior randomized Kaczmarz analyses.
- Preliminaries and Notation: The column space R(A) is the set of all vectors Ax, while R(A)⊥ is its orthogonal complement.Any b can be decomposed into projections onto R(A) and R(A)⊥; the minimum-norm least-squares solution is A†b.
- 3.1 Randomized Approximate Orthogonal Projection: Randomized Orthogonal Projection repeatedly selects columns of A with probability proportional to squared column norms and projects onto their orthogonal complements.The output approximates the projection of b onto R(A) through orthogonality.
- 3.1 Randomized Approximate Orthogonal Projection: Theorem 2 bounds the expected error after k iterations of Algorithm 1, while each iteration requires at most 5Cavg expected arithmetic operations.The algorithm therefore combines an expected convergence guarantee with a cost expressed through average column sparsity.
- 3.2 Randomized Kaczmarz: Randomized Kaczmarz selects rows randomly and projects iterates onto row-defined constraints; under a solvable system, it converges in expectation to xLS.The analysis relies on expected error reduction from random row selection and applies the Moore–Penrose solution xLS := A†b.
- 3.3 Randomized Kaczmarz Applied to Noisy Linear Systems: The Strohmer–Vershynin analysis assumes Ax = b has a solution, whereas Needell considers inconsistent systems with full column rank.For inconsistent systems, randomized Kaczmarz approaches an estimate within a fixed distance of the least-squares solution.
- 3.3 Randomized Kaczmarz Applied to Noisy Linear Systems: The paper removes the full-column-rank assumption from Needell’s setting and improves the additive term to ∥w∥2^2 / ∥A∥2_F.The stated technical difference is that full column rank is unnecessary in the revised analysis.
4 Randomized Extended Kaczmarz
Randomized Extended Kaczmarz combines randomized orthogonal projection with randomized Kaczmarz to approach the minimum Euclidean-norm least-squares solution in expectation. Its runtime guarantees provide high-probability termination and an expected cost proportional to the matrix sparsity and squared condition number.
- Algorithm: The algorithm reduces the component of b orthogonal to the column space, then applies randomized Kaczmarz to the resulting increasingly consistent system.The two components maintain an approximation b − z(k) to the column-space projection of b and solve using that approximation.
- Termination: The stopping rule certifies bounded forward error after termination, using conditions on the residual approximation and the iterate.The analysis shows that the termination criteria imply the forward-error guarantee in Equation (6).
- Iteration bounds: For accuracy 0 < ε < 2 and failure probability 0 < δ < 1, REK terminates after a bounded number of iterations with probability at least 1 − δ.Lemma 9 establishes the iteration bound under exact arithmetic.
- Time complexity: The worst-case operation bound scales with (m + n) rank(A) κ_F^2(A), while sparse implementation reduces expected running time to nnz(A) times the squared condition number.Each iteration uses sparse row and column operations in compressed sparse formats, yielding the tighter expected bound.
5 Implementation and Experimental Results
The implementation uses C-based REK variants with BLAS, transpose storage, sparse formats, and alias-method sampling, then evaluates them against direct and randomized solvers across sparse, dense, and ill-conditioned systems.
- Implementation: REK-C directly translates Algorithm 3, while REK-BLAS adds level-1 BLAS operations and explicit transpose storage for row and column access.REK-BLAS-PRECOND additionally supports preconditioning, although the supplied passage truncates further details.
- Implementation: Sparse matrices are handled with compressed column and row formats, and non-uniform sampling uses the alias method with constant-time sample generation after linear preprocessing.The alias method assumes constant-time access to a uniform random variable on [0, 1].
- Experimental setup: The experiments compare three REK implementations with LAPACK’s DGELSY and DGELSD, Blendenpik, and MATLAB’s backslash operator.DGELSY uses QR factorization with pivoting, whereas DGELSD uses the singular value decomposition.
- Experimental setup: The evaluation uses sparse, dense, and ill-conditioned random matrices in strongly rectangular overdetermined and underdetermined settings, with normalized column norms and Gaussian right-hand sides.Running times average 10 independent executions on a Pentium Dual-Core E5300 system with 5GB RAM.
- Sparse least squares: REK-BLAS is superior to the other methods for the tested random sparse overdetermined and underdetermined least-squares problems.The sparse matrices have density 0.25; the passage reports even better relative performance for still sparser matrices.
- Dense and well-conditioned least squares: For dense well-conditioned problems, Blendenpik wins in the overdetermined case, while REK-BLAS almost matches Blendenpik in the underdetermined case.In the overdetermined case, REK-BLAS is marginally superior to LAPACK’s routines, whereas REK-C is inferior.
- Dense and ill-conditioned least squares: For dense ill-conditioned matrices with condition number 10e6, REK-BLAS-PRECOND is inferior to LAPACK’s routines and Blendenpik, while Blendenpik is superior to LAPACK’s solvers.These tests use 500 columns and vary the number of rows from 1000 to 20000.
7 Appendix
The appendix proves supporting lemmas and Theorem 7 using randomized Kaczmarz projections onto affine hyperplanes, row-space invariance, conditional expectations, and induction.
- Supporting lemmas: Lemma 5 analyzes an update of the form x(k+1) = x(k) + α_i^k A(i_k) for an arbitrary selected row.The argument proceeds from the definition of the next iterate and bounds the resulting quantity row by row.
- Proof setup: The appendix introduces affine hyperplanes associated with rows of A and analyzes randomized Kaczmarz projections onto them.For a selected row, the iterate is described as the projection onto the corresponding hyperplane.
- Proof setup: The two relevant hyperplanes are parallel with a common normal, so one projection can be interpreted as a projection of the other iterate onto the target hyperplane.This geometric relation underlies the comparison between the auxiliary and Kaczmarz iterates.
- Invariant subspace: The proof uses that iterates remain in the row space of A when the initial vector does, and that the minimum-norm least-squares solution also lies there.This row-space property is invoked to apply the supporting lemmas.
- Theorem 7: Inequality (5) follows by taking expectations in Equation (14), bounding its right-hand side with Equations (15) and (16), and applying the inequality inductively.The derivation also uses that the initial iterate lies in the row space of A.