Page 406 - Applied Numerical Methods Using MATLAB

PROBLEMS 395
which implies that
−1
A = [H N−1 H N−2 ··· H 1 ] R = H −1 ·· · H −1 H −1 R
1 N−2 N−1
= H 1 ·· · H N−2 H N−1 R = QR (P8.4.14)

where the product of all the Householder matrices
Q = H 1 ··· H N−2 H N−1 (P8.4.15)

turns out to be not only symmetric, but also orthogonal like each H k :
T T
Q Q = [H 1 ·· · H N−2 H N−1 ] H 1 ·· · H N−2 H N−1
T T T
= H H ·· · H H 1 ·· · H N−2 H N−1 = I
N−1 N−2 1
This suggests a QR factorization method that is cast into the following
routine “qr_my()”. You can try it for a nonsingular 3 × 3matrixgener-
ated by the MATLAB command rand(3) and compare the result with
that of the MATLAB built-in routine “qr()”.

function [Q,R] = qr_my(A)
%QR factorization
N = size(A,1); R = A; Q = eye(N);
for k = 1:N - 1
H = Householder(R(:,k),k);
R = H*R; %Eq.(P8.4.13)
Q = Q*H; %Eq.(P8.4.15)
end

8.5 Hessenberg Form Using Householder Transform

function [Hs,HH] = Hessenberg(A)
%Transform into an almost upper triangular matrix
% having only zeros below lower subdiagonal
N = size(A,1); Hs = A; HH = eye(N); %HH*A*HH’ = Hs
fork=1:N-2
H = Householder(Hs(:,k), );
Hs = H*Hs*H; HH = H*HH;
end

We can make use of Householder transform (introduced in Problem 8.4) to
zero-out the elements below the lower subdiagonal of a matrix so that it
becomes an upper Hessenberg form which is almost upper-triangular matrix.
Complete the above routine “Hessenberg()” by filling in the second input
argument of the routine “Householder()” and try it for a 5 × 5matrix
generated by the MATLAB command rand(5) to check if it works.