A rate that is known in advance

The convergence bound for conjugate gradients is a bound, and the rate the condition number predicts measures how loose it is: on the model problem it permits four hundred iterations for a job that takes forty.

The stationary methods are different in a way that makes them worth an essay of their own even though nobody uses them as solvers any more. Their convergence rates are not bounds. On the model problem they are exact closed-form numbers, known before a single iteration runs, and what the iteration does can be compared against them digit for digit.

That makes this the one place in the whole field where a convergence claim is a prediction rather than a guarantee — and predictions can be wrong in a way guarantees cannot, which is what makes checking them worth the effort.

What a stationary method is #

One idea, written three ways, and writing it once makes the differences visible instead of buried in three similar loops.

Split the matrix as A = P − N, where P is something easy to invert. Then Ax = b becomes Px = Nx + b, and the iteration is x ← P⁻¹(Nx + b), or equivalently x ← x + P⁻¹r with r the current residual. Every stationary method is that, and the only difference between them is P.

Jacobi takes P to be the diagonal. Each update divides the residual entry by the diagonal entry, and every component can be computed independently — which is why the method survives on parallel hardware long after it stopped being competitive serially.

Gauss–Seidel takes P to be the lower triangle including the diagonal. The updates are sequential, each using values already updated in the current sweep, and the “solve” is a forward substitution.

SOR takes P to be D/ω + L: Gauss–Seidel with the step scaled by a relaxation parameter. At ω = 1 it is Gauss–Seidel, which the figure asserts rather than assumes — the two must produce the same iteration count at that position of the slider, and it is a check on the implementation rather than on the mathematics.

The error after one step is multiplied by the iteration matrix M = I − P⁻¹A, so after k steps it is multiplied by $M^k$, and the whole convergence question is: what does the spectral radius of M equal.

When they converge at all #

The rates above assume convergence, and convergence is not automatic. The iteration contracts exactly when the spectral radius of M is below 1, and for an arbitrary matrix it need not be — a stationary method can diverge on a perfectly nonsingular system, which is a failure mode the Krylov methods do not have.

Two conditions are worth carrying because they cover most matrices anybody meets. Strict diagonal dominance — each diagonal entry larger in magnitude than the sum of the others in its row — guarantees that both Jacobi and Gauss–Seidel converge. Symmetric positive definiteness guarantees Gauss–Seidel and SOR for any ω in (0, 2), and says nothing about Jacobi, which can diverge on a symmetric positive definite matrix.

The model problem satisfies both, which is convenient and also slightly misleading: it is the best-behaved matrix in the subject, and every number in this essay is measured on it. The closed forms are properties of that matrix rather than of the methods, and on anything else the rate has to be measured — which is what the measured-contraction route in the last section is for.

The rates, in closed form #

For the one-dimensional discrete Laplacian, all of them are known.

These are formulas, not estimates. No eigensolver, no iteration, no rounding beyond evaluating a trigonometric function — which is what makes them usable as ground truth in the way an answer that is known uses the exact Hilbert inverse.

The measured values, taken as the geometric mean of the contraction over two hundred steps, are 0.99211 for Jacobi and 0.98429 for Gauss–Seidel. Five significant figures, on both.

What the numbers mean once they are believed #

Jacobi’s rate is 0.992. That is not a slow method; it is barely a method at all. Reducing the error by ten requires log(0.1)/log(0.99211) ≈ 291 iterations, and the site’s own run does not reach a relative residual of 10⁻¹³ within three thousand.

Gauss–Seidel at 0.984 is exactly twice as fast in the sense that matters — its rate is Jacobi’s square, so it needs half as many steps for the same reduction. It converges in 1,574 iterations.

Optimally relaxed SOR converges in 128.

The last of those is not a better constant. It is a different asymptotic behaviour, and the closed forms say so plainly. As n grows, cos(π/(n+1)) approaches 1 like 1 − O(h²) where h is the mesh spacing; ω_opt − 1 approaches 1 like 1 − O(h). One power of h, and it is the difference between needing O(n²) iterations and O(n).

That is why SOR mattered historically, and it is why the figure’s slider is the relaxation parameter rather than anything else.

Why Gauss–Seidel is exactly Jacobi squared, and why that is fragile #

The relationship in the table is unusually clean — one rate is the square of the other — and it is worth knowing that this is a theorem about a class of matrices rather than a general fact.

The class is the consistently ordered matrices with a particular two-cyclic structure, which the discrete Laplacian with its natural ordering belongs to. Young’s theory gives, for such matrices, exactly the relations in the table: ρ_GS = ρ_J², the formula for ω_opt in terms of ρ_J, and ρ_SOR = ω_opt − 1 at the optimum. All three come from one quadratic relating the eigenvalues of the SOR iteration matrix to those of the Jacobi one.

The fragility is that the property depends on the ordering of the unknowns. Renumber the variables and the matrix is symmetrically permuted; its eigenvalues are unchanged, so Jacobi’s rate is unchanged — Jacobi does not care about ordering at all. But Gauss–Seidel sweeps in index order, so its iteration matrix is a different matrix under a different numbering, and the exact squaring relation can fail.

That is a real dependency and not a technicality. It is the same variable numbering that the order decides the memory shows deciding whether a factorisation fits in memory — one relabelling, two entirely different consequences, in two different fields. Red-black ordering, chosen precisely to make Gauss–Seidel parallelisable, changes the rate as well as the parallelism, and a code that adopts it for the second reason inherits the first.

The optimum is sharp, and its location is not guessable #

Dragging ω is the part worth doing slowly.

At ω = 1 the method is Gauss–Seidel, at 1,574 iterations. At ω_opt = 1.7773 it is 128. In between it falls steeply, and past the optimum it rises again — and the curve is not symmetric about its minimum, because the two sides are governed by different branches. Below ω_opt the iteration matrix has two real eigenvalues and the rate is the larger of them; above it, the two have collided into a complex conjugate pair of modulus exactly ω − 1, so the rate rises linearly and predictably.

Along the slider the SOR count reads 900, 900, 900, 793, 637, 500, 376, 257, 152, 123, 183 and 549 at ω = 1, 1.095, 1.19, 1.285, 1.38, 1.475, 1.57, 1.665, 1.76, 1.777, 1.855 and 1.95. The descent takes eight stops and the ascent takes two.

The upper branch has a closed form and it can be checked against the count. Above ω_opt the spectral radius is exactly ω − 1, so the iterations to a relative residual of 10⁻¹³ should be ln(10⁻¹³)/ln(ω − 1). At ω = 1.777, 1.855 and 1.95 that predicts 118, 191 and 583 against a measured 123, 183 and 549 — within seven per cent at all three, with the residue being the transient before the asymptotic rate takes hold.

Which makes the asymmetry a statement about what a mis-set ω costs, and it is not the shape the word ‘sharp’ suggests. Seventeen thousandths below the optimum costs 24 per cent (152 against 123); seventy-eight thousandths above costs 49 per cent (183). Per unit of ω the descent into the optimum is steeper than the climb out of it — about 1,200 iterations a unit against 770. What makes the upper side dangerous is not its slope at the optimum but its curvature: by ω = 1.95 the climb is 3,850 iterations a unit and still rising, because ln(ω − 1) is heading for zero.

So the practical rule is the opposite of the one the picture suggests at a glance. Err low. A little under the optimum is on the branch whose rate is a real eigenvalue moving smoothly; a little over is on a branch whose cost diverges as ω approaches 2, and the divergence is a logarithm rather than a polynomial.

The practical content is that ω_opt must be computed rather than tuned. Being wrong by 0.05 costs more than the entire difference between Jacobi and Gauss–Seidel, and the formula for it involves the spectral radius of the Jacobi iteration matrix — which for a problem without a closed form is not available in advance. That is the reason SOR is a beautiful method that is difficult to use, and a large part of why Krylov methods displaced it: conjugate gradients have no parameter to get wrong.

The rate that cannot be measured the way the others were #

There is a complication at the optimum, and it is the most interesting thing in this essay because it was found by an assertion failing.

Jacobi’s and Gauss–Seidel’s measured contraction factors agree with their closed forms to five figures. SOR’s, at the optimum, came back as 0.7706 against a predicted 0.77725 — off by 0.9%, which is far outside the tolerance the other two met and far too close to be a coding error.

The cause is the reason ω_opt is optimal. Below the optimum the iteration matrix has two real eigenvalues; ω_opt is exactly where they collide. At the collision the matrix is defective — it does not have a full set of eigenvectors — and a defective matrix’s powers do not decay like $\rho^k$. They decay like $k\,\rho^k$, and the residual history acquires an oscillating envelope rather than a clean geometric one.

Measured over the last thirty steps, successive contraction ratios span 0.317 to 0.874 around a rate of 0.777. A single ratio is meaningless; a geometric mean over a long window recovers the rate to within a few per cent and no better.

So the assertion was rewritten to say what is true: the rate holds in the mean, within 4%, and the oscillation is real, with the spread of successive ratios asserted to exceed 1.5. The second half is the more informative one, and it would have been averaged away by a check that only looked at the mean.

This is the same object the spectra field is about from the other direction. Symmetry is worth more than precision measures what defectiveness does to eigenvalue sensitivity; here it is what defectiveness does to the observability of a convergence rate. A defective matrix is not merely one whose eigenvalues are delicate — it is one whose powers behave differently from what its eigenvalues alone suggest.

Two routes to the spectral radius #

The site’s habit, applied to the quantity the whole essay is about.

Jacobi’s iteration matrix on this problem is symmetric — the diagonal is constant, so I − D⁻¹A is a scalar multiple of a symmetric matrix — and its spectral radius is therefore available from the site’s Jacobi eigensolver directly. That is one route. The closed form cos(π/(n+1)) is the other. They agree to 10⁻¹⁰, checked on every build.

Gauss–Seidel’s iteration matrix is not symmetric, so that route is closed, and the only available measurement is to watch the residual contract. Which is the honest situation for any problem without a closed form, and is why the measured-contraction route is the one the figure leads with even where a cheaper one exists.

There is a trap in the symmetric route worth recording. The Jacobi iteration matrix for this problem has eigenvalues cos(kπ/(n+1)) for k = 1…n, which includes both +cos(π/(n+1)) and its negative — the spectrum is symmetric about zero. So a power iteration, the obvious cheap way to find a spectral radius, does not converge here at all: it oscillates between the two extreme eigenvectors forever. The same failure, for the same reason, as the algorithm the libraries actually run exhibits on a matrix with eigenvalues ±1.

What the rate does not say #

A closed-form spectral radius is a strong thing to have and it is worth being precise about what it covers, because it is an asymptotic statement and the runs it describes are finite.

The spectral radius governs the error’s decay eventually. It says nothing about the first several steps, which can be much better or much worse — the transient behaviour is governed by the norm of the iteration matrix rather than its spectral radius, and for a non-symmetric iteration matrix those two can be far apart — which is the spectrum that predicts nothing’s whole subject. A method with ρ < 1 can have ‖M‖ > 1 and the error can grow substantially before it starts to fall.

It is tempting to set that aside as a caveat about exotic matrices. Two of the three methods here do behave as the rate says, and the third does not — on this same model problem. Taking ‖M^k‖ maximised over starting vectors rather than measured from one:

Jacobi’s iteration matrix is symmetric, so its norm is its spectral radius to five digits, and Gauss–Seidel’s differ in the third. For those two the rate describes the whole run and there is no transient to worry about.

Optimal SOR’s norm is 2.2 against a spectral radius of 0.79. The error can more than double in a single step, on the very problem the closed form is derived for, and by a factor of 2.8 between the two quantities — a matrix nobody would call exotic.

The qualification worth keeping is that the excursion is one step. The worst ‖M^k‖ is at k = 1 and falls monotonically after it, so this is a single doubling rather than the sustained growth a strongly non-normal iteration produces. A run that measures its residual at every step sees one bad step; a run that checks every tenth sees nothing at all.

And it belongs to the method whose rate is predicted best. SOR at the optimal ω is where the closed form is sharpest and where the norm is furthest from it, because the parameter that pushes every eigenvalue onto a small circle is the parameter that makes the eigenvectors most nearly parallel — so the sharpest rate in the essay comes with the largest gap between what it promises asymptotically and what the first step does.

It is the same caveat conjugate gradients carry from the other direction — there the bound is asymptotically pessimistic and the early steps beat it — and the shared moral is that a rate is a statement about the tail of a run and a user cares about its length.

assertTheModelProblemHasATransientAfterAll measures both quantities for all three methods and requires the two symmetric-ish ones to have no transient and SOR to have one, at a single step.

The second thing the rate does not say is anything about cost. ρ counts iterations, and an iteration of Jacobi and an iteration of SOR cost the same only because the model problem’s P is trivial to invert in both cases. Comparing methods with different splittings on rate alone would be comparing iteration counts of unequal price, which is the error changing the condition number on purpose warns about for preconditioners and is the same error here.

What is asserted here #

Jacobi contracts at cos(π/(n+1)) and Gauss–Seidel at its square, both to within 3·10⁻³, measured over two hundred steps — an answer known in advance, in the sense an answer that is known means it.

Jacobi’s spectral radius from the eigensolver equals the closed form to 10⁻¹⁰ — the second route.

SOR at ω = 1 is Gauss–Seidel, identical iteration counts, checked at that slider position.

Optimal SOR beats Gauss–Seidel by more than a constant, asserted as a factor of at least eight, which is the claim about the exponent rather than about the constant.

And the oscillation at the optimum is real, with successive ratios spanning more than 1.5 — the assertion added when the clean version failed.

The refusal: Jacobi’s measured rate, compared against Gauss–Seidel’s closed form, must throw. Without it, a rate check that matched anything in the right ballpark would pass on both.

Why these are still worth knowing #

They are not competitive as solvers and this essay has not tried to argue otherwise. Three reasons they remain in every course and every library.

They are the smoothers inside multigrid. A few Jacobi or Gauss–Seidel sweeps do not solve a problem, but they annihilate the high-frequency components of the error very effectively — which is exactly what a multigrid cycle needs before it coarsens. The slow convergence is entirely in the low-frequency modes, and the coarse grid handles those. The method’s defect and multigrid’s design are the same fact.

They parallelise. Jacobi’s updates are independent, which conjugate gradients’ inner products are not: a global reduction per step is a synchronisation point, and on a large enough machine that dominates. This is why a method with a rate of 0.992 has not disappeared.

And they are where the vocabulary comes from. Spectral radius, splitting, relaxation, smoothing — the whole language of iterative methods was built here, and the Krylov methods that displaced them are still described with it.

A rate that is not known and is measured instead #

The model problem’s stationary rates are closed forms. Where a method’s step count has no closed form it is worth fitting one, and an inverse-free orthogonalisation has a step count of 5.67 per decade of κ against a predicted log 10 / log 1.5 = 5.68.