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