A1 Linear Algebra (II)

MT 2024

corrections to: [Link]@[Link]

With grateful acknowledgement to Professor Chris Martin for permission to

use his course materials from previous years.
References
• Strang(2016). Introduction to Linear Algebra (5th edn). Wellesley-Cambridge
Press.

• Boyd & Vandenberghe (2018). Introduction to Applied Linear Algebra

(1st edn). Cambridge Univ. Press.

• Kreyszig (2006). Advanced engineering mathematics (9th edn). Wiley.

Version: Thursday 10th October, 2024 at 12:12.

i
Contents

5 The singular value decomposition (SVD) 68

5.1 Eigenvectors and eigenvalues of symmetric matrices . . . . . 69
5.2 The singular value decomposition . . . . . . . . . . . . . . . 75
5.3 SVD and the four subspaces . . . . . . . . . . . . . . . . . . 80
5.4 Matrix approximation and data compression . . . . . . . . . 81

6 Norms and conditioning 84

6.1 Vector norms . . . . . . . . . . . . . . . . . . . . . . . . . . 84
6.2 Matrix norms . . . . . . . . . . . . . . . . . . . . . . . . . . 87
6.3 Conditioning of linear systems . . . . . . . . . . . . . . . . . 95

7 Iterative methods for Ax = b 104

7.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . 104
7.2 Jacobi method (‘simultaneous displacements’) . . . . . . . . 105
7.3 Gauss-Seidel method (‘successive displacements’) . . . . . . . 107
7.4 Compact Formulations . . . . . . . . . . . . . . . . . . . . . 110
7.5 Convergence criteria . . . . . . . . . . . . . . . . . . . . . . 116

8 Computing Eigenvalues and Eigenvectors 124

8.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . 124
8.2 Power method . . . . . . . . . . . . . . . . . . . . . . . . . 125
8.3 Rayleigh method . . . . . . . . . . . . . . . . . . . . . . . . 131
8.4 Matrix deflation . . . . . . . . . . . . . . . . . . . . . . . . 135

ii
5 The singular value decomposition (SVD)

Let us look again at the eigenvalue-eigenvalue equation, i.e.

Ax = λx.

We found that it is sometimes possible to diagonalise the matrix A by using

its eigenvectors and eigenvalues. If the eigenvectors of A are the vectors si
with corresponding eigenvalues λi , then we could form the matrices

 
h i λ1
S = s1 . . . sn , Λ= ...
.
 

λn

We could then write down all of the eigenvector-eigenvalue pairs together in

the single equation AS = SΛ. If the matrix is defective, i.e. if it does not
have a full set of eigenvectors, then this was as far as we could go. Otherwise,
we could write
A = SΛS−1 .

In other words, the matrix S could be used to diagonalise A, and we find this
to be very useful in a wide variety of applications. Ideally we would like to be
able to diagonalise A even if it doesn’t have a full set of eigenvectors. This
could be because it is defective, but it can also happen because A might not
even be square. If it is not square then there are no eigenvectors to talk about;
the equation Ax = λx would not make any sense in the first place.

68
5.1 Eigenvectors and eigenvalues of symmetric matrices
On the other hand, there is one case where the eigenvectors are guaranteed
to exist, and are even mutually orthonormal1 . This is the case when A is
symmetric, i.e. A = AT . In that case we can write

A = VΛVT ,

because we can be sure that V−1 will exist and that V−1 = VT . We will use
this idea to help us diagonalise any matrix, whether or not it is defective or
even square.

Let’s first prove the two nice properties of a symmetric matrix. The first was
this one:

All of the eigenvalues of a symmetric matrix are real.

To prove this, assume that v is an eigenvector of a symmetric matrix A and

λ is its eigenvalue. We can’t be sure initially if the eigenvalue might be a
complex number λ = a + bȷ, so we can also define its complex conjugate
λ̄ = a − bȷ. The same is true for the vector v; it might also have one or more
entries that are complex numbers. Define the complex conjugate of v to be
v̄, which we form by taking the complex conjugate of each of its elements.

Our eigenvalue equation tells us that

Av = λv, and v̄T A = v̄T λ̄.

To obtain the second of those equations from the first we just take the trans-
pose of both sides, followed by the complex conjugate of both sides. Re-
member that A has all real values, so its complex conjugate (again taken
element-by-element) is just Ā = A.
1
We will show that they are all mutually orthogonal. To get orthonormal vectors, we just make them all
unit length. This is why we can write VT V = I.

69
Now take the product of the left hand equation with v̄T (multiplying from the
left), and the product of the right hand equation with v (multiplying from the
right) to get
v̄T Av = v̄T λv and v̄T Av = v̄T λ̄v

Looking at the right-hand side of each of these equations, we have λv̄T v =

λ̄v̄T v. The term v̄T v is just the length of the (possibly complex-valued)
vector though, i.e.
Xn
T
v̄ v = |vi |2 = ∥v∥2 ,
i=1

and it can’t be zero because v is not zero. We therefore have λ = λ̄, which
means that λ must be real.

The corresponding eigenvector v then comes from solving (A − λI)v = 0,

i.e. a system of linear equations with real coefficients, so the solution v will
also be real. However, the more important fact is that the eigenvectors are all
orthogonal. This is the second of our nice properties of symmetric matrices:

The eigenvectors of a symmetric matrix are mutually orthogonal.

To prove this one, assume that v and w are eigenvectors corresponding to two
different eigenvalues λ and ω, with λ ̸= ω. They will satisfy the two equations

Av = λv and wT AT = ωwT .

where in the second equation we have transposed both sides of Aw = ωw.

Multiply the first equation from the left by wT and the second equation from
the right by v to get

wT Av = λwT v and wT AT v = ωwT v.

Since A = AT , the left-hand terms are the same. This means that we can

70
equate the two right-hand terms to get

λwT v = ωwT v or (λ − ω)wT v = 0.

Since we assumed that ω ̸= λ, we can conclude that wT v, i.e. that the vectors
are orthogonal.
5.1.1 Positive-definite matrices
If a matrix is symmetric, then its eigenvalues are all real and its eigenvectors
are mutually orthogonal. If, in addition, all of the eigenvalues are positive,
then we say that the matrix is positive definite. If the eigenvalues are only
non-negative (i.e. we allow either positive numbers or zeros to be eigenvalues),
then we say that the matrix is positive semidefinite.

Definition (Positive Definiteness) A symmetric matrix is positive def-

inite if all of its eigenvalues are positive. It is positive semidefinite if its
eigenvalues are all non-negative.

Positive definite matrices play a role similar to positive scalars. If a is a positive

number, then ax2 > 0 for any number x. If A is a positive definite matrix,
then xT Ax > 0 is for any vector x.

To see why this is the case, suppose we take any non-zero vector x and write
it as a weighted sum of the eigenvectors vi of a symmetric matrix A with
positive eigenvalues, i.e.

x = c1 v1 + · · · + cn vn .

If we multiply by A, we get

Ax = c1 Av1 + · · · + cn Avn = c1 λ1 v1 + · · · + cn λn vn .

Recalling that all of the eigenvectors can be chosen to be orthonormal, i.e.

71
that viT vj = 0 and viT vi = 1, multiply Ax from the left by xT to get

xT Ax = (c1 v1 + · · · + cn vn )T Ax
= (c1 v1 + · · · + cn vn )T (c1 λ1 v1 + · · · + cn λn vn )
= c21 λ1 + · · · + c2n λn > 0.

We can actually use this property as an alternative test (or even a definition)
for positive definiteness, as in the following:

(Positive Definiteness) A matrix is positive definite if xT Ax > 0 for

every x. It is positive semidefinite if xT Ax ≥ 0 for every x.

There is a reason that we have gone to so much effort to understand the

eigenvalues and eigenvectors of symmetric matrices, and in particular to discuss
when the eigenvalues will be positive. It has led us to the following statement,
which we will use to help diagonalise matrices of arbitrary size and shape, i.e.
even non-square ones. Here it is :

For any matrix A, the matrices AT A and AAT are both symmetric and
positive semidefinite.

The proof of this is easy. Take any vector x and form the product xT (AT A)x =
(Ax)T (Ax) = ∥Ax∥2 ≥ 0. The same argument holds for the product AAT .

Example Find the eigenvectors and eigenvalues of the matrices AT A and

AAT , where  
1 1
A = 1 2
 

2 1

72
Solution First just compute the two products :
 
" # 2 3 3
6 5
AT A = , AAT = 3 5 4 .
 
5 6
3 4 5

Let’s take the smaller of these first, i.e. the 2 × 2 matrix AT A. Its eigen-
values are the roots of the equation

6−λ 5
det(A − λI) =
5 6−λ
= (6 − λ)2 − 25
= λ2 − 12λ + 11 = (λ − 11)(λ − 1) = 0.

The eigenvalues are therefore λ1 = 11 and λ2 = 1.

What about the eigenvectors? They are the solutions to

" # " #
6 − λ1 5 6 − λ2 5
v = 0, and v = 0.
5 6 − λ1 5 6 − λ2

The solutions are

" # " #
1 1 1 1
v1 = √ and v2 = √ ,
2 1 2 −1

where we have scaled the eigenvectors so that they are unit length. We
note also that they are orthogonal (v1T v2 = 0) as expected.

Now take the 3 × 3 matrix AAT . Its eigenvalues are the roots of the

73
 equation

2−λ 3 3
det(A − λI) = 3 5−λ 4
3 4 5−λ

5−λ 4 3 4 3 5−λ
= (2 − λ) −3 +3
4 5−λ 3 5−λ 3 4
= −λ3 + 12λ2 − 11λ
= −λ(λ − 11)(λ − 1),

or (λ1 , λ2 , λ3 ) = (11, 1, 0).

Now compute the corresponding eigenvectors. They are the solutions to

 
2−λ 3 3
 3 5−λ 4  ui = 0,
 

3 4 5−λ

or
   
2 0 3
1   1   1  
u1 = √ 3 and u2 = √  1 and u3 = √ −1 ,
22 2 11
3 −1 −1

where we have once again scaled the solutions to be unit length. We notice
again that they are all orthogonal to each other.

There are two final points to make about the preceding example.

The first is that the eigenvalues turned out to be the same, with an additional
zero eigenvalue for the bigger matrix. This was not a coincidence, as will
become apparent in the next section.

The second point is that the rank of the matrix A is the same as the rank of
the matrix AT A, which is the same as the rank of the matrix AAT . This is

74
a general property of any matrix A:

If A has rank r, then the matrices AT A and AAT also have rank r.

5.2 The singular value decomposition

We can now go back to the question of how to diagonalise a general matrix
A. The nicest case we have so far is the one where A is symmetric, in which
case we could write its eigenvector equations as

Asi = λi si . (5.1)

and the symmetry of A ensures that the vectors si are all orthogonal. We can
then diagonalise A by writing A = SΛS−1 . Since A is symmetric, we can do
even better by choosing the eigenvectors to be orthonormal and writing

A = SΛST . (5.2)

The key point is that, for a square and symmetric matrix, we can diagonalise
the matrix by multiplying from the left and the right by (the same!) orthonor-
mal matrix.

If A isn’t square then this is obviously not going to work. Even if it is square,
the solutions to the eigenvector equation won’t all be orthogonal in general,
so S−1 ̸= ST . If the matrix is defective, then S−1 won’t even exist.

Let’s look instead for two different sets of orthogonal vectors to go on the
left and right-hand sides of an equation like (5.1). In other words, look for
solutions like
Avi = σi ui , i = 1, . . . , s

where the vectors ui are unit length and orthogonal to each other, and the
vectors vi are also unit length and orthogonal to each other. Now instead of
eigenvalues λi , we have the singular values σi .

In matrix form we are asking to again do something like (5.2), but by multi-

75
plying on the left and the right by different square matrices with orthonormal
columns. In other words, we want to write A as

A = UΣVT .

If we remember that A has m rows and n columns, then the square matrix U
of left singular vectors should have size m × m. Since it is orthonormal it
should also satisfy
UT U = UUT = I.

We should expect the matrix V of right singular vectors to instead have

size n × n, and also to satisfying

VT V = VVT = I.

The size of the matrix Σ of singular values will be m × n, but we said also
that it should be diagonal. This means that when m > n (more rows than
columns) or when n > m (more columns than rows), we will need to pad out
the matrix Σ with zeros to make it the right shape. It will be either

 
σ1

 ... 

Σ=  when m > n

 σs 

0
or
 
σ1
Σ=
 ... 0 when n > m.

σs

We can see immediately that the number s is just s = min(m, n).

76
How can we find all of these matrices? The answer is to make a symmetric
matrix out of A and use eigenvalue decomposition. First take AT A – we get

AT A = (UΣVT )T (UΣVT )
= VΣT UT UΣVT
= VΣT ΣVT .

The matrix ΣT Σ is diagonal and contains the positive numbers σi2 , along with
some additional zeros in the lower right-hand corner if m > n. What we have
written is the eigenvalue decomposition of the matrix AT A! We know that
this matrix will be positive semidefinite, since its eigenvalues λi will all satisfy
λi = σi2 ≥ 0. The matrix AT A will also have some extra zero eigenvalues if
m > n.

Let’s now take the product AAT . We get

AAT = (UΣVT )(UΣVT )T

= UΣT VT VΣT UT
= UΣΣT UT .

The matrix ΣΣT is diagonal and contains the positive numbers σi2 , along with
some additional zeros in the lower right-hand corner if n > m. This time,
what we have written is the eigenvalue decomposition of the matrix AAT !
We know that this matrix will be positive semidefinite, since its eigenvalues λi
will all satisfy λi = σi2 ≥ 0. The matrix AAT will also have some extra zero
eigenvalues if n > m.

We have allowed for up to s = min(m, n) singular values, but in truth only

r ≤ min(m, n) can be strictly positive. The number r is just the rank of A,
which we have said is the same as the rank of AT A and also of AAT . These
two matrices share the same non-zero eigenvalues, and there are r of them.
Their square roots give us the singular values of A.

We now have everything we need to compute the singular value decomposition

77
of a matrix A. Our procedure is as follows:

(Computing the SVD) We can compute the singular value decompo-

sition of the m × n matrix A as follows:

1) Find the eigenvalue decomposition of AT A, with eigenvectors vi and

corresponding eigenvalues λi . There will be r positive eigenvalues,
where r is the rank of A. The rest will be zero.
√
2) For each of the r non-zero eigenvalues λi , set σi = λi . Compute the
corresponding vector ui by setting

1
Avi = σi ui =⇒ ui = Avi .
σi

3) The remaining min(m, n) − r singular values are zero.

4) Choose n − r more vectors vi as an orthogonal basis for the null space

of AT A. They are the eigenvectors of AT A with zero eigenvalues.

5) Choose m − r more vectors ui as an orthogonal basis for the null space

of AAT . They are the eigenvectors of AAT with zero eigenvalues.

Let’s look at an example matrix and compute its singular value decomposition.
Example Compute the singular value decomposition of the matrix
 
1 1
A = 1 2 .
 

2 1

Solution We have seen this example already in the previous section, and
have computed the eigenvectors and eigenvalues of AT A. They were
" # " #
1 1 1 1
v1 = √ and v2 = √ .
2 1 2 −1

78
 These are our right singular vectors, i.e. we set V = [v1 v2 ]. The
corresponding eigenvalues λ1 = 11 and λ1 = 1. The singular values are
√ √
therefore σ1 = 11 and σ2 = 1 = 1.

Now we compute the left singular vectors ui . They should turn out to be
eigenvectors of AAT , but we instead compute them as:
   
1 1 " #! 2
1 1  1 1 1  
u1 = Av1 = √ 1 2 √ = √ 3

σ1 11 2 1 22
2 1 3

and
  
1 1 " #! 0
1 1  1 1 1  
u2 = Av2 = √ 1 2 √ = √ −1

σ2 1 2 −1 2
2 1 1

We still need a third left singular vector ui , which should be orthogonal

to the others. We already have it though; it is the remaining eigenvector
of AAT , which corresponded to the eigenvalue 0. Our full SVD for A is
therefore   √
√2 √3

0 11 0 " #
 322 −1 −1 11  √1 √1
 2 2
A= √ √ √ 
 22 2 11   0 1 √1 √

−1
3 1 −1
√
22
√
2
√
11
0 0 | {z 2 } 2
| {z } | {z } VT
U Σ

There is a small detail that should not be overlooked when computing the vec-
tors ui that we use to form the matrix of left singular vectors U = [u1 u2 u3 ].
They are all eigenvectors of AAT , with the last one corresponding to the
zero eigenvalue. We couldn’t just use the vectors we computed in the pre-
vious example though, because for eigenvectors the direction doesn’t matter.
When we compute our left singular vectors ui the sign does matter and we
need to make sure we get it the right way around, since we need it to satisfy
Avi = σi ui .

79
5.3 SVD and the four subspaces
The singular value decomposition can be very useful in applications, but it also
reveals something fundamental about the matrix A. If our matrix A is m × n
and has rank r, then it will have r singular values strictly greater than zero.

We can write the singular value decomposition as

 
v1T
  .. 
σ1  . 
...  
h i
 0  vrT 
A = u1 · · · ur ur+1 · · · um   (5.3)

 T


 σr vr+1

. 
0 0  ..


vnT

The singular vectors can be divided into four groups. They are the left and
right singular vectors that match up with the non-zero singular values, and
the remaining left and right vectors (that match up with zeros). They have
the following property:

Theorem (SVD and the four subspaces) Suppose that A is m × n

with rank r. Then

• The first r left singular vectors are an orthonormal basis for the column
space of A, i.e. im(A).

• The last m − r left singular vectors are an orthonormal basis for the left
nullspace of A, i.e. ker(AT ).

• The first r right singular vectors are an orthonormal basis for the row
space of A, i.e. im(AT ).

• The last n − r left singular vectors are an orthonormal basis for the
nullspace of A, i.e. ker(A).

80
5.4 Matrix approximation and data compression
If we look again at (5.3), we see that the zero blocks can be removed to get
the ‘economy size‘ SVD of A. We can write this as

T
A = U1:r Σ1:r V1:r

where U1:r is the matrix formed by taking the first r columns of U (and
likewise for V1:r ). The matrix Σ1:r is just the upper r × r block of Σ.

We can also write this out explicitly as the sum

A = u1 σ1 v1T + · · · + ur σr vrT ,

and we typically order the singular values such that σ1 > σ2 > · · · > σr > 0.
Each of the products ui σi viT is a matrix with rank one. All of the terms in this
summation are linearly independent, so each rank one matrix in the summation
increases the overall rank of the sum by 1, until we get back A with rank r.

If A is very large, then it might be difficult or expensive to store. This leads to

a natural question: What is the best low rank approximation of A? By ‘best‘
we mean the low rank matrix that is closest to A in a least-squares sense.

The singular value decomposition of A gives us the answer. The best rank k
approximation of A is just the sum of the first k terms in the singular value
decomposition, i.e. our best approximation is

A ≈ u1 σ1 v1T + · · · + ur σr vkT .

81
Example The image show below can be represented as a matrix of values
between 0 (black pixels) and 1 (white pixels). The size of the image is
2767 × 2006 pixels.

The plot next to the image shows the singular values, ordered by size. It is
clear from this figure that the first 100 or so singular values are dominant
(i.e. much bigger than the others). We should therefore expect that a low
rank approximation of our matrix-valued image data might give a good
representation.

10 4

10 3

10 2

10 1

0
10

-1
10

10 -2
0 500 1000 1500 2000 2500

We first compute the singular value decomposition of this image data

A = UΣVT . We then compute low-rank approximations using the first
k singular values and left/right vectors pairs. We show the reconstructed
images for k = (1, 2, 3, 5, 10, 15, 20, 30, 50) in the figure below.

82
 k=1 k=2 k=3

k=5 k = 10 k = 15

k = 20 k = 30 k = 50

In the example above, the full matrix A had 5550602 elements. The recon-
structed image with k = 50 requires us to store only the first 50 left and right
singular vectors, plus the 50 largest singular values. This amounts to about
238700 elements to store, or about 4% of the original space requirement.

83
6 Norms and conditioning

6.1 Vector norms

When we talk about the ‘size’ or ‘magnitude’ of a real number, we mean its
absolute value. Although the number −6.2 is algebraically smaller than the
number 2.8, most of us would agree that −6.2 is ‘bigger’ than 2.8, in the sense
that it lies further away from zero. The idea can be expressed more formally:
the length of the vector running from the origin to −6.2 is greater than that
of the vector running from the origin to 2.8.

What about higher dimensions? We know how to calculate the length of

vectors in R2 and R3 using Pythagoras’s theorem, and there is a natural
extension to Rn :
h iT q
2
in R : length of x1 x2 = x21 + x22
h iT q
3
in R : length of x1 x2 x3 = x21 + x22 + x23
h iT q
n
in R : length of x1 x2 . . . xn = x21 + x22 + · · · + x2n .

This general definition also works when n = 1, since

h i q
length of x1 = x21 = |x1 |.

Going back to the example in the opening paragraph, we get the desired
outcome: the length of the single-component vector [−6.2] is greater than
that of [2.8].

We are used to calculating the length of a vector in this way, but in fact a

84
whole family of ‘length’ measures can be defined as follows:
h iT 1
2
in R : length of x1 x2 = (|x1 |p + |x2 |p ) p
h iT 1
in R3 : length of x1 x2 x3 = (|x1 |p + |x2 |p + |x3 |p ) p
h iT 1
in Rn : length of x1 x2 . . . xn = (|x1 |p + |x2 |p + · · · + |xn |p ) p .

This generalized measure of length is called the p–norm (or Lp norm). To

denote the p–norm of a vector x we use the special notation ∥x∥p . We never
write |x|p , though the similarity is of course intentional; taking the norm of
a vector gives us a scalar measure of its magnitude, and thus is analogous to
taking the absolute value of a real (or complex) number. The p–norm can be
defined compactly as follows:

Definition (Vector p–norm) Suppose that x ∈ Rn and p ≥ 1 is a real

scalar. Then the p–norm of x is

n
! p1
X
∥x∥p = |xi |p
i=1

It is not essential that p be an integer, but in practice it almost always is. In

fact, there are only a few norms in common use, namely the ones that are
easy to calculate: those for p = 1, p = 2 and p = ∞.

The familiar Pythagorean measure of length is just the p–norm with p = 2. It

is commonly referred to as the Euclidean norm, or the Euclidean metric (note
that in this context, metric means “measure of size” – nothing to do with SI
units). The 2–norm can also be expressed in a compact vector form, since it
is the square root of the dot product of x with itself, i.e.
√ √
∥x∥2 = x·x= xT x.

It is sometimes simpler (e.g. when minimizing the length of a vector, see Q7

85
of the tutorial sheet) to work with the square of the 2-norm:

∥x∥22 = x · x = xT x.

If we set p = 1 in our general definition, we get

X
∥x∥1 = |xi |
i

i.e. the sum of the absolute values of the components. The 1–norm is some-
times referred to as the ‘taxicab’ or ‘Manhattan’ norm, since a hypothetical
taxi driving from (0; 0) to a destination with coordinates (x1 ; x2 ) has to cover
a distance of at least |x1 | + |x2 | if it is constrained to drive on a rectangular
grid of streets in R2 .

Finally, if we let p approach infinity in the general definition, we get

∥x∥∞ = max |xi |,

i

i.e. the maximum of the absolute values of the components. You are asked to
justify this in Q3 of the tutorial sheet.

It is also possible to define vector norms other than the p-norm, but any pro-
posed norm must possess certain properties in order to qualify as a meaningful
measure of vector size. These are similar to the properties of the absolute
value function. Specifically, we define a vector norm as being any function
that satisfies the following conditions:

Definition (Vector norms) The function ∥·∥ is a norm in Rn if it satisfies

the following properties for all x and y in Rn and all α ∈ R:

i) ∥x∥ = 0 if and only if x = 0

ii) ∥αx∥ = |α|∥x∥ (homogeneity)

iii) ∥x + y∥ ≤ ∥x∥ + ∥y∥ (triangle inequality)

86
It is easy to prove that these basic properties imply another property that we
would expect any measure of vector size to exhibit:

∥x∥ ≥ 0 for all x.

P P
For example, the 1-norm i |xi | is a valid vector norm, but i xi is not, since
the sum of the terms in a vector could be negative.

A good way to visualize vector norms is to consider what the unit circle (or
unit ball) looks like in R2 looks under the three common p–norms:

1 1 1

1 1 1

kxk1 = 1
<latexit sha1_base64="xzedbfybVbe3G7Ke0nCsXNgqpAI=">AAAHdnichVXdbts2FFbbLW69v7S9HFAQE7rtwjUkx1mCAgHa7KIbhmJZtqQFTM2gpCObC0VpJJXa5fQsu90eaW+yyx3K8hLLRUNA0NH5vvPLIzIuBdcmCP65dfvOBx/u9O7e63/08SeffrZ7/8G5LiqVwFlSiEK9jpkGwSWcGW4EvC4VsDwW8Cq++Nbhry5BaV7IX8yyhChnM8kznjCDqunuQyoLlVuaxxlZ1NOQHJFwuusHw6BZZFsIW8H32nUyvb/zA02LpMpBmkQwrSdhUJrIMmV4IqDu00pDyZILNoMJipLloAfpJS91I0a2KaQmjxFMSVYofKQhjfa6sWW51ss8RmbOzFx3Mad8FzapTHYYWS7LyoBMVoGyShBTENcVknIFiRFLFFiiOKZNkjlTLDHYuw1Pb3hq5kfj4T6XkZ1BkYNRy40KLcgq5wbyut+nKWS4M00hNl8eiwpqe/riuLYHB4PD/cFo75u6Q2JGsPiFApAtMxiEh8EgwKRdY5Iid23eLL1kSl/wEgMSqkDCG8diMrU0YzkXSwzAKmFqS3W2ljG568QLWMrCwHElBDiiYlxDXCxsMNyDhbPEvpRG87dADSxMvCLSAZ1rTANsiCznVINxcztxHUByZLEcEEfdCN34MXdRY5hxaVtT5GxQoKEAileETYbBHcMGUzcARWlpUYJiplBuyqxRtfO42Z+0KfadfAd1U8iTrKUnTHSxeI2hQLpYvMbirRxOm5Tj2J5uVfw77lZtnx5t2WDWv1JTlB0DUxbaeTMdPZeXaGCfhN0IuFv1ZBRZKiAz1BI/pKSR/yCU+CMcJz6bm2H7pl1zd3jUk9DZn4MyaK3ceyu6BFWqIm1DMTkTYP2wHlh/VOOouW/cy85u/liaE1XgSYYh9qNJZPuEtBPCBJ9JSGvaDJdzhdgVqhRb1laJldYfUcHxj9RT6+/VX64H9kmAE+uP6eXV/FK6cuMG3Ooq/g0PBTwjnl4zcrSNWG0mExM1agy3v8LdoK7R/vq7yay/idXdwl9y+X/hYyy8f70XE9+1O+d4PGDzXEXWH2/9Cd/L7CYXXGbvdfGSLW7Mgi3e6+LnqrzJha7Kjgu8fsLuZbMtnI+GYTAMfxr7z563F9Fd73PvC+9rL/QOvGfed96Jd+Yl3tL70/vL+3vn396j3uPeVyvq7VutzUNvY/WC/wD+yaLw</latexit>
kxk2 = 1
<latexit sha1_base64="XHQS0svNJ7monaKPBz1J2v7B7So=">AAAHdnichVXdbts2FFbbLW69v7S9HFAQE7rtwjUkx1mCAgHa7KIbhmJZtqQFTM2gpCObC0VpJJXa5fQsu90eaW+yyx3K8hLLRUNA0NH5vvPLIzIuBdcmCP65dfvOBx/u9O7e63/08SeffrZ7/8G5LiqVwFlSiEK9jpkGwSWcGW4EvC4VsDwW8Cq++Nbhry5BaV7IX8yyhChnM8kznjCDqunuQyoLlVuaxxlZ1NMROSLhdNcPhkGzyLYQtoLvtetken/nB5oWSZWDNIlgWk/CoDSRZcrwREDdp5WGkiUXbAYTFCXLQQ/SS17qRoxsU0hNHiOYkqxQ+EhDGu11Y8tyrZd5jMycmbnuYk75LmxSmewwslyWlQGZrAJllSCmIK4rJOUKEiOWKLBEcUybJHOmWGKwdxue3vDUzI/Gw30uIzuDIgejlhsVWpBVzg3kdb9PU8hwZ5pCbL48FhXU9vTFcW0PDgaH+4PR3jd1h8SMYPELBSBbZjAID4NBgEm7xiRF7tq8WXrJlL7gJQYkVIGEN47FZGppxnIulhiAVcLUlupsLWNy14kXsJSFgeNKCHBExbiGuFjYYLgHC2eJfSmN5m+BGliYeEWkAzrXmAbYEFnOqQbj5nbiOoDkyGI5II66EbrxY+6ixjDj0ramyNmgQEMBFK8ImwyDO4YNpm4AitLSogTFTKHclFmjaudxsz9pU+w7+Q7qppAnWUtPmOhi8RpDgXSxeI3FWzmcNinHsT3dqvh33K3aPj3assGsf6WmKDsGpiy082Y6ei4v0cA+CbsRcLfqySiyVEBmqCV+SEkj/0Eo8Uc4Tnw2N8P2Tbvm7vCoJ6GzPwdl0Fq591Z0CapURdqGYnImwPphPbD+qMZRc9+4l53d/LE0J6rAkwxD7EeTyPYJaSeECT6TkNa0GS7nCrErVCm2rK0SK60/ooLjH6mn1t+rv1wP7JMAJ9Yf08ur+aV05cYNuNVV/BseCnhGPL1m5GgbsdpMJiZq1Bhuf4W7QV2j/fV3k1l/E6u7hb/k8v/Cx1h4/3ovJr5rd87xeMDmuYqsP976E76X2U0uuMze6+IlW9yYBVu818XPVXmTC12VHRd4/YTdy2ZbOB8Nw2AY/jT2nz1vL6K73ufeF97XXugdeM+877wT78xLvKX3p/eX9/fOv71Hvce9r1bU27dam4fexuoF/wEF2qLx</latexit>
kxk1 = 1
<latexit sha1_base64="5VTNUW2YGJ2bKQzqr8OzWdzPnQE=">AAAHe3ichVXdjttEFHYLbEr425YbJG4sLBBCaWRnk+6q0kpluSgIVSwLu62UcaOxfZwMO54xM+NtwmCehlt4Hx4GiTOOw26cqjuS5eP5vvM7x2eSkjNtwvCfO3ffevudvd69d/vvvf/Bhx/t339woWWlUjhPJZfqRUI1cCbg3DDD4UWpgBYJh+fJ5TcOf34FSjMpfjarEuKCzgXLWUoNbs32PyFCqsKSIsn9ZT0jTORm5R/70Ww/CIdhs/xdIWqFwGvX6ez+3vckk2lVgDApp1pPo7A0saXKsJRD3SeVhpKml3QOUxQFLUAPsitW6kaMbZNN7X+OYObnUuEjjN/s3lS2tNB6VSTILKhZ6C7mNl+HTSuTH8WWibIyINK1o7zivpG+K42fMQWp4SsUaKoYhu2nC6poarCAW5ZescwsjsfDCROxnYMswKjVVoYWRFUwA0Xd75MMcjyeJhFbrE54BbU9e3pS28PDwdFkMDp4VHdI1HCaPFUAomWGg+goHIQYtCtMKgtX5u3US6r0JSvRoU8UCHjlWFRkluS0YHyFDmjFTW2JzjcyBneTeAkrIQ2cVJyDIyrKNCRyacPhASydJtalNJr9BsTA0iRrIhmQhcYwwEbIckY1GNe8U1cBJMcW0wF+3PXQ9Z8w5zWBORO2VUXOFgUaCqB4TdhmGDwxLDBxDSBLS2QJihqpXJdZo2pncbs+WZPsa/kO6oZQpHlLTynvYskGQ8HvYskGS3ZiOGtCThJ7tpPxr3hatX18vKODUb8kRpYdBVNK7ayZzj4TV6hgH0ZdD3ha9XQUW8IhN8T6QUT8Rv7dJ34wwnZi84UZtm/SVXcTpJ5GTv8ClEFt5d473gWoUsmsdUXFnIMNonpgg1GNrea+8Sw7p/lDaU6VxHGGLibxNLZ93287hHI2F5DVpGkuZwqxa1Qpuqqt4uvdYEQ4wz9Sz2xwUH+xadiHIXZsMCZX1/1LyNqMa3Crq+QXHAo4Ix7fUHK0LV9tJFMTN9vobrLGXaNu0P7mu4msv43V3cSfMfF/4mNMvH+zFtPAlbtgOB6weC4jG4x3/oTvRH6bCZz3bzTxjC5vjYIu32jip6q8zYSuyo4JvH6i7mWzK1yMhlE4jH4cB0++bi+ie96n3mfel17kHXpPvG+9U+/cS70/vD+9v7y/9/7tBb2veoM19e6dVudjb2v1Jv8BTqyldw==</latexit>

Example Find the 1-, 2-, 5-, and ∞–norms of the vector x = (3; −5; 4; −1).

Solution

∥x∥1 = |3| + | − 5| + |4| + | − 1| = 13

1 √
∥x∥2 = |3|2 + | − 5|2 + |4|2 + | − 1|2 2 = 51 = 7.14

1
∥x∥5 = |3|5 + | − 5|5 + |4|5 + | − 1|5 5 = 5.35
∥x∥∞ = max(|3|, | − 5|, |4|, | − 1|) = 5.

6.2 Matrix norms

We can also extend the notion of a norm from vectors to matrices. Intuition
suggests that a matrix such as
" #
100 −200
75 400

87
must somehow be ‘bigger’ than a matrix such as
" #
1 3
.
−1 2

One obvious approach would be to place all the entries of the matrix into
a single vector, then apply the vector p–norm introduced in the preceding
section. If we took p = 2, this would be equivalent to squaring each entry of
the matrix, summing and taking the square root. This is indeed a commonly
used matrix norm known as the Frobenius norm, which we can write as
sX X
∥A∥F = A2i,j .
i j

Another approach, and one that turns out to be much more useful, is to think
about the action of a matrix on an arbitrary vector. As we noted in the first
lecture, an m × n matrix A can be viewed as a linear operator that maps a
vector x ∈ Rn to a vector Ax ∈ Rm . If we take the norm of the ‘output’
vector Ax and divide it by the norm of the ‘input’ vector x, we should get
some measure of the ‘magnifying power’ of the matrix A (and thus its inherent
‘size’ or ‘magnitude’). We are therefore motivated to calculate the ratio

∥Ax∥p
∥x∥p

assuming of course that x is not the zero vector. The problem is that, unless
A happens to be the identity matrix (or some multiple thereof), this ratio will
depend on the direction of the input vector x. To use a very simple example,
consider the diagonal matrix " #
5 0
.
0 2
The vector x = (1; 0) is mapped to Ax = (5; 0) (magnification = 5), but
the vector (0; 1) is mapped to (0; 2) (magnification = 2). Vectors in other

88
directions undergo intermediate levels of magnification. This can be clarified
by looking at the image of the unit circle under the action of A. It is mapped
to an ellipse with semi-axes of length 5 and 2:
3

-1

-2

-3
-6 -5 -4 -3 -2 -1 0 1 2 3 4 5 6

Given that a general matrix A is likely to exert widely varying ‘magnifying

powers’ on different input vectors x, we take a worst-case approach and define
the matrix p–norm of A as

∥Ax∥p
∥A∥p = max .
x̸=0 ∥x∥p

With this definition, the matrix p–norm provides an upper bound on the ‘mag-
nifying power’ of A (as measured by applying the vector p–norm to both the
input vector x and output vector Ax). It means that we are assured of sat-
isfying the following important inequality for any vector x, regardless of its
direction:

Definition (Induced matrix norms) Let A be a matrix with n columns,

and let x be a vector in Rn . If the inequality

∥Ax∥ ≤ ∥A∥∥x∥

is satisfied for all vectors x ∈ Rn , then the matrix norm ∥A∥ is said to be

89
 compatible with, or induced by, the vector norm used as a metric for x
and Ax.

There is no general formula for calculating the p–norm of a matrix A in terms

of p and the entries Ai,j . For selected values of p, however, it is possible
to work out formulae for ∥A∥p that guarantee the satisfaction of the above
‘compatibility inequality’ for an arbitrary vector x.
6.2.1 Induced matrix norm for p = 1
For p = 1, the compatible matrix norm is
!
X
∥A∥1 = max |Ai,j |
j
i

which is the maximum absolute column sum of A (remember that i counts

along the rows, and j counts along the columns). In order to prove this we
must verify that the inequality

∥Ax∥1 ≤ ∥A∥1 ∥x∥1

90
is actually satisfied for every vector x. Recalling that the indices i and j are
understood to run from 1 . . . m and 1 . . . n respectively, we can write

X X
∥Ax∥1 = Ai,j xj
i j
XX
≤ |Ai,j xj |
i j
XX
= |xj | |Ai,j |
i j
!
X X
= |xj | |Ai,j |
j i
!
X X
≤ |xj | · max |Ai,j |
j
j i

= ∥A∥1 ∥x∥1 ,

where the second step holds due to the scalar version of the triangle inequality.
6.2.2 Induced matrix norm for p = ∞
For p = ∞, the compatible matrix norm is
!
X
∥A∥∞ = max |Ai,j |
i
j

which is the maximum absolute row sum of A. To verify the validity of this
definition, we need to show that the inequality

∥Ax∥∞ ≤ ∥A∥∞ ∥x∥∞

91
holds for all x. In this case the proof of our claim is a bit shorter, and we can
write
X
∥Ax∥∞ = max Ai,j xj
i
j
X
≤ max |Ai,j | |xj |
i
j
!
X
≤ max |Ai,j | max |xj |
i j
j

= ∥A∥∞ ∥x∥∞ ,
where we have one again applied the scalar triangle inequality in the second
step.

Note that you should also be able to verify that both of the above definitions
are tight, in the sense that, for any matrix A, you can construct a correspond-
ing vector x for which the norm compatibility condition is actually satisfied
with equality.
6.2.3 Induced matrix norm for p = 2
The matrix 2–norm is widely used and important because it is compatible with
the 2–norm for vectors, i.e. it ensures that

∥Ax∥2 ≤ ∥A∥2 ∥x∥2

for any vector x. The 2-norm of A is the square root of the largest eigenvalue
of the symmetric matrix AT A, i.e.
q
∥A∥2 = λmax (AT A).

In other words, the 2–norm of a matrix A is just its largest singular value. If
A happens to be square and symmetric to begin with, then the 2–norm is just
the maximum of the absolute values of the eigenvalues:

∥A∥2 = max(|λi |) (if A is symmetric)

i

92
You might like to consider this for the unit circle to ellipse mapping discussed
above.

Example Find the 1–norm, ∞–norm and Frobenius norm of the matrix
 
3 −2 6
A =  5 −8 4 .
 

−4 1 −3

Solution The 1–norm is the maximum absolute column sum:

∥A∥1 = max(12, 11, 13) = 13.

The ∞–norm is the maximum absolute row sum:

∥A∥∞ = max(11, 17, 8) = 17.

The Frobenius-norm is the square root of the sum of the squares of the
elements:

12 √
∥A∥F = 32 + (−2)2 + · · · + (−3)2 = 180 = 13.42.

Example Given the matrix and vector

" # " #
1.8 1.4 0.3
A= , x=
−0.1 0.6 0.2

investigate the satisfaction of the ‘compatibility inequality’

∥Ax∥p ≤ ∥A∥q ∥x∥p

when p and q take values 1 and ∞ in all combinations.

93
 Solution We have
" #" # " #
1.8 1.4 0.3 0.82
Ax = = .
−0.1 0.6 0.2 0.09

We therefore have

∥x∥1 = 0.5, ∥x∥∞ = 0.3 (vector norms of x)

∥A∥1 = 2.0, ∥A∥∞ = 3.2 (matrix norms of A)
∥Ax∥1 = 0.91, ∥Ax∥∞ = 0.82 (vector norms of Ax)

Combination ∥Ax∥p ≤ ∥A∥q ∥x∥p Compatibility satisfied?

p = 1, q = 1 0.91 ≤ 2.0 · 0.5 YES (expected)
p = ∞, q = ∞ 0.82 ≤ 3.2 · 0.3 YES (expected)
p = 1, q = ∞ 0.91 ≤ 3.2 · 0.5 YES (lucky)
p = ∞, q = 1 0.82 ≤ 2.0 · 0.3 NO (expected)

The conclusion to be drawn from the preceding example is that one should
not mix and match incompatible matrix and vector norms.

Example Find the 2–norm of the matrix

" #
1 2
A= .
2 −3

Solution This matrix is square and symmetric, so its 2–norm can be cal-
culated by finding the eigenvalues. Solve the matrix characteristic equation

0 = det(A − λI)
= (1 − λ)(−3 − λ) − 4
= λ2 + 2λ − 7.

94
 √
The roots of this quadratic equation are λ = −1 ± 8, hence
 √ √ 
∥A∥2 = max −1 + 8 , −1 − 8 = 3.83.

6.3 Conditioning of linear systems

Consider the 2 × 2 linear system of equations
" #" # " #
1 1 x1 2
= ,
1 −1 x2 0

which as the (unique) solution

" # " #
x1 1
= .
x2 1

We would now like to consider how sensitive this solution is to small changes
in the problem data in either the matrix A or the right-hand side b. If we
perturb one of the left-hand term’s coefficients by 1% to
" #" # " #
1 1.01 x1 2
= ,
1 −1 x2 0

then the solution changes to

" # " #
x1 0.995
= .
x2 0.995

This seems reasonable – both components of x have changed by 0.5%, which

is of the same order of magnitude as the perturbation.

Similarly, if we perturb the first term of the right hand side vector by 1% to
get " #" # " #
1 1 x1 2.02
= ,
1 −1 x2 0

95
then the solution changes to
" # " #
x1 1.01
= ,
x2 1.01

which is again a change of roughly the same size as the perturbation to our
problem data. This type of ‘nicely-behaved’ linear system said to be well-
conditioned, meaning that its solution is not strongly sensitive to small varia-
tions in the data.

Consider now a different 2 × 2 system

" #" # " #
1 1 x1 2
= ,
0.99 1.01 x2 2

which also happens to have the solution

" # " #
x1 1
= .
x2 1

In this case, when we perturb one of the terms on the left-hand side by 1% to
get the new system " #" # " #
1 1.01 x1 2
= ,
0.99 1.01 x2 2
we get a drastic change in the solution, which becomes
" # " #
x1 0
= .
x2 1.98

Likewise, if we perturb the right-hand side of the same system by 1%, to

" #" # " #
1 1 x1 2.02
= ,
0.99 1.01 x2 2

96
we get another drastic change in the solution, which becomes
" # " #
x1 2.01
= .
x2 0.01

This type of behaviour can be surprising and is clearly undesirable: the solution
x is highly sensitive to small perturbations in the problem data A and b
(such as those that might arise from noisy data, imperfect measurements, or
numerical roundoff errors during a computation). This type of ‘badly behaved’
linear system is said to be ill-conditioned.

We can gain some insight into what is happening if we plot the solutions of
our two scalar equations in each of the two situations (well-conditioned on the
left, ill-conditioned on the right):

0 0

-0.25 -0.25

-0.5 -0.5

-0.75 -0.75

-1 -1

-1.25 -1.25

-1.5 -1.5

-1.75 -1.75

-2 -2
0 0.25 0.5 0.75 1 1.25 1.5 1.75 2 0 0.25 0.5 0.75 1 1.25 1.5 1.75 2

The reason for the ill-conditioned behaviour is now clear. The lines in the
right-hand plot are so close to parallel that even small changes in slope and/or
intercept cause the intersection point to move far away from the nominal
solution point (1; 1). Note that the lines being nearly parallel is equivalent to
the rows of A (and hence also the columns!) being nearly linearly dependent,
i.e. the matrix A is nearly singular (i.e. nearly not full rank).

If the lines were parallel, then the determinant of the coefficient matrix would

97
be zero, which might suggest using det(A) as an indicator of possible condi-
tioning trouble, i.e. as an indicator that the matrix is nearly singular. In the
preceding example the ill-conditioned example has
" #
1 1
A= ,
0.99 1.01

with det(A) = 0.02, which is close to zero. Unfortunately, the determinant is

not a useful measure of ill-conditioning because a simple rescaling of the prob-
lem leads to a change in the determinant. If we multiply all of the equations
of the ill-conditioned system by 10, we get
" #" # " #
10 10 x1 20
= ,
9.9 10.1 x2 20

so that det(A) = 2. However, the conditioning has not changed, and the
solution remains just as sensitive to variations in A or b as it was before.

Fortunately, there is a better approach for quantifying the sensitivity of a

linear system Ax = b, thus allowing ill-conditioned behaviour to be predicted
in advance (i.e. without resorting to trial and error perturbations of the sort we
tried above). The perturbations to A and b are characterized by the relative
error ratios
∥δA∥ ∥δb∥
and
∥A∥ ∥b∥.
We would like to be able to predict the effect on the solution via the corre-
sponding relative error ratio
∥δx∥
.
∥x∥
Any convenient norm can be used, but it is essential that the matrix norm and
vector norm be ‘compatible’ in the sense defined in Section 6.2.

Suppose we have a linear system defined by a matrix A, a right-hand side

vector b, and an exact solution x, such that Ax = b. We will assume that

98
the matrix A is square and invertible.

If the matrix is perturbed from A to A + δA, the solution will be perturbed

from x to x + δx, such that

(A + δA)(x + δx) = b.

Subtract Ax from the left and b from the right to obtain

Aδx + δA(x + δx) = 0,

then left-multiply by A−1 and rearrange to obtain

δx = −A−1 δA(x + δx).

Now take the norm of both sides and apply the compatibility inequality twice
to obtain
∥δx∥ = ∥A−1 δA(x + δx)∥
≤ ∥A−1 ∥∥δA(x + δx)∥
≤ ∥A−1 ∥∥δA∥∥x + δx∥.
Hence
∥δx∥
≤ ∥A−1 ∥∥δA∥,
∥x + δx∥
or, equivalently,
∥δx∥
∥δA∥
≤ ∥A∥∥A−1 ∥ . (6.1)
∥x + δx∥ ∥A∥

Similarly, if the right-hand side vector is perturbed from b to b + δb, the

solution will be perturbed from x to x + δx, such that

A(x + δx) = b + δb.

Subtract Ax from the left and b from the right to obtain

Aδx = δb,

99
then left-multiply by A−1 to obtain

δx = A−1 δb.

Now take the norm of both sides and apply the compatibility inequality to
obtain
∥δx∥ = ∥A−1 δb∥
≤ ∥A−1 ∥∥δb∥.
Multiplying by ∥b∥ on the left and ∥Ax∥ on the right, then applying the
compatibility inequality once more, we obtain

∥b∥∥δx∥ ≤ ∥A−1 ∥∥δb∥∥A∥∥x∥,

then finally rearrange to

∥δx∥
∥δb∥
≤ ∥A∥∥A−1 ∥ . (6.2)
∥x∥ ∥b∥

Note that the bracketed term ∥A∥∥A−1 ∥ appears in both (6.1) and (6.2). It is
called the condition number of A, and is usually denoted κ(A). The condition
number is always at least 1, and very large values indicate an ill-conditioned
system.

Definition (Matrix condition number) Let A be a square invertible

matrix. The condition number of A is defined as

κ(A) = ∥A∥∥A−1 ∥.

For a linear system Ax = b, the condition number provides an upper

bound on the relative error in x caused by relative errors in either A or b,
i.e.

∥δx∥ ∥δA∥
≤ κ(A)
∥x + δx∥ ∥A∥

100
 and

∥δx∥ ∥δb∥
≤ κ(A) .
∥x∥ ∥b∥

Note that in the preceding definition we have not specified which matrix norm
we are using to define the condition number. The most common choice by
far is the matrix 2–norm, which is the default, for example, in the matlab
condition number function cond. However, any matrix norm can be used,
and in this course we will typically use either the matrix 1– or ∞–norm since
calculation of the general matrix 2–norm is generally more difficult. For the
sake of clarity, we can denote the condition number defined in terms of the
matrix p–norm as κp (·).
Example Find the condition numbers for the two systems considered at
the start of this section. Also, for the ill-conditioned system, verify that the
actual relative errors caused by the perturbations in A and b are within the
bounds predicted by the condition number. Use the matrix infinity norm
when computing the condition number.

Solution For the well-conditioned system,

" # " #
1 1 0.5 0.5
A= and A−1 = ,
1 −1 0.5 −0.5

hence
κ∞ (A) = ∥A∥∞ ∥A−1 ∥∞ = 2 · 1 = 2.

For the ill condition system, we have instead

" # " #
1 1 50.5 −50
A= and A−1 = ,
0.99 1.01 −49.5 50

hence
κ∞ (A) = ∥A∥∞ ∥A−1 ∥∞ = 2 · 100.5 = 201.

101
In the ill-condition case, we perturbed the matrix from
" # " #
1 1 1 1.01
A= to A + δA = ,
0.99 1.01 0.99 1.01

so that " #
0 0.01
δA = .
0 0
This caused the solution to change from
" # " # " #
1 0 −1
x= to x + δx = , with δx = .
1 1.98 0.98

We therefore expect

1 ∥δx∥∞ ∥δA∥∞ 0.01

= ≤ κ∞ (A) = 201 · ,
1.98 ∥x + δx∥∞ ∥A∥∞ 2

which turns out to hold.

Considering instead the change to the right-hand side, we changed

" # "#
2 2.02
b= to b + δb = ,
2 2

so that " #
0.02
δb = .
0
This caused the solution to change from
" # " # " #
1 2.01 1.01
x= to x + δx = , with δx = .
1 0.01 −0.99

102
We therefore expect

1.01 ∥δx∥∞ ∥δb∥∞ 0.02

= ≤ κ∞ (A) = 201 · ,
1 ∥x∥∞ ∥b∥∞ 2

which again turns out to hold.

103
7 Iterative methods for Ax = b

7.1 Introduction
Earlier in the course we saw how to reduce the linear system Ax = b to echelon
form using elementary row operations. Solution methods that rely on this
strategy (e.g. LU factorization) are robust and efficient, and are fundamental
tools for solving the systems of linear equations that arise in practice. These
are known as direct methods, since the solution x is obtained following a single
pass through the relevant algorithm. We will next consider some alternative
approaches that fall into the category of iterative or indirect methods. These
techniques can only be applied to square linear systems (n equations in n
unknowns), but this is of course a common and important case.

Iterative methods for Ax = b begin with a (possibly very poor) approximation

to the solution, x0 , then seek to provide a series of improved approximations
x1 , x2 , . . . that converge to the exact solution. For the engineer, this approach
is appealing because it can be stopped as soon as the approximations xi have
converged to an acceptable precision, which might be something as crude as
10−3 . With a direct method, early termination of the algorithm is not an
option; the process of elimination and back-substitution has to be carried out
to completion, or else abandoned altogether. By far the main attraction of
iterative methods, however, is that for certain problems – particularly those
where the matrix A is large and sparse (containing mostly zeros) – they can
be much faster than direct methods. On the other hand, iterative methods
can be unreliable; for some problems they may exhibit very slow convergence,
or they may not converge at all.

104
7.2 Jacobi method (‘simultaneous displacements’)
The Jacobi method is the simplest iterative method for solving a (square) linear
system Ax = b. Before developing a general formulation of the algorithm, it
is instructive to explain the basic workings of the method with reference to a
small example. Consider the system
    
4 2 3 x 8
 3 −5 2 y  = −14 . (7.1)
    

−2 3 8 z 27

Noting that there are no zeros on the leading diagonal, we can solve the first
row for x, the second for y and the third for z:

x = (8 − 2y − 3z)/4
y = (−14 − 3x − 2z)/(−5) (7.2)
z = (27 + 2x − 3y)/8.

This rearrangement is the key to the method, since it puts the equations in a
form that is amenable to iteration. At each stage of the process, new values of
x, y and z will be obtained by substituting the old values into the expressions
on the right-hand side of (7.2). In other words, we will apply the simple
iterative scheme
xn+1 = (8 − 2yn − 3zn )/4
yn+1 = (−14 − 3xn − 2zn )/(−5) (7.3)
zn+1 = (27 + 2xn − 3yn )/8.

In the absence of better information, we will begin with x0 = y0 = z0 = 0 as

105
our initial approximation. This makes the first iteration easy:

x1 = (8 − 2 · 0 − 3 · 0)/4 =2
y1 = (−14 − 3 · 0 − 2 · 0)/(−5) = 2.8
z1 = (27 + 2 · 0 − 3 · 0)/8 = 3.375.

The second iteration gives

x2 = (8 − 2 · 2.8 − 3 · 3.375)/4 = −1.931

y2 = (−14 − 3 · 2 − 2 · 3.375)/(−5) = 5.350
z2 = (27 + 2 · 2 − 3 · 2.8)/8 = 2.825.

The third iteration gives

x3 = (8 − 2 · 5.350 − 3 · 2.825)/4 = −2.794

y3 = (−14 − 3 · (−1.931) − 2 · 2.825)/(−5) = 2.771
z3 = (27 + 2 · (−1.931) − 3 · 5.350)/8 = 0.886

And so on. Although the early iterations don’t look highly promising in terms
of convergence, things do eventually settle down:

Iteration x y z Iteration x y z
1 2.000 2.800 3.375 26 -1.003 2.984 1.989
2 -1.931 5.350 2.825 27 -0.984 2.994 2.005
3 -2.794 2.771 0.886 28 -1.001 3.012 2.006
4 -0.050 1.478 1.637 29 -1.011 3.002 1.995
5 0.033 3.425 2.808 30 -0.997 2.992 1.997
6 -1.819 3.943 2.099 31 -0.993 3.000 2.004
7 -1.546 2.548 1.442 32 -1.003 3.006 2.002
8 -0.355 2.449 2.033 33 -1.004 2.999 1.997
9 -0.749 3.400 2.368 34 -0.997 2.996 1.999
10 -1.476 3.297 1.913 35 -0.998 3.001 2.002

106
 Iteration x y z Iteration x y z
11 -1.083 2.680 1.770 36 -1.002 3.002 2.000
12 -0.667 2.858 2.099 37 -1.001 2.999 1.999
13 -1.003 3.240 2.137 38 -0.998 2.999 2.000
14 -1.222 3.053 1.909 39 -1.000 3.001 2.001
15 -0.958 2.830 1.925 40 -1.001 3.001 2.000
16 -0.859 2.995 2.074 41 -1.000 2.999 1.999
17 -1.053 3.114 2.037 42 -0.999 3.000 2.000
18 -1.085 2.983 1.944 43 -1.000 3.001 2.000
19 -0.949 2.926 1.985 44 -1.001 3.000 2.000
20 -0.952 3.024 2.040 45 -1.000 3.000 2.000
21 -1.042 3.045 2.003 46 -1.000 3.000 2.000
22 -1.025 2.976 1.973 47 -1.000 3.000 2.000
23 -0.967 2.974 2.003 48 -1.000 3.000 2.000
24 -0.989 3.021 2.018 49 -1.000 3.000 2.000
25 -1.024 3.014 1.995 50 -1.000 3.000 2.000

It can be verified by substitution in (7.1) that the solution is indeed

   
x −1
y  =  3 .
   

z 2

7.3 Gauss-Seidel method (‘successive displacements’)

The Gauss-Seidel method is a variant of the Jacobi method that usually im-
proves the rate of convergence. Recall that in the previous section, the pro-
cedure for conducting a Jacobi iteration was

xn+1 = (8 − 2yn − 3zn )/4

yn+1 = (−14 − 3xn − 2zn )/(−5) (7.3)
zn+1 = (27 + 2xn − 3yn )/8.

107
Note that in the second line, when working out yn+1 , we used the old value
of x (i.e. xn ) even though a new and presumably more accurate value (i.e.
xn+1 ) had just been worked out in the first line. Similarly, in the third line,
when working out zn+1 , we used the old values of x and y (i.e. xn and yn )
even though updated values (i.e. xn+1 and yn+1 ) had just been worked out in
previous lines. The Gauss-Seidel method implements the strategy of always
using the latest available value of a particular variable. For our small example
it therefore gives the iterative scheme

xn+1 = (8 − 2yn − 3zn )/4

yn+1 = (−14 − 3xn+1 − 2zn )/(−5) (7.4)
zn+1 = (27 + 2xn+1 − 3yn+1 )/8.

These expressions should be compared very carefully with their Jacobi coun-
terparts in (7.3), where only ‘old’ variable values (from iteration n) appear on
the right-hand sides. In the Gauss-Seidel method, we use ‘new’ variable values
(from iteration n + 1) wherever possible.

To clarify the operation of the Gauss-Seidel method, we will go through the

first few iterations of the example, again starting from x0 = y0 = z0 = 0 as
the initial approximation. The first iteration gives

x1 = (8 − 2 · 0 − 3 · 0)/4 =2
y1 = (−14 − 3 · 2 − 2 · 0)/(−5) = 4
z1 = (27 + 2 · 2 − 3 · 4)/8 = 2.375.

which is already different from before. The second iteration gives

x2 = (8 − 2 · 4 − 3 · 2.375)/4 = −1.781
y2 = (−14 − 3 · (−1.781) − 2 · 2.375)/(−5) = 2.681
z2 = (27 + 2 · (−1.781) − 3 · 2.681)/8 = 1.924.

108
The third iteration gives

x3 = (8 − 2 · 2.681 − 3 · 1.924)/4 = −0.784

y3 = (−14 − 3 · (−0.784) − 2 · 1.924)/(−5) = 3.099
z3 = (27 + 2 · (−0.784) − 3 · 3.099)/8 = 2.017.

Bearing in mind that the exact solution is x = −1, y = 3, z = 2, things are

looking a great deal better than they were at the same stage of the Jacobi
iterations, where we had x3 = −2.794, y3 = 2.771 and z3 = 0.886. If we
carry on with the Gauss-Seidel iterations we obtain:
Iteration x y z
1 2.000 4.000 2.375
2 -1.781 2.681 1.924
3 -0.784 3.099 2.017
4 -1.062 2.969 1.996
5 -0.982 3.009 2.001
6 -1.006 2.997 2.000
7 -0.998 3.001 2.000
8 -1.001 3.000 2.000
9 -1.000 3.000 2.000
10 -1.000 3.000 2.000

Convergence to three decimal places in all components is achieved after just 9

iterations, whereas the Jacobi method took 45 iterations to achieve the same
level of accuracy. It would be wrong, however, to assume that Gauss-Seidel is
always superior to Jacobi; occasionally it is worse!

The Gauss-Seidel method can seem a more natural approach than the Jacobi
method. If the solution is converging and updated information is available for
some of the variables, surely it makes sense to use that information. From
a programming point of view, the Gauss-Seidel method can also seem more
convenient since the old value of a variable can be overwritten as soon as

109
a new value becomes available. With the Jacobi method, the values of all
variables from the previous iteration must be retained throughout the current
iteration, which means that twice as much storage is needed.

On the other hand, the Jacobi method is perfectly suited to parallel compu-
tation, whereas the Gauss-Seidel method is not. Because the Jacobi method
updates or ‘displaces’ all of the variables at the same time (at the end of each
iteration) it is often called the method of simultaneous displacements. The
Gauss-Seidel method updates the variables one by one (during each iteration)
so its corresponding name is the method of successive displacements. This
means that, although the Jacobi method may require more total computations
(because it needs more iterations), it may actually converge faster in terms of
computation time, because every variable can be updated simultaneously at
each iteration.

7.4 Compact Formulations

It is useful to be able to represent the Jacobi and Gauss-Seidel methods in
terms of matrices and vectors. In fact, it turns out that both iterative schemes
can be written (at least analytically – this is not necessarily how they are
programmed!) in the compact form

xn+1 = Pxn + q (7.5)

where P is a constant matrix and q is a constant vector. Our task is to derive

expressions for P and q in terms of the input data: the coefficient matrix A
and the right-hand side vector b.

To help establish these general formulae, we will organize our thoughts by

looking at a 3 × 3 system Ax = b:
    
A11 A12 A13 x b1
A21 A22 A23  y  = b2  .
    

A31 A32 A33 z b3

110
Solving the first row for x, the second for y and the third for z gives the basic
template
x = (b1 − A12 y − A13 z)/A11
y = (b2 − A21 x − A23 z)/A22
z = (b3 − A31 x − A32 y)/A33 .
Each iteration of the Jacobi method (cf. (7.3)) takes the form

xn+1 = (b1 − A12 yn − A13 zn )/A11

yn+1 = (b2 − A21 xn − A23 zn )/A22
zn+1 = (b3 − A31 xn − A32 yn )/A33 .

This can be rearranged as

A11 xn+1 = b1 − A12 yn − A13 zn

A22 yn+1 = b2 − A21 xn − A23 zn
A33 zn+1 = b3 − A31 xn − A32 yn

or, in matrix from,

       
A11 xn+1 0 A12 A13 xn b1
A22   yn+1  = − A21 0 A23   yn  + b2  .
       

A33 zn+1 A31 A32 0 zn b3

Thus, in general, the Jacobi method involves iterations of the form

Dxn+1 = −(L + U)xn + b (7.6)

where D is a 3 × 3 matrix containing the diagonal part of A, and L and U are

3 × 3 matrices containing the (strictly) lower and (strictly) upper triangular
parts of A. (Note that the L and U matrices here have nothing to do with
those in the factorization A = LU.) Pre-multiplying by D−1 gives

xn+1 = −D−1 (L + U)xn + D−1 b (7.7)

111
which is in the standard form of (7.5).

Each iteration of the Gauss-Seidel method (cf. (7.4)) takes the form

xn+1 = (b1 − A12 yn − A13 zn )/A11

yn+1 = (b2 − A21 xn+1 − A23 zn )/A22
zn+1 = (b3 − A31 xn+1 − A32 yn+1 )/A33 .

This can be rearranged as

A11 xn+1 = b1 − A12 yn − A13 zn

A22 yn+1 = b2 − A21 xn+1 − A23 zn
A33 zn+1 = b3 − A31 xn+1 − A32 yn+1

or, in matrix form,

          
A11 xn+1 0 0 0 xn+1 0 A12 A13 xn b1
A22   yn+1  = − A21 0 0  yn+1  − 0 0 A23   yn  + b2  .
          

A33 zn+1 A31 A32 0 zn+1 0 0 0 zn b3

Thus, in general, the Gauss-Seidel method involves iterations of the form

Dxn+1 = −Lxn+1 − Uxn + b (7.8)

where D, L and U are as defined above. Pre-multiplying by D−1 gives

xn+1 = −D−1 Lxn+1 − D−1 Uxn + D−1 b.

At first sight it looks strange to have xn+1 appearing on the right-hand side,
until we realize that the matrix −D−1 L is (strictly) lower triangular. Thus,
only part of the vector xn+1 is needed to process any given row, and that
part is always available (from the results of previous rows). To represent the
Gauss-Seidel method more compactly, we need to go back to (7.8) and bring

112
all terms related to xn+1 to the left-hand side:

(D + L)xn+1 = −Uxn + b. (7.9)

This expression should be compared with its Jacobi counterpart in eqn (7.6).
Pre-multiplying by (D + L)−1 gives

xn+1 = −(D + L)−1 Uxn + (D + L)−1 b (7.10)

which is in the standard form of (7.5).

We can summarize the key equations (7.6, 7.7, 7.9, 7.10) as follows:

Jacobi Gauss-Seidel
Dxn+1 = −(L + U)xn + b (D + L)xn+1 = −Uxn + b
xn+1 = Pxn + q xn+1 = Pxn + q
where where
P = −D−1 (L + U) P = −(D + L)−1 U
q = D−1 b q = (D + L)−1 b

In fact, both methods can be derived more quickly (though without such
obvious connections to the basic algorithms as explained in Sections 7.1 and
7.2) by starting directly from the additive decomposition

A = D + L + U.

The system Ax = b becomes

(D + L + U)x = b,

and it is then simply a question of how much we leave on the left-hand side
when setting up our iterative scheme.

113
If we leave just Dx on the left-hand side, we get

Dx = −(L + U)x + b

=⇒ Dxn+1 = −(L + U)xn + b,

which is the Jacobi method.

If we leave (D + L)x on the left-hand side, we get

(D + L)x = −Ux + b

=⇒ (D + L)xn+1 = −Uxn + b,

which is the Gauss-Seidel method.

Both iterative schemes are special cases of the more general

(D + ωL)xn+1 = − ((1 − ω)L + U) xn + b. (7.11)

With ω = 0 we recover the Jacobi method. With ω = 1 we recover the Gauss-

Seidel method. We could, if we prefer, take ω = 0.5 (say) to obtain a method
that lies ‘between’ Jacobi and Gauss-Seidel. We could even take ω > 1,
resulting in a method that goes ‘beyond’ Gauss-Seidel. This latter strategy
takes us into the realm of over-relaxation, and for certain problems it turns
out to be highly effective. It is possible to recast (7.11) in such a way that
the matrices on the left- and right-hand sides are lower and upper triangular
respectively. This allows over-relaxation to be implemented in a manner similar
to the Gauss-Seidel method, with new variable values overwriting old ones as
soon as they become available. This type of iterative method is known as
successive over-relaxation (SOR).

Example Derive the Jacobi iteration matrix P and iteration vector q for
the example used in Section 7.2.

114
Solution From (7.1) we have Ax = b with
   
4 2 3 8
A =  3 −5 2 , b = −14 .
   

−2 3 8 27

We first split A into

   
4 0 2 3
D =  −5 and L + U =  3 0 2 .
   

8 −2 3 0

Noting that D−1 is trivial to calculate from D (we just take the reciprocal
of each term on the diagonal) we can easily work out

P = −D−1 (L + U)
1    
4 0 2 3 0 −0.5 −0.75
= − − 51   3 0 2 =  0.6 0 0.4
    
1
8 −2 3 0 0.25 −0.375 0

and

q = D−1 b
1    
4 8 2
= − 15  −14 =  2.8  .
    
1
8 27 3.375

The Jacobi iterations are given by xn+1 = Pxn + q, i.e.

      
xn+1 0 −0.5 −0.75 xn 2
 yn+1  =  0.6 0 0.4  yn  +  2.8  .
      

zn+1 0.25 −0.375 0 zn 3.375

115
As expected, the preceding solution is identical to the iterative scheme derived
directly in Section 7.2 (by making x the subject of the first equation, etc.).
Recall that previously we had

xn+1 = (8 − 2yn − 3zn )/4

yn+1 = (−14 − 3xn − 2zn )/(−5) (7.3)
zn+1 = (27 + 2xn − 3yn )/8,

which is exactly the same iteration.

In our present and more compact form, the first few Jacobi iterations look like
(assuming that we start from all zeros):
        
x1 0 −0.5 −0.75 0 2 2
 y1  =  0.6 0 0.4 0 +  2.8  =  2.8 
        

z1 0.25 −0.375 0 0 3.375 3.375

        
x2 0 −0.5 −0.75 2 2 −1.913
 y2  =  0.6 0 0.4  2.8  +  2.8  =  5.35 
        

z2 0.25 −0.375 0 3.375 3.375 2.825

        
x3 0 −0.5 −0.75 −1.913 2 −2.794
 y3  =  0.6 0 0.4  5.35  +  2.8  =  2.771  .
        

z3 0.25 −0.375 0 2.825 3.375 0.886

The results are, of course, the same as those obtained previously.

7.5 Convergence criteria

Suppose we have a (square) linear system Ax = b with the exact solution
x = s. If we apply an iterative scheme in the standard form

xn+1 = Pxn + q (7.12)

and the iterations converge, we will – by definition – eventually reach a steady

state or fixed point in which each iteration simply regenerates the exact solu-

116
tion, i.e.
s = Ps + q. (7.13)

If we subtract (7.13) from (7.12) the q terms disappear and we obtain

(xn+1 − s) = P(xn − s)

or equivalently
en+1 = Pen (7.14)

where we have introduced the notation en to denote the error in the nth iterate
xn with respect to the exact solution:

en = xn − s.

We can now analyze how the solution evolves from one iteration to the next.
For convergence, we require the magnitude of the error vector e (as measured
by some vector norm) to approach zero as the iterations proceed. Taking
vector norms in (7.14) gives

∥en+1 ∥ = ∥Pen ∥ (7.15)

and applying the ‘compatibility inequality‘ (see §6.2) gives

∥en+1 ∥ ≤ ∥P∥∥en ∥. (7.16)

This is a very important relationship. If we can establish that the iteration

matrix P has the property
∥P∥ < 1

for any convenient induced matrix norm, then we can be certain that the
iterative scheme will converge. The reason is that if ∥P∥ is strictly less than
1, then (7.16) implies
∥en+1 ∥ < ∥en ∥

117
for all n. (Note: here we tacitly assume the use of a vector norm that is
compatible with the matrix norm used to establish ∥P∥ < 1). Even if our
initial approximation x0 has a large error (as measured by ∥e0 ∥ = ∥x0 − s∥,
i.e. we start very far away from the solution), we can use the compatibility
condition multiple times to establish that

∥en ∥ ≤ ∥P∥n ∥e0 ∥,

so that ∥en ∥ → 0 as n → ∞, and consequently

lim xn = s.
n→∞

If possible, we would like to find a matrix norm ∥P∥ that is not just less than 1,
but a good deal less than 1. If ∥P∥ = 0.1, for example, (7.16) shows that
each iteration will reduce ∥e∥ by a factor of at least 10. We know in advance
that the iterations will converge very quickly.

If we choose a particular matrix norm, say the ∞–norm, and find that ∥P∥
is greater than 1, this does not necessarily indicate that the iterative scheme
will fail to converge. For a start, there may be some other matrix norm (such
as the 1–norm) that is strictly less than 1, in which case convergence is still
guaranteed. In any case, however, the condition ∥P∥ < 1 is only a sufficient
condition for convergence, not a necessary one, for arbitrarily chosen norms.

However, the spectral radius of P (the maximum of the absolute values of the
eigenvalues) provides us with both a sufficient and necessary condition:

Theorem (Convergence of linear fixed-point iteration) The iterative

scheme xn+1 = Pxn + q converges to a fixed point if and only if every
eigenvalue λ of P satisfies |λ| < 1, i.e. if and only if the spectral radius
ρ(P) < 1.

118
 Caution: the spectral radius of a general square matrix is not a norm,
although it happens to agree with the matrix 2–norm for symmetric ma-
trices.

The idea of over-relaxation (see §7.4) is to choose the factor ω in such a way
that the spectral radius ρ(P) is made as small as possible, thereby giving the
fastest possible rate of convergence. A similar idea is explored in Q6 of the
tutorial sheet, in the context of a symmetric matrix (such that the matrix
2-norm ∥P∥ and the spectral radius ρ(P) are the same).

Of course, matrix norms such as the ∞–norm and 1–norm are much easier
to calculate than eigenvalues, so in practice we might first try to establish
convergence using ∥P∥∞ or ∥P∥1 .

Sometimes a guarantee of convergence can be established by direct inspection

of the coefficient matrix A (i.e. without needing to compute the iteration
matrix P). In particular, if A has the following property, then the Jacobi and
Gauss-Seidel schemes are both certain to converge:

Definition (Diagonal dominance) The matrix A is diagonally domi-

nant if, in each row, the absolute value of the entry on the diagonal is
greater than the sum of the absolute values of the other entries. More
compactly, A is diagonally dominant if
X
Ai,i > |Ai,j |, ∀i
j̸=i

In the case of the Jacobi scheme, it is not too hard to show that if A is
diagonally dominant then the iteration matrix P has an ∞–norm that is strictly
less than 1 (therefore guaranteeing convergence). From the formula

P = −D−1 (L + U)

119
it follows that 
−A /A when j ̸= i
i,j i,i
Pi,j =
0 when j = i.
To understand where these expressions come from, it is instructive to write
out P = −D−1 (L + U) in full:
   1  
P1,1 P1,2 . . . P1,n A1,1 0 A1,2 ... A1,n
1 .. .. 
 P2,1 P2,2 . . . P2,n 
   
 A2,2
  A2,1 0 . . 
 .
 . .. .. ..  = − 
..   ..

.. ..

 . . . . 
 
 .  . . . An−1,n 

1
Pn,1 Pn,2 . . . Pn,n An,n An,1 ... An,n−1 0

 A1,2 A1,n

0 A1,1
... A1,1
A
 2,1 ..
0 ... .
= −  A2,2
 
 ... ... ... An−1,n 

 An−1,n−1 
An,1 An,n−1
An,n An,n
... 0

The proof itself is then straightforward. Since A is assumed diagonally domi-

nant, we have X
|Ai,i | > |Ai,j |, ∀i
j̸=i
 
X |Ai,j |
=⇒ 1 > max   = ∥P∥∞ .
i |Ai,i |
j̸=i

Example Show that for each of the following matrices A, the system
Ax = b can be solved by Jacobi iteration with guaranteed convergence.
     
5 −1 3 −2 0 4 4 2 −2
(a) :  2 −8 1 (b) :  2 −8 1 (c) : 0 4 2
     

−2 0 4 5 −1 3 1 0 4

120
Solution a) This matrix is diagonally dominant since

|5| > | − 1| + |3| (first row)

| − 8| > |2| + |1| (second row)
|4| > | − 2| + |0| (third row)

From the above theorem we therefore know that the Jacobi iteration matrix
P must have an ∞–norm that is strictly less than 1. This is also easy to
verify directly by writing

P = −D−1 (L + U)
    
0.2 0 −1 3 0 0.2 −0.6
= − −0.125   2 0 1 = 0.25 0 0.125
    

0.25 −2 0 0 0.5 0 0

Hence ∥P∥∞ = 0.8, which (as expected) is < 1.

b) This matrix is not diagonally dominant since the entry in the top left-
hand corner (= −2) is not bigger than the rest of its row (with a similar
problem in the third row). However, the first and third equations of Ax =
b can be swapped to give a new system which is diagonally dominant; in
fact it is the matrix of part (a). The reordered system is suitable for Jacobi
iteration.

c) This matrix is not diagonally dominant since the upper left-hand term
is equal to the sum of the remaining terms in the row (we need strict
satisfaction of the inequality in each row), and reordering the rows will not
help.

121
We therefore compute the Jacobi iteration matrix

P = −D−1 (L + U)
    
0.25 0 2 −2 0 −0.5 0.5
= − 0.25  0 0 2  =  0 0 −0.5
    

0.25 1 0 0 −0.25 0 0

and see what its norms look like. The row-sum norm ∥P∥∞ and the
column-sum norm ∥P∥1 are both 1. We are still looking for an induced
matrix norm of P that is strictly less than 1, so we turn to the 2–norm.
Direct calculation of the 2-norm of an asymmetric matrix would require us
to compute the eigenvalues of AT A, but we can instead exploit a theorem
that says ∥M∥2 ≤ ∥M∥F for any matrix M. It is easy to compute ∥P∥F =
0.901 (check this for yourself) and it follows that ∥P∥2 ≤ 0.901 < 1,
so convergence is certain. An alternative approach is to calculate the
eigenvalues of P and show that ρ(P) < 1 but this involves a good deal
more work.

122
Example (A1 2002 Q3) (a) Define the Lp norm of a vector x. Deter-
mine the L1 , L2 and L∞ norms of the vector (2; 0.2; −5). Give an
example of an application where the L1 norm would be more applica-
ble than the L2 norm.

(b) Find ∥A∥∞ for the matrix

" #
17 2
A=
10 −16

(c) Describe the Jacobi algorithm for solution of the system of equations
Ax = b. State a condition on the norm of a matrix (not necessarily
A), which must be satisfied for the method to converge.

(d) On the basis of the appropriate L∞ matrix norm, choose which of the
equations below is more suitable for solution using Jacobi iteration.
Obtain an approximate solution using three iterations starting from
x = (0.5; 0.5) and comment on the solution.

" #" # " # " #" # " #

3 5 x1 2 3 2 x1 2
= , or = .
3 −4 x2 1 3 −4 x2 1

123
8 Computing Eigenvalues and Eigenvectors

8.1 Introduction
Let us first recall some terminology from an earlier lecture. If A is a square
(n × n) matrix, then a nonzero vector x ∈ Rn is an eigenvector of A if

Ax = λx (8.1)

for some scalar λ. The scalar λ is an eigenvalue of A, and we say that

x is an eigenvector corresponding to λ (or alternatively that x and λ are
an eigenvector-eigenvalue pair ). Geometrically, the n × n matrix A defines
a linear transformation from Rn → Rn . If the ‘output’ vector Ax points
in the same direction as the ‘input‘ vector x, then x is an eigenvector of
the transformation. The eigenvalue λ is the factor (which could be positive,
negative or zero) by which x is scaled to give Ax.

As you know, eigenvalues can be found by hand provided the matrix A is small,
e.g. n = 3. We rearrange (8.1) as (A − λI)x = 0 and, for non-trivial x, we
require the determinant of A − λI to be zero. This leads to the characteristic
polynomial (an nth degree polynomial in λ, e.g. a cubic when n = 3), the
roots of which are the eigenvalues.

If A is large, however, this approach is no longer an option. Given that

in practice we could be dealing with huge matrices (n = 100, 000 or more),
even finding the coefficients of the characteristic polynomial is computationally
very expensive. Instead we move to numerical methods that can find a small
number of eigenvalues cheaply. Typically it is only the largest few eigenvalues,

124
or else the smallest few, that are of interest in any given problem.

In this brief treatment, we will deal only with matrices that

• are real and symmetric;

• have a dominant eigenvalue.

Definition (Dominant Eigenvalue) An eigenvalue of an n × n matrix

A is called the dominant eigenvalue of A if its absolute value is strictly
greater than the absolute values of all remaining n − 1 eigenvalues. The
eigenvector corresponding to the dominant eigenvalue is called the domi-
nant eigenvector.

8.2 Power method

The power method is a simple but effective way to find the dominant eigen-
value and eigenvector of a matrix A. It starts with a (nonzero) estimate x0
for the dominant eigenvector, then iterates according to

xn+1 = Axn .

Let’s consider how well this works by way of a simple example. Suppose that
" #
−1 2
A= .
2 2

In order to assess the numerical results, we first find the exact solution. We
can set
0 = det(A − λI)
" #
−1 − λ 2
= det
2 2−λ
= (−1 − λ)(2 − λ) − 4
= (λ + 2)(λ − 3).

125
Hence A has a dominant eigenvalue equal to 3 since that eigenvalue is largest
inn magnitude. We can find the corresponding eigenvector by solving
" #" # " #" # " #
(−1 − 3) 2 x1 −4 2 x1 0
= = ,
2 (2 − 3) x2 2 −1 x2 0

which has (non-normalized) solution (x1 , x2 ) = (1; 2). (Note that we could
treat any non-zero scaling of this vector as the eigenvector corresponding to
the same eigenvalue.).

Our goal, however, is to see how the power method does in identifying this
same vector. We choose an arbitrary non-zero point to get started, say x =
(1; 1), then start iterating to obtain:

Iteration x Ax Magnitude ratio

" # " #" # " #
1 −1 2 1 1 4
1 = 1
= 4.000
1 2 2 1 4
" # " #" # " #
1 −1 2 1 7 10
2 = 4
= 2.500
4 2 2 4 10
" # " #" # " #
7 −1 2 7 13 34
3 = 10
= 3.400
10 2 2 10 34
" # " #" # " #
13 −1 2 13 55 94
4 = 34
= 2.765
34 2 2 34 94
" # " #" # " #
55 −1 2 55 133 298
5 = 94
= 3.170
94 2 2 94 298
" # " #" # " #
133 −1 2 133 463 862
6 = 298
= 2.893
298 2 2 298 862
" # " #" # " #
463 −1 2 463 1261 2650
7 = 862
= 3.074
862 2 2 862 2650

The method clearly works, in the sense that the magnitude ratio (which we

126
have defined as just the ratio of the largest component of x between subse-
quent iterations) appears to be converging. The magnitude of x is growing
steadily (this can easily be fixed if required, see below) but its direction is
converging to that of the dominant eigenvector (1; 2). This can be seen if we
look at the ratios of the components in the last few iterates:
" #
133
x5 = , 298/133 = 2.241
298
" #
463
x6 = , 862/463 = 1.862
862
" #
1261
x7 = , 2650/1261 = 2.102.
2650

Remember that eigenvectors are only unique up to a scalar multiple – it is

the direction that counts. Even if we normalize the eigenvector, there are still
two possibilities (x̂ and −x̂). Strictly speaking, we should not talk about ‘the’
dominant eigenvector, but ‘a’ dominant eigenvector, since there are infinitely
many of different lengths.

Note: we have not yet explained why x converges to the direction of the
dominant eigenvector. This is discussed later.

As for the dominant eigenvalue, the final column of the table tracks the evo-
lution of the ‘magnification ratio’, defined for convenience as the factor by
which the largest (in this case the second) component of x is scaled by an
iteration1 . As we work through the iterations, we see that the magnification
ratio is approaching the dominant eigenvalue (= 3).

The above results are reasonable since we know that A scales its dominant
eigenvector (1; 2) by a factor of exactly 3, so as x approaches the direction
(1; 2) in the iterations of the power method we expect A to scale it by a
1
Here the other component is of course magnified by a similar ratio, but in general it is advisable to
monitor the magnification of the largest element (in absolute magnitude terms). This ensures that we don’t
inadvertently monitor the noisy fluctuations of a component that is converging towards zero.

127
factor that gradually approaches 3. Note, however, that in any given iteration
prior to convergence, the individual components of x are not scaled by exactly
the same factor under the action of A (e.g. in the final row of the table,
2650/862 = 3.074, but 1261/463 = 2.724). If the components were all
scaled by the same factor, this would imply that x was an eigenvector of A,
and the iterations would have already converged.

The power method is fine for hand calculations on small matrices, but in
practice the steady growth of the components of x cannot be allowed because
it will eventually lead to floating-point overflow (or underflow, if the dominant
eigenvalue happens to be less than 1). Using IEEE 64-bit arithmetic, we can
make it to " #
307
3.3217 × 10
x645 = ,
6.6434 × 10307
before overflow. This many iterations may seem more than sufficient, and for
this problem it is (note how the ratio of the ratio of the second component to
the first component has reached 2), but in real problems the power method
may be very slow to converge. We will see why shortly, but for now what we
need is a method for keeping the magnitude of x under control.
The solution is simple. At each iteration, before computing Ax, we first
normalize the vector x by its ∞–norm (or some other convenient vector norm).
This is the power method with scaling. It makes sure that as the iterations
proceed, no component of x acquires an absolute value greater than 1. When
the method is applied to the above example with the same starting vector x0 ,
we get:

Iteration x Ax Magnitude ratio

" # " # " #" # " #
1 1 −1 2 1 1
1 → = 4.000
1 1 2 2 1 4
" # " # " #" # " #
1 0.25 −1 2 0.25 1.75
2 → = 2.500
4 1 2 2 1 2.5

128
 Iteration x Ax Magnitude ratio
" # " # " #" # " #
1.75 0.7 −1 2 0.7 1.3
3 → = 3.400
2.5 1 2 2 1 3.4
" # " # " #" # " #
1.3 0.382 −1 2 0.382 1.618
4 → = 2.765
3.4 1 2 2 1 2.765
" # " # " #" # " #
1.618 0.585 −1 2 0.585 1.415
5 → = 3.170
2.765 1 2 2 1 3.170
" # " # " #" # " #
1.415 0.446 −1 2 0.446 1.554
6 → = 2.893
3.170 1 2 2 1 2.893
" # " # " #" # " #
1.554 0.537 −1 2 0.537 1.463
7 → = 3.074
2.893 1 2 2 1 3.074

This involves slightly more work but is numerically better conditioned because
the size of the iterates is well regulated. The components of x remain bounded,
and at each iteration the magnification ratio can be read off directly from the
largest component of Ax (note that the magnification ratios are identical
to those that appeared in the unscaled power method). Re-normalizing at
each iteration also makes it more obvious that x is converging to a particular
direction: that of the dominant eigenvector (1; 2) ∝ (0.5; 1). Other scaling
strategies could have been implemented, e.g. normalizing x by its 2–norm to
obtain an eigenvector of unit length. But this would have involved extra work
to find the 2–norm, and we are trying to keep each iteration as cheap as
possible.

We now need to address two important questions: why does the power method
work, and what controls its rate of convergence? Consider an n × n matrix A
with a dominant eigenvalue λ1 , such that

|λ1 | > |λ2 | ≥ |λ3 | ≥ · · · ≥ |λn |. (8.2)

129
Note that the leftmost inequality is strict by assumption. Since A is real and
symmetric, we know that we can always find n corresponding eigenvectors

{v1 , v2 , v3 , . . . , vn }

that are linearly independent (in fact, as you may recall from P1, we can find
n corresponding eigenvectors that are mutually orthogonal). This means that
any vector x0 that is used to start the power method can be expressed as a
linear combination of the eigenvectors of A:

x0 = c1 v1 + c2 v2 + c3 v3 + · · · + cn vn . (8.3)

The first iteration gives

x1 = Ax0
= c1 Av1 + c2 Av2 + c3 Av3 + · · · + cn Avn
= c1 λ1 v1 + c2 λ2 v2 + c3 λ3 v3 + · · · + cn λn vn .

The second iteration gives

x2 = Ax1
= c1 λ1 Av1 + c2 λ2 Av2 + c3 λ3 Av3 + · · · + cn λn Avn
= c1 λ21 v1 + c2 λ22 v2 + c3 λ23 v3 + · · · + cn λ2n vn .

After the k th iteration we will have

xk = c1 λk1 v1 + c2 λk2 v2 + c3 λk3 v3 + · · · + cn λkn vn . (8.4)

Since the eigenvalue λ1 has been assumed dominant in (8.2), it is clear that
as k increases the first term on the right-hand side of eqn (8.4) will begin to
dominate all of the remaining terms. As this happens, x will gradually ap-
proach the direction of the dominant eigenvector v1 (and in the limit, become

130
parallel to it). Taking another view, we can write (8.4) as
"  k  k  k #
λ2 λ3 λn
x2 = λk1 c1 v1 + c2 v2 + c3 v3 + · · · + cn vn .
λ1 λ1 λ1
(8.5)
Since λ1 is dominant, the ratios λi /λ1 have absolute values that are strictly
less than one for all i ≥ 2. The corresponding powers (λi /λ1 )k will therefore
all tend to zero as k increases. This again explains why the power method
gradually ‘turns’ the vector x until it becomes parallel to the dominant eigen-
vector v1 .

Furthermore, (8.5) reveals that the rate of convergence of the power method
will depend on the ratio of the two largest eigenvalues, λ2 /λ1 . If |λ2 /λ1 | is
small, say 0.1, then the influence of the ‘second most dominant’ eigenvector
v2 (and all subsequent eigenvectors) will quickly tail off as k increases; con-
vergence will be very fast. If is only a little smaller than 1, e.g. 0.95, then
convergence will be much slower. In practice this difficulty can be partially
overcome by using the shifted power method, which temporarily offsets the
whole spectrum of eigenvalues in order to make the ratio |λ2 /λ1 | smaller.

A final important observation is that the power method can fail. If we are
unlucky enough to choose a starting vector x0 that has no component in the
direction of the dominant eigenvector v1 (i.e. x0 is orthogonal to v1 ) then the
coefficient c1 in the initial decomposition (8.4) will be zero. If c1 = 0 in (8.4)
and (8.5), it is clear that the iterations of the power method will cause x to
approach the direction of v2 , the second most dominant eigenvector, rather
than v1 . Also, the ‘magnification ratio’ in each iteration will converge to λ2
rather than λ1 . A reasonable precaution in practice is to try two different
starting vectors x0 that are not parallel to each other.

8.3 Rayleigh method

There is a very effective way to improve the estimate of the dominant eigen-
value λ1 at each iteration of the power method. So far we have effectively

131
been making a ‘direct’ estimate of λ1 by monitoring the factor by which the
(absolutely) largest component of x is magnified under the action of A. As x
gradually turns towards the dominant eigenvector v1 , this magnification ratio
approaches λ1 .

Another approach is to revisit the fundamental eigenvector/eigenvalue equa-

tion
Ax = λx. (8.1)

We rearrange this to read

Ax − λx = 0 (8.6)

noting that the 0 on the right-hand side is a vector of zeros.

When the power method eventually converges, (8.6) will be satisfied exactly
(with x = v1 and λ = λ1 ). Prior to convergence, however, x is only an
approximation of the dominant eigenvector, so the best we can hope for is to
find an approximation µ ≈ λ1 of the dominant eigenvalue such that

Ax − µx ≈ 0.

To do this rigorously, we should try to choose µ in such a way that some

norm of the ‘residual’ or error vector Ax − µx is minimized. It is most natural
to do this minimization with respect to the ordinary Euclidean norm (i.e. the
2-norm) since it is differentiable, in which case the optimum value of µ turns
out to be the Rayleigh quotient.

Definition (Rayleigh quotient) Let A be a real symmetric matrix with

a dominant eigenvalue λ1 , and let x be a given approximation of the
dominant eigenvector v1 . The Rayleigh quotient, defined as

xT Ax
r(x) =
xT x

provides an approximation of λ1 that is optimal in a least-squares sense,

132
 i.e. the residual error ∥Ax − µx∥22 is minimized by taking µ = r(x).

This is an important result, and you are asked to prove it in Q7 of the tutorial
sheet. A hint is provided there to get you started, and here is another one:
use the fact that
(Ax − µx)T = xT AT − xT µ,

then expand out the product with the other bracketed term. Don’t forget that
A is assumed symmetric, i.e. A = AT .

If we apply the power method as before, but in addition evaluate the Rayleigh
quotient at each iteration, we have the Rayleigh method for finding the dom-
inant eigenvalue (and the corresponding eigenvector). We will not go into
the details, but it can be shown that the power method provides only linear
convergence to the dominant eigenvalue, whereas the Rayleigh method pro-
vides quadratic convergence (so it always gives better performance). With
a more sophisticated scheme known as the shifted inverse power method, it
is even possible to obtain cubic convergence. Note however that the use of
the Rayleigh quotient does not improve the rate of convergence of x to the
dominant eigenvector v1 ; this remains the same as in the standard power
method.

It is worth pointing out that the Rayleigh quotient can also be expressed in
terms of dot products:

xT Ax x · (Ax)
r(x) = T = .
x x x·x

The normalization by x · x ensures that we obtain a correct estimate of the

eigenvalue even if x is not a unit vector.

To show the effectiveness of using the Rayleigh quotient, here is the first
example (no scaling) again, with an extra column appended:

133
 xT Ax
Iteration x Ax Magnitude ratio r(x) = xT x
" # " #" # " #
1 −1 2 1 1 4 5
1 = 1
= 4.000 2
= 2.500
1 2 2 1 4
" # " #" # " #
1 −1 2 1 7 10 47
2 = 4
= 2.500 17
= 2.765
4 2 2 4 10
" # " #" # " #
7 −1 2 7 13 34 431
3 = 10
= 3.400 139
= 2.893
10 2 2 10 34
" # " #" # " #
13 −1 2 13 55 94 3911
4 = 34
= 2.765 1325
= 2.952
34 2 2 34 94
" # " #" # " #
55 −1 2 55 133 298 35327
5 = 94
= 3.170 11861
= 2.978
94 2 2 94 298
" # " #" # " #
133 −1 2 133 463 862 318455
6 = 298
= 2.893 106493
= 2.990
298 2 2 298 862
" # " #" # " #
463 −1 2 463 1261 2650 2868143
7 = 862
= 3.074 957413
= 2.996
862 2 2 862 2650

As expected, the Rayleigh quotient provides a far superior rate of convergence

to the dominant eigenvalue, λ1 = 3.

134
 35
Power
Rayleigh

30

Percent error (absolute) 25

20

15

10

0
1 2 3 4 5 6 7
iteration number

8.4 Matrix deflation

We have already seen a result saying that a real symmetric matrix A can
always be factored as
A = MDMT , (8.7)

where  
λ1
 
 λ2 
D=
 ...


 
λn
is a matrix containing the eigenvalues of A on its diagonal, and
h i
M = v̂1 v̂2 . . . v̂n

is a matrix containing an orthonormal set of corresponding eigenvectors, ar-

ranged columnwise. The matrix M is orthogonal and has rank n (i.e. it satisfies
MMT = MT M = I).

135
We can expand (8.7) in an alternative form

A = λ1 v̂1 v̂1T + λ2 v̂2 v̂1T + · · · + λn v̂n v̂nT ,

that will be more useful for the derivations to come. We called this the spectral
decomposition or eigendecomposition of A. It shows that a real symmetric
matrix in Rn×n can be expressed as the sum of n component matrices, each of
which is formed by taking the outer product of a normalized eigenvector with
itself, then applying a ‘weighting factor’ equal to the corresponding eigen-
value. Note that all of the component matrices are themselves symmetric
(and, incidentally, have rank 1 since every column is a multiple of the same
eigenvector).

We need to establish one more important result. Suppose we have a real

symmetric n × n matrix A and form the spectral decomposition

A = λ1 v̂1 v̂1T + λ2 v̂2 v̂2T + · · · + λn v̂n v̂nT .

If we move the first component matrix λ1 v̂1 v̂1T to the left-hand side, we obtain
a new matrix
B = A − λ1 v̂1 v̂1T
(8.8)
= 0 · v̂1 v̂1T + λ2 v̂2 v̂2T + · · · + λn v̂n v̂nT .
From the second line, we expect the ‘deflated’ matrix B to have the same
eigenvectors and eigenvalues as A, but with zero replacing λ1 in the list of
eigenvalues. That is, we expect (v̂1 , v̂2 , . . . , v̂n ) to be eigenvectors of B
corresponding to eigenvalues (0, λ2 , . . . , λn ) respectively.

This is not too difficult to prove. If we right-multiply both sides of (8.8) by v̂1
we get
Bv̂1 = Av̂1 − λ1 v̂1 v̂1T v̂1
= λ1 v̂1 − λ1 v̂1 (v̂1T v̂1 )
= λ1 v̂1 − λ1 v̂1 = 0,

136
noting that v̂1T v̂1 = 1. This shows that v̂1 is an eigenvector of B with
eigenvalue 0 (whereas it was an eigenvector of A with eigenvalue λ1 ), or
equivalently that v1 is in the kernel of B.

If we post-multiply both sides of (8.8) by any other v̂i with i ̸= 1, we get

Bv̂i = Av̂i − λ1 v̂1 v̂1T v̂i

= λi v̂i − λ1 v̂1 (v̂1T v̂i )
= λi v̂i

where we have used the fact that v̂iT v̂j = 0 when i ̸= j (i.e. the eigenvectors
are mutually orthogonal). This shows that v̂i is an eigenvector of B with
eigenvalue λi (just as it was an eigenvector of A) with eigenvalue λi , and this
holds for all i = 2, . . . , n.

Finally we are ready to outline the technique of matrix deflation. Suppose we

have a real symmetric matrix A ∈ Rn×n with eigenvalues

|λ1 | > |λ2 | > |λ3 | ≥ · · · ≥ |λn |.

Note that we have now assumed that the leftmost two inequalities are strict.
Let us further assume that the dominant eigenvalue λ1 and dominant eigen-
vector v1 have been found using the power or Rayleigh method. The dominant
eigenvector v1 can be normalized to give a unit dominant eigenvector v̂1 (this
step is crucial!) and the ‘deflated’ matrix

B = A − λ1 v̂1 v̂1T

can then be formed. As we know from above, the eigenvalues of B are the
same as those of A, but with 0 replacing λ1 , i.e.

|λ2 | > |λ3 | ≥ · · · ≥ |λn | > 0.

The power method or Rayleigh method can now be applied to B. This will

137
give the dominant eigenvalue and eigenvector of B, namely λ2 and v2 , which
are the ‘second most dominant’ eigenvalue and eigenvector of the original
matrix A. The deflation process can, in principle, be repeated by forming

C = B − λ2 v̂2 v̂2T

and applying the power method or Rayleigh method to C, and so on, until
all n eigenvalues and eigenvectors of the original matrix A have been found
(assuming that they are all distinct). In practice, however, numerical roundoff
errors build up, and the accuracy of the computed results starts to deteriorate.
The deflation technique is therefore typically only suitable for obtaining the
largest few eigenvalues of A and their corresponding eigenvectors.

If the smallest few eigenvalues are required (e.g. the fundamental natural
frequency and first few harmonics of a dynamic system) then it is much more
efficient to use a technique known as the inverse power method. This method
computes the (absolutely) smallest eigenvalue of A by finding the dominant
eigenvalue of A−1 , but does so without requiring the inverse to be formed
explicitly. If all eigenvalues of a matrix are required, one popular choice is
the QR algorithm, which is based on the matrix QR decomposition. These
methods are beyond the scope of the present course.

We will finish this section by demonstrating the application of deflation to our

example problem involving the matrix
" #
−1 2
A=
2 2

We assume that the power method or Rayleigh method has been applied
to A and continued to convergence, such that its dominant eigenvalue and
eigenvector are known exactly:
" #
1
λ1 = 3, v1 = .
2

138
We then normalize v1 to make it a unit vector (crucial step!):
" #
1 1
v̂1 = √ ·
5 2

and form the ‘deflated’ matrix

B = A − λ1 v̂1 v̂1T
" # " #h i
−1 2 1 1 1 2
= −3·
2 2 5 2
" #
−1.6 0.8
= .
0.8 −0.4

The eigenvalues and eigenvectors of B are easy to find by setting

0 = det(B − λI)
" #
−1.6 − λ 0.8
= det
0.8 −0.4 − λ
= (1.6 + λ)(0.4 + λ) − 0.64
= λ(λ + 2),

so that the eigenvalues are (0, −2). This is exactly as expected. We happen
to know (from §8.2) that the eigenvalues of A are 3 and −2, so the dominant
eigenvalue (i.e. λ1 = 3) has been correctly ‘altered’ to zero in the list of eigen-
values of the deflated matrix B. The dominant eigenvalue of B is λ2 = −2,
and the corresponding eigenvector can be found in the usual way by solving
" #" # " #" # " #
(−1.6 + 2) 0.8 x1 0.4 0.8 x1 0
= =
0.8 (−0.4 + 2) x2 0.8 1.6 x2 0

139
which gives the (non-normalized) dominant eigenvector of B as
" # " #
x1 −2
= .
x2 1

We can therefore deduce that the ‘second most dominant’ eigenvalue and
eigenvector of the original matrix A must be
" #
−2
λ2 = −2, v2 = .
1

This can readily be verified, since

" #" # " # " #
−1 2 −2 4 −2
Av2 = = = −2 · = λ2 v̂2 .
2 2 1 −2 1

In practice, of course, only the dominant eigenvalue and eigenvector of A

would be known following the initial application of the power method or
Rayleigh method to A. Also, the dominant eigenvalue and eigenvector of
the deflated matrix B would need to be found using the power method or
Rayleigh method, rather than analytically. Incidentally, if we had done this
we would have seen the ‘magnification ratio’ changing sign at each iteration.
This simply indicates that the dominant eigenvalue, to which the method is
converging, is negative (in this case −2).

140