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

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)
  t2 = @elapsed  begin
    H, α = householder!(A2)
    x2 = solve_householder!(b2, H, α)
  end
  @show norm(A' * A * x1 .- A' * b)
  @show norm(A' * A * x2 .- A' * b)
  @show t2 / t1
end