Account for pivoting in sparse Cholesky decompositions (#175)

* Account for pivoting in sparse Cholesky decompositions

* Update Project.toml

Project.toml

@@ -1,6 +1,6 @@
 name = "PDMats"
 uuid = "90014a1f-27ba-587c-ab20-58faa44d9150"
-version = "0.11.19"
+version = "0.11.20"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

src/chol.jl

@@ -15,9 +15,12 @@ chol_lower(a::Matrix) = cholesky(Symmetric(a, :L)).L
 # NOTE: Formally, the line above should use Hermitian() instead of Symmetric(),
 # but this currently has an AutoDiff issue in Zygote.jl, and PDMat is
 # type-restricted to be Real, so they are equivalent.
+chol_upper(a::Matrix) = cholesky(Symmetric(a, :U)).U
 
 if HAVE_CHOLMOD
     CholTypeSparse{T} = SuiteSparse.CHOLMOD.Factor{T}
 
-    chol_lower(cf::CholTypeSparse) = cf.L
+    # Take into account pivoting!
+    chol_lower(cf::CholTypeSparse) = cf.PtL
+    chol_upper(cf::CholTypeSparse) = cf.UP
 end

src/pdsparsemat.jl

@@ -70,13 +70,17 @@ LinearAlgebra.sqrt(A::PDSparseMat) = PDMat(sqrt(Hermitian(Matrix(A))))
 ### whiten and unwhiten
 
 function whiten!(r::AbstractVecOrMat, a::PDSparseMat, x::AbstractVecOrMat)
-    r[:] = sparse(chol_lower(a.chol)) \ x
-    return r
+    # Can't use `ldiv!` due to missing support in SparseArrays
+    return copyto!(r, chol_lower(a.chol) \ x)
 end
 
 function unwhiten!(r::AbstractVecOrMat, a::PDSparseMat, x::AbstractVecOrMat)
-    r[:] = sparse(chol_lower(a.chol)) * x
-    return r
+    # `*` is not defined for `PtL` factor components,
+    # so we can't use `chol_lower(a.chol) * x`
+    C = a.chol
+    PtL = sparse(C.L)[C.p, :]
+    # Can't use `lmul!` due to missing support in SparseArrays
+    return copyto!(r, PtL * x)
 end
 
 
@@ -105,14 +109,23 @@ end
 ### tri products
 
 function X_A_Xt(a::PDSparseMat, x::AbstractMatrix)
-    z = x * sparse(chol_lower(a.chol))
+    # `*` is not defined for `PtL` factor components,
+    # so we can't use `x * chol_lower(a.chol)`
+    C = a.chol
+    PtL = sparse(C.L)[C.p, :]
+    z = x * PtL
     z * transpose(z)
 end
 
 
 function Xt_A_X(a::PDSparseMat, x::AbstractMatrix)
-    z = transpose(x) * sparse(chol_lower(a.chol))
-    z * transpose(z)
+    # `*` is not defined for `UP` factor components,
+    # so we can't use `chol_upper(a.chol) * x`
+    # Moreover, `sparse` is only defined for `L` factor components
+    C = a.chol
+    UP = transpose(sparse(C.L))[:, C.p]
+    z = UP * x
+    transpose(z) * z
 end
 
 

test/chol.jl

@@ -1,20 +1,48 @@
 using LinearAlgebra, PDMats
 using PDMats: chol_lower, chol_upper
 
-@testset "chol_lower" begin
-    A = rand(100, 100)
-    C = A'A
-    size_of_one_copy = sizeof(C)
-    @assert size_of_one_copy > 100  # ensure the matrix is large enough that few-byte allocations don't matter
-
-    chol_lower(C)
-    @test (@allocated chol_lower(C)) < 1.05 * size_of_one_copy  # allow 5% overhead
-
-    for uplo in (:L, :U)
-        ch = cholesky(Symmetric(C, uplo))
-        chol_lower(ch)
-        @test (@allocated chol_lower(ch)) < 33  # allow small overhead for wrapper types
-        chol_upper(ch)
-        @test (@allocated chol_upper(ch)) < 33  # allow small overhead for wrapper types
+@testset "chol_lower and chol_upper" begin
+    @testset "allocations" begin
+        A = rand(100, 100)
+        C = A'A
+        size_of_one_copy = sizeof(C)
+        @assert size_of_one_copy > 100  # ensure the matrix is large enough that few-byte allocations don't matter
+
+        @test chol_lower(C) ≈ chol_upper(C)'
+        @test (@allocated chol_lower(C)) < 1.05 * size_of_one_copy  # allow 5% overhead
+        @test (@allocated chol_upper(C)) < 1.05 * size_of_one_copy
+
+        for uplo in (:L, :U)
+            ch = cholesky(Symmetric(C, uplo))
+            @test chol_lower(ch) ≈ chol_upper(ch)'
+            @test (@allocated chol_lower(ch)) < 33  # allow small overhead for wrapper types
+            @test (@allocated chol_upper(ch)) < 33  # allow small overhead for wrapper types
+        end
+    end
+
+    # issue #120
+    @testset "correctness with pivoting" begin
+        A = [2 1 1; 1 2 0; 1 0 2]
+        x = randn(3)
+
+        # Compute `invquad` without explicit factorization
+        b = x' * (A \ x)
+
+        @test sum(abs2, PDMats.chol_lower(A) \ x) ≈ b
+        @test sum(abs2, PDMats.chol_upper(A)' \ x) ≈ b
+
+        for uplo in (:L, :U)
+            # dense version
+            ch_dense = cholesky(Symmetric(A, uplo))
+            @test sum(abs2, PDMats.chol_lower(ch_dense) \ x) ≈ b
+            @test sum(abs2, PDMats.chol_upper(ch_dense)' \ x) ≈ b
+
+            # sparse version
+            if PDMats.HAVE_CHOLMOD
+                ch_sparse = cholesky(Symmetric(sparse(A), uplo))
+                @test sum(abs2, PDMats.chol_lower(ch_sparse) \ x) ≈ b
+                @test sum(abs2, PDMats.chol_upper(ch_sparse)' \ x) ≈ b
+            end
+        end
     end
 end