2.8. Iterative Methods for Simultaneous Linear Equations#
Last revised on September 23-24, 2025, modifying the example matrix T in (2.17), adding code and testing for a version of the Jacobi method, and revising the exercises.
import numpy as np
import numpy.linalg as la
from numpy import Inf
References:
[Sauer, 2022] Section 2.5, Iterative Methods, sub-sections 2.5.1 to 2.5.3.
[Burden et al., 2016] Chapter 7, Iterative Techniques in Linear Algebra, sections 1 to 3.
[Chenney and Kincaid, 2013] Section 8.4, Iterative Solution of Linear Systems.
[Kincaid and Chenney, 1990] Section 4.6, Solution of Equations by Iterative Methods.
[Dionne, 2023] Chapter 5, Iterative Methods to Solve Systems of Nonlinear Equations.
[Dahlquist and Björck, 2009] Section 10.1, Classical Iterative Methods.
2.8.1. Introduction#
This topic is a huge area, with lots of ongoing research; this section just explores the first few methods in the field:
The Jacobi Method.
The Gauss-Seidel Method.
The next three major topics for further study are:
The Method of Succesive Over-Relaxation (“SOR”). This is usually done as a modification of the Gauss-Seidel method, though the strategy of “over-relaxation” can also be applied to other iterative methods such as the Jacobi method.
The Conjugate Gradient Method (“CG”). This is beyond the scope of this course; I mention it because in the realm of solving linear systems that arise in the solution of differential equations, CG and SOR are the basis of many of the most modern, advanced methods.
Preconditioning.
2.8.2. The Jacobi method#
The basis of the Jacobi method for solving \(Ax = b\) is splitting \(A\) as \(D + S\) where \(D\) is the diagonal of \(A\):
so that \(S = A-D\) has
Visually
It is easy to solve \(Dx = b\): the equations are just \(a_{ii} x_i = b_i\) with solution \(x_i = b_i/a_{ii}\).
Thus we rewrite the equation \(Ax = Dx + Sx = b\) in the fixed point form
and then use the familiar fixed point iteration strategy of inserting the currect approximation at right and solving for the new approximation at left:
Note: We could make this look closer to the standard fixed-point iteration form \(x_k = g(x_{k-1})\) by dividing out \(D\) to get
but — as is often the case — it will be better to avoid matrix inverses by instead solving this easy system. This “inverse avoidance” becomes far more important when we get to the Gauss-Seidel method!
See Exercise 2.17.
2.8.3. The underlying strategy#
To analyse the Jacobi method — answering questions like for which matrices it works, and how quickly it converges — and also to improve on it, it helps to described a key strategy underlying it, which is this: approximate the matrix \(A\) by another one \(E\) one that is easier to solve with, chosen so that the discrepacy \(R = A-E\) is small enough. Thus, repeatedly solving the new easier equations \(Ex^{(k)} = b^{(k)}\) plays a similar role to repeatedly solving tangent line approximations in Newton’s method.
Of course to be of any use, \(E\) must be somewhat close to \(A\); the remainder \(R\) must be small enough. We can make this requirement precise with the use of matrix norms introduced in Error bounds for linear algebra, matrix norms, condition numbers, etc. and an upgrade of the contraction mapping theorem seen in section Solving Equations by Fixed Point Iteration (of Contraction Mappings).
Thus consider a general splitting of \(A\) as \(A = E + R\). As above, we rewrite \(Ax = Ex + Rx = b\) as \(Ex = b - Rx\) and thence as \(x = E^{-1}b - (E^{-1}R)x\). (It is alright to use the matrix inverse here, since we are not actually computing it; only using it for a theoretical argument!) The fixed point iteration form is thus
where \(c = E^{-1}b\) and \(S = E^{-1}R\).
For vector-valued functions we extend the previous Definition 1.2 in Section Solving Equations by Fixed Point Iteration (of Contraction Mappings) as:
Definition 2.12 (Vector-valued contraction mapping)
For a set \(D\) of vectors in \(\mathbb{R}^n\), a mapping \(g:D \to D\) is called a contraction or contraction mapping if there is a constant \(C < 1\) such that
for any \(x\) and \(y\) in \(D\). We then call \(C\) a contraction constant.
Next, the contraction mapping theorem Theorem 1.2 extends to
Theorem 2.13 (Contraction mapping theorem for vector-valued functions)
Any contraction mapping \(g\) on a closed, bounded set \(D \in \mathbb{R}^n\) has exactly one fixed point \(p\) in \(D\).
This can be calculated as the limit \(\displaystyle p = \lim_{k \to \infty} x^{(k)}\) of the iteration sequence given by \(x^{(k)} = g(x^{(k-1)})\) for any choice of the starting point \(x^{(0)} \in D\).
The errors decrease at a guaranteed minimum speed: \(\| x^{(k)} - p \| \leq C \| x^{(k-1)} - p \|\), so \(\| x^{(k)} - p \| \leq C^k \| x^{(0)} - p \|\).
With this, it turns out that the above iteration converges if \(S\) is “small enough” in the sense that \(\|S\| = C < 1\) — and it is enough that this works for any choice of matrix norm!
Theorem 2.14
If \(S := E^{-1}R = E^{-1}A - I\) has \(\|S\| = C < 1\) for any choice of matrix norm, then the iterative scheme \(x^{(k)} = c - S x^{(k-1)}\) with \(c = E^{-1}b\) converges to the solution of \(Ax = b\) for any choice of the initial approximation \(x^{(0)}\). (Aside: the zero vector is an obvious and popular choice for \(x^{(0)}\).)
Incidentally, since this condition guarantees that there exists a unique solution to \(Ax=b\), it also shows that \(A\) is non-singular.
Proof. (sketch)
The main idea is that for \(g(x) = c - S x\),
so with \(C < 1\), it is a contraction.
(The omitted more “technical” detail is to find a suitable bounded domain \(D\) that all the iterates x^{(k)} stay inside it.)
2.8.3.1. What does this say about the Jacobi method?#
For the Jacobi method, \(E = D\) so \(E^{-1}\) is the diagonal matrix with elements \(1/a_{i,i}\) on the main diagonal, zero elsewhere. The product \(E^{-1}A\) then multiplies each row \(i\) of \(A\) by \(1/a_{i,i}\), giving
so that subtracting the identity matrix to get \(S\) cancels the ones on the main diagonal:
Here is one of many places that using the maximum-norm, a.k.a. \(\infty\)-norm, makes life much easier! Recalling that this is given by
First, sum the absolute values of elements in each row \(i\); with the common factor \(1/|a_{i,i}|\), this gives \( \left( |a_{i,1}| + |a_{i,2}| + \cdots |a_{i,i-1}| + |a_{i,i+1}| + \cdots |a_{i,n}| \right)/|a_{i,i}| \).
Such a sum, skipping index \(j=i\), can be abbreviated as
Then get the norm as the maximum of these:
and the contraction condition \(C < 1\) becomes the requirement that each of these \(n\) “row sums” is less than 1:
Multiplying each of the inequalities by the denominator \(|a_{i,i}|\) gives \(n\) conditions
This is strict diagonal dominance, as in Definition 2.2 in the section Row Reduction (Gaussian Elimination), and as discussed there, one way to think of this is that such a matrix \(A\) is close to its main diagonal \(D\), which is the intuitive condition that the approximation of \(A\) by \(D\) as done in the Jacobi method is “good enough”.
And indeed, combining this result with Theorem 2.14 gives:
Theorem 2.15 (Convergence of the Jacobi method)
If \(A\) is strictly diagonally dominant or positive defininte, the Jacobi Method converges for any initial approximation \(x^{(0)}\).
Further, the error goes down by at least a factor of \(\| I - D^{-1}A \|\) at each iteration.
Remark 2.27
Note: convergence is not guaranteed for all positive definite matrices, but can be achieved with modification to weighted Jacobi method: see the Wikipedia page on the Jacpbi Method
By the way, other matrix norms give other conditions guaranteeing convergence; perhaps the most useful of these others is that it is also sufficient for \(A\) to be column-wise strictly diagonally dominant as in Definition 2.3.
2.8.4. The Gauss-Seidel method#
To recap, two key ingredients for a splitting \(A = E + R\) to be useful are that
the matrix \(E\) is “easy” to solve with, and
it is not too far from \(A\).
The Jacobi method choice of \(E\) being the main diagonal of \(A\) strongly emphasizes the “easy” part, but we have seen another larger class of matrices for which it is fairly quick and easy to solve \(Ex = b\): triangular matrices, which can be solved with forward or backward substitution, not needing row reduction.
The Gauss-Seidel Method takes \(E\) be the lower triangular part of \(A\), which intuitively leaves more of its entries closer to \(A\) and makes the remainder \(R = A-E\) “smaller”.
To discuss this and other splittings, we write the matrix as \(A = L + D + U\) where:
\(D\) is the diagonal of \(A\), as for Jacobi
\(L\) is the strictly lower diagonal part of \(A\) (just the elements with \(i > j\))
\(U\) is the strictly upper diagonal part of \(A\) (just the elements with \(i < j\))
That is,
Thus for the Jacobi method, \(E = D\) and \(S = L+U\).
So now we use \(E = L + D\), the lower triangular part of \(A\), which will be called \(A_L\), and the remainder is \(S = U\). The fixed point form becomes
giving the fixed point iteration
Here we definitely do not calculate the inverse of \(A_L\) explicitly; instead, we can solve with forward substitution.
However to analyse convergence, the mathematical form
is useful: the iteration map is now \(g(x) = c - G x\) with \(c = (L + D)^{-1} b\) and \(G = (L + D)^{-1} U\).
Arguing as above, we see that convergence is guaranteed if \(\| (L + D)^{-1} U\| < 1\). However it is not so easy in general to get a formula for \(\| (L + D)^{-1} U\|\) as used to prove cnvergence fot the Jacobi method.
One can get a slightly disappointing result: despite the \(S=U\) here being in some sense “smaller” than the \(E=L+U\) for the Jacobi method, the general convergence guarantee is no better:
Theorem 2.16 (Convergence of the Gasuss-Seidel method)
The Gauss-Seidel method converges if \(A\) is strictly diagonally dominant, for any initial approximation \(x^{(0)}\).
There is some good news: in practice, the convergence rate as given by \(C = C_{GS} := \| (L + D)^{-1} U\|\) is often better than for \(C = C_J := \| D^{-1} (L+U) \|\) for the Jacobi method.
Also, when one writes out the steps in detail, the two methods involve the same amount of calculation in each iteration, so a smaller contraction constant (\(C_{GS} < C_{J}\)) gives a clear “cost” advantage.
See Exercise 2.18.
2.8.5. Stopping conditions#
To turn these iterative schemes into practical algorithms that terminate in a finite number of iterations, we as usual need an error estimate, and a natural one is the backward error, or more specifically the dimension-less “relative backward error” as in Theorem 2.11.
That is, set an error tolerance \(\epsilon\), and stop when \(\displaystyle \frac{\| b - A x^{(k)} \|}{\| b \|} < \epsilon\).
def jacobi(A, b, error_tolerance, maximum_iterations=100, demo_mode=False):
"""Splitting A = D + S, where D is its diagonal part,
and S is "L + U", the off diagonal part.
D is stored as just the vector of its values.
"""
n = np.size(b)
D = np.zeros(n)
S = A.copy() # but the diagonal elements must be removed:
for i in range(n):
D[i] = A[i,i]
S[i,i] = 0.0
if demo_mode:
print("D=", D)
print("S=", S)
# Initial approximation and its residual:
x = np.zeros_like(b)
residual = b # b-Ax, initial value just b
norm_b = la.norm(b, Inf)
for iteration in range(maximum_iterations):
if demo_mode:
print(f"iteration {iteration+1}")
x = (b - S@x)/D
residual = b - A@x
error_estimate = la.norm(residual)/norm_b
if demo_mode:
print()
print(f"x={x}")
print(f"r={residual}")
print(f"{error_estimate=}")
if error_estimate < error_tolerance: # We are done here
if demo_mode:
print()
print("Error tolerance achieved.")
return x, error_estimate, iteration+1
if demo_mode:
print()
print("Iteration limit reached; error tolerance not achieved.")
return x, error_estimate, iteration+1
2.8.6. Test Case 1: a very diagonally dominant matrix#
The exact solution is specified, so that we can check the true error, not just the backward error.
A_1 = np.array([[2, -0.1, -0.1],[-0.1, 2, -0.1],[-0.1, -0.1, 2]], dtype=float)
print("A_1=\n", A_1)
print(f"with condition number {la.cond(A_1, Inf)}")
# Create a known solution:
x_exact = np.ones(3) # an array of all ones.
b_1 = A_1 @ x_exact
print(f"b={b_1}")
A_1=
[[ 2. -0.1 -0.1]
[-0.1 2. -0.1]
[-0.1 -0.1 2. ]]
with condition number 1.2222222222222223
b=[1.8 1.8 1.8]
[x_1, error_estimate_1, iterations_1] = jacobi(A_1, b_1, error_tolerance = 1e-8, demo_mode=True)
print()
print(f"The approximate solution is")
print(f"x={x_1},")
print(f"with relative backward error {error_estimate_1}")
print(f"and true relative error {la.norm(x_1-x_exact)/la.norm(x_exact)}")
print(f"This required {iterations_1} iterations.")
D= [2. 2. 2.]
S= [[ 0. -0.1 -0.1]
[-0.1 0. -0.1]
[-0.1 -0.1 0. ]]
iteration 1
x=[0.9 0.9 0.9]
r=[0.18 0.18 0.18]
error_estimate=0.1732050807568878
iteration 2
x=[0.99 0.99 0.99]
r=[0.018 0.018 0.018]
error_estimate=0.017320508075688787
iteration 3
x=[0.999 0.999 0.999]
r=[0.0018 0.0018 0.0018]
error_estimate=0.0017320508075689713
iteration 4
x=[0.9999 0.9999 0.9999]
r=[0.00018 0.00018 0.00018]
error_estimate=0.0001732050807569541
iteration 5
x=[0.99999 0.99999 0.99999]
r=[1.8e-05 1.8e-05 1.8e-05]
error_estimate=1.7320508075617067e-05
iteration 6
x=[0.999999 0.999999 0.999999]
r=[1.8e-06 1.8e-06 1.8e-06]
error_estimate=1.7320508076542938e-06
iteration 7
x=[0.9999999 0.9999999 0.9999999]
r=[1.8e-07 1.8e-07 1.8e-07]
error_estimate=1.732050806301098e-07
iteration 8
x=[0.99999999 0.99999999 0.99999999]
r=[1.79999999e-08 1.79999999e-08 1.79999999e-08]
error_estimate=1.7320508013156392e-08
iteration 9
x=[1. 1. 1.]
r=[1.80000015e-09 1.80000015e-09 1.79999993e-09]
error_estimate=1.7320508796585665e-09
Error tolerance achieved.
The approximate solution is
x=[1. 1. 1.],
with relative backward error 1.7320508796585665e-09
and true relative error 1.0000000457329382e-09
This required 9 iterations.
2.8.7. Test Case 2: a just barely diagonally dominant matrix#
It must still converge, but as you might suspect it will be slower.
A_2 = np.array([[2.1, -1, -1],[-1, 2.1, -1],[-1, -1, 2.1]], dtype=float)
print("A=\n", A_2)
print(f"with condition number {la.cond(A_2, Inf)}")
# Create a known solution:
x_exact = np.ones(3) # an array of all ones.
b_2 = A_2 @ x_exact
print(f"b={b_2}")
A=
[[ 2.1 -1. -1. ]
[-1. 2.1 -1. ]
[-1. -1. 2.1]]
with condition number 41.00000000000003
b=[0.1 0.1 0.1]
Compare the condition numbers: they once again give a hint that this problem is “harder” than the previous one.
[x_2, error_estimate_2, iterations_2] = jacobi(A_2, b_2, error_tolerance = 1e-3, maximum_iterations=20, demo_mode=True)
print()
print(f"The approximate solution is")
print(f"x={x_2},")
print(f"with relative backward error {error_estimate_2}")
print(f"and true relative error {la.norm(x_2-x_exact)/la.norm(x_exact)}")
print(f"This required {iterations_2} iterations.")
D= [2.1 2.1 2.1]
S= [[ 0. -1. -1.]
[-1. 0. -1.]
[-1. -1. 0.]]
iteration 1
x=[0.04761905 0.04761905 0.04761905]
r=[0.0952381 0.0952381 0.0952381]
error_estimate=1.649572197684645
iteration 2
x=[0.09297052 0.09297052 0.09297052]
r=[0.09070295 0.09070295 0.09070295]
error_estimate=1.5710211406520427
iteration 3
x=[0.1361624 0.1361624 0.1361624]
r=[0.08638376 0.08638376 0.08638376]
error_estimate=1.4962106101448025
iteration 4
x=[0.17729753 0.17729753 0.17729753]
r=[0.08227025 0.08227025 0.08227025]
error_estimate=1.424962485852193
iteration 5
x=[0.21647383 0.21647383 0.21647383]
r=[0.07835262 0.07835262 0.07835262]
error_estimate=1.3571071293830408
iteration 6
x=[0.2537846 0.2537846 0.2537846]
r=[0.07462154 0.07462154 0.07462154]
error_estimate=1.2924829803648006
iteration 7
x=[0.28931867 0.28931867 0.28931867]
r=[0.07106813 0.07106813 0.07106813]
error_estimate=1.2309361717760003
iteration 8
x=[0.32316064 0.32316064 0.32316064]
r=[0.06768394 0.06768394 0.06768394]
error_estimate=1.1723201635961908
iteration 9
x=[0.35539108 0.35539108 0.35539108]
r=[0.06446089 0.06446089 0.06446089]
error_estimate=1.1164953939011344
iteration 10
x=[0.38608675 0.38608675 0.38608675]
r=[0.06139133 0.06139133 0.06139133]
error_estimate=1.063328946572509
iteration 11
x=[0.41532071 0.41532071 0.41532071]
r=[0.05846793 0.05846793 0.05846793]
error_estimate=1.0126942348309609
iteration 12
x=[0.44316258 0.44316258 0.44316258]
r=[0.05568374 0.05568374 0.05568374]
error_estimate=0.96447069983901
iteration 13
x=[0.46967865 0.46967865 0.46967865]
r=[0.05303214 0.05303214 0.05303214]
error_estimate=0.9185435236562005
iteration 14
x=[0.49493205 0.49493205 0.49493205]
r=[0.0505068 0.0505068 0.0505068]
error_estimate=0.8748033558630481
iteration 15
x=[0.5189829 0.5189829 0.5189829]
r=[0.04810171 0.04810171 0.04810171]
error_estimate=0.8331460532029019
iteration 16
x=[0.54188848 0.54188848 0.54188848]
r=[0.04581115 0.04581115 0.04581115]
error_estimate=0.7934724316218114
iteration 17
x=[0.56370331 0.56370331 0.56370331]
r=[0.04362967 0.04362967 0.04362967]
error_estimate=0.7556880301160114
iteration 18
x=[0.58447935 0.58447935 0.58447935]
r=[0.04155207 0.04155207 0.04155207]
error_estimate=0.7197028858247726
iteration 19
x=[0.60426604 0.60426604 0.60426604]
r=[0.0395734 0.0395734 0.0395734]
error_estimate=0.6854313198331167
iteration 20
x=[0.62311052 0.62311052 0.62311052]
r=[0.03768895 0.03768895 0.03768895]
error_estimate=0.6527917331743975
Iteration limit reached; error tolerance not achieved.
The approximate solution is
x=[0.62311052 0.62311052 0.62311052],
with relative backward error 0.6527917331743975
and true relative error 0.3768894828730004
This required 20 iterations.
2.8.8. Test Case 2 should converge, so let it run#
Test case 2 again, allowing enough iterations that it succeeds — without the huge amount of “demo” output.
[x3, error_estimate_3, iterations_3] = jacobi(A_2, b_2, error_tolerance = 1e-3, maximum_iterations=1000)
print()
print(f"The approximate solution is")
print(f"x={x3},")
print(f"with relative backward error {error_estimate_3}")
print(f"and true relative error {la.norm(x3-x_exact)/la.norm(x_exact)}")
print(f"This required {iterations_3} iterations.")
The approximate solution is
x=[0.99942715 0.99942715 0.99942715],
with relative backward error 0.0009921969088145376
and true relative error 0.0005728451523927536
This required 153 iterations.
2.8.9. A family of test cases, arising from boundary value problems for differential equations#
The family of tri-diagonal matrices \(T\) of the form
and variants of this arise in the solutions of boundary value problems for ODEs like
and also in related problems for partial differential equations.
This uses the approximation
to be seen in Example 4.3.
Thus these provide useful test cases — usually with \(h = (b-a)/n\).
2.8.10. Exercises#
Exercise 2.17 (Test the Jacobi method)
Test the above implementation jacobi of the Jacobi method with matrices of form \(T\) above with \(h= 1/n\) (as for \(a=0\), \(b=1\));
first \(K=1\) and then \(K=0\).
The latter case might be more challenging.
Do this for several values of \(n\), increasingly geometrically, like \(n=4, 8, \dots\).
Also, check the condition number for each matrix, and look for a connection between that and how well or badly the algorithm performs.
Exercise 2.18 (Implement and test the Gauss-Seidel method, and compare to Jacobi)
Implement the Gauss-Siedel method in the format
[x, error_estimate, iterations] = gauss_seidel(A, b, error_tolerance, maximum_iterations, demo_mode)
as for jacobi above, and use the same test cases.
(It might help to test your code first with the “easy” matrix \(A_1\) in Test Case 1 above.)
Then compare the speed/cost of the two methods.
One way to do this is by using Python’s “stop watch”, function time.time:
see the description of Python module time in the Python tutorial.