Page 405 - Applied Numerical Methods Using MATLAB

394 MATRICES AND EIGENVALUES
into

y = x 1 ·· · x k−1 −g k 0 ·· · 0 (P8.4.8)
where g k is fixed in such a way that the norms of these two vectors are
the same:

N

g k = x n 2 (P8.4.9)
n=k
First, we find the difference vector of unit norm as
1
w k = (x − y)
c

1
= 0 ··· 0 x k + g k x k+1 ··· x N (P8.4.10)
c
with

2
c =||x − y|| 2 = (x k + g k ) + x 2 +· · · + x 2 (P8.4.11)
k+1 N
Then, one more thing we should do is to substitute this difference vector
into Eq. (P8.4.1).
T
H k = [I − 2w k w ] (P8.4.12)
k
Complete the following routine “Householder()” by permuting the
statements and try it with k = 1, 2, 3, and 4 for a four-dimensional
vector generated by the MATLAB command rand(5,1) to check if it
works fine.

>> x = rand(5,1), for k = 1:4, householder(x,k)*x, end

function H = Householder(x,k)
%Householder transform to zero out tail part starting fromk+1
H = eye(N) - 2*w*w’; %Householder matrix
N = length(x);
w = zeros(N,1);
w(k) =(x(k) + g)/c; w(k + 1:N) = x(k + 1:N)/c; %Eq.(P8.4.10)
tmp = sum(x(k + 1:N).^ 2);
c = sqrt((x(k) + g)^2 + tmp); %Eq.(P8.4.11)
g = sqrt(x(k)^2 + tmp); %Eq.(P8.4.9)

(c) QR Factorization Using Householder Transform
We can use Householder transform to zero out the part under the main
diagonal of each column of an N × N matrix A successively and then
make it a lower triangular matrix R in (N − 1) iterations. The necessary
operations are collectively written as
H N−1 H N−2 ··· H 1 A = R (P8.4.13)