Rearranging into a fixed point
The previous guide argued WHY we iterate on a huge sparse A x = b instead of factoring it: factoring fills in the zeros and chokes on memory, while iteration only ever multiplies A by a vector. Now we build the two oldest iterative methods, and they come from one childlike move. Take the i-th equation of A x = b — it ties together all the unknowns — and solve it for just the i-th unknown x_i, pretending for a moment that every OTHER unknown is already known. You get a formula for x_i in terms of the others. That formula is not the answer; it is a recipe for IMPROVING a guess, and feeding a guess through it over and over is a stationary iteration.
Make this exact by splitting the matrix. Write A = M - N, where M is some easy-to-invert piece and N is the leftover. Then A x = b becomes M x = N x + b, which suggests the loop x_{n+1} = M^{-1} (N x_n + b): plug your current guess x_n into the right, solve the easy M-system, and out comes the next guess x_{n+1}. This is a matrix splitting, and every choice of M gives a different classical method. The art is picking an M that is cheap to solve against yet close enough to A that the loop actually converges. Notice we never invert A itself — only the easy M — which is exactly the bargain iteration offers.
Jacobi: update everyone at once
The simplest splitting picks M to be just the diagonal of A. Solving a diagonal system is trivial — divide each entry by its diagonal — so this is the cheapest possible M. That choice is the Jacobi method: to get x_i's new value, take its own equation, move every other term to the right using the OLD values of the other unknowns, and divide by the diagonal a_ii. Crucially, every component is updated from the same old guess, so the order you sweep the equations does not matter and you could compute all components in parallel. It is the most democratic update imaginable: nobody hears anybody else's news until the round is over.
Solve A x = b, A = [[10, -1, 2], b = [ 6, start x = (0, 0, 0)
[-1, 11, -1], 25,
[ 2, -1, 10]] -11]
Jacobi update (all from OLD x):
x1 = ( 6 + x2 - 2 x3) / 10
x2 = (25 + x1 + x3 ) / 11
x3 = (-11 - 2 x1 + x2) / 10
iter 0: x = (0.000, 0.000, 0.000)
iter 1: x = (0.600, 2.273, -1.100)
iter 2: x = (1.047, 2.227, -0.993) ... -> (1.043, 2.269, -1.082)Gauss-Seidel: use the news as it arrives
Jacobi wastes information. By the time you are updating x_2 you have ALREADY computed a fresh, better x_1 this round — yet Jacobi insists on using the stale x_1. The Gauss-Seidel method fixes this with one tiny change: as you sweep down the equations, immediately use each new component the moment it is available. So x_2 is computed from the brand-new x_1; x_3 from the brand-new x_1 and x_2; and so on. In splitting language M is now the lower-triangular part of A (diagonal plus everything below), solved by forward substitution. Using fresher numbers almost always converges faster — for the nice systems in the example, roughly twice as fast as Jacobi, halving the iteration count.
There is a cost to that speed: because each component now leans on the ones already updated this sweep, the result depends on the sweep ORDER, and you cannot trivially parallelize a single sweep the way you can with Jacobi. Gauss-Seidel can still be made parallel with clever colourings of the unknowns, but the plain version is inherently sequential. So the two methods trade off honestly — Jacobi is embarrassingly parallel but slower per step in serial; Gauss-Seidel is faster per step but order-dependent. Neither is the universal winner, and as we will see, both are usually too slow on their own for serious problems and serve mainly as building blocks.
The spectral radius: does the loop close in?
Now the real question. Run the loop x_{n+1} = G x_n + c, where G = M^{-1} N is the iteration matrix baked into your splitting. Why should the guesses settle on the true x rather than blow up or oscillate forever? Subtract the fixed-point equation x = G x + c from the iteration and you find the error obeys e_{n+1} = G e_n. So after n steps the error is e_n = G^n e_0: the iteration matrix is applied to the starting error again and again. The whole question of convergence is therefore the question of whether the powers G^n shrink to zero — and that is governed by a single number.
That number is the spectral radius of G, written rho(G): the largest absolute value among G's eigenvalues. The clean theorem is the spectral radius convergence criterion: the stationary iteration converges for EVERY starting guess if and only if rho(G) < 1, and it diverges if rho(G) > 1. The intuition is exact — along each eigen-direction the error is multiplied by that eigenvalue every step, so the error dies only if every eigenvalue lives strictly inside the unit circle. Better still, rho(G) sets the SPEED: each step shrinks the error by roughly the factor rho(G), so rho = 0.9 crawls (you need about 22 steps to gain one decimal digit) while rho = 0.1 sprints (one digit per step).
SOR, residuals, and when to stop
If Gauss-Seidel always nudges the guess in the right direction, why not push a little FURTHER each step? That is successive over-relaxation: instead of fully replacing x_i with its Gauss-Seidel value, take a weighted blend that overshoots, x_i_new = (1 - w) x_i_old + w (Gauss-Seidel value), with a relaxation factor w between 1 and 2. Tuned to the magic best w, SOR can turn a method that needs thousands of sweeps into one that needs dozens — a genuinely dramatic win. The honest catch: the optimal w depends on the spectral radius you usually do not know in advance, a poor w can be no better than Gauss-Seidel, and w outside (0, 2) diverges outright.
How do you know when to quit? You cannot watch the error e_n = x_n - x because you do not have the true x. What you CAN compute is the residual r_n = b - A x_n, which measures how badly the current guess fails the equations. The residual is your only honest, measurable feedback, so a practical stopping criterion watches it: halt when the relative residual ||b - A x_n||_2 / ||b||_2 drops below a tolerance you set, say 1e-8. One caveat to keep you honest: a small residual does NOT guarantee a small error. On an ill-conditioned A the error can be far larger than the residual suggests — the condition number is exactly the worst-case amplification from residual to error.
Step back and the limit of this whole family comes into view. Jacobi, Gauss-Seidel, and SOR are simple, low-memory, and matrix-free-friendly, but their spectral radius typically creeps toward 1 as the problem grows — refine a grid and rho might be 1 - O(h^2), so the iteration count explodes for fine meshes. That is precisely why the next guides move on: Krylov methods like the conjugate gradient method extract far more from each matrix-times-vector product, and these old stationary iterations survive mostly as the smoothers inside multigrid. Master them anyway — they are the clearest place to meet the splitting idea and the spectral radius that every later method still answers to.