Add QR factorization with modified Gram-Schmidt method

- Implement ModifiedGramSchmidt QR factorization method
- Add new QR factorization variant to QR module
- Update test suite to include modified Gram-Schmidt implementation
- Rename substitution method test files for consistency
- Add OMEinsum dependency to Project.toml

Author: yangjunjie0320

Project.toml

@@ -5,10 +5,12 @@ version = "1.0.0-DEV"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
+OMEinsum = "ebe7aa44-baf0-506c-a96f-8464559b3922"
 Printf = "de0858da-6303-5e67-8744-51eddeeeb8d7"
 
 [compat]
 LinearAlgebra = "1.10.0"
+OMEinsum = "0.8.4"
 Printf = "1.11.0"
 julia = "1.10"
 

src/SimpleLinearAlgebra.jl

@@ -1,7 +1,7 @@
 module SimpleLinearAlgebra
 
     using LinearAlgebra, Printf
-    import Printf.@sprintf  # Fixed macro import
+    using OMEinsum
 
     abstract type ProblemMixin end
     abstract type SolutionMixin end
@@ -14,8 +14,8 @@ module SimpleLinearAlgebra
         message::String
     end
 
-    include("forward-substitution.jl")
-    include("back-substitution.jl")
+    include("forward.jl")
+    include("back.jl")
     include("lu.jl")
     include("qr.jl")
 

src/back-substitution.jl renamed to src/back.jl

@@ -24,14 +24,17 @@ function kernel(prob::BackSubstitutionProblem)
     x = zeros(n)
 
     for i in n:-1:1
-        @assert abs(u[i, i]) > tol "matrix is singular"
-        x[i] += b[i] / u[i, i]
+        uii = u[i, i]
+        @assert abs(uii) > tol "matrix is singular"
+        x[i] += b[i] / uii
 
         if i == n
             continue
         end
 
-        x[i] -= u[i, i+1:n]' * x[i+1:n] / u[i, i]
+        ui = u[i, i+1:n] / uii # shape (n-i, )
+        xi = x[i+1:n]          # shape (n-i, )
+        x[i] -= ui' * xi 
     end
 
     return BackSubstitutionSolution(x)

src/forward-substitution.jl renamed to src/forward.jl

@@ -16,21 +16,26 @@ end
 
 function kernel(prob::ForwardSubstitutionProblem)
     l = prob.l
+    @assert istril(l) "matrix must be lower triangular"
+
     b = copy(prob.b)
     n = length(b)
     tol = prob.tol
 
     x = zeros(n)
 
     for i in 1:n
-        @assert abs(l[i, i]) > tol "matrix is singular"
-        x[i] += b[i] / l[i, i]
+        lii = l[i, i]
+        @assert abs(lii) > tol "matrix is singular"
+        x[i] += b[i] / lii
 
         if i == 1
             continue
         end
 
-        x[i] -= l[i, 1:i-1]' * x[1:i-1] / l[i, i]
+        li = l[i, 1:i-1] / lii # shape (i-1, )
+        xi = x[1:i-1]          # shape (i-1, )
+        x[i] -= li' * xi 
     end
     return ForwardSubstitutionSolution(x)
 end

src/lu.jl

@@ -16,57 +16,72 @@ module LU
         end
     end
 
-    struct Version1 <: LUFactorizationProblemMixin
+    struct GaussianEliminationV1 <: LUFactorizationProblemMixin
         a::AbstractMatrix
         tol::Real
-        function Version1(a, tol=1e-10)
+        function GaussianEliminationV1(a, tol=1e-10)
             prob = new(a, tol)
             return prob
         end
     end
 
-    function kernel(prob::Version1)
+    function kernel(prob::GaussianEliminationV1)
         tol = prob.tol
         u = copy(prob.a)
         n = size(u, 1)
 
         l = Matrix{Float64}(I, n, n)
-
-        for k in 1:n-1
-            @assert abs(u[k, k]) > tol "Gaussian elimination failed"
-
-            m_k = Matrix{Float64}(I, n, n)
-            m_k[k+1:n, k] = -u[k+1:n, k] / u[k, k]
-
-            l = l * inv(m_k)
-            u = m_k * u
+        # l * u = a now, we insert resolution of the identity matrix
+        # l * inv(lk) * lk * u to make l * inv(lk) lower triangular,
+        # and lk * u to make it upper triangular
+        # as lk is constructed from gaussian elimination, lk * lu
+        # will eliminate the k-th column of u
+        # and as inv(lk) is also lower triangular, the product of
+        # l * inv(lk) is lower triangular
+
+        for i in 1:n-1
+            uii = u[i, i]
+            @assert abs(uii) > tol "Gaussian elimination failed"
+            
+            # lk is lower triangular matrix
+            li = Matrix{Float64}(I, n, n)
+            li[i+1:n, i] -= u[i+1:n, i] / uii
+            li_inv = inv(li)
+
+            l = l * li_inv
+            u = li * u
         end
 
         soln = LUFactorizationSolution(tril(l), triu(u))
         return soln
     end
 
-    struct Version2 <: LUFactorizationProblemMixin
+    struct GaussianEliminationV2 <: LUFactorizationProblemMixin
         a::AbstractMatrix
         tol::Real
-        function Version2(a, tol=1e-10)
+        function GaussianEliminationV2(a, tol=1e-10)
             prob = new(a, tol)
             return prob
         end
     end
 
-    function kernel(prob::Version2)
+    function kernel(prob::GaussianEliminationV2)
         tol = prob.tol
         u = copy(prob.a)
         n = size(u, 1)
 
         # initialize l to be the identity matrix
         l = Matrix{Float64}(I, n, n)
 
-        for k in 1:n-1
-            @assert abs(u[k, k]) > tol "Gaussian elimination failed"
-            l[k+1:n, k] = u[k+1:n, k] / u[k, k]
-            u[k+1:n, k+1:n] -= l[k+1:n, k] * u[k, k+1:n]'
+        for i in 1:n-1
+            uii = u[i, i]
+            @assert abs(uii) > tol "Gaussian elimination failed"
+
+            ci = u[i+1:n, i] / uii # shape (n-i, ), the i-th column of U
+            ri = u[i, i+1:n] # shape (1, n-i), the i-th row of U
+            
+            l[i+1:n, i] = ci
+            u[i+1:n, i+1:n] -= ci * ri'
         end
 
         soln = LUFactorizationSolution(tril(l), triu(u))
@@ -99,29 +114,34 @@ module LU
         tol = prob.tol
         u = copy(prob.a)
         n = size(u, 1)
-    
+        
+        # initialize l to be the identity matrix
         l = zeros(n, n)
         p = collect(1:n)
 
-        for k in 1:n-1
-            v, x = findmax(i->abs(u[i, k]), k:n)
-            x += k - 1
+        for i in 1:n-1
+            # v is the max element, x is the index
+            v, j = findmax(u[i:n, i])
+            j += i - 1
 
-            if x != k
-                u[k, :], u[x, :] = u[x, :], u[k, :]
-                l[k, :], l[x, :] = l[x, :], l[k, :]
-                p[k], p[x] = p[x], p[k]
+            if i != j # if the max element is not on the diagonal
+                u[i, :], u[j, :] = u[j, :], u[i, :]
+                l[i, :], l[j, :] = l[j, :], l[i, :]
+                p[i], p[j] = p[j], p[i]
             end
-    
-            if abs(u[k, k]) < tol
+            
+            uii = u[i, i]
+            if abs(uii) < tol
                 continue
             end
-    
-            l[k, k] = 1.0
-            l[k+1:n, k] = u[k+1:n, k] / u[k, k]
-            u[k+1:n, k+1:n] -= l[k+1:n, k] * u[k, k+1:n]'
+            
+            l[i, i] = 1.0
+            ri = u[i, i+1:n]
+            ci = u[i+1:n, i] / uii
+            l[i+1:n, i] = ci
+            u[i+1:n, i+1:n] -= ci * ri'
         end
-    
+        
         l[n, n] = 1.0
         soln = PartialPivotingLUFactorizationSolution(tril(l), triu(u), p)
         return soln
@@ -130,5 +150,5 @@ end
 
 kernel(prob::LUFactorizationProblemMixin) = LU.kernel(prob)
 
-LUFactorization = LU.Version2
+LUFactorization = LU.GaussianEliminationV2
 export LUFactorization

src/qr.jl

@@ -30,18 +30,18 @@ module QR
         q = Matrix{Float64}(I, n, n)
         r = copy(prob.a)
 
-        for k in 1:n    
+        for i in 1:n    
             # Get the column vector below and including the diagonal
-            v = copy(r[k:n, k])
+            vi = copy(r[i:n, i])
 
-            h = Matrix{Float64}(I, n, n)
-            beta = -sign(v[1]) * norm(v)
-            v[1] -= beta
-            h[k:n, k:n] -= 2.0 * v * v' / dot(v, v)
+            hi = Matrix{Float64}(I, n, n)
+            beta = -sign(vi[1]) * norm(vi)
+            vi[1] -= beta
+            hi[i:n, i:n] -= 2.0 * vi * vi' / dot(vi, vi)
 
             # Apply reflection to R and update Q
-            r = h * r
-            q = q * h
+            r = hi * r
+            q = q * hi
         end
 
         soln = QRFactorizationSolution(q, triu(r))
@@ -62,22 +62,23 @@ module QR
         q = Matrix{Float64}(I, n, n)
         r = copy(prob.a)
 
-        for i in 1:n
-            for j in i+1:n
-                x, y = r[i, i], r[j, i]
-                c = x / sqrt(x^2 + y^2)
-                s = y / sqrt(x^2 + y^2)
+        ij = [(i, j) for i in 1:n for j in i+1:n]
+        for (i, j) in ij
+            x, y = r[i, i], r[j, i]
+            t = atan(y, x)
 
-                if abs(s) < prob.tol
-                    continue
-                end
+            if abs(y) < prob.tol
+                continue
+            end
 
-                h = Matrix{Float64}(I, n, n)
-                h[[i, j], [i, j]] = [c -s; s c]
+            # Manually calculate hij' * r and q * hij
+            ri =  cos(t) * r[i, :] + sin(t) * r[j, :]
+            rj = -sin(t) * r[i, :] + cos(t) * r[j, :]
+            r[i, :], r[j, :] = ri, rj
 
-                r = h' * r
-                q = q * h
-            end
+            qi =  cos(t) * q[:, i] + sin(t) * q[:, j]
+            qj = -sin(t) * q[:, i] + cos(t) * q[:, j]
+            q[:, i], q[:, j] = qi, qj
         end
         return QRFactorizationSolution(q, triu(r))
     end
@@ -105,12 +106,13 @@ module QR
         q[:, 1] /= r[1, 1]
 
         for i in 2:n
-            # q[:, i] = a[:, i]
-            
-            for j in 1:i-1
-                r[j, i] = dot(q[:, j], a[:, i])
-                q[:, i] -= q[:, j] * r[j, i]
-            end
+            # Method 1:
+            qi = q[:, 1:i-1] # shape (n, i-1)
+            ai = a[:, i] # shape (n, )
+            ri = qi' * ai # shape (i-1, )
+
+            r[1:i-1, i] = ri
+            q[:, i] -= qi * ri
 
             r[i, i] = norm(q[:, i])
             @assert r[i, i] > tol "Singular matrix"
@@ -120,6 +122,47 @@ module QR
 
         return QRFactorizationSolution(q, triu(r))
     end
+
+    struct ModifiedGramSchmidt <: QRFactorizationProblem
+        a::AbstractMatrix
+        tol::Real
+        function ModifiedGramSchmidt(a, tol=1e-10)
+            prob = new(a, tol)
+            return prob
+        end
+    end
+
+    function kernel(prob::ModifiedGramSchmidt)
+        a = copy(prob.a)
+        tol = prob.tol
+
+        n = size(a, 1)
+        q = zeros(n, n)
+        r = zeros(n, n)
+    
+        for i in 1:n
+            # vi is normalized vector of a[:, i]
+            # is orthogonal to all the vectors
+            # before i
+            ni = norm(a[:, i])
+            vi = a[:, i] / ni # shape (n, )
+            @assert size(vi) == (n, )
+
+            # ai contains all the vectors
+            # after i
+            ai = a[:, i+1:n] # shape (n, n-i)
+            # ri is the projection of vi into the 
+            # basis spanned by ai
+            ri = vi' * ai # shape (n-i, )
+            @assert size(ri) == (1, n-i)
+
+            r[i, i] = ni
+            r[i, i+1:n] = ri
+            q[:, i] = vi
+            a[:, i+1:n] -= vi * ri
+        end
+        return QRFactorizationSolution(q, triu(r))
+    end
 end
 
 kernel(prob::QRFactorizationProblem) = QR.kernel(prob)

test/back-substitution.jl renamed to test/back.jl

@@ -1,12 +1,15 @@
+using LinearAlgebra
+
 @testset "backward substitution" begin
     n = 10
-    u = LinearAlgebra.triu(rand(n, n))
+    u = triu(rand(n, n))
     b = rand(n)
 
-    p = BackSubstitution(u, b)
-    s = kernel(p)
-    x = s.x
-    tol = p.tol
+    prob = BackSubstitution(u, b)
+    tol = prob.tol
+
+    soln = kernel(prob)
+    x = soln.x
 
     @test isapprox(u * x, b, atol=tol)
 end