<< . .

. 31
( : 95)



. . >>


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

<< . .

. 31
( : 95)



. . >>