v = x ± x 2 em , (5.39)

5.6 The QR Method for Matrices in Hessenberg Form 205

em being the m-th unit vector of Rn . The matrix P computed by (5.38)

depends on the vector x itself, and it can be checked that

® T

Px = °0, 0, . . . , ± x 2 , 0, . . . , 0» . (5.40)

m

Example 5.5 Let x = [1, 1, 1, 1]T and m = 3; then

® ® ®

5 ’1 ’3 ’1

1 0

1 ’3 ’1 0

1 ’1 5

v= , P = , Px = .

3 6 ’3 ’3 ’3 ’3 ’2

° » ° » ° »

’1 ’1 ’3 5

1 0

•

If, for some k ≥ 1, the ¬rst k components of x must remain unaltered,

while the components from k + 2 on are to be set to zero, the Householder

matrix P = P(k) takes the following form

®

Ik 0 (k) (k) T

» , Rn’k = In’k ’ 2 w (w ) .

°

P(k) = (5.41)

w(k) 22

0 Rn’k

As usual, Ik is the identity matrix of order k, while Rn’k is the elementary

Householder matrix of order n ’ k associated with the re¬‚ection across the

hyperplane orthogonal to the vector w(k) ∈ Rn’k . According to (5.39), the

Householder vector is given by

(n’k)

w(k) = x(n’k) ± x(n’k) 2 e1 , (5.42)

where x(n’k) ∈ Rn’k is the vector formed by the last n ’ k components

(n’k)

is the ¬rst unit vector of the canonical basis of Rn’k . We

of x and e1

notice that P(k) is a function of x through w(k) . The criterion for ¬xing

the sign in the de¬nition of w(k) will be discussed in Section 5.6.5.

The components of the transformed vector y = P(k) x read

±

j = 1, · · · , k,

yj = xj

j = k + 2, · · · , n,

yj = 0

yk+1 = ± x(n’k) 2 .

The Householder matrices will be employed in Section 5.6.2 to carry out the

reduction of a given matrix A to a matrix H(0) in upper Hessenberg form.

This is the ¬rst step for an e¬cient implementation of the QR iteration

(5.32) with T(0) = H(0) (see Section 5.6).

206 5. Approximation of Eigenvalues and Eigenvectors

Example 5.6 Let x=[1, 2, 3, 4, 5]T and k = 1 (this means that we want to set to

zero the components xj , with j = 3, 4, 5). The matrix P(1) and the transformed

vector y=P(1) x are given by

® ®

1 0 0 0 0 1

0 0.6804 7.3485

0.2722 0.4082 0.5443

= 0 ’0.3816 , y= 0 .

’0.3053

0.4082 0.7710

P(1)

0 ’0.5089 0

’0.3053

0.5443 0.5929

° » ° »

’0.3816 ’0.5089

0 0.6804 0.3639 0

•

The Givens elementary matrices are orthogonal rotation matrices that al-

low for setting to zero in a selective way the entries of a vector or matrix.

For a given pair of indices i and k, and a given angle θ, these matrices are

de¬ned as

G(i, k, θ) = In ’ Y (5.43)

where Y∈ Rn—n is a null matrix except for the following entries: yii =

ykk = 1 ’ cos(θ), yik = ’ sin(θ) = ’yki . A Givens matrix is of the form

i k

®

0

1

1

..

.

cos(θ) sin(θ)

i

..

.

G(i, k, θ) = .

’ sin(θ) cos(θ) k

..

.

1

° »

0 1

T

For a given vector x ∈ Rn , the product y = (G(i, k, θ)) x is equivalent to

rotating x counterclockwise by an angle θ in the coordinate plane (xi , xk )

(see Figure 5.3, right). After letting c = cos θ, s = sin θ, it follows that

±

xj j = i, k,

cxi ’ sxk j = i,

yj = (5.44)

sxi + cxk j = k.

5.6 The QR Method for Matrices in Hessenberg Form 207

Let ±ik = x2 + x2 and notice that if c and s satisfy c = xi /±ik , s =

i k

’xk /±ik (in such a case, θ = arctan(’xk /xi )), we get yk = 0, yi = ±ik

and yj = xj for j = i, k. Similarly, if c = xk /±ik , s = xi /±ik (that is,

θ = arctan(xi /xk )), then yi = 0, yk = ±ik and yj = xj for j = i, k.

The Givens rotation matrices will be employed in Section 5.6.3 to carry

out the QR factorization step in the algorithm (5.32) and in Section 5.10.1

where the Jacobi method for symmetric matrices is considered.

Remark 5.3 (Householder de¬‚ation for power iterations) The ele-

mentary Householder tranformations can be conveniently employed to com-

pute the ¬rst (largest or smallest) eigenvalues of a given matrix A ∈ Rn—n .

Assume that the eigenvalues of A are ordered as in (5.16) and suppose

that the eigenvalue/eigenvector pair (»1 , x1 ) has been computed using the

power method. Then the matrix A can be transformed into the following

block form (see for the proof [Dat95], Theorem 8.5.4, p. 418)

bT

»1

A1 = HAH =

0 A2

where b ∈ Rn’1 , H is the Householder matrix such that Hx1 = ±x1 for

some ± ∈ R, the matrix A2 ∈ R(n’1)—(n’1) and the eigenvalues of A2 are

the same as those of A except for »1 . The matrix H can be computed using

(5.38) with v = x1 ± x1 2 e1 .

The de¬‚ation procedure consists of computing the second dominant (sub-

dominant) eigenvalue of A by applying the power method to A2 provided

that |»2 | = |»3 |. Once »2 is available, the corresponding eigenvector x2

can be computed by applying the inverse power iteration to the matrix A

taking µ = »2 (see Section 5.3.2) and proceeding in the same manner with

the remaining eigenvalue/eigenvector pairs. An example of de¬‚ation will

be presented in Section 5.12.2.

5.6.2 Reducing a Matrix in Hessenberg Form

A given matrix A∈ Rn—n can be transformed by similarity transforma-

tions into upper Hessenberg form with a cost of the order of n3 ¬‚ops. The

algorithm takes n ’ 2 steps and the similarity transformation Q can be

computed as the product of Householder matrices P(1) · · · P(n’2) . For this,

the reduction procedure is commonly known as the Householder method.

Precisely, the k-th step consists of a similarity transformation of A through

the Householder matrix P(k) which aims at setting to zero the elements in

positions k + 2, . . . , n of the k-th column of A, for k = 1, . . . , (n ’ 2) (see

208 5. Approximation of Eigenvalues and Eigenvectors

Section 5.6.1). For example, in the case n = 4 the reduction process yields

® ® ®

• ••• • ••• ••••

• • • • ’’ • • • • ’’ • • • •

P(1) P(2) ,

• • • • 0 • • • 0 • • •

° » ° » ° »

• • • • • • • • •

0 0 0

having denoted by • the entries of the matrices that are a priori non zero.

Given A(0) = A, the method generates a sequence of matrices A(k) that

are orthogonally similar to A

A(k) = PT A(k’1) P(k) = (P(k) · · · P(1) )T A(P(k) · · · P(1) )

(k)

(5.45)

k ≥ 1.

= QT AQ(k) ,

(k)

For any k ≥ 1 the matrix P(k) is given by (5.41), where x is substituted

by the k-th column vector in matrix A(k’1) . From the de¬nition (5.41) it

is easy to check that the operation PT A(k’1) leaves the ¬rst k rows of

(k)

(k’1) T (k’1)

P(k) = A(k) does the same on the

A unchanged, whilst P(k) A

¬rst k columns. After n ’ 2 steps of the Householder reduction, we obtain

a matrix H = A(n’2) in upper Hessenberg form.

Remark 5.4 (The symmetric case) If A is symmetric, the transforma-

tion (5.45) maintains such a property. Indeed

∀k ≥ 1,

(A(k) )T = (QT AQ(k) )T = A(k) ,

(k)

so that H must be tridiagonal. Its eigenvalues can be e¬ciently computed

using the method of Sturm sequences with a cost of the order of n ¬‚ops, as

will be addressed in Section 5.10.2.

A coding of the Householder reduction method is provided in Program

29. To compute the Householder vector, Program 32 is employed. In output,

the two matrices H and Q, respectively in Hessenberg form and orthogonal,

are such that H = QT AQ.

Program 29 - houshess : Hessenberg-Householder method

function [H,Q]=houshess(A)

n=max(size(A)); Q=eye(n); H=A;

for k=1:(n-2),

[v,beta]=vhouse(H(k+1:n,k)); I=eye(k); N=zeros(k,n-k);

m=length(v); R=eye(m)-beta*v*v™; H(k+1:n,k:n)=R*H(k+1:n,k:n);

H(1:n,k+1:n)=H(1:n,k+1:n)*R; P=[I, N; N™, R]; Q=Q*P;

end

5.6 The QR Method for Matrices in Hessenberg Form 209

The algorithm coded in Program 29 requires a cost of 10n3 /3 ¬‚ops and

is well-conditioned with respect to rounding errors. Indeed, the following

estimate holds (see [Wil65], p. 351)

¤ cn2 u A

H = QT (A + E) Q, E (5.46)

F F

where H is the Hessenberg matrix computed by Program 29, Q is an or-

thogonal matrix, c is a constant, u is the roundo¬ unit and · F is the

Frobenius norm (see (1.18)).

Example 5.7 Consider the reduction in upper Hessenberg form of the Hilbert

matrix H4 ∈ R4—4 . Since H4 is symmetric, its Hessenberg form should be a

triadigonal symmetric matrix. Program 29 yields the following results

® ®

1.00 0 0 0 1.00 0.65 0 0

0 0.20 0.65 0

’0.61

0.77 0.65 0.06

Q= , H= .

0 ’0.76 0 0.001

0.51 0.40 0.06 0.02

° » ° »

0 0.38 0.69 0.61 0 0 0.001 0.0003

The accuracy of the transformation procedure (5.45) can be measured by com-

puting the · F norm of the di¬erence between H and QT H4 Q. This yields

H ’ QT H4 Q F = 3.38 · 10’17 , which con¬rms the stability estimate (5.46). •

5.6.3 QR Factorization of a Matrix in Hessenberg Form

In this section we explain how to e¬ciently implement the generic step of

the QR iteration, starting from a matrix T(0) = H(0) in upper Hessenberg

form.

For any k ≥ 1, the ¬rst phase consists of computing the QR factorization

of H(k’1) by means of n ’ 1 Givens rotations

T T T

(k) (k)

(k) (k’1)

H(k’1) = R(k) ,

Q H = Gn’1 ... G1 (5.47)

(k)

where, for any j = 1, . . . , n ’ 1, Gj = G(j, j + 1, θj )(k) is, for any k ≥ 1,

the j-th Givens rotation matrix (5.43) in which θj is chosen according

to (5.44) in such a way that the entry of indices (j + 1, j) of the matrix

T T

(k) (k)

· · · G1 H(k’1) is set equal to zero. The product (5.47) requires

Gj

a computational cost of the order of 3n2 ¬‚ops.

The next step consists of completing the orthogonal similarity transfor-

mation

(k) (k)

H(k) = R(k) Q(k) = R(k) G1 . . . Gn’1 . (5.48)

210 5. Approximation of Eigenvalues and Eigenvectors

(k) (k)

The orthogonal matrix Q(k) = G1 . . . Gn’1 is in upper Hessenberg

form. Indeed, taking for instance n = 3, and recalling Section 5.6.1, we get

® ® ®

•• •••

0 100

Q(k) = G1 G2 = • • 0 0 • • = • • • .

(k) (k)

° »° »° »

0•• 0••

00 1

Also (5.48) requires a cost of the order of 3n2 operations, for an overall e¬ort

of the order of 6n2 ¬‚ops. In conclusion, performing the QR factorization

with elementary Givens rotations on a starting matrix in upper Hessenberg

form yields a reduction of the operation count of one order of magnitude

with respect to the corresponding factorization with the modi¬ed Gram-

Schmidt procedure of Section 5.5.

5.6.4 The Basic QR Iteration starting from Upper Hessenberg

Form

A basic implementation of the QR iteration to generate the real Schur

decomposition of a matrix A is given in Program 30.

This program uses Program 29 to reduce A in upper Hessenberg form;

then each QR factorization step in (5.32) is carried out with Program 31

which utilizes Givens rotations. The overall e¬ciency of the algorithm is

ensured by pre- and post-multiplying with Givens matrices as explained in

(k) (k)

Section 5.6.5, and by constructing the matrix Q(k) = G1 . . . Gn’1 in the

function prodgiv, with a cost of n2 ’ 2 ¬‚ops and without explicitly forming

(k)

the Givens matrices Gj , for j = 1, . . . , n ’ 1.

As for the stability of the QR iteration with respect to rounding er-

ror propagation, it can be shown that the computed real Schur form T is

orthogonally similar to a matrix “close” to A, i.e.

T = QT (A + E)Q

where Q is orthogonal and E u A 2 , u being the machine roundo¬

2

unit.

Program 30 returns in output, after niter iterations of the QR proce-

dure, the matrices T, Q and R in (5.32).

Program 30 - hessqr : Hessenberg-QR method

function [T,Q,R]=hessqr(A,niter)

n=max(size(A));

[T,Qhess]=houshess(A);

for j=1:niter

[Q,R,c,s]= qrgivens(T);

5.6 The QR Method for Matrices in Hessenberg Form 211

T=R;

for k=1:n-1,

T=gacol(T,c(k),s(k),1,k+1,k,k+1);

end

end

Program 31 - givensqr : QR factorization with Givens rotations

function [Q,R,c,s]= qrgivens(H)

[m,n]=size(H);

for k=1:n-1

[c(k),s(k)]=givcos(H(k,k),H(k+1,k));

H=garow(H,c(k),s(k),k,k+1,k,n);

end

R=H; Q=prodgiv(c,s,n);

function Q=prodgiv(c,s,n)

n1=n-1; n2=n-2;

Q=eye(n); Q(n1,n1)=c(n1); Q(n,n)=c(n1);

Q(n1,n)=s(n1); Q(n,n1)=-s(n1);

for k=n2:-1:1,

k1=k+1; Q(k,k)=c(k); Q(k1,k)=-s(k);

q=Q(k1,k1:n); Q(k,k1:n)=s(k)*q;

Q(k1,k1:n)=c(k)*q;

end

Example 5.8 Consider the matrix A (already in Hessenberg form)

®

3 17 ’37 18 ’40

1 0 0