Add symmetric QR factorization method

- Implement SymmetricQRFactorizationProblem for symmetric matrices
- Create kernel method for symmetric QR decomposition
- Add test case for symmetric QR factorization
- Remove unnecessary Printf and Unicode dependencies

Author: yangjunjie0320

examples/Manifest.toml

@@ -27,21 +27,12 @@ deps = ["Artifacts", "CompilerSupportLibraries_jll", "Libdl"]
 uuid = "4536629a-c528-5b80-bd46-f80d51c5b363"
 version = "0.3.27+1"
 
-[[deps.Printf]]
-deps = ["Unicode"]
-uuid = "de0858da-6303-5e67-8744-51eddeeeb8d7"
-version = "1.11.0"
-
 [[deps.SimpleLinearAlgebra]]
-deps = ["LinearAlgebra", "Printf"]
+deps = ["LinearAlgebra"]
 path = "/Users/yangjunjie/.julia/dev/SimpleLinearAlgebra"
 uuid = "555e9691-8231-4cde-a0d7-6dc20204a59a"
 version = "1.0.0-DEV"
 
-[[deps.Unicode]]
-uuid = "4ec0a83e-493e-50e2-b9ac-8f72acf5a8f5"
-version = "1.11.0"
-
 [[deps.libblastrampoline_jll]]
 deps = ["Artifacts", "Libdl"]
 uuid = "8e850b90-86db-534c-a0d3-1478176c7d93"

src/qr.jl

@@ -167,3 +167,53 @@ end
 
 kernel(prob::QRFactorizationProblem) = QR.kernel(prob)
 export QR
+
+struct SymmetricQRFactorizationProblem <: ProblemMixin
+    a::AbstractMatrix
+    tol::Real
+    function SymmetricQRFactorizationProblem(a, tol=1e-10)
+        prob = new(a, tol)
+        return prob
+    end
+end
+
+struct SymmetricQRFactorizationSolution <: SolutionMixin
+    q::AbstractMatrix
+    r::AbstractMatrix
+    function SymmetricQRFactorizationSolution(q, r)
+        r = triu(r, -2)
+        r = tril(r,  2)
+        @assert isapprox(q' * q, I) "matrix must be orthogonal"
+        soln = new(q, r)
+        return soln
+    end
+end
+
+function kernel(prob::SymmetricQRFactorizationProblem)
+    tol = prob.tol
+    n = size(prob.a, 1)
+    q = Matrix{Float64}(I, n, n)
+    r = copy(prob.a)
+    
+    for i in 1:n-1
+        # Get the column vector below and including the diagonal
+        vi = copy(r[(i+1):n, i])
+        vi[1] += sign(vi[1]) * norm(vi)
+        beta = 2.0 / dot(vi, vi)
+        hi = I - beta * vi * vi' # shape (n-i-1, n-i-1)
+
+        r[(i+1):n, i:n] = hi * r[(i+1):n, i:n]
+        r[i:n, (i+1):n] = r[i:n, (i+1):n] * hi'
+        q[:, (i+1):n] = q[:, (i+1):n] * hi'
+    end
+
+    @assert isapprox(q' * q, I, atol=tol) "q is not orthogonal"
+    @assert isapprox(q * r * q', prob.a, atol=tol) "q * r * q' is not equal to a"
+    
+    soln = SymmetricQRFactorizationSolution(q, r)
+    return soln
+end
+
+SymmetricQRFactorization = SymmetricQRFactorizationProblem
+SymmetricQRDecomposition = SymmetricQRFactorizationProblem
+export SymmetricQRFactorization, SymmetricQRDecomposition

test/qr.jl

@@ -64,3 +64,21 @@ end
     @test isapprox(q' * q, I, atol=tol)
     @test isapprox(a, q * r, atol=tol)
 end
+
+@testset "Symmetric QR factorization" begin
+    n = 4
+    a = rand(n, n)
+    a = a * a'
+
+    prob = SymmetricQRFactorization(a, 1e-6)
+    tol = prob.tol
+
+    soln = kernel(prob)
+    q = soln.q
+    r = soln.r
+    
+    @test isapprox(r, triu(r, -2), atol=tol)
+    @test isapprox(r, tril(r,  2), atol=tol)
+    @test isapprox(q' * q, I, atol=tol)
+    @test isapprox(a, q * r * q', atol=tol)
+end