The Kaczmarz Method

Imagine you have a billion equations and a billion unknowns. Gaussian elimination — the method you learned in school — would choke: storing even the matrix takes more memory than any computer has. Yet scientists solve systems like this every day. How?

One powerful answer is the Kaczmarz method, invented by the Polish mathematician Stefan Kaczmarz in 1937. The idea is almost absurdly simple: take one equation at a time and project your current guess onto it. Repeat, cycling through all the equations, and the guess bounces its way toward the solution.

Each equation $a_1 x + a_2 y = b$ describes a line in 2D (or a hyperplane in higher dimensions). Projection onto that line means moving your current point to its nearest neighbor on the line — a single step that perfectly satisfies just that one equation. Cycling through all equations, the sequence of projections zigzags, but if the system is consistent the zigzag converges to the unique solution.

This is not just a curiosity. The Kaczmarz method — especially its modern randomized variant — is the engine behind CT-scan reconstruction, large-scale machine learning, and signal recovery. Its magic is that it never needs to see all the equations at once.

The canvas below shows a 2D linear system — two lines whose intersection is the solution. The algorithm starts from a random point and alternately projects onto each line, bouncing between them.

<!-- {{c_html_intro}} -->
<div class="controls">
  <button id="btn-step" type="button">{{btn_step}}</button>
  <button id="btn-run" type="button">{{btn_run}}</button>
  <button id="btn-reset" type="button" class="ghost">{{btn_reset}}</button>
  <span class="iter-label">{{label_iter}} <span id="iter-count">0</span></span>
</div>
<canvas id="canvas" width="460" height="340"></canvas>
<div id="status" class="status"></div>

/* {{c_css_intro}} */
* { box-sizing: border-box; }
body { font-family: system-ui, sans-serif; margin: 0; color: #222; }
.controls { display: flex; align-items: center; gap: .5rem; flex-wrap: wrap; margin-bottom: .5rem; }
button { font: 600 13px system-ui; padding: .4rem .85rem; border: 1px solid #1d3557;
         background: #1d3557; color: #fff; border-radius: 7px; cursor: pointer; }
button.ghost { background: #fff; color: #1d3557; }
button:disabled { opacity: .45; cursor: default; }
canvas { display: block; border: 1px solid #d0d8e0; border-radius: 8px; max-width: 100%; }
.iter-label { font-size: .85rem; color: #555; margin-left: .3rem; }
.status { font-size: .9rem; font-weight: 600; margin-top: .4rem; min-height: 1.3em; color: #0a7d33; }

// Code not found

Notice how each projection lands exactly on one of the lines, satisfying that equation perfectly — but potentially violating the other. Cycling back and forth, the iterates spiral inward. For two equations the path is a classic zigzag; add more equations and the bouncing becomes a spiral that still converges to the same point.

Why does zigzagging between hyperplanes actually work?

Cyclic Kaczmarz. The classical version cycles through the equations in order. Each projection strictly reduces the distance to the solution (when the system is consistent), so the iterates converge. The rate depends on the geometry — how "nearly parallel" the hyperplanes are. If two equations are almost identical, the zigzag is nearly perpendicular and converges quickly; if they are nearly parallel, the angle is shallow and convergence crawls.

Randomized Kaczmarz (Strohmer & Vershynin, 2009). A breakthrough: pick each equation at random, with probability proportional to the squared norm of its row $\|a_i\|^2$. This single change gives an exponential convergence rate in expectation:

$$\mathbb{E}\,\|x^{(k)} - x^*\|^2 \;\le\; \left(1 - \frac{\sigma_{\min}^2}{\|A\|_F^2}\right)^k \|x^{(0)} - x^*\|^2$$

where $\sigma_{\min}$ is the smallest singular value of $A$ and $\|A\|_F$ is its Frobenius norm. The factor in the exponent is a condition-number-like quantity: the better-conditioned the system, the faster the convergence.

Per-step cost is just $O(n)$ — a dot product and an update. No matrix factorization, no matrix-vector products involving the full matrix. For a system with $m$ equations and $n$ unknowns, one pass through all equations costs $O(mn)$, the same as a single matrix-vector multiply. This is exponentially cheaper than Gaussian elimination's $O(n^3)$ for the sparse, streaming settings where Kaczmarz thrives.

The method is related to dimensionality reduction and shares the spirit of gradient descent: both are iterative algorithms that make cheap local steps and trust the geometry to guide them to a global answer.

"Solve a huge linear system without ever loading it all into memory" is the exact requirement of several critical real-world problems:

• CT and MRI reconstruction: a medical scanner measures line integrals through the body — each measurement is one linear equation. The image (millions of unknowns) is reconstructed by the algebraic reconstruction technique (ART), which is the Kaczmarz method. Kaczmarz himself could not have foreseen that his 1937 algorithm would one day produce medical images.
• Compressed sensing and signal recovery: recovering a sparse signal from random linear measurements is a linear system where the equations arrive one at a time. Randomized Kaczmarz is provably efficient here.
• Machine learning: stochastic gradient descent on a squared loss is mathematically equivalent to randomized Kaczmarz. Every time you train a neural network with mini-batches, you are running a variant of the 1937 algorithm.
• Radio astronomy: aperture synthesis telescopes (like the Event Horizon Telescope that photographed a black hole) reconstruct images from sparse Fourier measurements — again a linear system solved iteratively.
• Large-scale data fitting: systems with billions of observations and thousands of features are standard in recommendation engines and genomics; Kaczmarz-style passes are the practical tool.

The method's elegance is its frugality: no more data in RAM than one equation at a time, yet guaranteed progress toward the answer with every step.

Stefan Kaczmarz published his projection method in 1937, was killed in the Second World War two years later, and never knew that his algorithm would become a cornerstone of medical imaging, signal processing, and deep learning.

The idea is almost embarrassingly simple: satisfy one equation, then the next, then the next, and repeat. Yet the geometry of projection guarantees convergence, the randomized variant provides exponential speed, and the per-step cost is so low that the method scales to systems that no other technique can touch.

The next time you look at a CT scan or train a neural network, remember: behind the complexity lies a sequence of right angles, bouncing patiently toward the truth — one equation at a time. For more on the iterative spirit, see dimensionality reduction.