Gauss-Seidel Iteration

Solve a small diagonally dominant system by sweeping updates in place, using each new value immediately.

Problem

Gauss-Seidel is an iterative method for linear systems that updates each unknown in turn using the most recent values of the others. For diagonally dominant matrices it converges, and because it uses updated values immediately it typically beats the simpler Jacobi iteration.

Iteration

For each equation, solve for its diagonal unknown given the current values of the rest: x_i = (b_i - sum of a_ij x_j for j not equal to i) / a_ii. Sweeping through all rows once is one iteration.

python

import numpy as np
A=np.array([[4.,1.,0.],[1.,4.,1.],[0.,1.,3.]]); b=np.array([5.,6.,4.])
x=np.zeros(3)
for it in range(8):
    for i in range(3):
        s=sum(A[i,j]*x[j] for j in range(3) if j!=i)
        x[i]=(b[i]-s)/A[i,i]
    print(it,np.round(x,5))
print('residual',round(np.linalg.norm(A@x-b),3e-1*0+8))

Result

The iterates converge within a handful of sweeps because the matrix is strongly diagonally dominant, and the residual drops toward zero. Gauss-Seidel needs no matrix factorization and little storage, which is why it, and its accelerated cousin successive over-relaxation, are used as smoothers inside multigrid solvers for very large PDE systems.

• Convergence is guaranteed for diagonally dominant or symmetric positive definite matrices.
• Using updated values in place gives faster convergence than Jacobi but makes the sweep sequential.
• Gauss-Seidel smoothing is a component of the multigrid solvers used on Kronos field and transport grids.