Research Radar — 2026-11-08

Methodology Review: Randomized Kaczmarz: solve a large system one equation at a time

English · About 10 minutes

Primary author full texts and publisher metadata were inspected. Verified applications here are numerical signal reconstruction and linear-system benchmarks, not omics validation. The original worked examples were independently calculated; published experiments were not reproduced.

The intuition

Imagine solving a huge collection of linear equations when looking at one equation is cheap but processing the whole matrix is expensive. Each equation describes a flat surface containing its acceptable answers. Rather than solve everything at once, move the current guess to the nearest point satisfying one equation, then repeat with another. The useful surprise is that a carefully randomized choice can yield a clean convergence guarantee. Randomness organizes computational work; it does not create new measurements or quantify biological uncertainty. [R1]

This is randomized Kaczmarz, a row-action method. Its distinctive object is a sequence of projections onto individual equations, not a low-rank approximation to the matrix. Think of repeatedly correcting a map using one reliable landmark at a time. A correction can spoil agreement with a previous landmark, yet still bring the overall position closer to a location satisfying all of them. That last qualification matters: incompatible measurements need a different target, not merely more iterations. [R2] [R3]

A small worked example

For an original two-dimensional example, let the unknown coordinates be \(u\) and \(v\). Require \(u=1\), \(u+v=3\), and \(v=2\). Their common solution is \((1,2)\). Start at \((0,0)\). The three row vectors are \((1,0)\), \((1,1)\), and \((0,1)\), with squared lengths one, two, and one. Sampling in proportion to these squared lengths assigns probabilities one quarter, one half, and one quarter. The following deliberately chosen sequence illustrates possible draws; it is not a simulated average or a guaranteed trajectory.

First select the middle equation. Its current left side is zero, so the shortfall is three. Moving equally in both coordinates by one and a half gives \((1.5,1.5)\), the closest point on the line \(u+v=3\). Next select the first equation, giving \((1,1.5)\). Select the middle equation again: the shortfall is now one half, so add one quarter to each coordinate, reaching \((1.25,1.75)\). Finally select the third equation and obtain \((1.25,2)\). Each update satisfies the selected equation exactly.

The squared distance to \((1,2)\) falls from five initially to one half, one quarter, one eighth, and one sixteenth. Notice the second move breaks the middle equation even while improving the answer. Thus checking only the last selected equation would produce a misleading stopping rule. Conversely, selecting an already satisfied equation makes no progress. Try reversing the first two draws, then compute the distances again. This separates the geometric operation, which is deterministic once a row is chosen, from the random schedule that decides which correction comes next.

A useful counterexample uses just one unknown, \(w\), with incompatible equations \(w=0\) and \(w=2\). Every exact row projection sets the estimate to zero or two. The least-squares compromise is one, because it minimizes the sum of the two squared residuals. Ordinary row projection never lands there after a step. This hand-checkable example explains why the consistent-system theorem must not be silently rebranded as an exact least-squares theorem. [R3]

The necessary mathematics

Use real quantities. Let \(A\) have \(m\) rows and \(n\) columns, with \(m\geq n\); let \(b\) contain the \(m\) right-hand sides. Write row \(i\) as \(a_i^\top\), where \(a_i\) is a column vector, \(i\) ranges from one to \(m\), and the superscript denotes transpose. The scalar \(b_i\) is entry \(i\) of \(b\). The estimate after \(k\) projections is \(x_k\), with \(x_0\) the starting vector and \(k\) a nonnegative integer. The Euclidean norm \(\lVert a_i\rVert_2\) is its ordinary length. Define \(\lVert A\rVert_F^2\) as the sum of squared entries of \(A\). Independently draw row \(i\) with probability \(p_i\) and update:

Zero rows are excluded from sampling. The numerator is the selected equation's residual; division by squared row length gives precisely the orthogonal projection. [R1] [R4]

Here is the geometric accounting behind the result. Suppose \(x_*\) satisfies every equation and define the error \(e_k=x_k-x_*\). The correction is perpendicular to the selected hyperplane, so the Pythagorean identity subtracts the squared length of that correction. Averaging over the next row gives

Here \(\mathbb{E}\) means expectation over algorithmic randomness, and the conditioning bar means the current error is fixed. The row probabilities cancel the projection denominators. This identity explains the sampling rule without claiming it is always optimal. [R1] [R4]

Assume \(A\) has independent columns, so its smallest singular value \(\sigma_{\min}(A)\) is positive. This value is the smallest factor by which \(A\) stretches any unit-length vector. It implies \(\lVert Ae_k\rVert_2\geq\sigma_{\min}(A)\lVert e_k\rVert_2\). Iterating the previous identity yields

For our three-row example the squared Frobenius norm is four and the squared smallest singular value is one. The bound is therefore five times three quarters raised to the number of steps, much looser than our chosen trajectory. It bounds expected squared error, not every run by that same factor. [R1] [R3]

Assumptions and failure modes

The clean guarantee assumes a consistent fixed system, independent draws with the stated probabilities, full column rank, and exact arithmetic. Nearly dependent columns make the smallest singular value small and the bound slow. Needell shows that additive right-hand-side noise introduces a noise-dependent error threshold. Convergence of a solver cannot remove measurement noise, repair a wrong forward model, or identify parameters that the measurements do not distinguish. A numerical tolerance and a scientifically acceptable error are separate decisions. [R1] [R2]

Scaling is part of the algorithm. Replace the first toy equation by \(100u=100\). The solution and its projection surface are unchanged, but its squared row length becomes ten thousand. It now receives probability \(10000/10003\), leaving little computational effort for the other constraints. Censor, Herman and Jiang explain this representation dependence. Row length is therefore not an intrinsic measure of scientific importance. Normalization changes sampling; in inconsistent least squares, changing equation scales also changes the relative weighting of residuals. Document both effects rather than calling rescaling harmless preprocessing. [R4]

For genuinely inconsistent data, Zouzias and Freris add a column-projection process that estimates the part of the right-hand side outside the column space. Row projections then use a progressively corrected right-hand side. Their extended algorithm converges in mean square to the minimum-norm least-squares solution, including rank-deficient cases when initialized with a zero solution estimate and a residual estimate equal to the right-hand side. That extra machinery needs column access as well as row access. Do not promise the extension when a data interface only supplies isolated rows, and do not confuse a least-squares optimum with a validated biological model. [R3]

Verified applications

Strohmer and Vershynin demonstrate reconstruction of bandlimited signals from nonuniformly spaced samples using a finite trigonometric model. This is an inspected numerical signal-reconstruction application, not an omics result. Its weighting and representation matter, as the related scaling critique emphasizes. The evidence supports an example of useful numerical behavior, not a universal claim that randomized Kaczmarz beats every competing solver on every matrix. [R1] [R4]

Needell tests noisy systems with Gaussian, partial Fourier, and Bernoulli matrices; these experiments probe the error floor rather than validate a biological assay. Zouzias and Freris benchmark their extension on numerical least-squares problems and identify favorable sparse, well-conditioned, strongly rectangular cases, alongside limitations on ill-conditioned tests. These are computational evaluations. No omics benchmark was established by the four readings inspected here, and none of their runtime results was reproduced for this column. [R2] [R3]

Possible transfers — speculative

A speculative omics experiment would use a fixed, calibrated linear spectral-unmixing model with many measurement channels and a smaller set of unknown component abundances. Begin with synthetic mixtures whose true abundances are known. Compare a trusted least-squares baseline with the appropriate Kaczmarz variant at equal residual accuracy, including time spent sampling and accessing data. Repeat across random seeds, noise levels, row scalings, and nearly overlapping component spectra. Measure abundance error as well as runtime. This proposal is not an application demonstrated by the cited papers.

Nonnegative abundances, count noise, calibration uncertainty, and missing channels would each require explicit treatment before a real analysis. Simply clipping negative coordinates after every projection changes the algorithm and forfeits the specific theorem stated above. A sensible first question is whether the linear inverse model is scientifically defensible; only then ask whether row-wise computation is useful. Preserve an exact small reference problem so that apparent speed does not hide convergence toward the wrong target. These are proposed validation checks, not a recommendation to replace an existing assay pipeline.

What to try and what to read next

Try the three-equation example before using a large matrix. Record the selected row, the full residual, and the distance to the known answer after every step. Then repeat the contradictory one-variable example and the rescaling experiment. Those three exercises expose the central choices: what target exists, how a correction is made, and which equations receive attention. Read the primary paper for the guarantee, the noise paper for its practical boundary, the extension for least squares, and the scaling note before assuming randomness makes representation irrelevant. [R1] [R2] [R3] [R4]

Reading list

• [R1] Primary reading: A Randomized Kaczmarz Algorithm with Exponential Convergence. Thomas Strohmer and Roman Vershynin (2008). Read Algorithm 1 and Theorem 2, then the signal-reconstruction experiment. Online in 2008, journal issue in 2009; the open author manuscript is arXiv:math/0702226v1. The 2006 conference precursor is grouped with this reading for overlap checks.

• [R2] Related reading: Randomized Kaczmarz solver for noisy linear systems. Deanna Needell (2010). Use this paper to separate convergence toward a noise-dependent error floor from recovery of an exact noiseless target. Open author version: arXiv:0902.0958v2.

• [R3] Related reading: Randomized Extended Kaczmarz for Solving Least Squares. Anastasios Zouzias and Nikolaos M. Freris (2013). Read section 4 for why inconsistent equations need more than repeated exact row projections. The open author version is arXiv:1205.5770v3; the journal appeared in June 2013.

• [R4] Related reading: A Note on the Behavior of the Randomized Kaczmarz Algorithm of Strohmer and Vershynin. Yair Censor, Gabor T. Herman and Ming Jiang (2009). Read this short caution before treating row length as intrinsic scientific importance. Equivalent equations can induce different sampling schedules.