8  Jacobi & Gauss-Seidel Methods

Lecture 8: Introduction to Chapter 2: Iterative Methods in Matrix Algebra

Published

February 18, 2026

Overview
  1. Reminders on matrix analysis,
  2. Jacobi & Gauss-Seidel
Note

I encourage you to play around with the juptyer notebook for this lecture - you can copy the code with the this notebook button on the side of this page.

Warning

These notes are mainly a record of what we discussed and are not a substitute for attending the lectures and reading books! If anything is unclear/wrong, let me know and I will update the notes.

Tip

This lecture is mostly based on Burden, Faires, and Burden (2015) sections 7.1-7.3.

8.1 Recap: Matrix Analysis

Burden, Faires, and Burden (2015) 7.1 - 7.2

We will define \|\cdot\|_{p} to be the induced matrix norm from the vector norm |\cdot|_p: that is,

\begin{align} |x|_p := \begin{cases} \big(\sum_{i=1}^n|x_i|^p\big)^{\frac1p}, &\text{if } 1\leq p < \infty\\ % \max_{1\leq i\leq n} |x_i|, &\text{if } p= \infty. \end{cases} \end{align}

and

\begin{align} \|A\|_p := \max_{|x|_p = 1} |Ax|_p % = \max_{ x \not= 0} \frac{|Ax|_p}{|x|_p}. \end{align}

Notice that |Ax|_p \leq \|A\|_p |x|_p.

You can use LinearAlgebra to work with matrices in Julia:

using LinearAlgebra

Exercise 8.1 Compute \| A \|_2 and \| A\|_{\infty} when

\begin{align} A = \begin{pmatrix} 1 & 0 \\ -1 & 1 \end{pmatrix}. \end{align}

A = [1 0; -1 1]
2×2 Matrix{Int64}:
  1  0
 -1  1
opnorm(A), opnorm(A,Inf)
(1.6180339887498951, 2.0)
Tip

You can show

\begin{align} \| A \|_\infty = \max_{1\leq i \leq n} \sum_{j=1}^n |A_{ij}| \end{align}

Warning

In Julia, norm(A) is the Frobenius norm:

\begin{align} \| A \|_{\mathrm{F}} := \sqrt{\sum_{ij} |A_{ij}|^2} \end{align}

which is not an induced matrix norm.

Exercise 8.2 Find the eigenvalues and eigenvectors of

\begin{align} A = \begin{pmatrix} 0 & 1 \\ -1 & 0 \end{pmatrix}. \end{align}

A = [0 1 ; -1 0]
2×2 Matrix{Int64}:
  0  1
 -1  0
λ, v = eigen(A)
Eigen{ComplexF64, ComplexF64, Matrix{ComplexF64}, Vector{ComplexF64}}
values:
2-element Vector{ComplexF64}:
 0.0 - 1.0im
 0.0 + 1.0im
vectors:
2×2 Matrix{ComplexF64}:
 0.707107-0.0im       0.707107+0.0im
      0.0-0.707107im       0.0+0.707107im
det(A - λ[1]I), det(A - λ[2]I)
(0.0 + 0.0im, 0.0 + 0.0im)
(A - λ[1]I) * v[:,1], (A - λ[2]I) * v[:,2]
(ComplexF64[0.0 + 0.0im, 0.0 + 0.0im], ComplexF64[0.0 + 0.0im, 0.0 + 0.0im])

We will write \sigma(A) := \{ \lambda : \exists x \not= 0 \text{ s.t. } A x = \lambda x \} for the set of eigenvalues of A. The spectral radius \rho(A) is defined as

\begin{align} \rho(A) := \max_{\lambda \in \sigma(A)} |\lambda|. \end{align}

Lemma 8.1 \rho(A) \leq \|A\| for every submultiplicative matrix norm.

Proof.

Here, we prove the result for the induced matrix norm - see Appendix: Linear Algebra for a more general proof.

Let |\cdot| be any vector norm and \|\cdot\| the corresponding induced matrix norm. For all \lambda \in \sigma(A), there exists a normalised eigenvector x (i.e. Ax = \lambda x and |x| = 1). Therefore,

\begin{align} |\lambda| &= \left|\lambda x \right| % = \left| Ax \right| % \leq \| A \|. \end{align}

Theorem 8.1 If A is normal (i.e. A^\star A = A A^\star where A^\star is the conjugate transpose of A), then \rho(A) = \|A\|_2.

For fixed A and \epsilon > 0, there exists a submultiplicative matrix norm \| \cdot \|_\epsilon such that

\begin{align} \rho(A) \leq \| A \|_{\epsilon} \leq \rho(A) + \epsilon. \end{align}

I’ve put the proof here Appendix: Linear Algebra for completeness.

8.2 Jacobi Method

7.3 Burden, Faires, and Burden (2015)

Implementation of Jacobi:

function Jacobi(A,b,x; Niter=100, tol=1e-10) 
    
    X = [x]

    L = -tril(A, -1)
    U = -triu(A, 1)
    D = Diagonal(diag(A))

    if (any(iszero, diag(D)))
        @warn "diagonal of A is not invertible"
        return X
    end

    for i ∈ 1:Niter
        push!( X, inv(D)*( b + (L+U)*X[end] ) ) 
        r = norm(X[end]-X[end-1],Inf)/norm(X[end],Inf)
        if (r < tol)
            return X
        end 
    end

    @warn "max iterations exceeded"
    return X
end
Jacobi (generic function with 1 method)

Example 8.1 (Example 1 of Burden, Faires, and Burden (2015)) Applying the Jacobi iteration to Ax = b with

b = [6, 25, -11, 15];
A = [10 -1 2 0; 
    -1 11 -1 3; 
    2 -1 10 -1; 
    0 3 -1 8]
4×4 Matrix{Int64}:
 10  -1   2   0
 -1  11  -1   3
  2  -1  10  -1
  0   3  -1   8

gives:

X = Jacobi(A,b,[0. ,0., 0., 0.])
X[end]
4-element Vector{Float64}:
  1.0000000000277243
  1.9999999999534517
 -0.9999999999639979
  0.9999999999487015

Exercise 8.3 Verify that

\begin{align} x = \begin{pmatrix} 1 \\ 2 \\ -1 \\ 1 \end{pmatrix} \end{align}

solves Ax = b with A and b given by Example 8.1.

Here is an implementation of the Gauss-Siedel iteration:

function GaussSiedel(A,b,x; Niter=100, tol=1e-10) 
    
    X = [x]

    L = -tril(A, -1)
    U = -triu(A, 1)
    D = Diagonal(diag(A))

    if (any(iszero, diag(D)))
        @warn "diagonal of A is not invertible"
        return X
    end

    for i ∈ 1:Niter
        push!( X, inv(D-L)*( b + U*X[end] ) )
        r = norm(X[end]-X[end-1],Inf)/norm(X[end],Inf)
        if (r < tol)
            return X
        end 
    end

    @warn "max iterations exceeded"
    return X
end
GaussSiedel (generic function with 1 method)

Example 8.2 Returning to Example 8.1, we see that Gauss-Seidel converges in fewer iterations:

Y = GaussSiedel(A,b,[0.,0., 0., 0.])
Y[end]
4-element Vector{Float64}:
  0.9999999999812351
  1.9999999999852736
 -0.9999999999921585
  1.0000000000065026
println( "# iterations for Jacobi = $(length(X)) \n# iterations for Gauss-Siedel = $(length(Y))")
# iterations for Jacobi = 29 
# iterations for Gauss-Siedel = 12

In fact, we have the following error curves:

ξ = [1.,2.,-1.,1.]

errJ = zeros(length(X))
for (i,x) ∈ enumerate(X)
    errJ[i] = norm( x - ξ, Inf ) 
end

errGS = zeros(length(Y))
for (i,y) ∈ enumerate(Y)
    errGS[i] = norm( y - ξ, Inf )
end

plot( errJ, label="Jacobi", m=3,
    xlabel="Iteration (k)", yaxis=(:log10, "Error"), 
    legend=:bottomleft, color=:blue ) 

plot!( errGS, label="Gauss-Siedel", m=3, color=:red )     

L = -tril(A, -1)
U = -triu(A, 1)
D = Diagonal(diag(A))

ρJ = maximum( abs.( eigen(inv(D)*(L+U)).values ) )
ρGS = maximum( abs.( eigen(inv(D-L)*U).values ) )

plot!( 1:length(X), (errJ[1]/ρJ)ρJ.^(1:length(X)), 
    primary=false, linestyle=:dash, color=:blue )
plot!( 1:length(Y), (errGS[1]/ρGS)ρGS.^(1:length(Y)), 
    primary=false, linestyle=:dash, color=:red )

Exercise 8.4 What is the rate of approximation in the error curves of Example 8.2?

8.3 Convergence

Both Jacobi and Gauss-Siedel can be written as

\begin{align} x^{(k+1)} = T x^{(k)} + c \end{align}

Exercise 8.5 Write down what T is for (i) Jacobi, (ii) Gauss-Siedel.

TipTO-DO
  • Assignment 3: Due next Wednesday 11:59 PM,
  • Make sure you are up to date with reading from chapter 1: Burden, Faires, and Burden (2015) sections 5.1 - 5.6, Driscoll and Braun (2018) zero stability
  • This lecture: sections 7.1 - 7.3 of Burden, Faires, and Burden (2015)
  • (Sign-up to office hours if you would like)

Next

  • Read: sections 7.4 - 7.5
Burden, Richard L., Douglas J. Faires, and Annette M. Burden. 2015. Numerical Analysis. 10th ed. CENGAGE Learning.
Driscoll, Tobin A, and Richard J Braun. 2018. Fundamentals of Numerical Computation. https://fncbook.com/.