Introduction

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 a1x+a2y=ba_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.

Try It

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.

<!-- Kaczmarz method demo — controls and canvas -->
<div class="controls">
  <button id="btn-step" type="button">Step</button>
  <button id="btn-run" type="button">Run</button>
  <button id="btn-reset" type="button" class="ghost">Reset</button>
  <span class="iter-label">Iteration: <span id="iter-count">0</span></span>
</div>
<canvas id="canvas" width="460" height="340"></canvas>
<div id="status" class="status"></div>

/* Styles for the Kaczmarz demo */
* { 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; }

const i18nFmt = (s, v) => s.replace(/%(\w+)%/g, (_, k) => String(v[k] ?? ''));
/* Kaczmarz method — project alternately onto each equation's hyperplane */

const canvas = document.getElementById('canvas');
const ctx = canvas.getContext('2d');
const W = canvas.width, H = canvas.height;

/* Coordinate transform from math space to canvas pixels */
const PAD = 50;
const XMIN = -3, XMAX = 5, YMIN = -2, YMAX = 6;
function toScreen(x, y) {
  const sx = PAD + (x - XMIN) / (XMAX - XMIN) * (W - 2 * PAD);
  const sy = H - PAD - (y - YMIN) / (YMAX - YMIN) * (H - 2 * PAD);
  return [sx, sy];
}

/* Linear system with unique solution at (1, 3):
   Line 1: 2x + y = 5  (solution: x=1, y=3)
   Line 2: x - y = -2  (solution: x=1, y=3) */
const EQS = [
  { a: 2, b: 1, c: 5 },
  { a: 1, b: -1, c: -2 }
];
const SOL = { x: 1, y: 3 };

/* Orthogonal projection of point (px,py) onto line ax+by=c */
function projectOntoLine(px, py, a, b, c) {
  const denom = a * a + b * b;
  const t = (c - a * px - b * py) / denom;
  return { x: px + a * t, y: py + b * t };
}

/* Algorithm state — current point, step counter, path history */
let cur = { x: -2, y: 5 };
let step = 0;
let path = [{ x: cur.x, y: cur.y }];
let running = false;
let timer = null;

function drawGrid() {
  ctx.clearRect(0, 0, W, H);
  /* Draw light grid lines */
  ctx.strokeStyle = '#dce4ec';
  ctx.lineWidth = 1;
  for (let xi = Math.ceil(XMIN); xi <= Math.floor(XMAX); xi++) {
    const [sx] = toScreen(xi, 0);
    ctx.beginPath(); ctx.moveTo(sx, PAD); ctx.lineTo(sx, H - PAD); ctx.stroke();
  }
  for (let yi = Math.ceil(YMIN); yi <= Math.floor(YMAX); yi++) {
    const [, sy] = toScreen(0, yi);
    ctx.beginPath(); ctx.moveTo(PAD, sy); ctx.lineTo(W - PAD, sy); ctx.stroke();
  }
  /* Draw coordinate axes */
  ctx.strokeStyle = '#aab6c2';
  ctx.lineWidth = 1.4;
  const [ox, oy] = toScreen(0, 0);
  ctx.beginPath(); ctx.moveTo(PAD, oy); ctx.lineTo(W - PAD, oy); ctx.stroke();
  ctx.beginPath(); ctx.moveTo(ox, PAD); ctx.lineTo(ox, H - PAD); ctx.stroke();
}

function drawLine(eq, color, label) {
  /* Draw the equation's line across the visible canvas region */
  const { a, b, c } = eq;
  let pts = [];
  if (Math.abs(b) > 1e-9) {
    pts.push([XMIN, (c - a * XMIN) / b]);
    pts.push([XMAX, (c - a * XMAX) / b]);
  } else {
    const xv = c / a;
    pts.push([xv, YMIN]);
    pts.push([xv, YMAX]);
  }
  ctx.strokeStyle = color;
  ctx.lineWidth = 2.2;
  ctx.beginPath();
  pts.forEach(([x, y], i) => {
    const [sx, sy] = toScreen(x, y);
    i === 0 ? ctx.moveTo(sx, sy) : ctx.lineTo(sx, sy);
  });
  ctx.stroke();
  const [lx, ly] = toScreen(pts[1][0], pts[1][1]);
  ctx.fillStyle = color;
  ctx.font = '13px system-ui';
  ctx.fillText(label, lx - 30, ly - 6);
}

function drawPath() {
  /* Draw the dashed red path connecting all iterates */
  if (path.length < 2) return;
  ctx.strokeStyle = '#e63946';
  ctx.lineWidth = 1.6;
  ctx.setLineDash([4, 3]);
  ctx.beginPath();
  path.forEach(({ x, y }, i) => {
    const [sx, sy] = toScreen(x, y);
    i === 0 ? ctx.moveTo(sx, sy) : ctx.lineTo(sx, sy);
  });
  ctx.stroke();
  ctx.setLineDash([]);
  /* Draw a dot at each iterate (larger and red for the current one) */
  path.forEach(({ x, y }, i) => {
    const [sx, sy] = toScreen(x, y);
    ctx.beginPath();
    ctx.arc(sx, sy, i === path.length - 1 ? 5 : 3, 0, 2 * Math.PI);
    ctx.fillStyle = i === 0 ? '#555' : (i === path.length - 1 ? '#e63946' : '#e6888e');
    ctx.fill();
  });
}

function drawSolution() {
  /* Draw the true solution as a green filled circle */
  const [sx, sy] = toScreen(SOL.x, SOL.y);
  ctx.beginPath();
  ctx.arc(sx, sy, 7, 0, 2 * Math.PI);
  ctx.fillStyle = '#0a7d33';
  ctx.fill();
  ctx.fillStyle = '#0a7d33';
  ctx.font = 'bold 12px system-ui';
  ctx.fillText('Solution', sx + 10, sy - 8);
}

function render() {
  drawGrid();
  drawLine(EQS[0], '#1d3557', '2x+y=5');
  drawLine(EQS[1], '#457b9d', 'x-y=-2');
  drawPath();
  drawSolution();
}

function dist() {
  const dx = cur.x - SOL.x, dy = cur.y - SOL.y;
  return Math.sqrt(dx * dx + dy * dy);
}

function doStep() {
  /* Project current point onto the next equation in round-robin order */
  const eq = EQS[step % EQS.length];
  const p = projectOntoLine(cur.x, cur.y, eq.a, eq.b, eq.c);
  cur = p;
  path.push({ x: cur.x, y: cur.y });
  step++;
  document.getElementById('iter-count').textContent = step;
  render();
  const d = dist();
  if (d < 0.01) {
    document.getElementById('status').textContent = i18nFmt('Converged in %n% iterations!', { n: step });
    stopRun();
  } else {
    document.getElementById('status').textContent = i18nFmt('Distance to solution: %d%', { d: d.toFixed(4) });
  }
}

function startRun() {
  running = true;
  document.getElementById('btn-run').textContent = 'Pause';
  document.getElementById('btn-step').disabled = true;
  timer = setInterval(() => {
    doStep();
    if (dist() < 0.01) stopRun();
  }, 260);
}

function stopRun() {
  running = false;
  clearInterval(timer);
  document.getElementById('btn-run').textContent = 'Run';
  document.getElementById('btn-step').disabled = false;
}

function resetDemo() {
  /* Return to the starting point and clear the path */
  cur = { x: -2, y: 5 };
  step = 0;
  path = [{ x: cur.x, y: cur.y }];
  stopRun();
  document.getElementById('iter-count').textContent = 0;
  document.getElementById('status').textContent = '';
  render();
}

document.getElementById('btn-step').addEventListener('click', doStep);
document.getElementById('btn-run').addEventListener('click', () => {
  running ? stopRun() : startRun();
});
document.getElementById('btn-reset').addEventListener('click', resetDemo);

render();
Result
Console

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.

The Real Complexity

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 ∥ai∥2\|a_i\|^2. This single change gives an exponential convergence rate in expectation:

E ∥x(k)−x∗∥2  ≤  (1−σmin⁡2∥A∥F2)k∥x(0)−x∗∥2\mathbb{E}\,\|x^{(k)} - x^*\|^2 \;\le\; \left(1 - \frac{\sigma_{\min}^2}{\|A\|_F^2}\right)^k \|x^{(0)} - x^*\|^2

where σmin⁡\sigma_{\min} is the smallest singular value of AA and ∥A∥F\|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)O(n) — a dot product and an update. No matrix factorization, no matrix-vector products involving the full matrix. For a system with mm equations and nn unknowns, one pass through all equations costs O(mn)O(mn), the same as a single matrix-vector multiply. This is exponentially cheaper than Gaussian elimination's O(n3)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.

Where It Matters

"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.

Conclusion

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.

Share this article

Pick a channel — or use your device's native share sheet.

Comments

Loading comments...

https://www.kipuhub.com/en/article/kaczmarz-method/Content licensed under CC BY-NC 4.0.