🔁Linear Systems: Iterative Methods

Jacobi and Gauss–Seidel

Splitting A = M − N, the Jacobi and Gauss–Seidel iterations, spectral radius and diagonal dominance

Göttingen, 26 December 1823

Carl Friedrich Gauss is directing the geodetic survey of the Kingdom of Hanover — years of triangulation, every angle measured repeatedly, every measurement slightly wrong. Fitting the network by least squares produces linear systems far too large to eliminate comfortably by hand. On the day after Christmas he writes to his former student Christian Ludwig Gerling, who is facing the same problem in his own survey. Gauss describes his method. Do not eliminate. Start from a guess. Find the equation that is most violated, correct the one unknown that fixes it, and repeat. He recommends it warmly. It almost never goes wrong, he says, and it can be carried out half asleep, or while thinking about other things. Gauss never published it. Jacobi described a similar scheme in 1845 and Philipp Ludwig von Seidel another in 1874. A century later, the descendants of Gauss's 'indirect method' solve the largest systems in engineering — the ones no factorisation can touch.

The idea: rewrite

A\mathbf x = \mathbf b

so that each equation gives one unknown in terms of the others. Equation

i

reads

a_{ii}x_i + \sum_{j \ne i} a_{ij}x_j = b_i

, so

x_i = \dfrac{1}{a_{ii}}\Big(b_i - \sum_{j \ne i} a_{ij}x_j\Big).

If we knew the other

x_j

we would know

x_i

. We do not — so use the current guess for them, compute a new

x_i

, and repeat. Each sweep through the equations costs one pass over the nonzeros of

A

O(n)

for a sparse matrix, and no fill-in, ever.

💡Jacobi and Gauss–Seidel

Jacobi: compute every new

x_i^{(k+1)}

from the old values

x_j^{(k)}

only. All updates are independent — ideal for parallel hardware. Gauss–Seidel: sweep through

i = 1, \dots, n

and use each new value the moment it is computed:

x_i^{(k+1)} = \dfrac{1}{a_{ii}}\Big(b_i - \sum_{j<i} a_{ij}x_j^{(k+1)} - \sum_{j>i} a_{ij}x_j^{(k)}\Big).

It needs only one copy of

\mathbf x

and usually converges faster — but each update depends on the previous one, so it is inherently sequential.

Splitting. Both methods write

A = M - N

with

M

easy to solve with, and iterate

M\mathbf x^{(k+1)} = N\mathbf x^{(k)} + \mathbf b.

Write

A = D + L + U

(diagonal, strictly lower, strictly upper). Jacobi takes

M = D

; Gauss–Seidel takes

M = D + L

. Subtract the exact equation

M\mathbf x = N\mathbf x + \mathbf b

and the error

\mathbf e^{(k)} = \mathbf x^{(k)} - \mathbf x

obeys

\mathbf e^{(k+1)} = T\,\mathbf e^{(k)},\qquad T = M^{-1}N,\qquad \text{converges for every start} \iff \rho(T) < 1

ρ(T), the spectral radius, is the largest |eigenvalue| of the iteration matrix. Asymptotically the error shrinks by a factor ρ per sweep — lesson 2-3's fixed-point iteration, now in n dimensions.

Derivation

Why the largest eigenvalue decides

What we need about eigenvalues (Unit 6 shows how to compute them). A nonzero vector

\mathbf v

with

T\mathbf v = \lambda\mathbf v

is an eigenvector:

T

only stretches it, by the eigenvalue

\lambda

. Applying

T

k

times stretches it by

\lambda^k

. For the matrices of this unit,

n

independent eigenvectors exist, so any vector can be written as a combination of them.

The error after $k$ sweeps. Write the starting error as

\mathbf e^{(0)} = c_1\mathbf v_1 + \dots + c_n\mathbf v_n

. Each sweep multiplies by

T

, so

\mathbf e^{(k)} = T^k\mathbf e^{(0)} = c_1\lambda_1^k\mathbf v_1 + c_2\lambda_2^k\mathbf v_2 + \dots + c_n\lambda_n^k\mathbf v_n.

Every term goes to zero exactly when its

|\lambda_i| < 1

. If even one

|\lambda_i| \ge 1

, a starting error along

\mathbf v_i

never dies — so convergence for every start needs all of them below 1:

\rho(T) = \max|\lambda_i| < 1

. After many sweeps the term with the largest

|\lambda_i|

is the last one left, and each further sweep multiplies the error by about

\rho

. One more correct digit then costs

\ln 10/(-\ln\rho)

sweeps.

💡When are they guaranteed to converge?

If $A$ is strictly diagonally dominant — $|a_{ii}| > \sum_{j\ne i}|a_{ij}|$ in every row — both Jacobi and Gauss–Seidel converge.
If $A$ is symmetric positive definite, Gauss–Seidel converges.
For the matrices of the Laplace and heat equations (weakly dominant, irreducible), both converge, and for the standard five-point grid $\rho_{GS} = \rho_J^2$ : Gauss–Seidel needs exactly half the sweeps of Jacobi.

Without dominance there is no guarantee. For

\begin{bmatrix} 1 & 2 \\ 3 & 1 \end{bmatrix}

\rho_J = 2.45

and

\rho_{GS} = 6

: both diverge.

Derivation

Diagonal dominance makes Jacobi a contraction

Jacobi's iteration matrix is

T_J = -D^{-1}(L + U)

: row

i

holds

-a_{ij}/a_{ii}

off the diagonal and 0 on it. Its

\infty

-norm is the largest absolute row sum (lesson 3-2):

\|T_J\|_\infty = \max_i \sum_{j \ne i}\frac{|a_{ij}|}{|a_{ii}|} < 1,

and the inequality is exactly strict diagonal dominance. Then every sweep shrinks the largest error component:

\|\mathbf e^{(k+1)}\|_\infty = \|T_J\mathbf e^{(k)}\|_\infty \le \|T_J\|_\infty\,\|\mathbf e^{(k)}\|_\infty,

so the error falls at least geometrically — the contraction mapping argument of lesson 2-3, in

n

dimensions. No eigenvalues are needed. (For Gauss–Seidel the proof is a little longer but the conclusion is the same.)

📐Temperature in a plate

A square plate, top edge held at 100 °C and the other three at 0 °C. The steady temperature satisfies Laplace's equation; with the five-point stencil each interior temperature is the average of its four neighbours — and a Gauss–Seidel sweep just replaces each value by that average. By symmetry the exact centre temperature is 25 °C. Stopping when no value changes by more than

10^{-6}

, Gauss–Seidel needs 151 sweeps on a

9 \times 9

grid, 554 on

19 \times 19

and 2005 on

39 \times 39

. Halve the spacing and the work per sweep goes up 4×, and the number of sweeps goes up nearly 4×. For the model problem

\rho_J = \cos(\pi h)

, so

1 - \rho \approx \pi^2h^2/2

: each sweep removes less of the error the finer the grid.

Derivation

The eigenvalues of the second-difference matrix, and

\rho_J = \cos\pi h

Take the 1-D matrix

A = \operatorname{tridiag}(-1, 2, -1)

with

n

unknowns at

x_i = ih

h = 1/(n+1)

, and guess that its eigenvectors are sampled sine waves,

v_i = \sin(j\pi x_i)

for

j = 1, \dots, n

. Row

i

A\mathbf v

-\sin(j\pi x_{i-1}) + 2\sin(j\pi x_i) - \sin(j\pi x_{i+1}).

The identity

\sin(a - b) + \sin(a + b) = 2\sin a\cos b

, with

a = j\pi x_i

and

b = j\pi h

, turns the two outer terms into

2\sin(j\pi x_i)\cos(j\pi h)

(A\mathbf v)_i = \big(2 - 2\cos(j\pi h)\big)\,v_i = 4\sin^2\!\Big(\frac{j\pi h}{2}\Big)\,v_i.

The boundary rows work too, because the sine vanishes at

x_0 = 0

and

x_{n+1} = 1

. So the eigenvalues are

\lambda_j = 2 - 2\cos(j\pi h)

, from about

\pi^2h^2

up to almost 4.

Jacobi. The diagonal is 2, so

T_J = I - A/2

, whose eigenvalues are

1 - \lambda_j/2 = \cos(j\pi h)

. The largest in magnitude is

\cos(\pi h)

(for

j = 1

, and

-\cos\pi h

for

j = n

). On the 2-D five-point grid the eigenvectors are products of sines, the eigenvalues add, the diagonal is 4, and the same calculation gives

\rho_J = \tfrac12(\cos\pi h + \cos\pi h) = \cos\pi h

The cost. By Taylor,

\cos\pi h \approx 1 - \pi^2h^2/2

, and

-\ln\rho_J \approx \pi^2h^2/2

, so a digit costs about

2\ln 10/(\pi^2h^2)

sweeps — proportional to

1/h^2

. Halving

h

quadruples the sweeps. That the slowest mode is the smoothest one,

j = 1

, is the key to multigrid in lesson 4-4. (David Young proved that

\rho_{GS} = \rho_J^2

for this whole class of matrices; we use his result without proof.)

the unique end of science is the honor of the human mind

– Carl Gustav Jacob Jacobi, letter to Legendre (1830)

🕰Relaxation by hand

In 1910 Lewis Fry Richardson computed the stresses in a masonry dam by iterating finite-difference equations on a grid — the first serious engineering use of the idea. In the 1930s and 1940s Richard Southwell's 'relaxation method' turned it into a craft: human computers scanned a grid for the largest residual and 'relaxed' it by hand, guided by experience. Gauss's 'half asleep' was optimistic, but the method needed no matrices at all — only the stencil.

💡Where you will use this

CFD and Hydraulics: Gauss–Seidel sweeps are the 'smoothers' inside every multigrid pressure solver. Geotechnical Engineering: seepage and consolidation on large grids. Power Systems: the earliest load-flow programs (before Newton) used Gauss–Seidel iteration on the bus voltages. Heat transfer: the five-point stencil and its 3-D seven-point cousin are the first thing to try for steady conduction.

📚Further reading

Yousef Saad, Iterative Methods for Sparse Linear Systems (2nd ed., SIAM, 2003), chapter 4, 'Basic Iterative Methods' — Jacobi, Gauss–Seidel and SOR as splittings, with convergence proofs.
Randall LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations (SIAM, 2007), chapter 4, 'Iterative Methods for Sparse Linear Systems' — the model Poisson problem and its spectral radii, worked out as here.
Richard Varga, Matrix Iterative Analysis (Prentice-Hall, 1962; 2nd ed., Springer, 2000) — the classical theory, including diagonal dominance and irreducibility.

↓But there is a problem

Gauss–Seidel works, needs almost no memory and never fills in. But on a fine grid it crawls: for a

320 \times 320

grid it needs about 24,000 sweeps for each additional decimal digit. In 1950 a Harvard doctoral student, David Young, asked a simple question. Gauss–Seidel corrects each value by exactly the amount its equation demands. What happens if you correct it by more than that — deliberately overshoot? The answer turned 24,000 sweeps into about 120.