From f187feb673e4ba8ee0ebe3efb55d53ad9901f6d4 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Wed, 22 Jun 2022 12:52:29 +0200 Subject: [PATCH 1/7] Define functions for `Cholesky` --- Project.toml | 2 +- README.md | 11 ++++++++ src/chol.jl | 64 +++++++++++++++++++++++++++++++++++++++++++ src/utils.jl | 7 +++++ test/chol.jl | 13 +++++++-- test/specialarrays.jl | 37 +++++++++++++++++++++---- test/testutils.jl | 2 +- 7 files changed, 127 insertions(+), 9 deletions(-) diff --git a/Project.toml b/Project.toml index 0a1c0df..79167bf 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "PDMats" uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" -version = "0.11.12" +version = "0.11.13" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" diff --git a/README.md b/README.md index 407e7a2..74b6d17 100644 --- a/README.md +++ b/README.md @@ -222,6 +222,17 @@ While in theory all of them can be defined, at present only the following subset PRs to implement more generic fallbacks are welcome. +### Fallbacks for `LinearAlgebra.Cholesky` + +For Cholesky decompositions of type `Cholesky` the following functions are defined as well: + + - `dim` + - `whiten`, `whiten!` + - `unwhiten`, `unwhiten!` + - `quad`, `quad!` + - `invquad`, `invquad!` + - `X_A_Xt`, `Xt_A_X`, `X_invA_Xt`, `Xt_invA_X` + ## Define Customized Subtypes In some situation, it is useful to define a customized subtype of `AbstractPDMat` to capture positive definite matrices with special structures. For this purpose, one has to define a subset of methods (as listed below), and other methods will be automatically provided. diff --git a/src/chol.jl b/src/chol.jl index f2455b5..d7ffab8 100644 --- a/src/chol.jl +++ b/src/chol.jl @@ -21,3 +21,67 @@ if HAVE_CHOLMOD chol_lower(cf::CholTypeSparse) = cf.L end + +# Interface for `Cholesky` + +dim(A::Cholesky) = LinearAlgebra.checksquare(A) + +# whiten +whiten(A::Cholesky, x::AbstractVecOrMat) = chol_lower(A) \ x +whiten!(A::Cholesky, x::AbstractVecOrMat) = ldiv!(chol_lower(A), x) + +# unwhiten +unwhiten(A::Cholesky, x::AbstractVecOrMat) = chol_lower(A) * x +unwhiten!(A::Cholesky, x::AbstractVecOrMat) = lmul!(chol_lower(A), x) + +# 3-argument whiten/unwhiten +for T in (:AbstractVector, :AbstractMatrix) + @eval begin + whiten!(r::$T, A::Cholesky, x::$T) = whiten!(A, copyto!(r, x)) + unwhiten!(r::$T, A::Cholesky, x::$T) = unwhiten!(A, copyto!(r, x)) + end +end + +# quad +quad(A::Cholesky, x::AbstractVector) = sum(abs2, chol_upper(A) * x) +function quad(A::Cholesky, X::AbstractMatrix) + Z = chol_upper(A) * X + return vec(sum(abs2, Z; dims=1)) +end +function quad!(r::AbstractArray, A::Cholesky, X::AbstractMatrix) + Z = chol_upper(A) * X + return map!(Base.Fix1(sum, abs2), r, _eachcol(Z)) +end + +# invquad +invquad(A::Cholesky, x::AbstractVector) = sum(abs2, chol_lower(A) \ x) +function invquad(A::Cholesky, X::AbstractMatrix) + Z = chol_lower(A) \ X + return vec(sum(abs2, Z; dims=1)) +end +function invquad!(r::AbstractArray, A::Cholesky, X::AbstractMatrix) + Z = chol_lower(A) * X + return map!(Base.Fix1(sum, abs2), r, _eachcol(Z)) +end + +# tri products + +function X_A_Xt(A::Cholesky, X::AbstractMatrix) + Z = X * chol_lower(A) + return Z * transpose(Z) +end + +function Xt_A_X(A::Cholesky, X::AbstractMatrix) + Z = chol_upper(A) * X + return transpose(Z) * Z +end + +function X_invA_Xt(A::Cholesky, X::AbstractMatrix) + Z = X / chol_upper(A) + return Z * transpose(Z) +end + +function Xt_invA_X(A::Cholesky, X::AbstractMatrix) + Z = chol_lower(A) \ X + return transpose(Z) * Z +end diff --git a/src/utils.jl b/src/utils.jl index 708440e..6f5d9d4 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -127,3 +127,10 @@ end else _ldiv!(Y::AbstractArray, s::Number, X::AbstractArray) = ldiv!(Y, s, X) end + +# eachcol was introduced in Julia 1.1 +@static if VERSION < v"1.1.0-DEV.792" + _eachcol(x::AbstractMatrix) = (view(x, :, i) for i in axes(x, 2)) +else + const _eachcol = eachcol +end diff --git a/test/chol.jl b/test/chol.jl index d4a6875..a617b94 100644 --- a/test/chol.jl +++ b/test/chol.jl @@ -2,19 +2,28 @@ using LinearAlgebra, PDMats using PDMats: chol_lower, chol_upper @testset "chol_lower" begin - A = rand(100, 100) + d = 100 + A = rand(d, d) C = A'A + invC = inv(C) 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 + @assert size_of_one_copy > d # 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 + X = randn(d, 10) 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 + + # Only test dim, `quad`/`invquad`, `whiten`/`unwhiten`, and tri products + @test dim(ch) == size(C, 1) + pdtest_quad(ch, C, invC, X, 0) + pdtest_triprod(ch, C, invC, X, 0) + pdtest_whiten(ch, C, 0) end end diff --git a/test/specialarrays.jl b/test/specialarrays.jl index b4b6b95..a09073b 100644 --- a/test/specialarrays.jl +++ b/test/specialarrays.jl @@ -23,12 +23,15 @@ using StaticArrays X = @SMatrix rand(10, 4) Y = @SMatrix rand(4, 10) - for A in (PDS, D) - @test A * x isa SVector{4, Float64} - @test A * x ≈ Matrix(A) * Vector(x) + for A in (PDS, D, C) + if !(A isa Cholesky) + # `*(::Cholesky, ::SArray)` is not defined + @test A * x isa SVector{4, Float64} + @test A * x ≈ Matrix(A) * Vector(x) - @test A * Y isa SMatrix{4, 10, Float64} - @test A * Y ≈ Matrix(A) * Matrix(Y) + @test A * Y isa SMatrix{4, 10, Float64} + @test A * Y ≈ Matrix(A) * Matrix(Y) + end @test X / A isa SMatrix{10, 4, Float64} @test X / A ≈ Matrix(X) / Matrix(A) @@ -50,6 +53,30 @@ using StaticArrays @test Xt_invA_X(A, Y) isa SMatrix{10, 10, Float64} @test Xt_invA_X(A, Y) ≈ Matrix(Y)' * (Matrix(A) \ Matrix(Y)) + + @test quad(A, x) isa Float64 + @test quad(A, x) ≈ quad(Matrix(A), Vector(x)) + + @test quad(A, Y) isa SVector{10, Float64} + @test quad(A, Y) ≈ quad(Matrix(A), Matrix(Y)) + + @test invquad(A, x) isa Float64 + @test invquad(A, x) ≈ invquad(Matrix(A), Vector(x)) + + @test invquad(A, Y) isa SVector{10, Float64} + @test invquad(A, Y) ≈ invquad(Matrix(A), Matrix(Y)) + + @test whiten(A, x) isa SVector{4, Float64} + @test whiten(A, x) ≈ whiten(Matrix(A), Vector(x)) + + @test whiten(A, Y) isa SMatrix{4, 10, Float64} + @test whiten(A, Y) ≈ whiten(Matrix(A), Matrix(Y)) + + @test unwhiten(A, x) isa SVector{4, Float64} + @test unwhiten(A, x) ≈ unwhiten(Matrix(A), Vector(x)) + + @test unwhiten(A, Y) isa SMatrix{4, 10, Float64} + @test unwhiten(A, Y) ≈ unwhiten(Matrix(A), Matrix(Y)) end end diff --git a/test/testutils.jl b/test/testutils.jl index 75e179e..30f6816 100644 --- a/test/testutils.jl +++ b/test/testutils.jl @@ -17,7 +17,7 @@ function test_pdmat(C, Cmat::Matrix; t_cholesky::Bool=true, # whether to test cholesky method t_scale::Bool=true, # whether to test scaling t_add::Bool=true, # whether to test pdadd - t_det::Bool=true, # whether to test det method + t_det::Bool=true, # whether to test det method t_logdet::Bool=true, # whether to test logdet method t_eig::Bool=true, # whether to test eigmax and eigmin t_mul::Bool=true, # whether to test multiplication From 5719f3a577f7ebea070aecf6e70d81d39b84e1bc Mon Sep 17 00:00:00 2001 From: David Widmann Date: Wed, 22 Jun 2022 13:12:30 +0200 Subject: [PATCH 2/7] Only test `Cholesky` --- test/specialarrays.jl | 3 +++ 1 file changed, 3 insertions(+) diff --git a/test/specialarrays.jl b/test/specialarrays.jl index a09073b..af207d1 100644 --- a/test/specialarrays.jl +++ b/test/specialarrays.jl @@ -53,7 +53,10 @@ using StaticArrays @test Xt_invA_X(A, Y) isa SMatrix{10, 10, Float64} @test Xt_invA_X(A, Y) ≈ Matrix(Y)' * (Matrix(A) \ Matrix(Y)) + end + # TODO: Fix for `PDMat` and `PDiagMat` + for A in (C,) @test quad(A, x) isa Float64 @test quad(A, x) ≈ quad(Matrix(A), Vector(x)) From 60b591c229d9d8d149c994dbd34445167892af40 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Sun, 18 Sep 2022 23:55:55 +0200 Subject: [PATCH 3/7] Fix numerical issues in tests --- test/specialarrays.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/specialarrays.jl b/test/specialarrays.jl index 46ab99a..aa18cb7 100644 --- a/test/specialarrays.jl +++ b/test/specialarrays.jl @@ -4,7 +4,7 @@ using StaticArrays @testset "Special matrix types" begin @testset "StaticArrays" begin # Full matrix - S = (x -> x * x')(@SMatrix(randn(4, 7))) + S = (x -> x * x' + I)(@SMatrix(randn(4, 7))) PDS = PDMat(S) @test PDS isa PDMat{Float64, <:SMatrix{4, 4, Float64}} @test isbits(PDS) From bd90e68ecae473f906f97a13520e6a1c36345ec4 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Sat, 7 Oct 2023 16:33:49 +0200 Subject: [PATCH 4/7] Update chol.jl --- test/chol.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/test/chol.jl b/test/chol.jl index 6f7b6e1..dc7be33 100644 --- a/test/chol.jl +++ b/test/chol.jl @@ -5,6 +5,7 @@ using PDMats: chol_lower, chol_upper @testset "allocations" begin A = rand(100, 100) C = A'A + invC = inv(C) 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 From 4fa2ad151fb2d9c940373157ec15ce3c8eb3b283 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Fri, 13 Oct 2023 13:28:48 +0200 Subject: [PATCH 5/7] Fixes --- src/chol.jl | 4 ++-- test/chol.jl | 1 + 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/src/chol.jl b/src/chol.jl index a15b528..313694d 100644 --- a/src/chol.jl +++ b/src/chol.jl @@ -53,7 +53,7 @@ function quad(A::Cholesky, X::AbstractMatrix) end function quad!(r::AbstractArray, A::Cholesky, X::AbstractMatrix) Z = chol_upper(A) * X - return map!(Base.Fix1(sum, abs2), r, _eachcol(Z)) + return map!(Base.Fix1(sum, abs2), r, eachcol(Z)) end # invquad @@ -64,7 +64,7 @@ function invquad(A::Cholesky, X::AbstractMatrix) end function invquad!(r::AbstractArray, A::Cholesky, X::AbstractMatrix) Z = chol_lower(A) * X - return map!(Base.Fix1(sum, abs2), r, _eachcol(Z)) + return map!(Base.Fix1(sum, abs2), r, eachcol(Z)) end # tri products diff --git a/test/chol.jl b/test/chol.jl index dc7be33..797087e 100644 --- a/test/chol.jl +++ b/test/chol.jl @@ -13,6 +13,7 @@ using PDMats: chol_lower, chol_upper @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 + X = randn(d, 10) for uplo in (:L, :U) ch = cholesky(Symmetric(C, uplo)) @test chol_lower(ch) ≈ chol_upper(ch)' From b4c0767f6720378c2c467c06d2bc97654cc926ac Mon Sep 17 00:00:00 2001 From: David Widmann Date: Fri, 13 Oct 2023 13:51:33 +0200 Subject: [PATCH 6/7] Fix tests (hopefully) --- test/chol.jl | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/test/chol.jl b/test/chol.jl index 797087e..ebdabb4 100644 --- a/test/chol.jl +++ b/test/chol.jl @@ -3,11 +3,12 @@ using PDMats: chol_lower, chol_upper @testset "chol_lower and chol_upper" begin @testset "allocations" begin - A = rand(100, 100) + d = 100 + A = rand(d, d) C = A'A invC = inv(C) 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 + @assert size_of_one_copy > d # 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 From 41ecda71a64af638e58abca5715a5f2331ed3232 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Fri, 13 Oct 2023 14:10:15 +0200 Subject: [PATCH 7/7] Remove old but now unnecessary tests --- test/specialarrays.jl | 27 --------------------------- 1 file changed, 27 deletions(-) diff --git a/test/specialarrays.jl b/test/specialarrays.jl index d1107f2..cc819bb 100644 --- a/test/specialarrays.jl +++ b/test/specialarrays.jl @@ -82,33 +82,6 @@ using StaticArrays @test Xt_invA_X(A, Y) isa SMatrix{10, 10, Float64} @test Xt_invA_X(A, Y) ≈ Matrix(Y)' * (Matrix(A) \ Matrix(Y)) end - - # TODO: Fix for `PDMat` and `PDiagMat` - for A in (C,) - @test quad(A, x) isa Float64 - @test quad(A, x) ≈ quad(Matrix(A), Vector(x)) - - @test quad(A, Y) isa SVector{10, Float64} - @test quad(A, Y) ≈ quad(Matrix(A), Matrix(Y)) - - @test invquad(A, x) isa Float64 - @test invquad(A, x) ≈ invquad(Matrix(A), Vector(x)) - - @test invquad(A, Y) isa SVector{10, Float64} - @test invquad(A, Y) ≈ invquad(Matrix(A), Matrix(Y)) - - @test whiten(A, x) isa SVector{4, Float64} - @test whiten(A, x) ≈ whiten(Matrix(A), Vector(x)) - - @test whiten(A, Y) isa SMatrix{4, 10, Float64} - @test whiten(A, Y) ≈ whiten(Matrix(A), Matrix(Y)) - - @test unwhiten(A, x) isa SVector{4, Float64} - @test unwhiten(A, x) ≈ unwhiten(Matrix(A), Vector(x)) - - @test unwhiten(A, Y) isa SMatrix{4, 10, Float64} - @test unwhiten(A, Y) ≈ unwhiten(Matrix(A), Matrix(Y)) - end end @testset "BandedMatrices" begin