Kaczmarz's Method

In the modern world, from creating detailed medical images to training complex machine learning models, science and engineering are often faced with a monumental challenge: solving systems of linear equations with millions, or even billions, of variables. Directly inverting such a massive matrix is often computationally impossible. This presents a critical problem: how can we find the solution hidden within this vast sea of data without getting lost in the computational complexity? The answer lies not in a brute-force attack, but in an elegant and surprisingly simple iterative approach proposed by the Polish mathematician Stefan Kaczmarz. His method provides a geometric path to the solution, one step at a time.

This article explores the power and beauty of Kaczmarz's method. In the first part, ​​Principles and Mechanisms​​, we will journey into the geometry of high-dimensional spaces to understand how the method works. We will uncover the "Kaczmarz dance" of successive projections, explore the mathematics behind its guaranteed convergence, and examine modern randomized variations that have supercharged its performance. In the second part, ​​Applications and Interdisciplinary Connections​​, we will witness this abstract mathematical tool in action, revealing its crucial role in diverse fields from computed tomography and signal processing to computational chemistry, often appearing under different names but always based on the same core idea.

Imagine you are standing in a vast, featureless space. Your task is to find a hidden treasure, a single point. You don't have a map, but you receive a series of clues. The first clue is: "The treasure lies on a certain giant, flat plane." This helps, but the plane is infinite. Then you get a second clue: "The treasure also lies on another specific plane." Now things are getting interesting. The treasure must be on the line where these two planes intersect. A third clue, a third plane, will likely pinpoint the treasure's exact location—the single point where all three planes cross.

This is the essence of solving a system of linear equations like $A\mathbf{x} = \mathbf{b}$. Each equation, such as $a_1 x_1 + a_2 x_2 + \dots + a_n x_n = b$, defines a ​​hyperplane​​ in an $n$-dimensional space. The solution vector $\mathbf{x}$ is the single point that lies on all these hyperplanes simultaneously. But how do we find this point without solving the whole system at once, which can be a monumental task for systems with millions of equations, as is common in fields like medical imaging or machine learning?

This is where the genius of the Polish mathematician Stefan Kaczmarz comes in. His method offers a disarmingly simple, yet powerful, strategy: deal with the clues one at a time.

Let's say you start with a wild guess, any point $\mathbf{x}_0$ in space. This guess is almost certainly not the treasure. Now, you consider the first clue—the first hyperplane. What is the most sensible thing to do? The most efficient way to satisfy this first clue is to move from your current guess $\mathbf{x}_0$ to the closest point on that hyperplane. This move is an ​​orthogonal projection​​: you travel along a line perpendicular to the hyperplane until you hit it. This new point, $\mathbf{x}_1$, is your improved guess. It now perfectly satisfies the first equation.

Now, you take the second clue. You are at $\mathbf{x}_1$ on the first hyperplane. You repeat the process: project $\mathbf{x}_1$ orthogonally onto the second hyperplane to get $\mathbf{x}_2$. Now, $\mathbf{x}_2$ satisfies the second equation. But wait, did this move mess up our first achievement? Most likely, yes. The new point $\mathbf{x}_2$ is probably not on the first hyperplane anymore.

This is the Kaczmarz dance. You cycle through the equations, projecting your current guess onto one hyperplane after another. You take a step toward satisfying a new constraint, while partially stepping away from others. It might seem like a chaotic two-steps-forward, one-step-back dance, but here's the magic: as long as a solution exists, every single step brings you closer to it.

The mathematical formula for a single projection step is as beautiful as the idea itself: $\mathbf{x}_{k+1} = \mathbf{x}_{k} + \frac{b_i - \mathbf{a}_i \cdot \mathbf{x}_{k}}{\|\mathbf{a}_i\|^2} \mathbf{a}_i$ Let's unpack this. The vector $\mathbf{a}_i$ is the normal vector to the $i$-th hyperplane; it defines its orientation. The term in the numerator, $b_i - \mathbf{a}_i \cdot \mathbf{x}_{k}$, measures the signed distance from your current point $\mathbf{x}_k$ to the hyperplane. If you're already on it, this term is zero, and you don't move. Otherwise, it tells you how far to go. You then move in the direction of the normal vector $\mathbf{a}_i$, scaled by this distance and the length of the normal vector itself, $\|\mathbf{a}_i\|^2$. It's the most direct route to satisfying the current clue.

Why are we guaranteed to get closer to the true solution $\mathbf{x}^*$? Let's invoke Pythagoras. The true solution $\mathbf{x}^*$ lies on the hyperplane you are projecting onto. Your step, the vector $(\mathbf{x}_{k+1} - \mathbf{x}_k)$, is perpendicular to the hyperplane. Your old error $(\mathbf{x}_k - \mathbf{x}^*)$ and your new error $(\mathbf{x}_{k+1} - \mathbf{x}^*)$ form a right-angled triangle. By the Pythagorean theorem, the hypotenuse (the old error) must be longer than the new one. The distance to the solution strictly decreases with every step, unless you've already arrived. The iterates spiral inwards, converging gracefully to the treasure.

Even more remarkably, if there are infinitely many solutions (an underdetermined system), the Kaczmarz method, starting from the origin, will converge to the one with the smallest length—the ​​minimum-norm solution​​. It doesn't just find any solution; it finds the most "economical" one.

The speed of this convergence depends critically on the geometry of the hyperplanes. Imagine a city grid where all streets are perpendicular. If you are told to be on 1st Avenue and then told to be on 2nd Street, you can move along 2nd Street to the intersection without leaving 1st Avenue. You satisfy each new constraint without violating the previous ones.

This is analogous to a linear system where the rows of the matrix $A$ (the normal vectors $\mathbf{a}_i$) are mutually orthogonal. In this ideal scenario, the Kaczmarz method is astonishingly efficient: it converges to the exact solution in just ​​one single cycle​​ through the equations. Each projection lands you on a new hyperplane without moving you off the ones you've already visited. After $m$ steps, you are at the intersection of all $m$ hyperplanes. The algorithm gives you the famous least-squares solution, $A^T(AA^T)^{-1}\mathbf{b}$, revealing a deep connection between iterative methods and fundamental matrix theory.

In contrast, if two hyperplanes are nearly parallel (like two roads meeting at a very sharp angle), projecting from one to the other results in a very small step towards their distant intersection. Convergence becomes excruciatingly slow. The geometry of the problem is everything.

Real-world problems are rarely so pristine. What happens if the clues are contradictory? This is an ​​inconsistent system​​, where the hyperplanes do not share a common point. The treasure doesn't exist! Kaczmarz's method doesn't give up. It doesn't converge to a single point but instead gracefully settles into a ​​limit cycle​​, a repeating dance pattern between the hyperplanes. The set of points in this limit cycle represents a "best compromise," minimizing a weighted sum of the squared distances to all the conflicting hyperplanes. A clever trick allows us to analyze this behavior by embedding the problem into a higher-dimensional space where it becomes consistent, letting us calculate the precise rate at which it settles into this compromise.

We can also actively guide the Kaczmarz dance using a ​​relaxation parameter​​, $\omega$. The standard step is scaled by $\omega$: $\mathbf{x}_{k+1} = \mathbf{x}_{k} + \omega \left( \frac{b_i - \mathbf{a}_i \cdot \mathbf{x}_{k}}{\|\mathbf{a}_i\|^2} \mathbf{a}_i \right)$ Setting $0 \lt \omega \lt 1$ (​​under-relaxation​​) makes the algorithm take smaller, more cautious steps. Setting $1 \lt \omega \lt 2$ (​​over-relaxation​​) makes it overshoot the target hyperplane. This can be a daring and effective strategy to accelerate convergence, like taking a running leap across a chasm instead of shuffling to the edge.

However, this boldness comes with a risk. The beautiful Pythagorean guarantee of convergence only holds for $0 \lt \omega \lt 2$. If you choose $\omega \ge 2$, the method can become unstable, with the iterates flying off to infinity. The stability of the method is governed by an "error shrinkage factor" for each full cycle, called the spectral radius of the iteration operator. As long as this factor is less than 1, the errors are guaranteed to shrink over time.

The traditional Kaczmarz method cycles through the equations in a fixed order (1, 2, ..., m, 1, 2, ...). But what if you get stuck in a neighborhood of nearly-parallel hyperplanes? Your progress could grind to a halt. The modern answer is surprisingly elegant: inject chaos.

The ​​Randomized Kaczmarz (RK)​​ method doesn't follow a fixed cycle. At each step, it picks an equation uniformly at random. This simple change has profound consequences. By randomly jumping around the system, the algorithm is far less likely to get stuck in "bad" geometric configurations. On average, its performance is often dramatically better.

The beauty of this approach is that its expected convergence rate can be described by a wonderfully compact formula: $\rho = 1 - \frac{\sigma_{\min}^2(A)}{\|A\|_F^2}$ Here, $\rho$ is the factor by which the squared error is expected to shrink at each step. $\|A\|_F^2$ is the Frobenius norm, which is simply the sum of squares of all entries in the matrix—a measure of the system's total "size." The crucial term is $\sigma_{\min}(A)$, the smallest singular value of the matrix. This value is a sophisticated measure of how "well-conditioned" the matrix is—how far its rows are from being linearly dependent. A larger $\sigma_{\min}(A)$ corresponds to a more favorable geometry. The formula tells us that convergence is fastest when the system is well-conditioned relative to its size. A matrix with orthogonal rows (like in System 2 of will fare much better than one with nearly-parallel rows. Randomness, guided by the deep structure of the matrix revealed by its singular values, finds the path to the solution.

Let's return to the real world, for instance, in medical imaging. When a CT scanner takes a measurement, that data point is never perfect; it's always corrupted by noise. This means the system of equations you are trying to solve is not just massive—it's also inconsistent. Your clues are slightly wrong.

What happens when we apply Kaczmarz's method to such a noisy system? Something fascinating occurs. In the initial iterations, the algorithm makes rapid progress. It's capturing the dominant, large-scale features of the underlying "true" image, effectively ignoring the small, random noise. The iterate gets closer and closer to the true picture, $x_{\text{true}}$.

However, if you let the algorithm run for too long, it will start to "fit the noise." It will contort the solution to try to satisfy the contradictory, noisy clues as best as possible. This process will inevitably pull the iterate away from the clean, true solution and towards a noisy, artifact-ridden compromise. This behavior is called ​​semi-convergence​​.

In a striking example, after just two steps, the Kaczmarz iterate can land exactly on the true solution, only to be pulled away by a third step that tries to satisfy a noisy data point. This reveals a deep and practical insight: for noisy problems, more is not always better. Stopping the iteration early is a form of ​​iterative regularization​​. The number of steps you take is a parameter that controls the trade-off between fitting the signal and overfitting to the noise. Kaczmarz's method, therefore, is not just a solver; it's a tool for extracting truth from a sea of uncertainty, and its power lies not just in its ability to converge, but in our wisdom of knowing when to stop.

In the previous chapter, we marveled at the almost magical simplicity of Kaczmarz’s method. We saw how a point in a high-dimensional space, nudged by a series of gentle projections, unerringly finds its way to the intersection of many hyperplanes. It’s a beautiful piece of geometry. But is it just that—a curiosity for mathematicians? Or does this elegant dance of projections play out in the world around us?

In this chapter, we will see that the answer is a resounding 'yes'. We will embark on a journey through different scientific disciplines and discover Kaczmarz's method in the most unexpected places, often in disguise, solving some of the most pressing problems in science and engineering. It is a testament to the profound unity of scientific thought, where a single, powerful idea can illuminate a vast landscape of challenges.

Let's begin our journey in a place that has revolutionized modern medicine: the hospital. Inside a Computed Tomography, or CT, scanner, a patient lies still as an X-ray source and detector circle around them. The machine doesn't take a single picture, like a camera. Instead, it measures the total X-ray absorption along thousands of different lines passing through the body. Each of these measurements is a "projection sum."

Now, think about what this means. Let's say we divide the cross-section of the body into a million tiny squares, or pixels, each with an unknown density value. A single X-ray measurement, $b_i$, is the sum of the densities of all the pixels it passes through. If we write this down, we get a simple linear equation: $a_{i1}x_1 + a_{i2}x_2 + \dots + a_{iN}x_N = b_i$, where $x_j$ is the unknown density of pixel $j$, and the coefficients $a_{ij}$ are one if the ray passes through pixel $j$ and zero otherwise. This is nothing but the equation of a hyperplane in a million-dimensional space! The CT scanner provides us with thousands of such equations, one for each projection. The solution we seek—the intersection of all these hyperplanes—is the complete map of pixel densities, which is the final image.

How do we find this intersection? A direct solution is computationally unthinkable. But here, Kaczmarz’s method, known in this field as the ​​Algebraic Reconstruction Technique (ART)​​, comes to the rescue. We start with a complete guess—say, a uniform gray image. Then, we take the first projection equation and project our guess onto its corresponding hyperplane. The image changes slightly. Then we take the second equation and project the new image onto its hyperplane. We continue this process, cycling through the projections one by one. With each projection, the image is subtly corrected to better match one piece of the data. It's like a sculptor chipping away at a block of stone, first making rough cuts based on one profile, then turning the stone and refining it from another angle. Miraculously, as we sweep through the projection angles again and again, the blurry gray square sharpens, and anatomical structures emerge from the mist. The final image quality depends on many factors, such as the completeness of our projection data (using more angles, like horizontal, vertical, and diagonal projections, gives a better result) and the number of iterative sweeps we perform.

The power of sequential projection is not limited to static images. Let's move to the dynamic world of signal processing. Imagine you are on a phone call, and you hear a faint, annoying echo of your own voice. The job of an echo cancellation system is to create an "anti-echo" signal that perfectly cancels the real echo. To do this, it must first build a mathematical model of the path the echo takes—a process called "system identification." The algorithm observes a short snippet of your speech and the resulting echo, and from this, it must deduce the parameters of the filter that created that echo. This is, once again, a problem of solving a system of equations.

In the 1970s, signal processing engineers developed a clever and efficient method for this called the ​​Affine Projection Algorithm (APA)​​. For decades, it was a workhorse of the field. And here is the beautiful discovery: the APA is, at its core, a form of Kaczmarz's method!

Instead of projecting onto one hyperplane at a time, APA projects onto the intersection of a small block of hyperplanes, corresponding to the most recent data samples. This is known as the ​​Block Kaczmarz​​ method. This small change makes the algorithm more robust and efficient for streaming data. Furthermore, the real world is noisy. The Block Kaczmarz update, if applied naively, can amplify this noise. APA includes a small "regularization" term that stabilizes the process, sacrificing a tiny amount of accuracy for a huge gain in stability—a classic engineering trade-off that makes the algorithm practical in the face of noisy data.

This idea of processing data in blocks, or "mini-batches," has become fundamental to modern machine learning. In fact, a powerful modern variant is the ​​Randomized Block Kaczmarz​​ method, where instead of cycling through data in a fixed order, we pick random blocks of equations at each step. This randomization helps avoid getting stuck in slow convergence patterns and is remarkably effective for the enormous, overdetermined systems found in big data applications.

At this point, you might be thinking that any iterative method that solves a problem "piece by piece" is a form of Kaczmarz's method. This is a good intuition, but it's worth taking a moment to draw a sharp distinction, just as a physicist must distinguish between momentum and energy. Kaczmarz's method is a row-action method; its "pieces" are the rows of the system matrix $A$, which correspond to the individual equations.

There is another family of methods, called column-action methods. The most intuitive of these is ​​coordinate descent​​. Imagine you are trying to find the lowest point in a vast, hilly landscape. Instead of moving in a complex diagonal direction, you decide to only move along the north-south axis until you find the lowest point on that line. Then, you lock your north-south position and move only along the east-west axis to find the lowest point. You repeat this, cycling through the cardinal directions. This is the essence of coordinate descent.

It turns out that this coordinate descent method for the least-squares problem is mathematically identical to applying another famous iterative solver, the Gauss-Seidel method, to a related system called the "normal equations" ($A^{\top} A \mathbf{x} = A^{\top} \mathbf{b}$). So, while Kaczmarz projects onto the solution spaces of the original equations ($A\mathbf{x}=\mathbf{b}$), this other method takes steps along the coordinate axes of the variable space. They are distinct, beautiful ideas. Often, there is a very practical reason to prefer Kaczmarz's direct approach. Forming the normal equations involves computing $A^{\top} A$, an operation that can "square" the ill-conditioning of a problem, making it much harder to solve accurately. Kaczmarz, by working with $A$ directly, gracefully sidesteps this numerical pitfall.

Our final stop takes us from the macroscopic world of images and signals down to the building blocks of matter itself. In computational chemistry, scientists simulate the behavior of complex molecules like proteins by calculating the motion of every single atom over tiny fractions of a second. A crucial aspect of these simulations is maintaining the structural integrity of the molecule. The bond lengths between atoms and the angles between bonds are not allowed to change. These are rigid "holonomic constraints."

When the simulation takes a step forward in time, numerical errors can cause these bonds to slightly stretch or compress, a violation of the laws of this simulated universe. An algorithm is needed to nudge the atoms back into their correct constrained positions. One of the most famous algorithms for this is called ​​SHAKE​​. And its operation should now sound eerily familiar. SHAKE looks at the first bond constraint and moves the two atoms involved just enough to restore their correct distance. This, however, might slightly disturb other bonds connected to those atoms. So, the algorithm moves to the second constraint and corrects it. It proceeds sequentially through all constraints, and then repeats this entire sweep multiple times until all bond lengths and angles are correct to within a tiny tolerance.

This is the Kaczmarz spirit in action! It is a sequential, iterative enforcement of local constraints to achieve a global solution. Interestingly, a deeper analysis reveals that SHAKE is equivalent to the Gauss-Seidel method applied to a system of equations for the "forces" of constraint (the Lagrange multipliers), showing a deep connection to the methods we discussed in the last section. What's more, the very same questions we ask about Kaczmarz apply here: Does the order in which we process the bonds matter? Yes, it can dramatically affect the speed of convergence. Is it sometimes better to randomize the order? Yes, for complex molecules, a randomized SHAKE can be more robust than a fixed order. From medical imaging to machine learning to the simulation of life's molecules, the same core principle of iterative projection provides a simple, powerful, and universally applicable tool.

What a remarkable journey! We started with a simple geometric game of projecting a point onto lines. We then saw this very game being played out to reconstruct images from deep inside the human body. We found it in disguise, named APA, working silently to clear up the echoes in our phone calls. We saw it evolve with randomization and block processing to tackle the massive datasets of modern machine learning. We learned to distinguish it from its cousins, like coordinate descent. And finally, we witnessed it at the smallest scales, meticulously holding molecules together in a computer simulation.

The story of Kaczmarz’s method is a beautiful illustration of a deep truth in science: the most powerful ideas are often the simplest. They possess a universality that transcends the boundaries of disciplines, revealing the interconnectedness of all our scientific endeavors. It is a simple algorithm, yet it gives us the power to see the invisible, to learn from data, and to build worlds, one projection at a time.