QR decomposition via Givens rotations

QRGivens.jl

using LinearAlgebra

"Compute Givens rotation parameters to zero out b"
function givens_rotation(a, b)
  T = promote_type(typeof(a), typeof(b))
  if iszero(b)
    return one(T), zero(T)
  elseif abs(b) > abs(a)
    t = a / b
    s = 1 / sqrt(1 + t^2)
    c = s * t
    return c, s
  else
    t = b / a
    c = 1 / sqrt(1 + t^2)
    s = c * t
    return c, s
  end
end


"""
A Givens rotation killing A[j,k] using rows i and j.
(i < j)
"""
struct GivensRotation{T}
  i::Int
  j::Int
  c::T
  s::T
end

struct MPIQRGivensStruct{M<:AbstractMatrix, T}
  R::M
  rotations::Vector{GivensRotation{T}}
end

"""
Compute Givens QR.

Reusable solver:
1. Factorize A once → R, rotations
2. For any b:
     y = Qᵀ b (using rotations)
     x = R⁻¹ y

Returns MPIQRGivensStruct containing:
  - R matrix
  - List of GivensRotation objects (to apply Qᵀ to new RHS vectors)
"""
function qr_factorize_givens!(A::AbstractMatrix)
  m, n = size(A)

  rotations = GivensRotation{eltype(A)}[]

  @inbounds for j in 1:min(m, n) # for all columns
    for i = m:-1:(j+1) # for all rows from bottom up to and including j
      ajj = A[j, j]
      aij = A[i, j]
      aij == 0.0 && continue
      c, s = givens_rotation(ajj, aij)
      push!(rotations, GivensRotation(i, j, c, s))
      # Apply rotation to matrix rows j and i
      for k in j:n # now apply to all columns j and to the right of j
        ajk = A[j, k]
        aik = A[i, k]
        A[j, k] =  c * ajk + s * aik
        A[i, k] = -s * ajk + c * aik
      end
    end
  end
  return MPIQRGivensStruct(A, rotations)
end
qr_factorize_givens(A_in::AbstractMatrix) = qr_factorize_givens!(deepcopy(A_in))



"""
Apply stored Givens rotations to a RHS vector to get y = Qᵀ b.
(rotations must be in the exact order they were generated)
"""
function apply_Qt!(rotations, b)
  @inbounds for rot in rotations
    for j in axes(b, 2)
      bj = b[rot.j, j]
      bi = b[rot.i, j]
      b[rot.j, j] =  rot.c * bj + rot.s * bi
      b[rot.i, j] = -rot.s * bj + rot.c * bi
    end
  end
  return b
end
apply_Qt(rotations, b) = apply_Qt!(rotations, deepcopy(b))

function back_substitute(R::AbstractMatrix, y)
  T = promote_type(eltype(R), eltype(y)) 
  n = size(R, 2)
  x = zeros(T, n, size(y, 2))
  s = zeros(T, size(y, 2))
  for i in n:-1:1
    fill!(s, 0)
    @inbounds for j = i+1:n, k = 1:size(y, 2)
      s[k] += R[i, j] * x[j, k]
    end
    @inbounds for k in 1:size(y, 2)
      x[i, k] = (y[i, k] - s[k]) / R[i, i]
    end
  end
  return x
end


function qr_solve(qrG, b)
  y = apply_Qt(qrG.rotations, b)
  n = size(qrG.R, 2)
  y_small = view(y, 1:n, :)         # only first n entries
  return back_substitute(view(qrG.R, 1:n, 1:n), y_small)
end

A = randn(ComplexF64, 6, 4)
x_true = randn(ComplexF64, 4, 2)
b = A * x_true

# Factor once
qrG = qr_factorize_givens(A)

# Solve many times cheaply
x_est = qr_solve(qrG, b)

println("residual = ", norm(A * x_est - b))