Chapter 7

Iterative Methods for Solving Linear

Systems

7.1 Convergence of Sequences of Vectors and Matrices

In Chapter 5 we have discussed some of the main methods

for solving systems of linear equations. These methods
are direct methods, in the sense that they yield exact
solutions (assuming infinite precision!).

Another class of methods for solving linear systems con-

sists in approximating solutions using iterative methods.

409
410 CHAPTER 7. ITERATIVE METHODS FOR SOLVING LINEAR SYSTEMS

The basic idea is this: Given a linear system Ax = b

(with A a square invertible matrix), find another matrix
B and a vector c, such that
1. The matrix I B is invertible
2. The unique solution x e of the system Ax = b is iden-
e of the system
tical to the unique solution u

u = Bu + c,

and then, starting from any vector u0, compute the se-
quence (uk ) given by

uk+1 = Buk + c, k 2 N.

Under certain conditions (to be clarified soon), the se-

e which is the unique
quence (uk ) converges to a limit u
solution of u = Bu + c, and thus of Ax = b.
7.1. CONVERGENCE OF SEQUENCES OF VECTORS AND MATRICES 411

Let (E, k k) be a normed vector space. Recall that a

sequence (uk ) of vectors uk 2 E converges to a limit
u 2 E, if for every ✏ > 0, there some natural number N
such that
kuk uk  ✏, for all k N.

We write
u = lim uk .
k7!1

If E is a finite-dimensional vector space and dim(E) =

n, we know from Theorem 6.3 that any two norms are
equivalent, and if we choose the norm k k1, we see that
the convergence of the sequence of vectors uk is equivalent
to the convergence of the n sequences of scalars formed
by the components of these vectors (over any basis).
412 CHAPTER 7. ITERATIVE METHODS FOR SOLVING LINEAR SYSTEMS

The same property applies to the finite-dimensional vec-

tor space Mm,n(K) of m ⇥ n matrices (with K = R or
K = C), which means that the convergence of a sequence
(k)
of matrices Ak = (aij ) is equivalent to the convergence
(k)
of the m ⇥ n sequences of scalars (aij ), with i, j fixed
(1  i  m, 1  j  n).

The first theorem below gives a necessary and sufficient

condition for the sequence (B k ) of powers of a matrix B
to converge to the zero matrix.

Recall that the spectral radius ⇢(B) of a matrix B is the

maximum of the moduli | i| of the eigenvalues of B.
7.1. CONVERGENCE OF SEQUENCES OF VECTORS AND MATRICES 413

Theorem 7.1. For any square matrix B, the follow-

ing conditions are equivalent:
(1) limk7!1 B k = 0,
(2) limk7!1 B k v = 0, for all vectors v,
(3) ⇢(B) < 1,
(4) kBk < 1, for some subordinate matrix norm k k.

The following proposition is needed to study the rate of

convergence of iterative methods.

Proposition 7.2. For every square matrix B and ev-

ery matrix norm k k, we have

lim kB k k1/k = ⇢(B).

k7!1

We now apply the above results to the convergence of

iterative methods.
414 CHAPTER 7. ITERATIVE METHODS FOR SOLVING LINEAR SYSTEMS

7.2 Convergence of Iterative Methods

Recall that iterative methods for solving a linear system

Ax = b (with A invertible) consists in finding some ma-
trix B and some vector c, such that I B is invertible,
e of Ax = b is equal to the unique
and the unique solution x
e of u = Bu + c.
solution u

Then, starting from any vector u0, compute the sequence

(uk ) given by

uk+1 = Buk + c, k 2 N,

and say that the iterative method is convergent i↵

e,
lim uk = u
k7!1

for every initial vector u0.

7.2. CONVERGENCE OF ITERATIVE METHODS 415

Here is a fundamental criterion for the convergence of any

iterative methods based on a matrix B, called the matrix
of the iterative method .

Theorem 7.3. Given a system u = Bu + c as above,

where I B is invertible, the following statements are
equivalent:
(1) The iterative method is convergent.
(2) ⇢(B) < 1.
(3) kBk < 1, for some subordinate matrix norm k k.

The next proposition is needed to compare the rate of

convergence of iterative methods.

It shows that asymptotically, the error vector

ek = B k e0 behaves at worst like (⇢(B))k .
416 CHAPTER 7. ITERATIVE METHODS FOR SOLVING LINEAR SYSTEMS

Proposition 7.4. Let k k be any vector norm, let B

e be
be a matrix such that I B is invertible, and let u
the unique solution of u = Bu + c.

(1) If (uk ) is any sequence defined iteratively by

uk+1 = Buk + c, k 2 N,
then 
lim sup kuk ek1/k = ⇢(B).
u
k7!1 ku0 u
ek=1
7.2. CONVERGENCE OF ITERATIVE METHODS 417

(2) Let B1 and B2 be two matrices such that I B1

and I B2 are invertibe, assume that both u = B1u+c1
and u = B2u + c2 have the same unique solution u e,
and consider any two sequences (uk ) and (vk ) defined
inductively by
uk+1 = B1uk + c1
vk+1 = B2vk + c2,
with u0 = v0. If ⇢(B1) < ⇢(B2), then for any ✏ > 0,
there is some integer N (✏), such that for all k N (✏),
we have
 1/k
kvk uek ⇢(B2)
sup .
ku0 uek=1 kuk ek
u ⇢(B1) + ✏

In light of the above, we see that when we investigate new

iterative methods, we have to deal with the following two
problems:
1. Given an iterative method with matrix B, determine
whether the method is convergent. This involves de-
termining whether ⇢(B) < 1, or equivalently whether
there is a subordinate matrix norm such that
kBk < 1.
418 CHAPTER 7. ITERATIVE METHODS FOR SOLVING LINEAR SYSTEMS

By Proposition 6.9, this implies that I B is invertible

(since k Bk = kBk, Proposition 6.9 applies).
2. Given two convergent iterative methods, compare them.
The iterative method which is faster is that whose ma-
trix has the smaller spectral radius.

We now discuss three iterative methods for solving linear

systems:
1. Jacobi’s method
2. Gauss-Seidel’s method
3. The relaxation method.
7.3. METHODS OF JACOBI, GAUSS-SEIDEL, AND RELAXATION 419

7.3 Description of the Methods of Jacobi,

Gauss-Seidel, and Relaxation

The methods described in this section are instances of the

following scheme: Given a linear system Ax = b, with A
invertible, suppose we can write A in the form

A=M N,

with M invertible, and “easy to invert,” which means

that M is close to being a diagonal or a triangular matrix
(perhaps by blocks).

Then, Au = b is equivalent to
M u = N u + b,
that is,
1 1
u=M Nu + M b.

Therefore, we are in the situation described in the previ-

ous sections with B = M 1N and c = M 1b.
420 CHAPTER 7. ITERATIVE METHODS FOR SOLVING LINEAR SYSTEMS

In fact, since A = M N , we have

1 1 1
B=M N =M (M A) = I M A,

1
which shows that I B=M A is invertible.

The iterative method associated with the matrix

B = M 1N is given by
1 1
uk+1 = M N uk + M b, k 0,

starting from any arbitrary vector u0.

From a practical point of view, we do not invert M , and

instead we solve iteratively the systems

M uk+1 = N uk + b, k 0.
7.3. METHODS OF JACOBI, GAUSS-SEIDEL, AND RELAXATION 421

Various methods correspond to various ways of choosing

M and N from A. The first two methods choose M
and N as disjoint submatrices of A, but the relaxation
method allows some overlapping of M and N .

To describe the various choices of M and N , it is conve-

nient to write A in terms of three submatrices D, E, F ,
as

A=D E F,

where the only nonzero entries in D are the diagonal en-

tries in A, the only nonzero entries in E are entries in A
below the the diagonal, and the only nonzero entries in
F are entries in A above the diagonal.
422 CHAPTER 7. ITERATIVE METHODS FOR SOLVING LINEAR SYSTEMS

More explicitly, if
0 1
a11 a12 a13 ··· a1n 1 a1n
B C
B a21 a22 a23 ··· a2n
1 C
a2n
B C
B a31 a32 a33 ··· a3n C
a3n
B 1 C
A=B .. .. .. ... .. C,
..
B C
B C
Ba a a · · · a a C
@ n 1 1 n 1 2 n 1 3 n 1 n 1 n 1 n A
an 1 an 2 an 3 · · · an n 1 an n

then
0 1
a11 0 0 ··· 0 0
B C
B0 a22 0 ··· 0 0 C
B C
B0 0 a33 ··· 0 0 C
B C
D=B . .. .. ... .. .. C ,
B . C
B C
B0 0 0 · · · an 1 n 0 C
@ 1 A
0 0 0 ··· 0 an n
7.3. METHODS OF JACOBI, GAUSS-SEIDEL, AND RELAXATION 423
0 1
0 0 0 ··· 0 0
B C
B a21 0 0 ··· 0 0C
B C
B a31 a32 0 ··· 0 0C
B C
E=B .. .. ... ... .. .. C ,
B C
B C
Ba ... 0C
@ n 1 1 an 1 2 an 1 3 0 A
an 1 an 2 an 3 · · · an n 1 0

0 1
0 a12 a13 · · · a1n 1 a1n
B C
B0 0 a23 · · · a2n 1 a2n C
B C
B0 0 0 . . . a3n 1 a3n C
B C
F = B. .. .. ... ... .. C.
B. C
B C
B0 0 0 ··· 0 an 1 n C
@ A
0 0 0 ··· 0 0
424 CHAPTER 7. ITERATIVE METHODS FOR SOLVING LINEAR SYSTEMS

In Jacobi’s method , we assume that all diagonal entries

in A are nonzero, and we pick

M =D
N = E + F,

so that
1
B=M N = D 1(E + F ) = I D 1A.

As a matter of notation, we let

J =I D 1A = D 1(E + F ),

which is called Jacobi’s matrix .

7.3. METHODS OF JACOBI, GAUSS-SEIDEL, AND RELAXATION 425

The corresponding method, Jacobi’s iterative method ,

computes the sequence (uk ) using the recurrence

uk+1 = D 1(E + F )uk + D 1b, k 0.

In practice, we iteratively solve the systems

Duk+1 = (E + F )uk + b, k 0.

If we write uk = (uk1 , . . . , ukn), we solve iteratively the

following system:
a11 uk+1
1 = a12 uk2 a13 uk3 ··· a1n ukn + b1
a22 uk+1
2 = a21 uk1 a23 uk3 ··· a2n ukn + b2
.. .. .. .
. . .
an 1 n 1 uk+1
n 1 = an 1 1 uk1 ··· an k
1 n 2 un 2 an k
1 n un + bn 1
an n uk+1
n = an 1 uk1 an 2 uk2 ··· an n 1 ukn 1 + bn

Observe that we can try to “speed up” the method by

using the new value uk+1 1 instead of uk1 in solving for
uk+2
2 using the second equations, and more generally, use
uk+1 k+1 k k k+1
1 , . . . , ui 1 instead of u1 , . . . , ui 1 in solving for ui
in the ith equation.
426 CHAPTER 7. ITERATIVE METHODS FOR SOLVING LINEAR SYSTEMS

This observation leads to the system

a11 uk+1
1 = a12 uk2 a13 uk3 ··· a1n ukn + b1
a22 uk+1
2 = a21 uk+1
1 a23 uk3 ··· a2n ukn + b2
.. .. ..
. . .
an 1 n 1 uk+1
n 1 = an 1 1 uk+1
1 ··· an k+1
1 n 2 un 2 an k
1 n un + bn 1
an n uk+1
n = an 1 uk+1
1 an 2 uk+1
2 ··· an n 1 uk+1
n 1 + bn ,

which, in matrix form, is written

Duk+1 = Euk+1 + F uk + b.

Because D is invertible and E is lower triangular, the ma-

trix D E is invertible, so the above equation is equiva-
lent to

uk+1 = (D E) 1F uk + (D E) 1b, k 0.

The above corresponds to choosing M and N to be

M =D E
N = F,
and the matrix B is given by
1
B=M N = (D E) 1F.
7.3. METHODS OF JACOBI, GAUSS-SEIDEL, AND RELAXATION 427

Since M = D E is invertible, we know that I B=

M 1A is also invertible.

The method that we just described is the iterative method

of Gauss-Seidel , and the matrix B is called the matrix
of Gauss-Seidel and denoted by L1, with

L1 = (D E) 1F.

One of the advantages of the method of Gauss-Seidel is

that is requires only half of the memory used by Jacobi’s
method, since we only need

uk+1 k+1 k k
1 , . . . , ui 1 , ui+1 , . . . , un

to compute uk+1
i .

We also show that in certain important cases (for exam-

ple, if A is a tridiagonal matrix), the method of Gauss-
Seidel converges faster than Jacobi’s method (in this case,
they both converge or diverge simultaneously).
428 CHAPTER 7. ITERATIVE METHODS FOR SOLVING LINEAR SYSTEMS

The new ingredient in the relaxation method is to incor-

porate part of the matrix D into N : we define M and N
by

D
M= E
!
1 !
N= D + F,
!

where ! 6= 0 is a real parameter to be suitably chosen.

Actually, we show in Section 7.4 that for the relaxation

method to converge, we must have ! 2 (0, 2).

Note that the case ! = 1 corresponds to the method of

Gauss-Seidel.

If we assume that all diagonal entries of D are nonzero,

the matrix M is invertible.
7.3. METHODS OF JACOBI, GAUSS-SEIDEL, AND RELAXATION 429

The matrix B is denoted by L! and called the matrix of

relaxation, with

✓ ◆ 1✓ ◆
D 1 !
L! = E D+F
! !
= (D !E) 1((1 !)D + !F ).

The number ! is called the parameter of relaxation.

When ! > 1, the relaxation method is known as succes-
sive overrelaxation, abbreviated as SOR.

At first glance, the relaxation matrix L! seems at lot more

complicated than the Gauss-Seidel matrix L1, but the
iterative system associated with the relaxation method is
very similar to the method of Gauss-Seidel, and is quite
simple.
430 CHAPTER 7. ITERATIVE METHODS FOR SOLVING LINEAR SYSTEMS

Indeed, the system associated with the relaxation method

is given by
✓ ◆ ✓ ◆
D 1 !
E uk+1 = D + F uk + b,
! !

which is equivalent to

(D !E)uk+1 = ((1 !)D + !F )uk + !b,

and can be written

Duk+1 = Duk !(Duk Euk+1 F uk b).

Explicitly, this is the system

a11 uk+1
1 = a11 uk1 !(a11 uk1 + a12 uk2 + a13 uk3 + · · · + a1n 2 ukn k k
2 + a1n 1 un 1 + a1n un b1 )
a22 uk+1
2 = a22 uk2 !(a21 uk+1
1 + a22 uk2 + a23 uk3 + ··· + a2n 2 ukn 2 + a2n 1 ukn 1 + a2n ukn b2 )
..
.
an n uk+1
n = an n ukn !(an 1 uk+1
1 + an 2 uk+1
2 + · · · + an n 2 uk+1 k+1 k
n 2 + an n 1 u n 1 + an n u n bn ).
7.3. METHODS OF JACOBI, GAUSS-SEIDEL, AND RELAXATION 431

What remains to be done is to find conditions that ensure

the convergence of the relaxation method (and the Gauss-
Seidel method), that is:
1. Find conditions on !, namely some interval I ✓ R
so that ! 2 I implies ⇢(L! ) < 1; we will prove that
! 2 (0, 2) is a necessary condition.
2. Find if there exist some optimal value !0 of ! 2 I,
so that
⇢(L!0 ) = inf ⇢(L! ).
!2I

We will give partial answers to the above questions in the

next section.

It is also possible to extend the methods of this section by

using block decompositions of the form A = D E F ,
where D, E, and F consist of blocks, and with D an
invertible block-diagonal matrix.
432 CHAPTER 7. ITERATIVE METHODS FOR SOLVING LINEAR SYSTEMS

7.4 Convergence of the Methods of Jacobi,

Gauss-Seidel, and Relaxation

We begin with a general criterion for the convergence of

an iterative method associated with a (complex) Hermi-
tian, positive, definite matrix, A = M N . Next, we
apply this result to the relaxation method.

Proposition 7.5. Let A be any Hermitian, positive,

definite matrix, written as
A=M N,
with M invertible. Then, M ⇤ + N is Hermitian, and
if it is positive, definite, then
1
⇢(M N ) < 1,
so that the iterative method converges.
7.4. CONVERGENCE OF THE METHODS 433

Now, as in the previous sections, we assume that A is

written as A = D E F , with D invertible, possibly
in block form.

The next theorem provides a sufficient condition (which

turns out to be also necessary) for the relaxation method
to converge (and thus, for the method of Gauss-Seidel to
converge).

This theorem is known as the Ostrowski-Reich theorem.

Theorem 7.6. If A = D E F is Hermitian, pos-

itive, definite, and if 0 < ! < 2, then the relaxation
method converges. This also holds for a block decom-
position of A.
434 CHAPTER 7. ITERATIVE METHODS FOR SOLVING LINEAR SYSTEMS

Remark: What if we allow the parameter ! to be a

nonzero complex number ! 2 C? In this case, the relax-
ation method also converges for ! 2 C, provided that
|! 1| < 1.
This condition reduces to 0 < ! < 2 if ! is real.

Unfortunately, Theorem 7.6 does not apply to Jacobi’s

method, but in special cases, Proposition 7.5 can be used
to prove its convergence.

On the positive side, if a matrix is strictly column (or

row) diagonally dominant, then it can be shown that the
method of Jacobi and the method of Gauss-Seidel both
converge.
7.4. CONVERGENCE OF THE METHODS 435

The relaxation method also converges if ! 2 (0, 1], but

this is not a very useful result because the speed-up of
convergence usually occurs for ! > 1.

We now prove that, without any assumption on

A = D E F , other than the fact that A and D are
invertible, in order for the relaxation method to converge,
we must have ! 2 (0, 2).

Proposition 7.7. Given any matrix A = D E F ,

with A and D invertible, for any ! 6= 0, we have
⇢(L! ) |! 1|.
Therefore, the relaxation method (possibly by blocks)
does not converge unless ! 2 (0, 2). If we allow ! to
be complex, then we must have
|! 1| < 1
for the relaxation method to converge.
436 CHAPTER 7. ITERATIVE METHODS FOR SOLVING LINEAR SYSTEMS

We now consider the case where A is a tridiagonal ma-

trix , possibly by blocks.

We begin with the case ! = 1, which is technically easier

to deal with.

The following proposition gives us the precise relationship

between the spectral radii ⇢(J) and ⇢(L1) of the Jacobi
matrix and the Gauss-Seidel matrix.

Proposition 7.8. Let A be a tridiagonal matrix (pos-

sibly by blocks). If ⇢(J) is the spectral radius of the
Jacobi matrix and ⇢(L1) is the spectral radius of the
Gauss-Seidel matrix, then we have
⇢(L1) = (⇢(J))2.
Consequently, the method of Jacobi and the method of
Gauss-Seidel both converge or both diverge simultane-
ously (even when A is tridiagonal by blocks);

when they converge, the method of Gauss-Seidel con-

verges faster than Jacobi’s method.
7.4. CONVERGENCE OF THE METHODS 437

We now consider the more general situation where ! is

any real in (0, 2).

Proposition 7.9. Let A be a tridiagonal matrix (pos-

sibly by blocks), and assume that the eigenvalues of
the Jacobi matrix are all real. If ! 2 (0, 2), then the
method of Jacobi and the method of relaxation both
converge or both diverge simultaneously (even when
A is tridiagonal by blocks).

When they converge, the function ! 7! ⇢(L! ) (for

! 2 (0, 2)) has a unique minimum equal to !0 1 for
2
!0 = p ,
1+ 1 (⇢(J))2
where 1 < !0 < 2 if ⇢(J) > 0. We also have
⇢(L1) = (⇢(J))2, as before.
438 CHAPTER 7. ITERATIVE METHODS FOR SOLVING LINEAR SYSTEMS

Combining the results of Theorem 7.6 and Proposition

7.9, we obtain the following result which gives precise
information about the spectral radii of the matrices J,
L1, and L! .

Proposition 7.10. Let A be a tridiagonal matrix (pos-

sibly by blocks) which is Hermitian, positive, definite.
Then, the methods of Jacobi, Gauss-Seidel, and relax-
ation, all converge for ! 2 (0, 2). There is a unique
optimal relaxation parameter
2
!0 = p ,
1 + 1 (⇢(J)) 2

such that
⇢(L!0 ) = inf ⇢(L! ) = !0 1.
0<!<2
7.4. CONVERGENCE OF THE METHODS 439

Furthermore, if ⇢(J) > 0, then

⇢(L!0 ) < ⇢(L1) = (⇢(J))2 < ⇢(J),
and if ⇢(J) = 0, then !0 = 1 and ⇢(L1) = ⇢(J) = 0.

Remark: It is preferable to overestimate rather than

underestimate the relaxation parameter when the opti-
mum relaxation parameter is not known exactly.
440 CHAPTER 7. ITERATIVE METHODS FOR SOLVING LINEAR SYSTEMS