QR decomposition via Givens rotations

GivensQR.jl

using LinearAlgebra
function qr_factorization(A)
  m, n = size(A)
  
  Q = zeros(float(eltype(A)), m, m)
  for i in 1:m Q[i, i] = 1 end  
  R = zeros(float(eltype(A)), m, n)
  R .= A
  G = zeros(float(eltype(A)), 2, 2)
  
  for j in 1:m, i in j:n-1 # the iteration over these (i,j) pairs can take place in any order!
    givens_rotation!(G, R[i,j], R[i+1,j])
    R[i:i+1, :] .= G' * R[i:i+1, :]
    Q[:, i:i+1] .= Q[:, i:i+1] * G
  end
  return Q, R
end

function givens_rotation!(G, a::T, b::T) where T
  (c, s) = if abs(b) < eps(float(real(T)))
    (1.0, 0.0)
  elseif abs(a) < abs(b)
    r = a / b
    s = 1 / sqrt(1 + r^2)
    (s * r, s)
  else
    r = b / a 
    c = 1 / sqrt(1 + r^2)
    (c, c * r)
  end
  G[1, 1] = G[2, 2] = c
  G[1, 2] = s
  G[2, 1] = -s
  return G
end

# Example usage
A = [1.0 2 3; 4 5 6; 7 8 9; 10 11 12]
Q, R = qr_factorization(A)
println("Q:", Q)
println("R:", R)
using Test
@testset "Givens" begin
@test norm(Q * R .- A) < 1e-10
end

qrA = qr(A)
@show Matrix(qrA.Q)
@show Matrix(qrA.R)