# 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