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