In [1]:
## Householder QR
In [2]:
## Julia's QR code
A = randn(5,5)
Q,R = qr(A)
Out[2]:
([-0.216703 0.116679 … -0.2447 -0.864023; 0.595895 -0.422689 … -0.0810492 -0.411142; … ; 0.451327 0.1462 … -0.784294 0.257678; -0.590535 -0.0375122 … -0.555466 0.0548423], [2.95018 0.173269 … -1.34675 2.64297; 0.0 1.98873 … 0.833141 0.514899; … ; 0.0 0.0 … -1.08929 -0.954043; 0.0 0.0 … 0.0 1.72833])
In [3]:
## A simple
"""
`householder_qr_simple`
A pedagogical code for QR via householder reflections
"""
function householder_qr_simple(Ain)
  A = copy(Ain)
  m,n = size(A)

  # transform A into R by working in-place
  Q = eye(m)
  for i=1:n-1 # we always have Ain = Q*A
    ai = A[i:end,i] # extract ith column
    w = -ai
    w[1] = w[1] + norm(ai)
    w = w/norm(w)
    H = eye(m-i+1)-2*w*w' # form reflector
    Q[:,i:end] = Q[:,i:end]*H'
    A[i:end,i:end] = H*A[i:end,i:end]
    @printf("%4i : norm(Q*A - Ain) = %.16e\n",
      i, norm(Q*A - Ain))
  end
  R = triu(A)
  return Q, R
end

Q1,R1 = householder_qr_simple(A)
@show norm(Q1'*Q1 - I)
@show norm(A - Q1*R1)
  
In [4]:
##
# Let's do a runtime analysis of this code.
In [5]:
## Here's a better version, but this version encodes A itself.
function householder_qr(Ain)
  A = copy(Ain)
  n = size(A,2)
  for i = 1:n
    # compute householder vector
    ai = A[i:end,i]
    ai1 = ai[1]
    v = copy(ai)
    v[1] = v[1] + sign(ai[1])*norm(ai)
    v = v/norm(v)
    A[i:end,i] = v
    w = A[i:end,(i+1):end]'*v

    # update the lower part of the matrix
    A[i:end,(i+1):end] = A[i:end,(i+1):end] - 2*v*w'
    A[i,i] = - norm(ai)*sign(ai1)
  end
  R = triu(A)
  return A,R
end

# In julia
A = randn(5,5)
Qdata,R = householder_qr(A)
qrfact!(A) # this is exactly what Julia can compute in it's super-optimal QR form.

display(Qdata)
display(A)
5×5 Array{Float64,2}:
  2.05679   -0.678643   0.412513  -0.522338   0.0914413
 -0.138869  -1.2369    -1.26101   -0.139775   0.872205 
  0.105297  -0.237646  -2.03053   -0.29942    0.965397 
  0.186558   0.504497  -0.546458  -0.181808  -0.412213 
 -0.150333  -0.337065   0.371005  -0.605724  -0.212322 
5×5 Array{Float64,2}:
  2.05679   -0.678643   0.412513  -0.522338   0.0914413
  0.145397  -1.2369    -1.26101   -0.139775   0.872205 
 -0.110247  -0.313292  -2.03053   -0.29942    0.965397 
 -0.195328   0.665084  -0.727809  -0.181808  -0.412213 
  0.1574    -0.444357   0.494129  -0.761271   0.212322