Householder QR¶
We will use Householder reflections to produce a QR factorization of the matrix
using LinearAlgebra
A = rand(1.:9.,6,4)
m,n = size(A)
(6, 4)
Our first step is to introduce zeros below the diagonal in column 1. Define the vector
z = A[:,1];
Applying the Householder definitions gives us
v = z - norm(z)*[1;zeros(m-1)]
P = I - 2/(v'*v)*(v*v') # reflector
6×6 Array{Float64,2}:
0.372678 0.0745356 0.521749 0.596285 0.298142 0.372678
0.0745356 0.991144 -0.0619919 -0.0708479 -0.0354239 -0.0442799
0.521749 -0.0619919 0.566057 -0.495935 -0.247968 -0.30996
0.596285 -0.0708479 -0.495935 0.433217 -0.283392 -0.354239
0.298142 -0.0354239 -0.247968 -0.283392 0.858304 -0.17712
0.372678 -0.0442799 -0.30996 -0.354239 -0.17712 0.7786
(Julia automatically fills in an identity of the correct size for the I above.) By design we can use the reflector to get the zero structure we seek:
P*z
6-element Array{Float64,1}:
13.416407864998739
1.6653345369377348e-16
4.440892098500626e-16
1.3322676295501878e-15
8.881784197001252e-16
4.440892098500626e-16
Now we let
A = P*A
6×4 Array{Float64,2}:
13.4164 12.8201 14.3854 11.2549
8.32667e-17 6.59559 3.24132 6.01919
7.77156e-16 -3.8309 2.68923 1.13435
1.11022e-15 -5.23531 0.93055 -2.84646
8.88178e-16 0.382343 4.96528 -1.92323
5.55112e-16 1.97793 -2.79341 0.0959628
We are set to put zeros into column 2. We must not use row 1 in any way, lest it destroy the zeros we just introduced. To put it another way, we can repeat the process we just did on the smaller submatrix
A[2:m,2:n]
5×3 Array{Float64,2}:
6.59559 3.24132 6.01919
-3.8309 2.68923 1.13435
-5.23531 0.93055 -2.84646
0.382343 4.96528 -1.92323
1.97793 -2.79341 0.0959628
z = A[2:m,2];
v = z - norm(z)*[1;zeros(m-2)];
P = I - 2/(v'*v)*(v*v');
We now apply the reflector to the submatrix.
A[2:m,2:n] = P*A[2:m,2:n]
A
6×4 Array{Float64,2}:
13.4164 12.8201 14.3854 11.2549
8.32667e-17 9.46808 0.27226 5.25039
7.77156e-16 5.6983e-16 -1.27046 0.109037
1.11022e-15 1.1591e-15 -4.48077 -4.24765
8.88178e-16 -8.83814e-17 5.36047 -1.8209
5.55112e-16 -5.0072e-16 -0.748982 0.62534
We need two more iterations of this process.
for j = 3:n
z = A[j:m,j]
v = z - norm(z)*[1;zeros(m-j)]
P = I - 2/(v'*v)*(v*v')
A[j:m,j:n] = P*A[j:m,j:n]
end
We have now reduced the original to an upper triangular matrix using four orthogonal Householder reflections:
R = A
6×4 Array{Float64,2}:
13.4164 12.8201 14.3854 11.2549
8.32667e-17 9.46808 0.27226 5.25039
7.77156e-16 5.6983e-16 7.14052 1.21349
1.11022e-15 1.1591e-15 -2.23307e-16 4.50429
8.88178e-16 -8.83814e-17 3.23365e-16 -5.05821e-16
5.55112e-16 -5.0072e-16 -6.93346e-17 1.73333e-16