From 5edb00fa754d951fc33a38380f3a7c2f24abed66 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Mon, 31 Jul 2023 23:46:44 +0200 Subject: [PATCH 1/2] Account for pivoting in sparse Cholesky decompositions --- Project.toml | 2 +- src/chol.jl | 5 +++- src/pdsparsemat.jl | 27 ++++++++++++++------ test/chol.jl | 58 ++++++++++++++++++++++++++++++++----------- test/pdmtypes.jl | 7 ++++-- test/specialarrays.jl | 18 +++++++------- 6 files changed, 82 insertions(+), 35 deletions(-) diff --git a/Project.toml b/Project.toml index edba39a..474dd9a 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "PDMats" uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" -version = "0.11.17" +version = "0.11.18" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/src/chol.jl b/src/chol.jl index f2455b5..56c5c7c 100644 --- a/src/chol.jl +++ b/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 diff --git a/src/pdsparsemat.jl b/src/pdsparsemat.jl index 95ac63a..4d88065 100644 --- a/src/pdsparsemat.jl +++ b/src/pdsparsemat.jl @@ -62,13 +62,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 @@ -97,14 +101,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 diff --git a/test/chol.jl b/test/chol.jl index d4a6875..06b43cf 100644 --- a/test/chol.jl +++ b/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 diff --git a/test/pdmtypes.jl b/test/pdmtypes.jl index 286288a..3c2a90f 100644 --- a/test/pdmtypes.jl +++ b/test/pdmtypes.jl @@ -88,8 +88,11 @@ using Test A = rand(1, 1) x = randn(1) y = x / A - @assert x / A isa Matrix{Float64} - @assert size(y) == (1, 1) + # See https://github.com/JuliaLang/julia/pull/44358 and https://github.com/JuliaLang/julia/pull/49915 + if !(v"1.9.0-DEV.17" <= VERSION < v"1.11.0-DEV.192") + @assert x / A isa Matrix{Float64} + @assert size(y) == (1, 1) + end for M in (PDiagMat(vec(A)), ScalMat(1, first(A))) @test x / M isa Matrix{Float64} diff --git a/test/specialarrays.jl b/test/specialarrays.jl index ba3cb5e..06621e8 100644 --- a/test/specialarrays.jl +++ b/test/specialarrays.jl @@ -63,17 +63,17 @@ using StaticArrays x = rand(5) X = rand(2, 5) Y = rand(5, 2) - @test P * x ≈ A * x - @test P * Y ≈ A * Y + @test P * x ≈ x + @test P * Y ≈ Y # Right division with Cholesky requires https://github.com/JuliaLang/julia/pull/32594 if VERSION >= v"1.3.0-DEV.562" - @test X / P ≈ X / A + @test X / P ≈ X end - @test P \ x ≈ A \ x - @test P \ Y ≈ A \ Y - @test X_A_Xt(P, X) ≈ X * A * X' - @test X_invA_Xt(P, X) ≈ X * (A \ X') - @test Xt_A_X(P, Y) ≈ Y' * A * Y - @test Xt_invA_X(P, Y) ≈ Y' * (A \ Y) + @test P \ x ≈ x + @test P \ Y ≈ Y + @test X_A_Xt(P, X) ≈ X * X' + @test X_invA_Xt(P, X) ≈ X * X' + @test Xt_A_X(P, Y) ≈ Y' * Y + @test Xt_invA_X(P, Y) ≈ Y' * Y end end From 045e49c6c5b77a2a5c70335df4d82660ca750c97 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Mon, 25 Sep 2023 13:41:59 +0200 Subject: [PATCH 2/2] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 474dd9a..ae30fd8 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "PDMats" uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" -version = "0.11.18" +version = "0.11.20" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"