Chapter 3

Eigenvalues and Eigenvectors

3.1 Power and Inverse Power Method

Power method is an iterative method for finding the dominant eigenvalue, i.e. one with the
largest magnitude, of an n × n matrix and the corresponding eigenvector.

The method assumes the following:

(i) There is a single eigenvalue of maximum magnitude
(ii) There is a linearly independent set of n eigenvectors.

The method proceeds as follows:

According to the first assumption, let the n × n matrix A have eigenvalues λ1 , λ2 , · · · , λn satisfying

|λ1 | > |λ2 | ≥ · · · ≥ |λn | (3.1)

The second assumption implies that there is a basis u(1) , u(2) , · · · , u(n) for Cn such that
 

Au(j) = λj u(j) , 1 ≤ j ≤ n. (3.2)

Let x(0) be any element of Cn such that when x(0) is expressed as a linear combination of the basis
elements u(1) , u(2) , · · · , u(n) , and we assume that the coefficient of u(1) is not zero. Thus

x(0) = a1 u(1) + a2 u(2) + · · · + an u(n) (a1 6= 0) (3.3)

and we form
x(1) = Ax(0) , x(2) = Ax(1) , · · · , x(k) = Ax(k−1) (3.4)
so that
x(k) = Ak x(0) . (3.5)
Without loss of generality, we will assume that all the coefficients aj of the vectors u(j) are absorbed,
i.e. we may write (3.5) as
x(0) = u(1) + u(2) + · · · + u(n) . (3.6)

47
48 NUM702S: Numerical Methods 2

Using Eqs. (3.5) and (3.6) we have

x(k) = A(k) u(1) + A(k) u(2) + · · · + A(k) u(n) (3.7)

Using Eq. (3.2) we obtain

x(k) = λk1 u(1) + λk2 u(2) + · · · + λkn u(n) (3.8)
and this can be rewritten as
"  k  k #
λ2 λn
x(k) = λk1 u(1) + u(2) + · · · + u(n) (3.9)
λ1 λ1

Since |λ1 | > |λj | for 2 ≤ j ≤ n, we see that the coefficients

 k
λj
→ 0 as k → ∞
λ1

and the vector within the brackets converges to u(1) . All the latter information can be written in
a compact form as
x(k) = λk1 u(1) + (k) with lim (k) = 0.
 
(3.10)
k→∞

In order to handle ratios, let assume that φ is an arbitrary linear mapping on Cn such that
φ(u(1) ) 6= 0. From Eq. (3.10) and by linearity of φ, we have

φ(x(k) ) = λk φ(u(1) ) + φ((k) )

 
(3.11)

As a result, the following ratio converges to λ1 as k → ∞, i.e.

φ(x(k+1) )
lim rk = lim = λ1 (3.12)
k→∞ k→∞ φ(x(k) )

Indeed
φ(x(k+1) ) φ(u(1) ) + φ((k+1) )
 
rk = = λ1 , (3.13)
φ(x(k) ) φ(u(1) ) + φ((k) )
and since φ is linear, φ((k) ) and φ((k+1) ) → 0 as k → ∞ and the term in brackets in Eq. (3.13)
converges to 1, and thus we obtain Eq. (3.12).

Definition 3.1.1 (Power Method)

The procedure consisting of computing λ1 through Eqs. (3.1)-(3.10) is called the power method
(for computing λ1 ).

The method can be written in algorithmic form as follows.

NUM702S: Numerical Methods 2 49

Algorithm 2 (Power Method)

Input: n, A, x, kmax
For j = 1, 2, . . . , kmax Do
y ← Ax
r ← φ(y)/φ(x)
x←y
Output: k, x, r
End For

where A is the target matrix, n is the matrix dimension, x is the eigenvalue and kmax is the
maximum number of iterations.

It is worth mentioning here that in practical implementation of the power method, the eigen-
vector x(k) may converge to zero or become unbounded. To avoid this, it is advisable to normalise
x(k) . Thus we obtain a modified Power method.

Algorithm 3 (Modified Power Method)

Input: n, A, x, kmax
For j = 1, 2, . . . , kmax Do
y ← Ax
r ← φ(y)/φ(x)
x ← y/kyk
Output: k, x, r
End For

Example 3.1.2 Find the largest eigenvalue of the following matrix

 
1 1 −1
A = 1 2 1
2 −1 1

using the initial vector x(0) = (1, 1, 1)T .

Solution

Using Eq. (3.4)

x(1) = Ax(0) , x(2) = Ax(1) , · · · , x(k) = Ax(k−1) (3.14)
we have  
1
x(1) = Ax(0) = 4
2
Since the second component of x(1) has the largest magnitude, we take k = 2 so that the first
50 NUM702S: Numerical Methods 2

approximation to λ1 is
(1)
(1) x2 4
λ1 = (0)
= = 4.
x2 1
Further iterations are given by
    
1 1 −1 1 3
x(2) = Ax(1) = 1 2 1  4 = 11 (k = 2)
2 −1 1 2 0
(2)
(2) x2 11
λ1 = (1)
= = 2.75.
x2 4

    
1 1 −1 3 14
x(3) = Ax(2) = 1 2
 1  11 = 25  (k = 2)
 
2 −1 1 0 −5
(3)
(3) x2 25
λ1 = (2)
= = 2.27.
x2 11

    
1 1 −1 14 44
x(4) = Ax(3) = 1 2 1   25  =  59  (k = 2)
2 −1 1 −5 −2
(4)
(4) x2 59
λ1 = (3)
= = 2.36.
x2 25

..
.  
1731
x(8) = Ax(7) = 3152 (k = 2)
108
(8)
(8) x2 3152
λ1 = (7)
= = 2.5857.
x2 1219

 
4775
x(9) = Ax(8) = 8143 (k = 2)
418
(9)
(9) x2 8143
λ1 = (8)
= = 2.5834.
x2 3152
..
.
NUM702S: Numerical Methods 2 51

We therefore conclude that the dominant eigenvalue is 2.6 after 8 iterations.

MATLAB DEMO

%pw_method_inf
function [sigma, x] = pw_method_inf(A,x_0,tol,max_N)
if (nargin <= 2)
disp(’tolerance undefined’)
return;
end

if (nargin == 3)
max_N = 20;
end

x_0 = x_0/norm(x_0, inf);

sigma_0 = 1;

x = A*x_0;

[x_max, i] = max(abs(x));
sigma = x(i);

if (abs(sigma) < 1e-16)

return;
end

x = x/sigma;

j = 1;
while ( (norm(x_0 - x, inf) > tol) | ...
(norm(A*x - sigma*x, inf) > tol) | ...
(abs(sigma - sigma_0) > tol) & ...
( j <= max_N ) )

x_0 = x;
sigma_0 = sigma;

x = A*x_0;
[x_max, i] = max(abs(x));
sigma = x(i);

if ( abs(sigma) < 1e-16 )

return;
52 NUM702S: Numerical Methods 2

end

x = x/sigma;
j = j+1;
end

if ( norm(x_0 - x, inf) > tol | ...

norm(A*x - sigma*x, inf) > tol | ...
abs(sigma - sigma_0) > tol & ...
( j > max_N ) )

disp(’divergence likely’)
end

Theorem 3.1.3 If A is a non-singular matrix and λ is an eigenvalue of A, then λ−1 is an eigen-

value of A−1 .
Proof 3.1.4 Let Ax = λx, x 6= 0. Then x = A−1 (λx) = λA−1 x.
Hence A−1 x = λ−1 x. Thus λ−1 is an eigenvalue of A−1 .


Remarks
The above theorem suggests a means of computing the smallest eigenvalue of A.
Let the eigenvalues of A satisfy the following

|λ1 | > |λ2 | ≥ · · · ≥ |λn−1 | > |λn | > 0 (3.15)

Eq. (3.15) implies that A is non-singular since 0 is not an eigenvalue. Thus the eigenvalues of A−1
are λ−1
j satisfying
|λ−1 −1 −1
n | > |λn−1 | ≥ · · · ≥ |λ1 | > 0. (3.16)
Now we can compute λ−1
n by applying the Power method to the inverse A
−1
from the procedure
of solving
Ax(k+1) = x(k) (3.17)
in contrast to Ax(k−1) = x(k) given by Eq. (3.14).
The above equation (3.17) is employed to avoid finding first the inverse A−1 and then using it in
x(k+1) = A−1 x(k) suggested by Eq. (3.14), which is not a very good computational approach.
Eq. (3.17) can be solved by the LU decomposition method.
Definition 3.1.5 The above method given by Eq. (3.17) and Theorem 3.1.3 is known as the in-
verse power method. it finds the smallest eigenvalue of A through the computation of the largest
eigenvalue of the inverse A−1 of A and then take the reciprocal.
Exercise
Use the Matlab code for the Power method to obtain the smallest eigenvalue of the matrix in
Example 3.1.2 with the same initial vector.
NUM702S: Numerical Methods 2 53

3.1.1 Recall of the LU-decomposition of a square matrix

Since this method tackled in Numerical Method 1 will be used in the inverse power method, it is
important to recall its definition and main features.
Definition 3.1.6 The LU-decomposition of a square matrix A is an expression of the form

A = LU

where L and U are lower and upper triangular matrices respectively.

• When all teh diagonal elements of L is 1, that is lkk = 1, 1 ≤ k ≤ n, the method is called
Doolittle decomposition/factorosation.
• When all the diagonal elements of U are 1, that is ukk = 1, 1 ≤ k ≤ n, the method is called
Crout’s decomposition/factorisation.
• When U = LT so that lkk = ukk , 1 ≤ k ≤ n, the method is called Cholesky’s factorisation.
In this method, A is required to be real, symmetric and positive definite.

Algorithm 4 (Doolittle LU-Decomposition Method)

Input: n, A = (ai,j )1≤i,j≤n
For k = 1 to n Do
lk,k ← 1
For j = k + 1 to n Do !
k−1
X
uk,j ← ak,j − lk,s us,j
s=1
End For
For i = k + 1 to n Do !
k−1
X
li,k ← ai,k − li,s us,k /uk,k
s=1
End For
End For
Output: (li,j ), (ui,j )

The following theorem gives us a sufficient condition for the existence of the LU -Decomposition of
a square matrix.

Theorem 3.1.7 (LU-Decomposition)

If all n leading principal minors of an n × n matrix A are non-singular, then A has an LU-
Decomposion.
Example 3.1.8 Find the largest and the smallest eigenvalues of the matrix
 
6 5 −5
A = 2 6 −2
2 5 −1
54 NUM702S: Numerical Methods 2

using the Power and Inverse Power methods respectively with the initiale vector x(0) = (−1, 1, 1)T .

In summary, we have the following results for the two variants of the power method

Method Equation Description

(k+1)
Power Method x = Ax(k) Solving for largest
eigenvalue
Inverse Power Method Ax(k+1) = x(k) Solving for smallest
eigenvalue

3.2 Shifted Inverse Power Method

Theorem 3.2.1 (Shifting Eigenvalues)
Suppose that λ, V is an eigenpair of A. If α is any constant, then λ − α is an eigenpair of the
matrix A − αI.

Theorem 3.2.2 ( Shifted-Inverse Power Method )

Assume that n × n matrix A has distinct eigenvalues λ1 , λ2 , . . . , λn and consider the eigenvalue
λj . Then a constant α can be chosen so that µ1 = 1/(λj − α) is the dominant eigenvalues of
(A − αI)−1 . Furtehrmore, if x(0) is chosen appropriately, then the sequences {x(k) } and {ck } are
generated recursively by
1
y(k) = (A − αI)−1 x(k) and x(k+1) = y(k) ,
ck+1
where
(k) (k) (k)
ck+1 = yj and yj = max {|yi |}
1≤n≤n

will converge to the dominant eigenpair µ1 , Vj of matrix (A − αI)−1 . Finally the corresponding
eigenvalue for the matrix A is given by the calculation
1
λj = + α.
µ1

Example
Use the shifted-inverse power method to find the eigenpairs of the matrix
 
0 11 −5
A = −2 17 −7 
−4 26 −10

Use the fact that the eigenvalues of A are λ1 = 4, λ2 = 2, and λ3 = 1 and select an appropriate
α and starting vector for each case.
NUM702S: Numerical Methods 2 55

Solution
Case-1: For the eigenvalue λ1 = 4, we select α = 4.2 and the starting vector x(0) = (1, 1, 1)T . We
then form the matrix A − 4.2I and compute the y(0) as the solution of the linear sytem
   
−4.2 11 −5 1
 −2 12.8 −7  y(0) = x(0) = 1
−4 26 −14.2 1

and we obtain y(0) = (−9.545454545, −14.09090909, −23.18181818)T .

(0) (0)
We then compute c1 = yj = max{|yi |} = −23.18181818 and x(1) =
T
(0.4117647059, 0.6078431373) . We proceed with the iteration and the results are given in
the table below

(A − αI)−1 x(k) = ck+1 x(k+1)

(A − αI)−1 x(0) = −23.18181818(0.4117647059, 0.6078431373, 1)T = c1 x(1)
(A − αI)−1 x(1) = −5.356506239(0.4009983361, 0.6006655574, 1)T = c2 x(2)
(A − αI)−1 x(2) = −5.030252609(0.4000902120, 0.6000601413, 1)T = c3 x(3)
(A − αI)−1 x(3) = −5.002733697(0.4000081966, 0.6000054644, 1)T = c4 x(4)
(A − αI)−1 x(4) = −500248382(0.4000007451, 0.6000004967, 1)T = c5 x(5)
(A − αI)−1 x(5) = −5.000022579(0.4000000677, 0.6000000452, 1)T = c6 x(6)
(A − αI)−1 x(6) = −5.000002053(0.4000000062, 0.6000000041, 1)T = c7 x(7)
(A − αI)−1 x(7) = −5.000000187(0.4000000006, 0.6000000004, 1)T = c8 x(8)
(A − αI)−1 x(8) = −5.000000017(0.4000000001, 0.6000000000, 1)T = c9 x(9)
We can observe that the sequence {ck } converges to µ1 = −5, which the dominant eigenvalue of
(A − 4.2I)−1 , and {x(k) } converges to V1 = (2/5, 3/5, 1)T . The eigenvalue λ1 of A is given by the
computation λ1 = 1/µ1 + α = 1/(−5) + 4.2 = −0.2 + 4.2 = 4.

Case-2: For the eigenvalue λ2 = 2, we select α = 2.1 and the starting vector x(0) = (1, 1, 1)T .
Proceeding in a manner similar to the above case, we obtain µ1 = −10 and V2 = (1/4, 1/2, 1)T .

Case-3: Similarly for the eigenvalue λ3 = 1, we select α = 0.875 and the starting vector x(0) =
(0, 1, 1)T . We obtain µ1 = 8 and V3 = (1/2, 1/2, 1)T .

3.3 QR Factorisation
Definition 3.3.1 Given a matrix A that can be written as

A = QR (3.18)

where Q is an orthogonal matrix, i.e Q−1 = QT , and R is an upper triangular matrix. Then the
above decomposition is called QR factorisation.

The QR algorithm is a factorisation technique employed to determine all the eigenvalue of a

56 NUM702S: Numerical Methods 2

symmetric matrix.
QR algorithm is applicable to a symmetric and tridiagonal matrix. If this is not the case, the first
step consists of using another method called Householder’s method to construct a (symmetric)
tridiagonal matrix having the same eigenvalues as the given/target matrix.

Let A be a symmetric tridiagonal matrix written as

 
a1 b 2 0 · · · 0
. . . .. 
 b 2 a2 b 3 .


A=
 . . .

(3.19)
 0 b 3 a 3 0 

. .
 .. . . . . . . . . bn 


0 . . . 0 b n an

Remark: If b2 = 0 or bn = 0, then the 1 × 1 matrix [a1 ] or [an ] immediately gives an eigenvalue

a1 0r [an ] of A.
When bj = 0 for some j, where 2 < j < n, the problem can be reduced to one considering some
other smaller matrices, instead of A.
Now if none of the bj are zero, the QR method proceeds by forming a sequence of matrices
A = A(1) , A(2) , A(3) , . . . as follows:
1. A(1) = A is factored as a product A(1) = Q(1) R(1) , where Q(1) is orthogonal and R(1) is upper
triangular.
2. A(2) is defined as A(2) = R(1) Q(1) .
Thus, A(i+1) is defined by the product of R(i) and Q(i) in the reverse direction

A(i+1) = R(i) Q(i) .

Since Q(i) is orthogonal, then R(i) = Q(i)T A(i) and we have

A(i+1) = R(i) Q(i) = Q(i)T A(i) Q(i) = Q(i)T A(i) Q(i)

 
(3.20)

So that A(i+1) is symetric (and moreover similar) with the same eigenvalues as A(i) .
Besides, following the definition of R(i) and Q(i) , we also ensure that A(i+1) is tridiagonal.
By induction, A(i+1) has the same eigenvalues as A and A(i+1) tends to a diagonal matrix with the
eigenvalues of Acalong the diagonal.

Theorem 3.3.2 If A is a symmetric matrix and D is a diagonal matrix whose diagonal entries
are the eigenvalues of A, then there exists an orthogonal matrix Q such that

D = Q−1 AQ = QT AQ. (3.21)

Example 3.3.3 Matlab has a built-in function for QR decomposition as seen in the following
example.
NUM702S: Numerical Methods 2 57

>> A = [3 1 0;1 3 1;0 1 3] % Create a matrix

A =
3 1 0
1 3 1
0 1 3
>>[Q R] = qr(A) % command to create matrices Q and R respectively
Q =
-0.9487 0.2941 0.1162
-0.3162 -0.8823 -0.3487
0 -0.3676 0.9300

R =
-3.1623 -1.8974 -0.3162
0 -2.7203 -1.9851
0 0 2.4412

>> QR % check that A = QR

ans =
3.0000 1.0000 -0.0000
1.0000 3.0000 1.0000
0 1.0000 3.0000

It is observed that Q is orthogonal

Example 3.3.4 Consider the following modification to the above example using the following com-
mands

>> B = [3 1 0;1 3 1]
B =
3 1 0
1 3 1

>> [q1 R1] = qr(B)

Q1 =
-0.9487 -0.3162
-0.3162 0.9487

R1 =
-3.1623 -1.8974 -0.3162
0 2.5298 0.9487

>> Q1*R1 % Check

ans =
3.0000 1.0000 0.0000
1.0000 3.0000 1.0000
58 NUM702S: Numerical Methods 2

The above example shows that for an m × n matrix B, the Matlab command

[Q R] = qr(B)

produces an m × n upper triangular matrix R and an m × n Unitary Matrix Q (i.e QT Q = I) so

that B = QR.
In what follows, we gives two method to (munually) compute the QR factorisation of a matrix A.

QR factorisation via Gram-Schmidt Process

We illustrate the method by applying it to the matrix

 
1 1 4
A = 1 4 2
1 1 0

which can be written as a row vector of column vectors

A = [x1 , x2 , x3 ]

The construction of Q is done via the Gram-Schmidt Process by determining an orthonormal

basis of V = span{x1 , x2 , x3 }.

• Orthogonal basis.

Let v1 = [1, 1, 1, 1]T

   
1 1
4 x2 · v1 1
v2 = x2 − projv1 x2 = 
4 − v1 · v1
  
1
1 1
     
1 1 −3/2
4 5 1  3/2 
= 
4 − 2 1 =  3/2 
    

1 1 −3/2

     
4 1 −3/2
2 x3 · v1 1 x3 · v2  3/2 
v3 = x3 − projv1 x3 − projv2 x3 = 
2 − v1 · v1
  −
1 v2 · v2
 
 3/2 
0 1 −3/2
     
4 1 2
2 1 0
2 − 2 1 − 0 =  0 
=      

0 1 −2
NUM702S: Numerical Methods 2 59

Thus       

 1 −3/2 2 

1 3/2 0
 
0 0 0
    
v1 =   , v2 = 
    , v =   (orthogonal basis).

 1 3/2  3  0  
1 −3/2 −2
 

• Orthogonal basis

√
kv10 k = 2, kv20 k = 3, kv30 k = 2 2
Defining vi = vi0 /kvi0 k for i = 1, 2, 3, the orthogonal matrix Q = [v1 , v2 , v3 ] is given by
 √ 
1/2 −1/2 1/ 2
1/2 1/2 0 
Q= 1/2 1/2

0√ 
1/2 −1/2 −1/ 2

• Computation of R.

From A = QR, we have R = QT A since QT Q = I. We find

 
2 5 4
R = QT A = 0 3 √ 0 
0 0 2 2

And therefore  √ 
1/2 −1/2 1/ 2  
1/2 1/2 2 5 4
0 

A=
1/2 1/2 0 3 √0 
0√ 
0 0 2 2
1/2 −1/2 −1/ 2

Following the same procedure as above, we have

 
1 2 4
B := 0 0
 5
0 3 6
  
1 0 0 1 2 4
=  0 0 1 0 3 6
0 1 0 0 0 5

QR factorisation by Givens rotations (plane rotation)

In an Rn , Givens rotation is defined as

60 NUM702S: Numerical Methods 2

where c = cos(θ), s = sin(θ).

Gi,j (θ)A only affects rows i and j of matrix A.
In R2 , it takes the form  
c −s
G(θ) =
s c
Let A be a general matrix  
a a
A = 11 12
a21 a22
A = QR ⇔ QT A = R
And  
r r
Q A = 11 12
T
0 r22
This transformation can be achieved using a Givens rotation that rotates
   
a11 r
7→ 11
a21 0
    
T c −s a11 a12 r11 r12
Q A= =
s c a21 a22 0 r22
Using the trigonometric circle and Pythagoras theorem, we have
a11
cos(−θ) = c = p ,
a211 + a221
a21
sin(−θ) = −s = p 2 .
a11 + a221

Thus    
T c −s 1 a11 a21
Q = =p 2
s c a11 + a221 −a21 a11
NUM702S: Numerical Methods 2 61

Example 3.3.5 Find the QR factorisation of the following matrix using Givens rotation method
 
4 −3 5
A = 3 4 5
0 12 0

The procedure is to eliminate the entries of A so as to obtain an upper triangular matrix R by

applying a sequence of Givens rotations.

• Multiply A by G12 to make a21 zero.

   
4 3 0 4 0 7
1
G12 = −3 4 0 then G12 A = 0 5 1
5
0 0 5 0 12 0

• Next, multiply the previous result by G23 to eliminate the 12 in the second column. There is no
need to multiply by G13 as there is already a zero in the (3,1) position. So
   
13 0 0 5 0 7
1 
G23 = 0 5 12 then G23 G12 A = 0 13 5/13  = R
13
0 −12 5 0 0 −12/13

Thus
G23 G12 A = R ⇒ A = (G23 G12 )−1 R = (GT12 GT23 )R = QR

where    
52 −15 36 5 0 7
1 
Q= 39 20 −48 , R = 0 13 5/13  .
65
0 60 25 0 0 −12/13

3.4 Householder’s Method

The Householder’s method is used to find a symmetric tridiagonal matrix B that is similar to
a given matrix A. After implementing the Householder’s method, efficient methods such as QR
algorithm can be used for accurate approximation of the eigenvalues of the resulting symmetric
tridiagonal matrix.

Definition 3.4.1 Let w ∈ Rn be an n × 1 vector such that

n
X
T
w w= wk2 = 1 (3.22)
k=1

Then the n × n matrix

H = I − 2wwT (3.23)
is called the Householder Transformation or Matrix.
62 NUM702S: Numerical Methods 2

Theorem 3.4.2 A Householder transformation (3.23) is symmetric and orthogonal, so that

H −1 = H

Proof 3.4.3 Since

(wwT )T = (wT )T wT = wwT
thus the transformation is symmetric.
Moreover, since wT w = 1, then we have

HH T = (I − 2wwT )(I − 2wwT )

= I − 2wwT − 2wwT + 4w (wT w) w
| {z }
=1
T T
= I − 4ww + 4ww = I (3.24)

Thus we have
HH T = I or H −1 = H T (3.25)
Hence, H is orthogonal. Therefore
H −1 = H T = H (3.26)


Householder’s method usually begins by determining a transformation H (1) with the property that
A(2) = H (1) AH (1) has
(2)
aj,1 = 0, for j = 3, 4, . . . , n (3.27)
(2)
and by symmetricity, a1,j = 0.
The vector w = (w1 , w2 , · · · , wn )T is chosen such that Eq. (3.22) is satisfied. Then the process
continues until we we get
H (i) = I − 2w(i) (w(i) )T (3.28)
and
A(i+1) = H (i) A(i) H (i) (3.29)
where A(i+1) is a symmetric tridiagonal matrix.
Coming back to the practical implementation of w for the construction of H (1) . As mentioned
above, we choose of w such that wT w = 1 and in the matrix

A(2) = H (1) AH (1) = (I − 2wwT )A(I − 2wwT ),

(2) (2)
we have a11 = a11 and aj1 = 0, for each j = 3, 4, . . . , n. This choice imposes n conditions on the
n unknowns w1 , w2 , . . . , wn .
NUM702S: Numerical Methods 2 63

Setting w1 = 0 ensures that a(2) = a11 . We want H (1) to satisfy

   
a11 a11
 a21   α 
H (1)  ...  =  0 
   
(3.30)
   
  . 
an−1,1   .. 

an1 0

where α will be determined below. To simplify the notation, let

ŵ = (w2 , . . . , wn )T ∈ Rn−1 , ŷ = (a21 , a31 . . . , an1 )T ∈ Rn−1

and Ĥ be the (n − 1) × (n − 1) Householder’s transformation

Ĥ = In−1 − 2ŵŵT .

Eq. (3.30) reads

 
       a 11
a11 1 0 ··· 0 a11 a11 −−
 
 a21   0   −−   α 
−−  

H (1)  ..  =  .. = =
      
ŷ   Ĥ ŷ   0 
 
 .   . Ĥ  .. 

an1 0  . 
0
with
Ĥ ŷ = (In−1 − 2ŵŵT )ŷ = ŷ − 2(ŵT ŷ)ŷ = (α, 0, · · · , 0)T (3.31)
Let r = ŵT ŷ. Then

(α, 0, · · · , 0)T = (a21 − 2rw2 , a31 − 2rw3 , · · · , an1 − 2rwn )T ,

and we can determine all the wi once we know α and r.
Equating components gives
α = a21 − 2rw2
and
0 = aj1 − 2rwj , for each j = 3, . . . , n.
Thus
2rw2 = a21 − α (3.32)
and
2rwj = aj1 , for each j = 3, . . . , n. (3.33)
64 NUM702S: Numerical Methods 2

Squaring both sides of each of the equations and adding corresponding terms gives
n
X n
X
2
4r wj2 = (a21 − α) + 2
a2j1 .
j=2 j=3

Since ŵT ŵ = 1 and w1 = 0, we have

n
X
2
4r = a2j1 − 2αa21 + α2 . (3.34)
j=2

Eq. (3.31) and the fact that P is orthogonal imply that

α2 = (α, · · · , 0)(α, · · · , 0)T = (Ĥ ŷ)T Ĥ ŷ = ŷT ĤT Ĥŷ = ŷT ŷ.

Thus n
X
2
α = a2j1 ,
j=2

which, when substituted into Eq. (3.34), gives

n
X
2r2 = a2j1 − αa21 .
j=2

To ensure that 2r2 = 0 only if a21 = a31 = · · · = an1 , we choose

n
!1/2
X
α = −sign(a21 ) a2j1 ,
j=2

which implies that

n n
!1/2
X X
2r2 = a2j1 + |a21 | a2j1 ,
j=2 j=2

with this change of α and 2r2 , we solve Eq. (3.32) and Eq. (3.33) to obtain
a21 − α
w2 =
2r
and
aj1
wj = , for each j = 3, . . . , n.
2r
NUM702S: Numerical Methods 2 65

To summarise the choice of H (1) , we have

n
!1/2 
X 
α = −sign(a21 ) a2j1 ,





j=2 


 1/2 
1 2 1



r= α − a21 α , 

2 2
w1 = 0, 



a21 − α


w2 = , and



2r 

aj1


wj = , for each j = 3, . . . , n.


2r
With this choice,  (2) (2) 
a11 a12 0 ··· 0
 (2) (2) (2) (2) 
a21 a22 a23 ··· a2n 
 (2) (2) (2) 
A(2) = H (1) AH (1)  0 a32 a33
= ··· a3n 
 .. .. .. .. 

 . . . . 
(2) (2)
0 an2 an3 ··· a(2)
nn

Having found H (1) and computed A(2) , the process is repeated for k = 2, 3, . . . , n − 2 as follows
!1/2 
X n 
(k) 2 
α = −sign(ak+1,k ) ajk , 



j=k+1




 1/2 
1 2 1 (k) 

r= α − ak+1,k α ,


2 2





(k) (k) (k) 
w1 = w2 = · · · = wk = 0, 

(k)
(k) ak+1,k − α 
wk+1 = , and


2r




(k) 
(k) ajk 

wj = , for each j = k + 2, k + 3, . . . , n,



2r 


H (k) = I − 2w(k) (w(k) )T , and






A(k+1) = H (k) A(k) H (k) ,


Continuing the process, the tridiagonal and symmetric matrix A(n−1) is formed, where

A(n−1) = H (n−2) H (n−3) · · · H (1) AH (1) · · · H (n−3) H (n−2) .

66 NUM702S: Numerical Methods 2

Example 3.4.4 Apply Householder’s transformations to the symmetric 4 × 4 matrix

 
4 1 −2 2
1 2 0 1
A= −2 0 3 −2


2 1 −2 −1

to produce a symmetric tridiagonal matrix that is similar to A.

Solution
For the determination of H (1) we first need to compute α, r and then w.

4
!1/2 1/2
√

X 1 1
α = −(1) a2j1 = −3, r = 2
(−3) − (1)(−3) = 6,
j=2
2 2

√ √ √ !T
6 6 6
w= 0, ,− , ,
3 6 6

   
1 0 0 0√ !2 1
0 1 0 0 6   2  · (0, 2, −1, 1)

H (1) = 
0
−2
0 1 06 −1
0 0 0 1 1
 
1 0 0 0
0 −1/3 2/3 −2/3
= 
0
,
2/3 2/3 1/3 
0 −2/3 1/3 2/3

then  
4 −3 0 0
−3 10/3 1 4/3 
A(2) = H (1) AH (1) =
0
.
1 5/3 −4/3
0 4/3 −4/3 −1

Continuing to the second iteration,

√ √ !T
5 2 5 √ 5
α=− , r= , w= 0, 0, 2 5, ,
3 3 5
 
1 0 0 0
0 1 0 0 
H (2) =
0
,
0 −3/5 −4/5
0 0 −4/5 3/5
NUM702S: Numerical Methods 2 67

and the symmetric tridiagonal matrix is

 
4 −3 0 0
−3 10/3 −5/3 0 
A(2) = 
 0 −5/3 −33/25 68/75  .


0 0 68/75 149/75

3.5 Jacobi Method for Symmetric Matrices

Introduction
• Jacobi method is an iterative method to determine the eigenvalues and eigen vectors of a
symmetric matrix.
• Using the method, a solution is guaranteed for all real symmetric matrices.
• It is based on series of rotations called Jacobi rotations which are similarity transformations.
• The rotations are chosen to discard the off-diagonal elements in such a way that eigenvalues are
preserved.

Consider the rotation matrix J(p, q, θ) similar to Givens matrix G(p, q, θ) with θ replaced with
−θ.

• The matrix J(p, q, θ) is applied to a symmetric A. It is a similarity transformation which rotates

rows and columns p and q of A through an angle θ so that its (p, q) and (q, p) entries become zero.
Let A0 = J(p, q, θ)T AJ(p, q, θ) be the similarity transformation of A.
Let off(A) and off(A0 ) be the square root of the sum of squares of all off-diagonal elements of A
and A0 respectively. Then
X n
off(A)2 = kAk2F − a2kk .
k=1

Since the Frobenius norm is invariant under orthogonal transformations, and only p and q columns
are reformed in matrix A0 , we have
n
X
0 2
off(A ) = kA0 k2F − a02
kk
k=1
n
X
= kAk2F − a02 02 02
kk − (app + aqq )
k6=p,q
X n
= kAk2F − a2kk − (a2pp + 2a2pq + a2qq )
k6=p,q
X n
= kAk2F − a2kk − 2a2pq
k=1
= off(A) − 2a2pq < off(A)2 .
2
68 NUM702S: Numerical Methods 2

Which shows that the size of the off-diagonal part decreases by applying the above similarity
transformation, that is, the off-diagonal elements of A0 are reducing to zero.

• The post multiplication of matrix A by J(p, q, θ) yields changes in columns p and q.

• In the same way, the pre multiplication of A by J(p, q, θ)T brings changes in rows p and q.
• Hence, the transformation A0 = J(p, q, θ)T AJ(p, q, θ) alters only rows p and q and columns p
and q of A.

The elements a0jk of A0 are given by the formulae

a0jp = cajp − sajq


when j 6= p and j 6= q, 

a0jq = sajp + cajq

when j 6= p and j 6= q, 



a0pp 2 2
= c app + s aqq − 2csapq (3.35)

a0qq = s2 app + c2 aqq + 2csapq





a0pq = (c2 − s2 )apq + cs(app − aqq )


and the rest of the elements are found by symmetry.

Zeroing out a0pq and a0qp

The goal at every step of Jacobi’s iteration is to make the off-diagonal elements a0pq and a0qp zero.
From c = cos(θ) and s = sin(θ), we have

c2 − s 2
φ = cot(2θ) = . (3.36)
2cs
If apq 6= 0 and we have to generate a0pq = 0, then using the last equation (i.e Eq. (3.36)) in
Eq. (3.35), we have
0 = (c2 − s2 )apq + cs(app − aqq ). (3.37)
Eq. (3.37) yields
c2 − s 2 app − aqq
= . (3.38)
cs apq
which is used in Eq. (3.36) to solve for φ:
aqq − app
φ= . (3.39)
2apq

Less round-off error is generated if we use tan(θ) in the later computations. So we define
s
t = tan(θ) = . (3.40)
c
NUM702S: Numerical Methods 2 69

From Eq. (3.36), dividing the numerator and denominator by c2 , we obtain

s2
1− 2 2
φ= c = 1−t
s 2t
2
c
which gives the equation
t2 + 2tφ − 1 = 0. (3.41)
The smaller root of Eq. (3.41) corresponds to the smaller angle of rotation |θ| ≤ π/4. Solving this
equation, we find that this root is

sign(φ)
t = −φ ± (φ2 + 1)1/2 = (3.42)
|φ| + (φ2 + 1)1/2

where (
1, if φ ≥ 0
sign(φ) =
−1, otherwise
Hence, c and s are found to be
1

c= 2 
(t + 1)1/2 (3.43)
s = ct


As application, let’s consider the following example.

Example 3.5.1 Use the Jacobi method to find the eigen-pairs of the matrix
 √ 
√1 2 √2
A =  2 √3 2
2 2 1

This is a real symmetric matrix.

To start, we select the off-diagonal element with the maximum modulus: a31 = a13 = 2.
We will perform a similarity transformation to make it zero. It means p = 1, q = 3, thus

app = 1, aqq = 1
apq = aqp = 2
 
1 −1 2apq π
θ = tan = (from Eq. (3.39))
2 aqq − app 4
1
c=s= √
2
70 NUM702S: Numerical Methods 2

We can write down J1 = J(1, 3, π/4)

 1 1 
√ 0 −√
 2 2
J1 =  0 1 0 ,
 
 1 1 
√ 0 √
2 2
 
3 2 0
A1 = J1T AJ1 = 2 3 0 
0 0 −1
We repeat the process with the new matrix A1 . This time

p = 1, q = 2
apq = aqp = 2
app = aqq = 3
θ = π/4.
 1 1 
√ −√ 0
 2 2 
 1 1
J2 =  √

√ 0 
2 2
 
0 0 1
 
5 0 0
A2 = J2T A1 J2 = 0 1 0 
0 0 −1
Therefore: λ1 = 5, λ2 = 1, λ3 = −1.
The eigenvectors are the columns of the matrix
 √ 
1/2
√ −1/2
√ −1/ 2
J = J1 J2 = 1/ 2 1/ 2
 0√ 
1/2 −1/2 1/ 2

Hence
√
v1 = [1/2 1/ 2 1/2]T is the eigenvector associated with λ1 = 5,
√
v2 = [−1/2 1/ 2 − 1/2]T is the eigenvector associated with λ2 = 1 and
√ √
v3 = [−1/ 2 0 1/ 2]T is the eigenvector associated with λ3 = −1.