Solve a square or least square linear system with the QR householder reflector algorithm

HouseholderQR.jl

using LinearAlgebra, Random, Base.Threads

Random.seed!(0)
alphafactor(x::Real) = -sign(x)
alphafactor(x::Complex) = -exp(im * angle(x))

function householder!(A::Matrix{T}) where T
  m, n = size(A)
  H = A
  α = zeros(T, min(m, n)) # Diagonal of R
  @inbounds @views for j in 1:n
    s = norm(H[j:m, j])
    α[j] = s * alphafactor(H[j, j])
    f = 1 / sqrt(s * (s + abs(H[j, j])))
    H[j,j] -= α[j]
    for i in j:m
      H[i, j] *= f
    end
    for jj in j+1:n
      s = dot(H[j:m, j], H[j:m, jj])
      for i in j:m
        H[i, jj] -= H[i, j] * s
      end
    end
  end
  return (H, α)
end

function solve_householder!(b, H, α)
  m, n = size(H)
  # multuply by Q' ...
  @inbounds @views for j in 1:n
    s = dot(H[j:m, j], b[j:m])
    for i in j:m
      b[i] -= H[i, j] * s
    end
  end
  # now that b holds the value of Q'b
  # we may back sub with R
  @inbounds @views for i in n:-1:1
    for j in i+1:n
      b[i] -= H[i, j] * b[j]
    end
    b[i] /= α[i]
  end
  return b[1:n]
end

function solve_householder!(x, b, H, α)
  m, n = size(H)
  # multuply by Q' ...
  x .= b[1:length(x)]
  c = b[n+1:m] # work array
  for j in 1:n
    s = dot(H[j:n, j], x[j:n])
    s += dot(H[n+1:m, j], c)
    for i in j:n
      x[i] -= H[i, j] * s
    end
    for i in n + 1:m
      c[i - n] -= H[i, j] * s
    end
  end
  # now that b holds the value of Q'b
  # we may back sub with R
  @inbounds @views for i in n:-1:1
    for j in i+1:n
      x[i] -= H[i, j] * x[j]
    end
    x[i] /= α[i]
  end
  return x
end


for T in (Float64, ComplexF64), mn in ((50, 50), (60, 50), (500, 500), (600, 500),
                                       (1200, 1000), (2400, 2000))
m, n = mn
@show T, m, n
A = rand(T, m,n)
A .= real.(A)
b = rand(T,m)
A1 = deepcopy(A)
b1 = deepcopy(b)
t1 = @elapsed x1 = qr!(A1, NoPivot()) \ b1
A2 = deepcopy(A)
b2 = deepcopy(b)
b3 = deepcopy(b)
t2 = @elapsed  begin
  H, α = householder!(A2)
  x2 = solve_householder!(b2, H, α)
  x3 = similar(x2) .* 0
  solve_householder!(x3, b3, H, α)
end
@show norm(A' * A * x1 .- A' * b)
@show norm(A' * A * x2 .- A' * b)
@show norm(A' * A * x3 .- A' * b)
@show t2 / t1
end