Add Singular Value Decomposition and Cholesky Factorization

- Implement SingularValueDecompositionProblem and SingularValueDecompositionSolution
- Create kernel method for SVD using eigen and QR decomposition
- Add CholeskyFactorizationProblem and CholeskyFactorizationSolution
- Implement Cholesky factorization with lower triangular matrix
- Update test suites to verify SVD and Cholesky decomposition
- Export new decomposition types for easy access

Author: yangjunjie0320

src/eig.jl

@@ -116,3 +116,52 @@ end
 
 kernel(prob::EigenProblemMixin) = Eigen.kernel(prob)
 export Eigen
+
+struct SingularValueDecompositionProblem <: ProblemMixin
+    a::AbstractMatrix
+    tol::Real    
+end
+
+struct SingularValueDecompositionSolution <: SolutionMixin
+    u::AbstractMatrix
+    s::AbstractMatrix
+    v::AbstractMatrix
+
+    function SingularValueDecompositionSolution(u, s, v)
+        s = Diagonal(s)
+        @assert isapprox(I, u' * u) "u must be orthogonal"
+        @assert isapprox(I, v' * v) "v must be orthogonal"
+
+        soln = new(u, s, v)
+        return soln
+    end
+end
+
+# we cheat by importing eigen and qr from LinearAlgebra
+import LinearAlgebra: eigen, qr
+function kernel(prob::SingularValueDecompositionProblem)
+    a = prob.a
+    tol = prob.tol
+   
+    n = size(a, 1)
+    e, c = eigen(a' * a)
+    q, r  = qr(a * c) |> (x -> (x.Q, x.R))
+
+    s = zeros(n)
+    rs = zeros(n, n)
+
+    for (i, ei) in enumerate(e)
+        @assert ei >= (tol * tol) "s2 must be non-negative"
+        s[i] = sqrt(ei)
+        rs[:, i] = r[:, i] / s[i]
+    end
+
+    u = Matrix(q)
+    v = c * rs
+
+    soln = SingularValueDecompositionSolution(u, s, v)
+    return soln
+end
+
+SingularValueDecomposition = SingularValueDecompositionProblem
+export SingularValueDecomposition

src/lu.jl

@@ -152,3 +152,43 @@ kernel(prob::LUFactorizationProblemMixin) = LU.kernel(prob)
 
 LUFactorization = LU.GaussianEliminationV2
 export LUFactorization
+
+struct CholeskyFactorizationProblem <: ProblemMixin
+    a::AbstractMatrix
+    tol::Real
+    function CholeskyFactorizationProblem(a, tol=1e-10)
+        prob = new(a, tol)
+        return prob
+    end
+end
+
+struct CholeskyFactorizationSolution <: SolutionMixin
+    l::AbstractMatrix
+    function CholeskyFactorizationSolution(l)
+        @assert istril(l) "matrix must be lower triangular"
+        soln = new(l)
+        return soln
+    end
+end
+
+function kernel(prob::CholeskyFactorizationProblem)
+    tol = prob.tol
+    l = copy(prob.a)
+    n = size(l, 1)
+
+    for i in 1:n
+        @assert l[i, i] > tol "Cholesky factorization failed"
+        l[i, i] = sqrt(l[i, i])
+        l[i+1:n, i] /= l[i, i]
+
+        li = l[i+1:n, i]
+        l[i+1:n, i+1:n] -= li * li'
+    end
+
+    soln = CholeskyFactorizationSolution(tril(l))
+    return soln
+end
+
+CholeskyFactorization = CholeskyFactorizationProblem
+CholeskyDecomposition = CholeskyFactorization
+export CholeskyFactorization, CholeskyDecomposition

test/eig.jl

@@ -7,8 +7,8 @@ import SimpleLinearAlgebra.Eigen
     # Make the first eigenvalue significantly larger
     a = a * a'  # Create a symmetric positive-definite matrix
 
-    prob = Eigen.PowerMethod(a, 1000, 1e-8)
-    tol = sqrt(prob.tol * n)
+    prob = Eigen.PowerMethod(a, 1000, 1e-10)
+    tol = sqrt(prob.tol) * n
 
     soln = kernel(prob)
     e = soln.e[1]
@@ -24,12 +24,28 @@ end
     # Make the first eigenvalue significantly larger
     a = a * a'  # Create a symmetric positive-definite matrix
 
-    prob = Eigen.RayleighQuotient(a, 1000, 1e-8)
-    tol = sqrt(prob.tol * n)
+    prob = Eigen.RayleighQuotient(a, 1000, 1e-10)
+    tol = sqrt(prob.tol) * n
 
     soln = kernel(prob)
     e = soln.e[1]
     c = soln.c[:, 1]
 
     @test isapprox(a * c, e * c, atol=tol)
 end
+
+@testset "Singular Value Decomposition" begin
+    n = 4
+    a = randn(n, n)
+    a = a * a'  # Create a symmetric positive-definite matrix
+
+    prob = SingularValueDecomposition(a, 1e-8)
+    tol = prob.tol
+
+    soln = kernel(prob)
+    u = soln.u
+    s = soln.s
+    v = soln.v
+
+    @test isapprox(a, u * s * v', atol=tol)
+end

test/lu.jl

@@ -78,3 +78,18 @@ end
     @test isapprox(pa, l * u, atol=tol)  # Verify PA = LU
     @test isapprox(p * a, l * u, atol=tol)
 end
+
+@testset "Cholesky factorization" begin
+    n = 10
+    q, r = rand(n, n) |> qr
+    a = q * Diagonal(rand(n)) * q'
+
+    prob = CholeskyFactorization(a)
+    tol = prob.tol
+
+    soln = kernel(prob)
+    l = soln.l
+    
+    @test istril(l)
+    @test isapprox(a, l * l', atol=tol) 
+end