Golub Kahan Bidiagonalization and Corresponding Tests

Project.toml

@@ -1,7 +1,17 @@
 name = "RidgeRegression"
 uuid = "739161c8-60e1-4c49-8f89-ff30998444b1"
-authors = ["Vivak Patel <vp314@users.noreply.github.com>"]
 version = "0.1.0"
+authors = ["Eton Tackett <etont@icloud.com>", "Vivak Patel <vp314@users.noreply.github.com>"]
+
+[deps]
+CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
+Downloads = "f43a241f-c20a-4ad4-852c-f6b1247861c6"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [compat]
+CSV = "0.10.15"
+DataFrames = "1.8.1"
+Downloads = "1.7.0"
+LinearAlgebra = "1.12.0"
 julia = "1.12.4"

docs/make.jl

@@ -14,6 +14,7 @@ makedocs(;
     ),
     pages=[
         "Home" => "index.md",
+        "Design" => "design.md",
     ],
 )
 

src/RidgeRegression.jl

@@ -1,5 +1,10 @@
 module RidgeRegression
+using CSV
+using DataFrames
+using LinearAlgebra
 
-# Write your package code here.
+include("bidiagonalization.jl")
+
+export compute_givens, rotate_rows!, rotate_cols!, bidiagonalize_with_H, apply_Ht_to_b
 
 end

src/bidiagonalization.jl

@@ -0,0 +1,165 @@
+"""
+    compute_givens(a, b)
+
+Compute Givens rotation coefficients for scalars `a` and `b`.
+
+# Arguments
+- `a::Number`: First scalar
+- `b::Number`: Second scalar
+
+# Returns
+- `c::Number`: Cosine coefficient of the Givens rotation
+- `s::Number`: Sine coefficient of the Givens rotation
+
+"""
+function compute_givens(a::Number, b::Number) # Compute Givens rotation coefficients for scalars a and b
+    if b == 0
+        return one(a), zero(a)
+    elseif a == 0
+        throw(ArgumentError("a is zero, cannot compute Givens rotation"))
+    else
+        r = hypot(a, b)
+        c = a / r
+        s = b / r
+        return c, s
+    end
+end
+
+"""
+    rotate_rows!(M::AbstractMatrix,i::Int,j::Int,c::Number,s::Number)
+
+Apply a Givens rotation to rows `i` and `j` of matrix `M`.
+
+# Arguments
+- `M::AbstractMatrix`: The matrix to be rotated
+- `i::Int`: First row index
+- `j::Int`: Second row index
+- `c::Number`: Cosine coefficient of the Givens rotation
+- `s::Number`: Sine coefficient of the Givens rotation
+
+# Returns
+- `M::AbstractMatrix`: The rotated matrix
+
+"""
+function rotate_rows!(M::AbstractMatrix,i::Int,j::Int,c::Number,s::Number)
+    for k in 1:size(M,2) # Loop over columns
+        temp = M[i,k] # Store the original value of M[i,k] before modification
+        M[i,k] = c*temp + s*M[j,k]
+        M[j,k] = -conj(s)*temp + c*M[j,k] #Apply the Givens rotation to the elements in rows i and j
+    end
+    return M
+end
+
+
+"""
+    rotate_cols!(M::AbstractMatrix,i::Int,j::Int,c::Number,s::Number)
+
+Apply a Givens rotation to columns `i` and `j` of matrix `M`.
+
+# Arguments
+- `M::AbstractMatrix`: The matrix to be rotated
+- `i::Int`: First column index
+- `j::Int`: Second column index
+- `c::Number`: Cosine coefficient of the Givens rotation
+- `s::Number`: Sine coefficient of the Givens rotation
+
+# Returns
+- `M::AbstractMatrix`: The rotated matrix
+
+"""
+function rotate_cols!(M::AbstractMatrix,i::Int,j::Int,c::Number,s::Number)
+    for k in 1:size(M,1) # Loop over rows
+        temp = M[k,i] # Store the original value of M[k,i] before modification
+        M[k,i] = c*temp + s*M[k,j]
+        M[k,j] = -conj(s)*temp + c*M[k,j] # Apply the Givens rotation to the elements in columns i and j
+    end
+    return M
+end
+"""
+    bidiagonalize_with_H(A::AbstractMatrix, L::AbstractMatrix, b::AbstractVector)
+
+Performs bidiagonalization of a matrix A using a sequence of Givens transformations while explicitly accumulating
+the orthogonal left factor `H` and right factor `K` such that
+
+    H' * A * K = B.
+
+# Arguments
+- `A::AbstractMatrix`: The matrix to be bidiagonalized
+- `L::AbstractMatrix`: The banded matrix to be updated in-place
+- `b::AbstractVector`: The vector to be transformed
+
+# Returns
+- `B::AbstractMatrix`: The bidiagonalized form of the input matrix A with dimension (n,n)
+- `C::AbstractMatrix`: The matrix resulting from the sequence of Givens transformations
+- `H::AbstractMatrix`: The orthogonal left factor
+- `K::AbstractMatrix`: The orthogonal right factor
+- `Ht::AbstractMatrix`: The transpose of the orthogonal left factor
+- `b::AbstractVector`: The vector resulting from applying the Givens transformations to the input vector b
+
+"""
+
+function bidiagonalize_with_H(A::AbstractMatrix, L::AbstractMatrix, b::AbstractVector)
+    m, n = size(A)
+
+    B = copy(A) #Will be transformed into bidiagonal form
+    C = copy(L)
+    bhat = copy(b) #Will recieve same left transformations applied to A
+
+    Ht = Matrix{eltype(A)}(I, m, m) #Ht will accumulate the left transformations, initialized as identity
+    K  = Matrix{eltype(A)}(I, n, n) #K will accumulate the right transformations, initialized as identity
+
+    imax = min(m, n)
+
+    for i in 1:imax
+        # Left Givens rotations: zero below diagonal in column i
+        for j in m:-1:(i + 1)
+            if B[j, i] != 0
+                c, s = compute_givens(B[i, i], B[j, i]) #Build Givens rotation to zero B[j, i]
+                rotate_rows!(B, i, j, c, s) #Apply the Givens rotation to rows i and j of B
+                rotate_rows!(Ht, i, j, c, s) #Accumulate the left transformations in Ht
+                B[j, i] = zero(eltype(B))
+            end
+        end
+
+        # Right Givens rotations: zero entries right of superdiagonal
+        if i <= n - 2
+            for k in n:-1:(i + 2)
+                if B[i, k] != 0
+                    c, s = compute_givens(B[i, i + 1], B[i, k]) #Build Givens rotation to zero B[i, k]
+                    #s = -s
+                    rotate_cols!(B, i + 1, k, c, s) #Apply the Givens rotation to columns i+1 and k of B
+                    rotate_cols!(C, i + 1, k, c, s) #Apply the same right rotation to C, since C is updated by the right transformations
+                    rotate_cols!(K, i + 1, k, c, s) #Accumulate the right transformations in K
+                    B[i, k] = zero(eltype(B))
+                end
+            end
+        end
+    end
+
+    H = transpose(Ht)
+    bhat = apply_Ht_to_b(Ht, b)
+    return B, C, H, K, Ht, bhat
+end
+
+"""
+    apply_Ht_to_b(Ht::AbstractMatrix, b::AbstractVector)
+
+Apply the accumulated left orthogonal transformation `H'` (stored as `Ht`)
+to the constant vector `b`.
+
+# Arguments
+- `Ht::AbstractMatrix`: The transpose of the orthogonal left factor `H`.
+- `b::AbstractVector`: The vector to be transformed.
+
+# Returns
+- `bhat::AbstractVector`: The transformed vector satisfying `bhat = Ht * b`.
+
+# Throws
+- `DimensionMismatch`: If the number of columns of `Ht` does not match the length of `b`.
+"""
+function apply_Ht_to_b(Ht::AbstractMatrix, b::AbstractVector)
+    size(Ht, 2) == length(b) || throw(DimensionMismatch(
+        "Ht has $(size(Ht, 2)) columns but b has length $(length(b))"
+    ))
+    return Ht * b
+end

test/Project.toml

@@ -1,2 +1,5 @@
 [deps]
+CSV = "336ed68f-0bac-5ca0-87d4-7b16caf5d00b"
+DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

test/apply_Ht_to_b_test.jl

@@ -0,0 +1,34 @@
+@testset "Testset 1" begin
+    Ht = Matrix{Float64}(I, 3, 3)
+    b  = [1.0, 2.0, 3.0]
+
+    @test apply_Ht_to_b(Ht, b) == b
+end
+
+@testset "Testset 2" begin
+    Ht =[1.0  0.0  0.0;
+        0.0  0.0  1.0;
+        0.0  1.0  0.0]
+
+    b = [4.0, 5.0, 6.0]
+
+    @test apply_Ht_to_b(Ht, b) == [4.0, 6.0, 5.0]
+end
+
+@testset "Testset 3" begin
+    c, s = 3/5, 4/5
+
+    Ht =[c   s;
+        -s   c]
+
+    b = [5.0, 0.0]
+
+    @test apply_Ht_to_b(Ht, b) ≈ [3.0, -4.0]
+end
+
+@testset "Testset 4" begin
+    Ht = Matrix{Float64}(I, 3, 3)
+    b  = [1.0, 2.0]
+
+    @test_throws DimensionMismatch apply_Ht_to_b(Ht, b)
+end

test/bidiagonalize_with_H_tests.jl

@@ -0,0 +1,101 @@
+@testset "Testset 4" begin
+    A = [1.0 2.0 3.0;
+         4.0 5.0 6.0;
+         7.0 8.0 10.0]
+
+    L = Matrix{Float64}(I, 3, 3)
+    b = [1.0, 2.0, 3.0]
+
+    B, C, H, K, Ht, bhat = bidiagonalize_with_H(A, L, b)
+
+    # Ht should be the transpose of H
+    @test Ht ≈ H'
+
+    # Orthogonality of H and K
+    @test H' * H ≈ Matrix{Float64}(I, 3, 3)
+    @test H * H' ≈ Matrix{Float64}(I, 3, 3)
+    @test K' * K ≈ Matrix{Float64}(I, 3, 3)
+    @test K * K' ≈ Matrix{Float64}(I, 3, 3)
+
+    # Main factorization identity
+    @test H' * A * K ≈ B
+
+    @test H' * b ≈ bhat
+
+    # Applies right transforms to L, implicitly assuming J = I
+    @test L * K ≈ C
+end
+
+@testset "Testset 5" begin
+    A = [2.0 1.0 3.0;
+         4.0 5.0 6.0;
+         7.0 8.0 9.0]
+
+    L = Matrix{Float64}(I, 3, 3)
+    b = [1.0, -1.0, 2.0]
+
+    B, C, H, K, Ht, bhat = bidiagonalize_with_H(A, L, b)
+
+    m, n = size(B)
+
+    for i in 1:m
+        for j in 1:n
+            # In an upper bidiagonal matrix, only diagonal and superdiagonal may be nonzero
+            if !(j == i || j == i + 1)
+                @test isapprox(B[i, j], 0.0; atol=1e-10, rtol=0)
+            end
+        end
+    end
+end
+
+@testset "Testset 6" begin
+    A = [1.0 2.0 3.0;
+         4.0 5.0 6.0;
+         7.0 8.0 9.0;
+         1.0 0.0 1.0]
+
+    L = Matrix{Float64}(I, 3, 3)
+    b = [1.0, 2.0, 3.0, 4.0]
+
+    B, C, H, K, Ht, bhat = bidiagonalize_with_H(A, L, b)
+
+    m, n = size(A)
+
+    @test size(B) == (m, n)
+    @test size(C) == size(L)
+    @test size(H) == (m, m)
+    @test size(K) == (n, n)
+    @test size(Ht) == (m, m)
+    @test length(bhat) == length(b)
+
+    @test H' * H ≈ Matrix{Float64}(I, m, m)
+    @test K' * K ≈ Matrix{Float64}(I, n, n)
+
+    @test H' * A * K ≈ B
+    @test H' * b ≈ bhat
+    @test L * K ≈ C
+
+    for i in 1:m
+        for j in 1:n
+            if !(j == i || j == i + 1)
+                @test isapprox(B[i, j], 0.0; atol=1e-10, rtol=0)
+            end
+        end
+    end
+end
+
+@testset "Testset 7" begin
+    A = [3.0 0.0;
+         4.0 5.0]
+
+    L = Matrix{Float64}(I, 2, 2)
+    b = [1.0, 2.0]
+
+    B, C, H, K, Ht, bhat = bidiagonalize_with_H(A, L, b)
+
+    @test H' * A * K ≈ B
+    @test H' * b ≈ bhat
+    @test L * K ≈ C
+
+    @test isapprox(B[2,1], 0.0; atol=1e-12, rtol=0)
+end

test/compute_givens_test.jl

@@ -0,0 +1,13 @@
+@testset "Testset 1" begin
+    c, s = compute_givens(3.0, 0.0)
+    @test c == 1.0
+    @test s == 0.0
+    a = 3.0
+    b = 4.0
+    c, s = compute_givens(a, b)
+    v1 = c*a + s*b
+    v2 = -s*a + c*b
+    @test isapprox(v2, 0.0; atol=1e-12, rtol=0)
+    @test isapprox(abs(v1), hypot(a, b); atol=1e-12, rtol=0)
+    @test_throws ArgumentError compute_givens(0.0, 2.0)
+end

test/rotate_cols_test.jl

@@ -0,0 +1,9 @@
+@testset "Testset 1" begin
+    M = [3.0 4.0;
+         1.0 2.0]
+
+    c, s = compute_givens(3.0, 4.0)
+    rotate_cols!(M, 1, 2, c, s)
+
+    @test isapprox(M[1, 2], 0.0; atol=1e-12, rtol=0)
+end

test/rotate_rows_test.jl

@@ -0,0 +1,9 @@
+@testset "rotate_rows!" begin
+    M = [1.0 2.0;
+         3.0 4.0]
+
+    c, s = compute_givens(1.0, 3.0)
+    rotate_rows!(M, 1, 2, c, s)
+
+    @test isapprox(M[2, 1], 0.0; atol=1e-12, rtol=0)
+end