Randomized Iterative Methods for Linear Systems

Robert M. Gower; Peter Richtárik

Source-linked AI summary

Randomized Iterative Methods for Linear Systems

Robert M. Gower, Peter Richtárik

arXiv:1506.03296v5Jun 10, 2015math.NA

TL;DR

The paper addresses the need for broadly applicable randomized methods for solving consistent linear systems. It develops a two-parameter method with six equivalent interpretations and unifies several established solvers and variants. The framework proves exponential convergence and exact expected-iterate evolution, while experiments report competitive Gaussian methods on sparse problems and Gauss-pd across tested problems.

Problem
Linear systems require efficient iterative solvers for large-scale applications, where direct methods are generally not competitive.
Method
A randomized iterative framework varies a positive definite geometry matrix B and an i.i.d. random matrix S, yielding six equivalent formulations and multiple known and new algorithms.
Results
The analysis proves exponential convergence and exactly characterizes expected iterates; experiments find new Gaussian methods competitive on sparse problems and Gauss-pd competitive with CD-pd on all tests.
Takeaways & Limitations
The framework provides a common basis for analyzing existing randomized linear solvers and designing block, Gaussian, importance-sampled, and structure-tailored variants.
Takeaways & Limitations
The results leave the non-full-column-rank case and accelerated convergence through preconditioning or optimized sampling for future work.

Abstract

from arXiv · show

We develop a novel, fundamental and surprisingly simple randomized iterative method for solving consistent linear systems. Our method has six different but equivalent interpretations: sketch-and-project, constrain-and-approximate, random intersect, random linear solve, random update and random fixed point. By varying its two parameters$-$a positive definite matrix (defining geometry), and a random matrix (sampled in an independently and identically distributed fashion in each iteration)$-$we recover a comprehensive array of well-known algorithms as special cases, including the randomized Kaczmarz method, randomized Newton method, randomized coordinate descent method and random Gaussian pursuit. We naturally also obtain variants of all these methods using blocks and importance sampling. However, our method allows for a much wider selection of these two parameters, which leads to a number of new specific methods. We prove exponential convergence of the expected norm of the error in a single theorem, from which existing complexity results for known variants can be obtained. However, we also give an exact formula for the evolution of the expected iterates, which allows us to give lower bounds on the convergence rate.

1. Introduction.

The paper introduces a simple randomized framework for consistent linear systems, motivated by broad applications and the limits of direct methods at large scale. Its flexible geometry and sampling recover established algorithms, support new variants, and yield unified convergence and complexity analyses.

Motivation: Linear systems underpin applications across optimization, machine learning, signal processing, engineering, and numerical computing.Large-scale optimization can use fast approximate linear-system solvers within inexact Newton-type methods.
Motivation: Randomized iterative methods are presented as an alternative to direct methods, which are generally uncompetitive in large-scale settings.The introduction contrasts deterministic classical iterations with randomized approaches that can be competitive or better in many situations.
Contributions: The paper develops a novel randomized iterative method with six equivalent views: sketch-and-project, constrain-and-approximate, random linear solve, random intersection, random update, and random fixed point.These formulations describe projection, affine constraints, simpler random systems, geometric intersection, closed-form updates, and fixed-point iteration.
Contributions: Two parameters—a positive definite geometry matrix B and an i.i.d. random matrix S—recover randomized Kaczmarz, coordinate descent, random pursuit, and stochastic Newton methods.The same framework also permits block methods, importance sampling, Gaussian variants, and other distributions tailored to problem structure.
Convergence analysis: For full-column-rank systems, invertibility of E[Z] implies ρ < 1 and convergence, while the analysis provides expected-error bounds and an exact evolution of expected iterates.The paper also derives iteration-complexity bounds and recovers known rates in special cases.
Future directions: The framework leaves extensions to rank-deficient systems, structure-specific stochasticity, preconditioning, and large-matrix inversion as future research directions.Theorem 5.2 is identified as suggesting a preconditioning strategy for higher-precision solutions on ill-conditioned problems.

2. Optimization Viewpoint: Constrain-and-Approximate.

The method constrains each iterate to a random affine space through the current point and selects the point that best approximates a solution. Equivalent geometric, linear-system, and update formulations characterize the same next iterate.

Constrain-and-Approximate: The next iterate is chosen from a random affine space containing xk to best approximate x∗.This constrain-and-approximate formulation does not depend on which solution x∗ is selected.
Geometric interpretation: The two defining affine spaces intersect at a unique point because they are orthogonal complements under the B-inner product.The construction therefore determines xk+1 uniquely.
Algebraic viewpoint: xk+1 can be characterized as the unique solution of a linear system involving x and y.The system is equivalent to the intersection formulation.
Random Linear Solve: The framework reduces solving the original linear system to a sequence of simpler random systems.This reduction is the central computational consequence of the algebraic formulation.
Random Update: Eliminating x yields a sketched system for y, whose minimum-norm solution gives an explicit update from xk to xk+1.The update uses the Moore–Penrose pseudoinverse of STAB−1AT S.

6. Analytic viewpoint: Random Fixed Point.

The random fixed-point viewpoint expresses the iteration through projection matrices induced by Z and the B-inner-product geometry. This structure supports contraction and convergence analysis, while encompassing the other equivalent formulations.

Equivalent formulations: The six viewpoints are formally equivalent and produce the same unique point xk+1.The equivalence links the method’s geometric, optimization, algebraic, and fixed-point descriptions.
Projection structure: B−1Z is an orthogonal projection onto a d-dimensional subspace in the B-inner-product geometry.Its complement I−B−1Z is also a projection matrix.
Convergence mechanism: The complementary operator I−B−1Z is a B-norm contraction and can become a strict contraction in expectation under reasonable assumptions on S.This provides the basis for convergence of the random fixed-point iteration.
Parameter flexibility: The framework recovers existing methods by selecting the positive definite geometry matrix B and random sampling matrix S.The cited discussion presents this parameterization as a flexible unifying framework for known algorithms.

3. Special Cases: Examples.

Specific choices of B and S recover randomized Kaczmarz, coordinate descent, randomized Newton, and block variants. Sampling distributions and subspaces determine their updates and convergence behavior.

General sketch: An invertible square sketch recovers the exact solution in one step, independently of B.The corresponding convergence factor is ρ = 0.
Randomized Kaczmarz: Choosing a random coordinate sketch with B = I yields randomized Kaczmarz, including importance sampling by row magnitude.The row-proportional distribution is pi = ∥Ai:∥2/∥A∥2F.
Randomized Kaczmarz: For randomized Kaczmarz, the norm of the expected error has a squared convergence rate, although this is weaker than expected-norm convergence.The paper notes that the stronger expected-norm guarantee also controls variance.
Coordinate descent: With appropriate positive-definite choices, the framework recovers randomized coordinate descent for positive-definite and least-squares problems.The methods perform exact minimization along a selected coordinate or coordinate subspace.
Block methods: Random subsets produce block Kaczmarz and randomized Newton variants, with convergence rates depending on the sampled block structure.The framework supports overlapping partitions and higher-dimensional random coordinate subspaces.

4. Convergence: General Theory.

The convergence theory analyzes both expected iterates and expected error norms through a common rate ρ. It establishes exponential convergence under positive-definiteness and relates faster possible convergence to larger search spaces.

Convergence measures: Convergence of the expected norm of the error is stronger than convergence of the norm of the expected error because it also controls coordinate variances.The stronger notion guarantees that Var(xi) converges to zero for each coordinate.
Rate bounds: The convergence factor satisfies 0 ≤ ρ ≤ 1, with a lower bound that decreases as the search-space dimension increases.Larger search spaces therefore permit faster convergence according to the bound.
Convergence conditions: If E[Z] is invertible, then ρ < 1, A has full column rank, and the consistent system has a unique solution.Invertibility of the expected projection-related matrix supplies the key condition for strict convergence.
Expected iterates: Theorem 4.4 exactly characterizes the expected iterates through a linear fixed-point iteration with rate ρ.The spectral radius and induced B-norm of I−B−1E[Z] both equal ρ.
Expected error norm: Theorem 4.6 gives exponential convergence of the expected norm of the error when E[Z] is positive definite.The paper distinguishes this stronger convergence guarantee from the norm-of-expectation result.
Discrete sampling: For discrete sampling, the framework can optimize probabilities and express convergence rates using the scaled condition number.This supplies a unified route to complexity results for the specialized methods.

5. Methods Based on Discrete Sampling.

For complete discrete samplings, the framework guarantees convergence and supports probability choices that optimize or simplify convergence-rate bounds across many methods.

Sampling assumptions: Complete discrete sampling uses finitely many sample matrices S_i with positive probabilities, full row-rank S_i^T A, and collectively full row-rank A^T S.These conditions define the sampling assumption used for the discrete analysis.
Convergence: Complete discrete sampling replaces the pseudoinverse with an inverse and guarantees convergence of the resulting method.The guarantee follows because the sampled systems have full row rank.
Convergence: E[Z] is positive definite under complete discrete sampling, enabling the general convergence theorems to apply.Positive definiteness is established propositionally and then used to invoke the convergence results.
Practical limitation: A practical implementation of the optimal-distribution approach requires an efficiently computable approximate solution to the semidefinite program.Solving the program may be harder than solving the original linear system, despite potential convergence benefits.
Scope: The resulting framework covers randomized Kaczmarz, randomized coordinate descent, their block variants, and additional algorithms defined by other sample-matrix sets.The authors note that an earlier work had characterized an optimal distribution for randomized Kaczmarz using an SDP.

Then the iterates (2.7) satisfy

The discrete-sampling theory yields interpretable convergence bounds, recovers known rates for several methods, and clarifies when block methods require sharper analysis.

Rate interpretation: The scaled condition-number rate specializes to equality for column-vector sampling but leaves more slack for block methods.For block methods, the exact expression for λ_min(D2) is needed to assess convergence accurately.
Block methods: The generic rate bound does not reveal an advantage for block methods, despite numerical experiments showing such an advantage.Matrix paving can be used to explore the sharper block convergence rate, although the resulting bound is harder to interpret.
Special cases: Choosing B and S appropriately recovers convergence results for randomized Kaczmarz, CD-LS, and CD-pd, including single-column and block sampling.The unified theorem supplies these special-case results through parameter selection.
Preconditioning: The convergence theorem suggests preconditioning by choosing S as an approximate inverse of B^-1/2 A^T.For randomized Kaczmarz with B = I, this means sampling rows of a preconditioner of A.

6. Methods Based on Gaussian Sampling.

Gaussian sampling produces randomized line-search and block variants within the same framework, with exponential convergence but higher per-iteration cost than discrete methods.

Gaussian sampling: Gaussian sampling draws S from a zero-mean multivariate normal distribution with positive definite covariance Σ.The transformed random matrix is Gaussian with covariance Ω = B^-1/2 A^T Σ A B^-1/2 when A has full column rank.
Computational trade-off: Dense Gaussian vectors substantially increase iteration cost because each iteration requires a matrix–vector product.Numerical tests indicate that faster convergence can nevertheless offset the higher iteration cost.
Gaussian Kaczmarz: Gaussian Kaczmarz forms a search direction A^Tη and performs exact line search along it to minimize Euclidean distance from the solution.This is the Gaussian analogue of randomized Kaczmarz in the discrete setting.
Gauss-LS: Gauss-LS takes a random Gaussian direction and performs exact line search to minimize the least-squares error, matching Random Pursuit on the least-squares function.The construction uses B = A^T A and covariance Σ = A A^T.
Gauss-pd: For positive definite A, Gauss-pd takes an exact line-search step in a random Gaussian direction and is equivalent to Random Pursuit on f(x).The method is introduced using B = A and S drawn from a standard Gaussian distribution.
Block methods: Block Gaussian methods use matrix-valued Gaussian samples; Block Gauss-pd updates with S(S^TAS)^-1S^T(Ax^k-b).The block construction uses i.i.d. Gaussian entries and B = A.
Experiments: In experiments, Gauss-pd and CD-pd had similar wall-clock performance on both ill-conditioned dense and sparse positive definite problems.These tests used relative error stopping criteria and a 300-second time limit.

7. Numerical Experiments.

The experiments compare randomized linear-system and ridge-regression methods across synthetic, benchmark, sparse, ill-conditioned, and block settings. Performance depends on matrix sparsity, conditioning, probability selection, and iteration cost.

Overdetermined linear systems: Gaussian methods perform poorly on dense synthetic systems because their iteration cost is high, but are more efficient on sparse systems while converging at similar rates to discrete methods.The comparison uses rand and sprandn systems with (m; n) = (1000, 500).
Overdetermined linear systems: The four methods have similar performance on well1033, while Gauss-LS and CD-LS converge faster than Gauss-Kaczmarz and Kaczmarz on illc1033.The systems come from Matrix Market and do not have full column rank, so Theorem 4.6 does not apply.
Overdetermined linear systems: Gaussian methods have wall-clock performance comparable to discrete methods on sparse problems despite requiring more flops because of expensive matrix-vector products.MATLAB automatic multithreading partly explains differences between flop counts and wall-clock time.
Bound for Gaussian convergence: The theoretical convergence bound is tight for a condition number of 1.94 but has substantially more slack when the condition number is 41.70.The comparison uses empirical convergence, the theoretical bound, and 95% and 5% quantiles from 100 runs.
Positive Definite: Block methods can substantially improve convergence: Block CD-pd is faster on both Matrix Market problems, while Block Gauss-pd converges in one iteration on covtype.binary.The bcsstk18 problem has condition number 4.3 · 10^4, and all methods converge slowly there.
Comparison between Optimized and Convenient probabilities: Optimized sampling probabilities usually accelerate CD-pd, saving 55 seconds on aloi after 7.401 seconds of probability computation, though gains vary by problem.Mushrooms is insensitive to the probabilities, while covtype.binary initially slows before CD-pd-opt catches up as CD-pd stagnates.
Conclusions from experiments: The framework unifies randomized Kaczmarz, randomized Newton, randomized coordinate descent, and random Gaussian pursuit under one analysis, including block variants.The paper also derives a new lower bound and recovers known discrete convergence rates.

8. Conclusion.

The conclusion emphasizes preconditioning, flexible stochastic design, and block methods as routes to improved randomized linear-system solvers. Experiments support competitiveness on sparse problems, while ill-conditioning and computationally difficult optimization remain boundaries.

Conclusion: Preconditioning is important for reaching higher precision on ill-conditioned problems.The experiments show randomized methods struggling to reach 10^-2 relative error when the matrix is ill-conditioned.
Conclusion: The framework supports new block and Gaussian variants and allows probability distributions to be tailored to particular problem classes.The authors mention extensions to other continuous and discrete distributions.
Conclusion: New Gaussian methods are competitive with Kaczmarz and CD-LS on sparse overdetermined problems, while Gauss-pd is competitive with CD-pd across tests.When applicable, Block CD-pd combined direct solving with iteration most efficiently.
Conclusion: Future work includes preconditioning strategies and approximate optimized probability distributions to accelerate convergence.The paper also identifies extensions to rank-deficient systems and randomized methods for computing large matrix inverses.

Appendix A. A Bound on the Expected Gaussian Projection Matrix.

The appendix derives properties of expected Gaussian projection matrices using rotational transformations, Gaussian symmetry, covariance calculations, and integration. These results support the paper’s Gaussian randomized methods and their experimental variants.

Expected Gaussian Projection Matrix: For a Gaussian vector with covariance Ω, the appendix analyzes the normalized direction S(ξ) = ξ/∥ξ∥2 and its expected projection matrix.The notation sets S(ξ) = 0 when ξ = 0.
Expected Gaussian Projection Matrix: Orthogonal transformations reduce the Gaussian projection analysis to diagonal covariance matrices while preserving the normalized-vector relation.The proof uses ξ = Uu and derives US(u) = S(ξ).
Expected Gaussian Projection Matrix: The proof uses the central symmetry of the Gaussian distribution to establish E[S(ξ)] = 0.This identity is used alongside the orthogonal-coordinate transformation in the expected projection argument.
Expected Gaussian Projection Matrix: For diagonal covariance, independence and central symmetry make the expected projection matrix diagonal.The diagonal entries are then computed directly and related to the covariance parameters.
Two-dimensional calculation: The appendix also invokes Jensen’s inequality for the convex map (x, y) → x^2/y in bounding the relevant expression.The calculation concludes after evaluating the singular integral limits around θ = π/2.
Two-dimensional calculation: The two-dimensional calculation handles unequal variances by evaluating covariance integrals in polar coordinates and applying a tangent substitution with partial fractions.The derivation separately treats equal and unequal standard deviations.