From 5ab95842b22380378882b7c7d025fa1c720b14bb Mon Sep 17 00:00:00 2001 From: David Widmann Date: Wed, 4 Oct 2023 09:08:11 +0200 Subject: [PATCH 01/14] Unify `PDMat` and `PDSparseMat` --- .github/workflows/ci.yml | 3 +- Project.toml | 15 +- README.md | 33 +---- .../PDMatsSparseArraysExt.jl | 15 ++ ext/PDMatsSparseArraysExt/chol.jl | 5 + ext/PDMatsSparseArraysExt/pdsparsemat.jl | 109 ++++++++++++++ src/PDMats.jl | 10 +- src/addition.jl | 38 ----- src/chol.jl | 8 - src/pdiagmat.jl | 7 +- src/pdmat.jl | 37 ++--- src/pdsparsemat.jl | 140 ------------------ src/scalmat.jl | 13 +- src/utils.jl | 40 ++--- test/addition.jl | 2 +- test/chol.jl | 2 +- test/pdmtypes.jl | 64 ++++---- test/specialarrays.jl | 5 +- test/testutils.jl | 34 +++-- 19 files changed, 231 insertions(+), 349 deletions(-) create mode 100644 ext/PDMatsSparseArraysExt/PDMatsSparseArraysExt.jl create mode 100644 ext/PDMatsSparseArraysExt/chol.jl create mode 100644 ext/PDMatsSparseArraysExt/pdsparsemat.jl delete mode 100644 src/pdsparsemat.jl diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 7abc2de..93b0cf0 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -12,8 +12,7 @@ jobs: fail-fast: false matrix: version: - - '1.0' - - '1.6' + - '1.9' - '1' - nightly os: diff --git a/Project.toml b/Project.toml index 6a2ba90..cc9242d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,19 +1,24 @@ name = "PDMats" uuid = "90014a1f-27ba-587c-ab20-58faa44d9150" -version = "0.11.22" +version = "0.12.0" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" -SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" [compat] -julia = "1" +julia = "1.9" + +[extensions] +PDMatsSparseArraysExt = "SparseArrays" [extras] BandedMatrices = "aae01518-5342-5314-be14-df237901396f" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["BandedMatrices", "StaticArrays", "Test"] +test = ["BandedMatrices", "SparseArrays", "StaticArrays", "Test"] + +[weakdeps] +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" diff --git a/README.md b/README.md index bea493d..82bc7d2 100644 --- a/README.md +++ b/README.md @@ -27,17 +27,17 @@ Elemenent types are in princple all Real types, but in practice this is limited * `PDMat`: full covariance matrix, defined as ```julia -struct PDMat{T<:Real,S<:AbstractMatrix} <: AbstractPDMat{T} - dim::Int # matrix dimension - mat::S # input matrix - chol::Cholesky{T,S} # Cholesky factorization of mat +struct PDMat{T<:Real,S<:AbstractMatrix,C<:Factorization} <: AbstractPDMat{T} + dim::Int # matrix dimension + mat::S # input matrix + chol::C # Cholesky factorization of mat end # Constructors PDMat(mat, chol) # with both the input matrix and a Cholesky factorization -PDMat(mat) # with the input matrix, of type Matrix or Symmetric +PDMat(mat) # with the input matrix # Remarks: the Cholesky factorization will be computed # upon construction. @@ -76,29 +76,6 @@ ScalMat(d, v) # with dimension d and diagonal value v ``` -* `PDSparseMat`: sparse covariance matrix, defined as - -```julia -struct PDSparseMat{T<:Real,S<:AbstractSparseMatrix} <: AbstractPDMat{T} - dim::Int # matrix dimension - mat::SparseMatrixCSC # input matrix - chol::CholTypeSparse # Cholesky factorization of mat -end - -# Constructors - -PDSparseMat(mat, chol) # with both the input matrix and a Cholesky factorization - -PDSparseMat(mat) # with the sparse input matrix, of type SparseMatrixCSC - # Remarks: the Cholesky factorization will be computed - # upon construction. - -PDSparseMat(chol) # with the Cholesky factorization - # Remarks: the sparse matrix 'mat' will be computed upon - # construction. -``` - - ## Common interface All subtypes of `AbstractPDMat` share the same API, *i.e.* with the same set of methods to operate on their instances. These methods are introduced below, where `a` is an instance of a subtype of `AbstractPDMat` to represent a positive definite matrix, `x` be a column vector or a matrix with `size(x,1) == size(a, 1) == size(a, 2)`, and `c` be a positive real value. diff --git a/ext/PDMatsSparseArraysExt/PDMatsSparseArraysExt.jl b/ext/PDMatsSparseArraysExt/PDMatsSparseArraysExt.jl new file mode 100644 index 0000000..7f684e0 --- /dev/null +++ b/ext/PDMatsSparseArraysExt/PDMatsSparseArraysExt.jl @@ -0,0 +1,15 @@ +module PDMatsSparseArraysExt + +using PDMats +using SparseArrays + +using PDMats.LinearAlgebra + +const HAVE_CHOLMOD = isdefined(SparseArrays, :CHOLMOD) + +if HAVE_CHOLMOD + include("chol.jl") + include("pdsparsemat.jl") +end + +end # module diff --git a/ext/PDMatsSparseArraysExt/chol.jl b/ext/PDMatsSparseArraysExt/chol.jl new file mode 100644 index 0000000..0c18ec7 --- /dev/null +++ b/ext/PDMatsSparseArraysExt/chol.jl @@ -0,0 +1,5 @@ +const CholTypeSparse{T} = SparseArrays.CHOLMOD.Factor{T} + +# Take into account pivoting! +PDMats.chol_lower(cf::CholTypeSparse) = cf.PtL +PDMats.chol_upper(cf::CholTypeSparse) = cf.UP diff --git a/ext/PDMatsSparseArraysExt/pdsparsemat.jl b/ext/PDMatsSparseArraysExt/pdsparsemat.jl new file mode 100644 index 0000000..50d43ec --- /dev/null +++ b/ext/PDMatsSparseArraysExt/pdsparsemat.jl @@ -0,0 +1,109 @@ +""" +Sparse positive definite matrix together with a Cholesky factorization object. +""" +const PDSparseMat{T<:Real,S<:AbstractSparseMatrix,C<:CholTypeSparse} = PDMat{T,S,C} + +function PDMats.PDMat(mat::AbstractSparseMatrix, chol::CholTypeSparse) + d = size(mat, 1) + size(chol, 1) == d || + throw(DimensionMismatch("Dimensions of mat and chol are inconsistent.")) + PDMat{eltype(mat),typeof(mat),typeof(chol)}(d, mat, chol) +end + +PDMats.PDMat(mat::SparseMatrixCSC) = PDMat(mat, cholesky(mat)) +PDMats.PDMat(fac::CholTypeSparse) = PDMat(sparse(fac), fac) + +PDMats.AbstractPDMat(A::CholTypeSparse) = PDMat(A) + +### Conversion +function Base.convert(::Type{PDMat{T}}, a::PDSparseMat) where {T<:Real} + # CholTypeSparse only supports Float64 and ComplexF64 type parameters! + # So there is no point in recomputing `cholesky(mat)` and we just reuse + # the existing Cholesky factorization + mat = convert(AbstractMatrix{T}, a.mat) + return PDMat{T,typeof(mat),typeof(a.chol)}(a.dim, mat, a.chol) +end + +### Arithmetics + +Base.:\(a::PDSparseMat{T}, x::AbstractVecOrMat{T}) where {T<:Real} = convert(Array{T},a.chol \ convert(Array{Float64},x)) #to avoid limitations in sparse factorization library CHOLMOD, see e.g., julia issue #14076 +Base.:/(x::AbstractVecOrMat{T}, a::PDSparseMat{T}) where {T<:Real} = convert(Array{T},convert(Array{Float64},x) / a.chol ) + +### Algebra + +Base.inv(a::PDSparseMat{T}) where {T<:Real} = PDMat(inv(a.mat)) +LinearAlgebra.det(a::PDSparseMat) = det(a.chol) +LinearAlgebra.logdet(a::PDSparseMat) = logdet(a.chol) +LinearAlgebra.sqrt(A::PDSparseMat) = PDMat(sqrt(Hermitian(Matrix(A)))) + +### whiten and unwhiten + +function PDMats.whiten!(r::AbstractVecOrMat, a::PDSparseMat, x::AbstractVecOrMat) + # Can't use `ldiv!` due to missing support in SparseArrays + return copyto!(r, PDMats.chol_lower(a.chol) \ x) +end + +function PDMats.unwhiten!(r::AbstractVecOrMat, a::PDSparseMat, x::AbstractVecOrMat) + # `*` 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 + + +### quadratic forms + +PDMats.quad(a::PDSparseMat, x::AbstractVector) = dot(x, a * x) +PDMats.invquad(a::PDSparseMat, x::AbstractVector) = dot(x, a \ x) + +function PDMats.quad!(r::AbstractArray, a::PDSparseMat, x::AbstractMatrix) + PDMats.@check_argdims eachindex(r) == axes(x, 2) + for i in axes(x, 2) + r[i] = quad(a, x[:,i]) + end + return r +end + +function PDMats.invquad!(r::AbstractArray, a::PDSparseMat, x::AbstractMatrix) + PDMats.@check_argdims eachindex(r) == axes(x, 2) + for i in axes(x, 2) + r[i] = invquad(a, x[:,i]) + end + return r +end + + +### tri products + +function PDMats.X_A_Xt(a::PDSparseMat, x::AbstractMatrix) + # `*` 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 PDMats.Xt_A_X(a::PDSparseMat, x::AbstractMatrix) + # `*` 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 + + +function PDMats.X_invA_Xt(a::PDSparseMat, x::AbstractMatrix) + z = a.chol \ collect(transpose(x)) + x * z +end + +function PDMats.Xt_invA_X(a::PDSparseMat, x::AbstractMatrix) + z = a.chol \ x + transpose(x) * z +end diff --git a/src/PDMats.jl b/src/PDMats.jl index 4cf988d..56f49ad 100644 --- a/src/PDMats.jl +++ b/src/PDMats.jl @@ -1,6 +1,6 @@ module PDMats - using LinearAlgebra, SparseArrays, SuiteSparse + using LinearAlgebra import Base: +, *, \, /, ==, convert, inv, Matrix, kron @@ -8,7 +8,6 @@ module PDMats # Types AbstractPDMat, PDMat, - PDSparseMat, PDiagMat, ScalMat, @@ -35,19 +34,12 @@ module PDMats """ abstract type AbstractPDMat{T<:Real} <: AbstractMatrix{T} end - const HAVE_CHOLMOD = isdefined(SuiteSparse, :CHOLMOD) - # source files include("chol.jl") include("utils.jl") include("pdmat.jl") - - if HAVE_CHOLMOD - include("pdsparsemat.jl") - end - include("pdiagmat.jl") include("scalmat.jl") diff --git a/src/addition.jl b/src/addition.jl index ab4cbc0..05f4b84 100644 --- a/src/addition.jl +++ b/src/addition.jl @@ -4,38 +4,18 @@ +(a::PDMat, b::AbstractPDMat) = PDMat(a.mat + Matrix(b)) +(a::PDiagMat, b::AbstractPDMat) = PDMat(_adddiag!(Matrix(b), a.diag, true)) +(a::ScalMat, b::AbstractPDMat) = PDMat(_adddiag!(Matrix(b), a.value)) -if HAVE_CHOLMOD - +(a::PDSparseMat, b::AbstractPDMat) = PDMat(a.mat + Matrix(b)) -end +(a::PDMat, b::PDMat) = PDMat(a.mat + b.mat) +(a::PDMat, b::PDiagMat) = PDMat(_adddiag(a.mat, b.diag)) +(a::PDMat, b::ScalMat) = PDMat(_adddiag(a.mat, b.value)) -if HAVE_CHOLMOD - +(a::PDMat, b::PDSparseMat) = PDMat(a.mat + b.mat) -end +(a::PDiagMat, b::PDMat) = PDMat(_adddiag(b.mat, a.diag)) +(a::PDiagMat, b::PDiagMat) = PDiagMat(a.diag + b.diag) +(a::PDiagMat, b::ScalMat) = PDiagMat(a.diag .+ b.value) -if HAVE_CHOLMOD - +(a::PDiagMat, b::PDSparseMat) = PDSparseMat(_adddiag(b.mat, a.diag)) -end +(a::ScalMat, b::PDMat) = PDMat(_adddiag(b.mat, a.value)) +(a::ScalMat, b::PDiagMat) = PDiagMat(a.value .+ b.diag) +(a::ScalMat, b::ScalMat) = ScalMat(a.dim, a.value + b.value) -if HAVE_CHOLMOD - +(a::ScalMat, b::PDSparseMat) = PDSparseMat(_adddiag(b.mat, a.value)) -end - -if HAVE_CHOLMOD - +(a::PDSparseMat, b::PDMat) = PDMat(Matrix(a) + b.mat) - +(a::PDSparseMat, b::PDiagMat) = PDSparseMat(_adddiag(a.mat, b.diag)) - +(a::PDSparseMat, b::ScalMat) = PDSparseMat(_adddiag(a.mat, b.value)) - +(a::PDSparseMat, b::PDSparseMat) = PDSparseMat(a.mat + b.mat) -end - # between pdmat and uniformscaling (multiple of identity) @@ -45,34 +25,16 @@ end pdadd(a::PDMat, b::AbstractPDMat, c::Real) = PDMat(a.mat + Matrix(b * c)) pdadd(a::PDiagMat, b::AbstractPDMat, c::Real) = PDMat(_adddiag!(Matrix(b * c), a.diag, one(c))) pdadd(a::ScalMat, b::AbstractPDMat, c::Real) = PDMat(_adddiag!(Matrix(b * c), a.value)) -if HAVE_CHOLMOD - pdadd(a::PDSparseMat, b::AbstractPDMat, c::Real) = PDMat(a.mat + Matrix(b * c)) -end pdadd(a::PDMat, b::PDMat, c::Real) = PDMat(a.mat + b.mat * c) pdadd(a::PDMat, b::PDiagMat, c::Real) = PDMat(_adddiag(a.mat, b.diag, c)) pdadd(a::PDMat, b::ScalMat, c::Real) = PDMat(_adddiag(a.mat, b.value * c)) -if HAVE_CHOLMOD - pdadd(a::PDMat, b::PDSparseMat, c::Real) = PDMat(a.mat + b.mat * c) -end pdadd(a::PDiagMat, b::PDMat, c::Real) = PDMat(_adddiag!(b.mat * c, a.diag, one(c))) pdadd(a::PDiagMat, b::PDiagMat, c::Real) = PDiagMat(a.diag + b.diag * c) pdadd(a::PDiagMat, b::ScalMat, c::Real) = PDiagMat(a.diag .+ b.value * c) -if HAVE_CHOLMOD - pdadd(a::PDiagMat, b::PDSparseMat, c::Real) = PDSparseMat(_adddiag!(b.mat * c, a.diag, one(c))) -end pdadd(a::ScalMat, b::PDMat, c::Real) = PDMat(_adddiag!(b.mat * c, a.value)) pdadd(a::ScalMat, b::PDiagMat, c::Real) = PDiagMat(a.value .+ b.diag * c) pdadd(a::ScalMat, b::ScalMat, c::Real) = ScalMat(a.dim, a.value + b.value * c) -if HAVE_CHOLMOD - pdadd(a::ScalMat, b::PDSparseMat, c::Real) = PDSparseMat(_adddiag!(b.mat * c, a.value)) -end -if HAVE_CHOLMOD - pdadd(a::PDSparseMat, b::PDMat, c::Real) = PDMat(a.mat + b.mat * c) - pdadd(a::PDSparseMat, b::PDiagMat, c::Real) = PDSparseMat(_adddiag(a.mat, b.diag, c)) - pdadd(a::PDSparseMat, b::ScalMat, c::Real) = PDSparseMat(_adddiag(a.mat, b.value * c)) - pdadd(a::PDSparseMat, b::PDSparseMat, c::Real) = PDSparseMat(a.mat + b.mat * c) -end diff --git a/src/chol.jl b/src/chol.jl index 56c5c7c..125a217 100644 --- a/src/chol.jl +++ b/src/chol.jl @@ -16,11 +16,3 @@ chol_lower(a::Matrix) = cholesky(Symmetric(a, :L)).L # 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} - - # Take into account pivoting! - chol_lower(cf::CholTypeSparse) = cf.PtL - chol_upper(cf::CholTypeSparse) = cf.UP -end diff --git a/src/pdiagmat.jl b/src/pdiagmat.jl index df6f74f..948a8d5 100644 --- a/src/pdiagmat.jl +++ b/src/pdiagmat.jl @@ -62,12 +62,7 @@ function \(a::PDiagMat, x::AbstractVecOrMat) end function /(x::AbstractVecOrMat, a::PDiagMat) @check_argdims a.dim == size(x, 2) - if VERSION < v"1.9-" - # return matrix for 1-element vectors `x`, consistent with LinearAlgebra < 1.9 - return reshape(x, Val(2)) ./ permutedims(a.diag) # = (a' \ x')' - else - return x ./ (x isa AbstractVector ? a.diag : a.diag') - end + return x ./ (x isa AbstractVector ? a.diag : a.diag') end Base.kron(A::PDiagMat, B::PDiagMat) = PDiagMat(vec(permutedims(A.diag) .* B.diag)) diff --git a/src/pdmat.jl b/src/pdmat.jl index b8fd077..f6f1972 100644 --- a/src/pdmat.jl +++ b/src/pdmat.jl @@ -1,19 +1,19 @@ """ Full positive definite matrix together with a Cholesky factorization object. """ -struct PDMat{T<:Real,S<:AbstractMatrix} <: AbstractPDMat{T} +struct PDMat{T<:Real,S<:AbstractMatrix,C<:Factorization} <: AbstractPDMat{T} dim::Int mat::S - chol::Cholesky{T,S} + chol::C - PDMat{T,S}(d::Int,m::AbstractMatrix{T},c::Cholesky{T,S}) where {T,S} = new{T,S}(d,m,c) + PDMat{T,S,C}(d::Int, m::AbstractMatrix{T}, c::Factorization) where {T,S,C} = new{T,S,C}(d,m,c) end function PDMat(mat::AbstractMatrix,chol::Cholesky{T,S}) where {T,S} d = size(mat, 1) size(chol, 1) == d || throw(DimensionMismatch("Dimensions of mat and chol are inconsistent.")) - PDMat{T,S}(d, convert(S, mat), chol) + PDMat{T,S,Cholesky{T,S}}(d, convert(S, mat), chol) end PDMat(mat::AbstractMatrix) = PDMat(mat, cholesky(mat)) @@ -24,17 +24,16 @@ AbstractPDMat(A::Cholesky) = PDMat(A) ### Conversion Base.convert(::Type{PDMat{T}}, a::PDMat{T}) where {T<:Real} = a function Base.convert(::Type{PDMat{T}}, a::PDMat) where {T<:Real} - chol = convert(Cholesky{T}, a.chol) - S = typeof(chol.factors) - mat = convert(S, a.mat) - return PDMat{T,S}(size(mat, 1), mat, chol) + chol = convert(Factorization{float(T)}, a.chol) + mat = convert(AbstractMatrix{T}, a.mat) + return PDMat{T,typeof(mat),typeof(chol)}(size(mat, 1), mat, chol) end Base.convert(::Type{AbstractPDMat{T}}, a::PDMat) where {T<:Real} = convert(PDMat{T}, a) ### Basics Base.size(a::PDMat) = (a.dim, a.dim) -Base.Matrix(a::PDMat) = copy(a.mat) +Base.Matrix(a::PDMat) = Matrix(a.mat) LinearAlgebra.diag(a::PDMat) = diag(a.mat) LinearAlgebra.cholesky(a::PDMat) = a.chol @@ -55,17 +54,7 @@ end *(a::PDMat, x::AbstractMatrix) = a.mat * x \(a::PDMat, x::AbstractVecOrMat) = a.chol \ x function /(x::AbstractVecOrMat, a::PDMat) - # /(::AbstractVector, ::Cholesky) is not defined - if VERSION < v"1.9-" - # return matrix for 1-element vectors `x`, consistent with LinearAlgebra - return reshape(x, Val(2)) / a.chol - else - if x isa AbstractVector - return vec(reshape(x, Val(2)) / a.chol) - else - return x / a.chol - end - end + return x isa AbstractVector ? vec(reshape(x, Val(2)) / a.chol) : x / a.chol end ### Algebra @@ -75,9 +64,15 @@ LinearAlgebra.det(a::PDMat) = det(a.chol) LinearAlgebra.logdet(a::PDMat) = logdet(a.chol) LinearAlgebra.eigmax(a::PDMat) = eigmax(a.mat) LinearAlgebra.eigmin(a::PDMat) = eigmin(a.mat) -Base.kron(A::PDMat, B::PDMat) = PDMat(kron(A.mat, B.mat), Cholesky(kron(A.chol.U, B.chol.U), 'U', A.chol.info)) LinearAlgebra.sqrt(A::PDMat) = PDMat(sqrt(Hermitian(A.mat))) +function Base.kron( + A::PDMat{<:Real,<:AbstractMatrix,<:Cholesky}, + B::PDMat{<:Real,<:AbstractMatrix,<:Cholesky}, +) + return PDMat(kron(A.mat, B.mat), Cholesky(kron(A.chol.U, B.chol.U), 'U', A.chol.info)) +end + ### tri products function X_A_Xt(a::PDMat, x::AbstractMatrix) diff --git a/src/pdsparsemat.jl b/src/pdsparsemat.jl deleted file mode 100644 index 2899b59..0000000 --- a/src/pdsparsemat.jl +++ /dev/null @@ -1,140 +0,0 @@ -""" -Sparse positive definite matrix together with a Cholesky factorization object. -""" -struct PDSparseMat{T<:Real,S<:AbstractSparseMatrix} <: AbstractPDMat{T} - dim::Int - mat::S - chol::CholTypeSparse - - PDSparseMat{T,S}(d::Int,m::AbstractSparseMatrix{T},c::CholTypeSparse) where {T,S} = - new{T,S}(d,m,c) #add {T} to CholTypeSparse argument once #14076 is implemented -end - -function PDSparseMat(mat::AbstractSparseMatrix,chol::CholTypeSparse) - d = size(mat, 1) - size(chol, 1) == d || - throw(DimensionMismatch("Dimensions of mat and chol are inconsistent.")) - PDSparseMat{eltype(mat),typeof(mat)}(d, mat, chol) -end - -PDSparseMat(mat::SparseMatrixCSC) = PDSparseMat(mat, cholesky(mat)) -PDSparseMat(fac::CholTypeSparse) = PDSparseMat(sparse(fac), fac) - -AbstractPDMat(A::SparseMatrixCSC) = PDSparseMat(A) -AbstractPDMat(A::CholTypeSparse) = PDSparseMat(A) - -### Conversion -Base.convert(::Type{PDSparseMat{T}}, a::PDSparseMat{T}) where {T<:Real} = a -function Base.convert(::Type{PDSparseMat{T}}, a::PDSparseMat) where {T<:Real} - # CholTypeSparse only supports Float64 and ComplexF64 type parameters! - # So there is no point in recomputing `cholesky(mat)` and we just reuse - # the existing Cholesky factorization - mat = convert(AbstractMatrix{T}, a.mat) - return PDSparseMat{T,typeof(mat)}(a.dim, mat, a.chol) -end -Base.convert(::Type{AbstractPDMat{T}}, a::PDSparseMat) where {T<:Real} = convert(PDSparseMat{T}, a) - -### Basics - -Base.size(a::PDSparseMat) = (a.dim, a.dim) -Base.Matrix(a::PDSparseMat) = Matrix(a.mat) -LinearAlgebra.diag(a::PDSparseMat) = diag(a.mat) -LinearAlgebra.cholesky(a::PDSparseMat) = a.chol - -### Inheriting from AbstractMatrix - -Base.getindex(a::PDSparseMat,i::Integer) = getindex(a.mat, i) -Base.getindex(a::PDSparseMat,I::Vararg{Int, N}) where {N} = getindex(a.mat, I...) - -### Arithmetics - -# add `a * c` to a dense matrix `m` of the same size inplace. -function pdadd!(r::Matrix, a::Matrix, b::PDSparseMat, c) - @check_argdims size(r) == size(a) == size(b) - _addscal!(r, a, b.mat, c) -end - -*(a::PDSparseMat, c::Real) = PDSparseMat(a.mat * c) -*(a::PDSparseMat, x::AbstractMatrix) = a.mat * x # defining these seperately to avoid -*(a::PDSparseMat, x::AbstractVector) = a.mat * x # ambiguity errors -\(a::PDSparseMat{T}, x::AbstractVecOrMat{T}) where {T<:Real} = convert(Array{T},a.chol \ convert(Array{Float64},x)) #to avoid limitations in sparse factorization library CHOLMOD, see e.g., julia issue #14076 -/(x::AbstractVecOrMat{T}, a::PDSparseMat{T}) where {T<:Real} = convert(Array{T},convert(Array{Float64},x) / a.chol ) - -### Algebra - -Base.inv(a::PDSparseMat{T}) where {T<:Real} = PDMat(inv(a.mat)) -LinearAlgebra.det(a::PDSparseMat) = det(a.chol) -LinearAlgebra.logdet(a::PDSparseMat) = logdet(a.chol) -LinearAlgebra.sqrt(A::PDSparseMat) = PDMat(sqrt(Hermitian(Matrix(A)))) - -### whiten and unwhiten - -function whiten!(r::AbstractVecOrMat, a::PDSparseMat, x::AbstractVecOrMat) - # 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) - # `*` 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 - - -### quadratic forms - -quad(a::PDSparseMat, x::AbstractVector) = dot(x, a * x) -invquad(a::PDSparseMat, x::AbstractVector) = dot(x, a \ x) - -function quad!(r::AbstractArray, a::PDSparseMat, x::AbstractMatrix) - @check_argdims eachindex(r) == axes(x, 2) - for i in axes(x, 2) - r[i] = quad(a, x[:,i]) - end - return r -end - -function invquad!(r::AbstractArray, a::PDSparseMat, x::AbstractMatrix) - @check_argdims eachindex(r) == axes(x, 2) - for i in axes(x, 2) - r[i] = invquad(a, x[:,i]) - end - return r -end - - -### tri products - -function X_A_Xt(a::PDSparseMat, x::AbstractMatrix) - # `*` 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) - # `*` 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 - - -function X_invA_Xt(a::PDSparseMat, x::AbstractMatrix) - z = a.chol \ collect(transpose(x)) - x * z -end - -function Xt_invA_X(a::PDSparseMat, x::AbstractMatrix) - z = a.chol \ x - transpose(x) * z -end diff --git a/src/scalmat.jl b/src/scalmat.jl index ea672ec..878c166 100644 --- a/src/scalmat.jl +++ b/src/scalmat.jl @@ -54,12 +54,7 @@ function \(a::ScalMat, x::AbstractVecOrMat) end function /(x::AbstractVecOrMat, a::ScalMat) @check_argdims a.dim == size(x, 2) - if VERSION < v"1.9-" - # return matrix for 1-element vectors `x`, consistent with LinearAlgebra < 1.9 - return reshape(x, Val(2)) / a.value - else - return x / a.value - end + return x / a.value end Base.kron(A::ScalMat, B::ScalMat) = ScalMat(A.dim * B.dim, A.value * B.value ) @@ -77,7 +72,7 @@ LinearAlgebra.sqrt(a::ScalMat) = ScalMat(a.dim, sqrt(a.value)) function whiten!(r::AbstractVecOrMat, a::ScalMat, x::AbstractVecOrMat) @check_argdims LinearAlgebra.checksquare(a) == size(x, 1) - _ldiv!(r, sqrt(a.value), x) + ldiv!(r, sqrt(a.value), x) end function unwhiten!(r::AbstractVecOrMat, a::ScalMat, x::AbstractVecOrMat) @@ -130,10 +125,10 @@ end function X_invA_Xt(a::ScalMat, x::Matrix) @check_argdims LinearAlgebra.checksquare(a) == size(x, 2) - _rdiv!(x * transpose(x), a.value) + rdiv!(x * transpose(x), a.value) end function Xt_invA_X(a::ScalMat, x::Matrix) @check_argdims LinearAlgebra.checksquare(a) == size(x, 1) - _rdiv!(transpose(x) * x, a.value) + rdiv!(transpose(x) * x, a.value) end diff --git a/src/utils.jl b/src/utils.jl index c85d938..b34fc50 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -10,27 +10,27 @@ end _rcopy!(r, x) = (r === x || copyto!(r, x); r) -function _addscal!(r::Matrix, a::Matrix, b::Union{Matrix, SparseMatrixCSC}, c::Real) - if c == one(c) - for i in eachindex(a) +function _addscal!(r::AbstractMatrix, a::AbstractMatrix, b::AbstractMatrix, c::Real) + if isone(c) + for i in eachindex(r, a, b) @inbounds r[i] = a[i] + b[i] end else - for i in eachindex(a) - @inbounds r[i] = a[i] + b[i] * c + for i in eachindex(r, a, b) + @inbounds r[i] = muladd(b[i], c, a[i]) end end return r end -function _adddiag!(a::Union{Matrix, SparseMatrixCSC}, v::Real) +function _adddiag!(a::AbstractMatrix, v::Real) for i in diagind(a) @inbounds a[i] += v end return a end -function _adddiag!(a::Union{Matrix, SparseMatrixCSC}, v::AbstractVector, c::Real) +function _adddiag!(a::AbstractMatrix, v::AbstractVector, c::Real) @check_argdims eachindex(v) == axes(a, 1) == axes(a, 2) if c == one(c) for i in eachindex(v) @@ -44,9 +44,9 @@ function _adddiag!(a::Union{Matrix, SparseMatrixCSC}, v::AbstractVector, c::Real return a end -_adddiag(a::Union{Matrix, SparseMatrixCSC}, v::Real) = _adddiag!(copy(a), v) -_adddiag(a::Union{Matrix, SparseMatrixCSC}, v::AbstractVector, c::Real) = _adddiag!(copy(a), v, c) -_adddiag(a::Union{Matrix, SparseMatrixCSC}, v::AbstractVector{T}) where {T<:Real} = _adddiag!(copy(a), v, one(T)) +_adddiag(a::AbstractMatrix, v::Real) = _adddiag!(copy(a), v) +_adddiag(a::AbstractMatrix, v::AbstractVector, c::Real) = _adddiag!(copy(a), v, c) +_adddiag(a::AbstractMatrix, v::AbstractVector{T}) where {T<:Real} = _adddiag!(copy(a), v, one(T)) function wsumsq(w::AbstractVector, a::AbstractVector) @@ -104,23 +104,3 @@ function colwise_sumsqinv!(r::AbstractArray, a::AbstractMatrix, c::Real) return r end -# `rdiv!(::AbstractArray, ::Number)` was introduced in Julia 1.2 -# https://github.com/JuliaLang/julia/pull/31179 -@static if VERSION < v"1.2.0-DEV.385" - function _rdiv!(X::AbstractArray, s::Number) - @simd for I in eachindex(X) - @inbounds X[I] /= s - end - X - end -else - _rdiv!(X::AbstractArray, s::Number) = rdiv!(X, s) -end - -# `ldiv!(::AbstractArray, ::Number, ::AbstractArray)` was introduced in Julia 1.4 -# https://github.com/JuliaLang/julia/pull/33806 -@static if VERSION < v"1.4.0-DEV.635" - _ldiv!(Y::AbstractArray, s::Number, X::AbstractArray) = Y .= s .\ X -else - _ldiv!(Y::AbstractArray, s::Number, X::AbstractArray) = ldiv!(Y, s, X) -end diff --git a/test/addition.jl b/test/addition.jl index 32b48e1..6902a5e 100644 --- a/test/addition.jl +++ b/test/addition.jl @@ -25,7 +25,7 @@ Base.getindex(a::ScalMat3D, i::Int, j::Int) = i == j ? a.value : zero(a.value) pm2 = PDiagMat(V) pm3 = PDiagMat(sparse(V)) pm4 = ScalMat(3, X) - pm5 = PDSparseMat(sparse(M)) + pm5 = PDMat(sparse(M)) pm6 = ScalMat3D(X) pmats = Any[pm1, pm2, pm3, pm4, pm5, pm6] diff --git a/test/chol.jl b/test/chol.jl index 06b43cf..b5a8abe 100644 --- a/test/chol.jl +++ b/test/chol.jl @@ -38,7 +38,7 @@ using PDMats: chol_lower, chol_upper @test sum(abs2, PDMats.chol_upper(ch_dense)' \ x) ≈ b # sparse version - if PDMats.HAVE_CHOLMOD + if 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 diff --git a/test/pdmtypes.jl b/test/pdmtypes.jl index 37c2e83..40ff000 100644 --- a/test/pdmtypes.jl +++ b/test/pdmtypes.jl @@ -1,4 +1,4 @@ -using LinearAlgebra, PDMats, SparseArrays, SuiteSparse +using LinearAlgebra, PDMats, SparseArrays using Test @testset "pd matrix types" begin @@ -10,8 +10,10 @@ using Test @test @test_deprecated(PDiagMat(d, d)) == PDiagMat(d) x = one(T) @test @test_deprecated(ScalMat(2, x, x)) == ScalMat(2, x) - s = SparseMatrixCSC{T}(I, 2, 2) - @test PDSparseMat(s, cholesky(s)).mat == PDSparseMat(s).mat == PDSparseMat(cholesky(s)).mat + if HAVE_CHOLMOD + s = SparseMatrixCSC{T}(I, 2, 2) + @test PDMat(s, cholesky(s)).mat == PDMat(s).mat == PDMat(cholesky(s)).mat + end end @testset "test the functionality" begin @@ -32,8 +34,10 @@ using Test @testset "ScalMat" begin test_pdmat(ScalMat(3,X), X*Matrix{T}(I, 3, 3), cmat_eq=true, verbose=1) end - @testset "PDSparseMat" begin - test_pdmat(PDSparseMat(sparse(M)), M, cmat_eq=true, verbose=1, t_eig=false) + if HAVE_CHOLMOD + @testset "PDMat from sparse matrix" begin + test_pdmat(PDMat(sparse(M)), M, cmat_eq=true, verbose=1, t_eig=false) + end end end end @@ -75,10 +79,10 @@ using Test end if HAVE_CHOLMOD - A = PDSparseMat(SparseMatrixCSC{T}(I, 2, 2)) - for R in (AbstractArray{S}, AbstractMatrix{S}, AbstractPDMat{S}, PDSparseMat{S}) + A = PDMat(SparseMatrixCSC{T}(I, 2, 2)) + for R in (AbstractArray{S}, AbstractMatrix{S}, AbstractPDMat{S}, PDMat{S}) B = @inferred(convert(R, A)) - @test B isa PDSparseMat{S} + @test B isa PDMat{S} @test B == A @test (B === A) === (S === T) @test (B.mat === A.mat) === (S === T) @@ -122,18 +126,15 @@ using Test @test z ≈ y end - # requires https://github.com/JuliaLang/julia/pull/32594 - if VERSION >= v"1.3.0-DEV.562" - z = x / PDMat(A) - @test typeof(z) === typeof(y) - @test size(z) == size(y) - @test z ≈ y - end + z = x / PDMat(A) + @test typeof(z) === typeof(y) + @test size(z) == size(y) + @test z ≈ y # right division not defined for CHOLMOD: - # `rdiv!(::Matrix{Float64}, ::SuiteSparse.CHOLMOD.Factor{Float64})` not defined + # `rdiv!(::Matrix{Float64}, ::SparseArrays.CHOLMOD.Factor{Float64})` not defined if !HAVE_CHOLMOD - z = x / PDSparseMat(sparse(first(A), 1, 1)) + z = x / PDMat(sparse(first(A), 1, 1)) @test typeof(z) === typeof(y) @test size(z) == size(y) @test z ≈ y @@ -141,15 +142,11 @@ using Test end @testset "PDMat from Cholesky decomposition of diagonal matrix (#137)" begin - # U'*U where U isa UpperTriangular etc. - # requires https://github.com/JuliaLang/julia/pull/33334 - if VERSION >= v"1.4.0-DEV.286" - x = rand(10, 10) - A = Diagonal(x' * x) - M = PDMat(cholesky(A)) - @test M isa PDMat{Float64, typeof(A)} - @test Matrix(M) ≈ A - end + x = rand(10, 10) + A = Diagonal(x' * x) + M = PDMat(cholesky(A)) + @test M isa PDMat{Float64, typeof(A)} + @test Matrix(M) ≈ A end @testset "AbstractPDMat constructors (#136)" begin @@ -158,10 +155,12 @@ using Test M = @inferred AbstractPDMat(A) @test M isa PDMat + @test cholesky(M) isa Cholesky @test Matrix(M) ≈ A M = @inferred AbstractPDMat(cholesky(A)) @test M isa PDMat + @test cholesky(M) isa Cholesky @test Matrix(M) ≈ A M = @inferred AbstractPDMat(Diagonal(A)) @@ -177,16 +176,13 @@ using Test @test Matrix(M) ≈ Diagonal(A) M = @inferred AbstractPDMat(sparse(A)) - @test M isa PDSparseMat + @test M isa PDMat + @test cholesky(M) isa SparseArrays.CHOLMOD.Factor @test Matrix(M) ≈ A - if VERSION < v"1.6" - # inference fails e.g. on Julia 1.0 - M = AbstractPDMat(cholesky(sparse(A))) - else - M = @inferred AbstractPDMat(cholesky(sparse(A))) - end - @test M isa PDSparseMat + M = @inferred AbstractPDMat(cholesky(sparse(A))) + @test M isa PDMat + @test cholesky(M) isa SparseArrays.CHOLMOD.Factor @test Matrix(M) ≈ A end end diff --git a/test/specialarrays.jl b/test/specialarrays.jl index b812a80..b02e1db 100644 --- a/test/specialarrays.jl +++ b/test/specialarrays.jl @@ -68,10 +68,7 @@ using StaticArrays Y = rand(5, 2) @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 - end + @test X / P ≈ X @test P \ x ≈ x @test P \ Y ≈ Y @test X_A_Xt(P, X) ≈ X * X' diff --git a/test/testutils.jl b/test/testutils.jl index 7571583..486a132 100644 --- a/test/testutils.jl +++ b/test/testutils.jl @@ -4,10 +4,16 @@ # the implementation of a subtype of AbstractPDMat # -using PDMats, SuiteSparse, Test +using PDMats, LinearAlgebra, SparseArrays, Test -const HAVE_CHOLMOD = isdefined(SuiteSparse, :CHOLMOD) -const PDMatType = HAVE_CHOLMOD ? Union{PDMat, PDSparseMat, PDiagMat} : Union{PDMat, PDiagMat} +const HAVE_CHOLMOD = isdefined(SparseArrays, :CHOLMOD) +const PDMatCholesky{T<:Real, S<:AbstractMatrix, C<:Cholesky} = PDMat{T, S, C} +if HAVE_CHOLMOD + const PDSparseMat{T<:Real, S<:AbstractSparseMatrix, C<:SparseArrays.CHOLMOD.Factor} = PDMat{T, S, C} + const PDMatType = Union{PDMatCholesky, PDSparseMat, PDiagMat, ScalMat} +else + const PDMatType = Union{PDMatCholesky, PDiagMat, ScalMat} +end ## driver function function test_pdmat(C, Cmat::Matrix; @@ -17,7 +23,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 @@ -120,7 +126,7 @@ function pdtest_diag(C, Cmat::Matrix, cmat_eq::Bool, verbose::Int) end end -function pdtest_cholesky(C::Union{PDMat, PDiagMat}, Cmat::Matrix, cmat_eq::Bool, verbose::Int) +function pdtest_cholesky(C::PDMatType, Cmat::Matrix, cmat_eq::Bool, verbose::Int) _pdt(verbose, "cholesky") if cmat_eq @test cholesky(C).U == cholesky(Cmat).U @@ -132,8 +138,8 @@ end if HAVE_CHOLMOD function pdtest_cholesky(C::PDSparseMat, Cmat::Matrix, cmat_eq::Bool, verbose::Int) _pdt(verbose, "cholesky") - # We special case PDSparseMat because we can't perform equality checks on - # `SuiteSparse.CHOLMOD.Factor`s and `SuiteSparse.CHOLMOD.FactorComponent`s + # We handle this case specially because we can't perform equality checks on + # `SparseArrays.CHOLMOD.Factor`s and `SparseArrays.CHOLMOD.FactorComponent`s @test diag(cholesky(C)) ≈ diag(cholesky(Cmat).U) # NOTE: `==` also doesn't work because `diag(cholesky(C))` will return `Vector{Float64}` # even if the inputs are `Float32`s. @@ -168,7 +174,7 @@ function pdtest_det(C, Cmat::Matrix, verbose::Int) @test det(C) ≈ det(Cmat) # generic fallback in LinearAlgebra performs LU decomposition - if C isa Union{PDMat,PDiagMat,ScalMat} + if C isa Union{PDMatCholesky,PDiagMat,ScalMat} @test iszero(@allocated det(C)) end end @@ -178,7 +184,7 @@ function pdtest_logdet(C, Cmat::Matrix, verbose::Int) @test logdet(C) ≈ logdet(Cmat) # generic fallback in LinearAlgebra performs LU decomposition - if C isa Union{PDMat,PDiagMat,ScalMat} + if C isa Union{PDMatCholesky,PDiagMat,ScalMat} @test iszero(@allocated logdet(C)) end end @@ -230,10 +236,10 @@ function pdtest_div(C, Imat::Matrix, X::Matrix, verbose::Int) @assert d == size(C, 1) == size(C, 2) @assert size(Imat) == size(C) @test C \ X ≈ Imat * X - # Right division with Choleskyrequires https://github.com/JuliaLang/julia/pull/32594 + # Right division with Cholesky requires https://github.com/JuliaLang/julia/pull/32594 # CHOLMOD throws error since no method is found for - # `rdiv!(::Matrix{Float64}, ::SuiteSparse.CHOLMOD.Factor{Float64})` - check_rdiv = !(C isa PDMat && VERSION < v"1.3.0-DEV.562") && !(C isa PDSparseMat && HAVE_CHOLMOD) + # `rdiv!(::Matrix{Float64}, ::SparseArrays.CHOLMOD.Factor{Float64})` + check_rdiv = !(HAVE_CHOLMOD && C isa PDSparseMat) check_rdiv && @test Matrix(X') / C ≈ (C \ X)' for i = 1:n @@ -347,7 +353,9 @@ end _randPDMat(T, n) = (X = randn(T, n, n); PDMat(X * X' + LinearAlgebra.I)) _randPDiagMat(T, n) = PDiagMat(rand(T, n)) _randScalMat(T, n) = ScalMat(n, rand(T)) -_randPDSparseMat(T, n) = (X = T.(sprand(n, 1, 0.5)); PDSparseMat(X * X' + LinearAlgebra.I)) +if HAVE_CHOLMOD + _randPDSparseMat(T, n) = (X = sprand(T, n, 1, 0.5); PDMat(X * X' + LinearAlgebra.I)) +end function _pd_compare(A::AbstractPDMat, B::AbstractPDMat) @test size(A) == size(B) From e4ed4ed78bbebeb55a2b880b83c64b328a0e7e21 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Wed, 4 Oct 2023 12:03:45 +0200 Subject: [PATCH 02/14] Replace tab with whitespace --- test/testutils.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/test/testutils.jl b/test/testutils.jl index 486a132..20c810b 100644 --- a/test/testutils.jl +++ b/test/testutils.jl @@ -23,7 +23,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 @@ -175,7 +175,7 @@ function pdtest_det(C, Cmat::Matrix, verbose::Int) # generic fallback in LinearAlgebra performs LU decomposition if C isa Union{PDMatCholesky,PDiagMat,ScalMat} - @test iszero(@allocated det(C)) + @test iszero(@allocated det(C)) end end @@ -185,7 +185,7 @@ function pdtest_logdet(C, Cmat::Matrix, verbose::Int) # generic fallback in LinearAlgebra performs LU decomposition if C isa Union{PDMatCholesky,PDiagMat,ScalMat} - @test iszero(@allocated logdet(C)) + @test iszero(@allocated logdet(C)) end end From 91f6d14074a8e96577a794509ca49fbe22dde4e7 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Wed, 4 Oct 2023 12:24:41 +0200 Subject: [PATCH 03/14] Try to add back support for Julia < 1.9 --- .github/workflows/ci.yml | 3 ++- Project.toml | 5 ++++- src/pdiagmat.jl | 7 ++++++- src/pdmat.jl | 8 +++++++- src/scalmat.jl | 13 +++++++++---- src/utils.jl | 20 ++++++++++++++++++++ test/pdmtypes.jl | 32 ++++++++++++++++++++++---------- test/specialarrays.jl | 5 ++++- test/testutils.jl | 2 +- 9 files changed, 75 insertions(+), 20 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 93b0cf0..7abc2de 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -12,7 +12,8 @@ jobs: fail-fast: false matrix: version: - - '1.9' + - '1.0' + - '1.6' - '1' - nightly os: diff --git a/Project.toml b/Project.toml index cc9242d..5dc5781 100644 --- a/Project.toml +++ b/Project.toml @@ -4,9 +4,11 @@ version = "0.12.0" [deps] LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" [compat] -julia = "1.9" +julia = "1" [extensions] PDMatsSparseArraysExt = "SparseArrays" @@ -22,3 +24,4 @@ test = ["BandedMatrices", "SparseArrays", "StaticArrays", "Test"] [weakdeps] SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" diff --git a/src/pdiagmat.jl b/src/pdiagmat.jl index 948a8d5..db075c7 100644 --- a/src/pdiagmat.jl +++ b/src/pdiagmat.jl @@ -62,7 +62,12 @@ function \(a::PDiagMat, x::AbstractVecOrMat) end function /(x::AbstractVecOrMat, a::PDiagMat) @check_argdims a.dim == size(x, 2) - return x ./ (x isa AbstractVector ? a.diag : a.diag') + if VERSION < v"1.9-" + # return matrix for 1-element vectors `x`, consistent with LinearAlgebra < 1.9 + return reshape(x, Val(2)) ./ permutedims(a.diag) # = (a' \ x') + else + return x ./ (x isa AbstractVector ? a.diag : a.diag') + end end Base.kron(A::PDiagMat, B::PDiagMat) = PDiagMat(vec(permutedims(A.diag) .* B.diag)) diff --git a/src/pdmat.jl b/src/pdmat.jl index f6f1972..efabf8a 100644 --- a/src/pdmat.jl +++ b/src/pdmat.jl @@ -54,7 +54,13 @@ end *(a::PDMat, x::AbstractMatrix) = a.mat * x \(a::PDMat, x::AbstractVecOrMat) = a.chol \ x function /(x::AbstractVecOrMat, a::PDMat) - return x isa AbstractVector ? vec(reshape(x, Val(2)) / a.chol) : x / a.chol + # /(::AbstractVector, ::Cholesky) is not defined + if VERSION < v"1.9-" + # return matrix for 1-element vectors `x`, consistent with LinearAlgebra + return reshape(x, Val(2)) / a.chol + else + return x isa AbstractVector ? vec(reshape(x, Val(2)) / a.chol) : x / a.chol + end end ### Algebra diff --git a/src/scalmat.jl b/src/scalmat.jl index 878c166..ea672ec 100644 --- a/src/scalmat.jl +++ b/src/scalmat.jl @@ -54,7 +54,12 @@ function \(a::ScalMat, x::AbstractVecOrMat) end function /(x::AbstractVecOrMat, a::ScalMat) @check_argdims a.dim == size(x, 2) - return x / a.value + if VERSION < v"1.9-" + # return matrix for 1-element vectors `x`, consistent with LinearAlgebra < 1.9 + return reshape(x, Val(2)) / a.value + else + return x / a.value + end end Base.kron(A::ScalMat, B::ScalMat) = ScalMat(A.dim * B.dim, A.value * B.value ) @@ -72,7 +77,7 @@ LinearAlgebra.sqrt(a::ScalMat) = ScalMat(a.dim, sqrt(a.value)) function whiten!(r::AbstractVecOrMat, a::ScalMat, x::AbstractVecOrMat) @check_argdims LinearAlgebra.checksquare(a) == size(x, 1) - ldiv!(r, sqrt(a.value), x) + _ldiv!(r, sqrt(a.value), x) end function unwhiten!(r::AbstractVecOrMat, a::ScalMat, x::AbstractVecOrMat) @@ -125,10 +130,10 @@ end function X_invA_Xt(a::ScalMat, x::Matrix) @check_argdims LinearAlgebra.checksquare(a) == size(x, 2) - rdiv!(x * transpose(x), a.value) + _rdiv!(x * transpose(x), a.value) end function Xt_invA_X(a::ScalMat, x::Matrix) @check_argdims LinearAlgebra.checksquare(a) == size(x, 1) - rdiv!(transpose(x) * x, a.value) + _rdiv!(transpose(x) * x, a.value) end diff --git a/src/utils.jl b/src/utils.jl index b34fc50..1da6c17 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -104,3 +104,23 @@ function colwise_sumsqinv!(r::AbstractArray, a::AbstractMatrix, c::Real) return r end +# `rdiv!(::AbstractArray, ::Number)` was introduced in Julia 1.2 +# https://github.com/JuliaLang/julia/pull/31179 +@static if VERSION < v"1.2.0-DEV.385" + function _rdiv!(X::AbstractArray, s::Number) + @simd for I in eachindex(X) + @inbounds X[I] /= s + end + X + end +else + _rdiv!(X::AbstractArray, s::Number) = rdiv!(X, s) +end + +# `ldiv!(::AbstractArray, ::Number, ::AbstractArray)` was introduced in Julia 1.4 +# https://github.com/JuliaLang/julia/pull/33806 +@static if VERSION < v"1.4.0-DEV.635" + _ldiv!(Y::AbstractArray, s::Number, X::AbstractArray) = Y .= s .\ X +else + _ldiv!(Y::AbstractArray, s::Number, X::AbstractArray) = ldiv!(Y, s, X) +end diff --git a/test/pdmtypes.jl b/test/pdmtypes.jl index 40ff000..961e592 100644 --- a/test/pdmtypes.jl +++ b/test/pdmtypes.jl @@ -126,10 +126,13 @@ using Test @test z ≈ y end - z = x / PDMat(A) - @test typeof(z) === typeof(y) - @test size(z) == size(y) - @test z ≈ y + # requires https://github.com/JuliaLang/julia/pull/32594 + if VERSION >= v"1.3.0-DEV.562" + z = x / PDMat(A) + @test typeof(z) === typeof(y) + @test size(z) == size(y) + @test z ≈ y + end # right division not defined for CHOLMOD: # `rdiv!(::Matrix{Float64}, ::SparseArrays.CHOLMOD.Factor{Float64})` not defined @@ -142,11 +145,15 @@ using Test end @testset "PDMat from Cholesky decomposition of diagonal matrix (#137)" begin - x = rand(10, 10) - A = Diagonal(x' * x) - M = PDMat(cholesky(A)) - @test M isa PDMat{Float64, typeof(A)} - @test Matrix(M) ≈ A + # U'*U where U isa UpperTriangular etc. + # requires https://github.com/JuliaLang/julia/pull/33334 + if VERSION >= v"1.4.0-DEV.286" + x = rand(10, 10) + A = Diagonal(x' * x) + M = PDMat(cholesky(A)) + @test M isa PDMat{Float64, typeof(A)} + @test Matrix(M) ≈ A + end end @testset "AbstractPDMat constructors (#136)" begin @@ -180,7 +187,12 @@ using Test @test cholesky(M) isa SparseArrays.CHOLMOD.Factor @test Matrix(M) ≈ A - M = @inferred AbstractPDMat(cholesky(sparse(A))) + if VERSION < v"1.6" + # inference fails e.g. on Julia 1.0 + M = AbstractPDMat(cholesky(sparse(A))) + else + M = @inferred AbstractPDMat(cholesky(sparse(A))) + end @test M isa PDMat @test cholesky(M) isa SparseArrays.CHOLMOD.Factor @test Matrix(M) ≈ A diff --git a/test/specialarrays.jl b/test/specialarrays.jl index b02e1db..b812a80 100644 --- a/test/specialarrays.jl +++ b/test/specialarrays.jl @@ -68,7 +68,10 @@ using StaticArrays Y = rand(5, 2) @test P * x ≈ x @test P * Y ≈ Y - @test X / P ≈ X + # Right division with Cholesky requires https://github.com/JuliaLang/julia/pull/32594 + if VERSION >= v"1.3.0-DEV.562" + @test X / P ≈ X + end @test P \ x ≈ x @test P \ Y ≈ Y @test X_A_Xt(P, X) ≈ X * X' diff --git a/test/testutils.jl b/test/testutils.jl index 20c810b..2a1d476 100644 --- a/test/testutils.jl +++ b/test/testutils.jl @@ -239,7 +239,7 @@ function pdtest_div(C, Imat::Matrix, X::Matrix, verbose::Int) # Right division with Cholesky requires https://github.com/JuliaLang/julia/pull/32594 # CHOLMOD throws error since no method is found for # `rdiv!(::Matrix{Float64}, ::SparseArrays.CHOLMOD.Factor{Float64})` - check_rdiv = !(HAVE_CHOLMOD && C isa PDSparseMat) + check_rdiv = !(C isa PDMatCholesky && VERSION < v"1.3.0-DEV.562") && !(HAVE_CHOLMOD && C isa PDSparseMat) check_rdiv && @test Matrix(X') / C ≈ (C \ X)' for i = 1:n From e817ceb910999b7ebf62660e554a7b5810e59ebb Mon Sep 17 00:00:00 2001 From: David Widmann Date: Wed, 4 Oct 2023 12:41:10 +0200 Subject: [PATCH 04/14] Fix support on older Julia versions --- ext/PDMatsSparseArraysExt/PDMatsSparseArraysExt.jl | 7 ++++++- ext/PDMatsSparseArraysExt/chol.jl | 7 ++++++- src/PDMats.jl | 4 ++++ test/pdmtypes.jl | 4 ++-- test/testutils.jl | 14 ++++++++++++-- 5 files changed, 30 insertions(+), 6 deletions(-) diff --git a/ext/PDMatsSparseArraysExt/PDMatsSparseArraysExt.jl b/ext/PDMatsSparseArraysExt/PDMatsSparseArraysExt.jl index 7f684e0..4c1f88b 100644 --- a/ext/PDMatsSparseArraysExt/PDMatsSparseArraysExt.jl +++ b/ext/PDMatsSparseArraysExt/PDMatsSparseArraysExt.jl @@ -5,7 +5,12 @@ using SparseArrays using PDMats.LinearAlgebra -const HAVE_CHOLMOD = isdefined(SparseArrays, :CHOLMOD) +if isdefined(Base, :get_extension) + const HAVE_CHOLMOD = isdefined(SparseArrays, :CHOLMOD) +else + import SuiteSparse + const HAVE_CHOLMOD = isdefined(SuiteSparse, :CHOLMOD) +end if HAVE_CHOLMOD include("chol.jl") diff --git a/ext/PDMatsSparseArraysExt/chol.jl b/ext/PDMatsSparseArraysExt/chol.jl index 0c18ec7..f04c084 100644 --- a/ext/PDMatsSparseArraysExt/chol.jl +++ b/ext/PDMatsSparseArraysExt/chol.jl @@ -1,4 +1,9 @@ -const CholTypeSparse{T} = SparseArrays.CHOLMOD.Factor{T} + +if isdefined(Base, :get_extension) + const CholTypeSparse{T} = SparseArrays.CHOLMOD.Factor{T} +else + const CholTypeSparse{T} = SuiteSparse.CHOLMOD.Factor{T} +end # Take into account pivoting! PDMats.chol_lower(cf::CholTypeSparse) = cf.PtL diff --git a/src/PDMats.jl b/src/PDMats.jl index 56f49ad..ff90ac5 100644 --- a/src/PDMats.jl +++ b/src/PDMats.jl @@ -48,4 +48,8 @@ module PDMats include("deprecates.jl") + # Support for sparse arrays + if !isdefined(Base, :get_extension) + include("../ext/PDMatsSparseArraysExt/PDMatsSparseArraysExt.jl") + end end # module diff --git a/test/pdmtypes.jl b/test/pdmtypes.jl index 961e592..2fa6e43 100644 --- a/test/pdmtypes.jl +++ b/test/pdmtypes.jl @@ -184,7 +184,7 @@ using Test M = @inferred AbstractPDMat(sparse(A)) @test M isa PDMat - @test cholesky(M) isa SparseArrays.CHOLMOD.Factor + @test cholesky(M) isa CHOLMOD.Factor @test Matrix(M) ≈ A if VERSION < v"1.6" @@ -194,7 +194,7 @@ using Test M = @inferred AbstractPDMat(cholesky(sparse(A))) end @test M isa PDMat - @test cholesky(M) isa SparseArrays.CHOLMOD.Factor + @test cholesky(M) isa CHOLMOD.Factor @test Matrix(M) ≈ A end end diff --git a/test/testutils.jl b/test/testutils.jl index 2a1d476..2a7d2b1 100644 --- a/test/testutils.jl +++ b/test/testutils.jl @@ -6,10 +6,20 @@ using PDMats, LinearAlgebra, SparseArrays, Test -const HAVE_CHOLMOD = isdefined(SparseArrays, :CHOLMOD) +if isdefined(Base, :get_extension) + const HAVE_CHOLMOD = isdefined(SparseArrays, :CHOLMOD) +else + import SuiteSparse + const HAVE_CHOLMOD = isdefined(SuiteSparse, :CHOLMOD) +end const PDMatCholesky{T<:Real, S<:AbstractMatrix, C<:Cholesky} = PDMat{T, S, C} if HAVE_CHOLMOD - const PDSparseMat{T<:Real, S<:AbstractSparseMatrix, C<:SparseArrays.CHOLMOD.Factor} = PDMat{T, S, C} + if isdefined(Base, :get_extension) + const CHOLMOD = SparseArrays.CHOLMOD + else + const CHOLMOD = SuiteSparse.CHOLMOD + end + const PDSparseMat{T<:Real, S<:AbstractSparseMatrix, C<:CHOLMOD.Factor} = PDMat{T, S, C} const PDMatType = Union{PDMatCholesky, PDSparseMat, PDiagMat, ScalMat} else const PDMatType = Union{PDMatCholesky, PDiagMat, ScalMat} From dd947d65dcd9c869606b9c5b8a811610294ceec7 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Tue, 24 Oct 2023 00:45:04 +0200 Subject: [PATCH 05/14] Fix test errors on Julia 1.0 --- Project.toml | 11 +++++------ ext/PDMatsSparseArraysExt/PDMatsSparseArraysExt.jl | 5 +++++ src/utils.jl | 4 ---- 3 files changed, 10 insertions(+), 10 deletions(-) diff --git a/Project.toml b/Project.toml index c6e8534..5dc5781 100644 --- a/Project.toml +++ b/Project.toml @@ -7,16 +7,12 @@ LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" -[weakdeps] -SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" -SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" +[compat] +julia = "1" [extensions] PDMatsSparseArraysExt = "SparseArrays" -[compat] -julia = "1" - [extras] BandedMatrices = "aae01518-5342-5314-be14-df237901396f" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" @@ -26,3 +22,6 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] test = ["BandedMatrices", "SparseArrays", "StaticArrays", "Test"] +[weakdeps] +SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9" diff --git a/ext/PDMatsSparseArraysExt/PDMatsSparseArraysExt.jl b/ext/PDMatsSparseArraysExt/PDMatsSparseArraysExt.jl index 4c1f88b..e3fa624 100644 --- a/ext/PDMatsSparseArraysExt/PDMatsSparseArraysExt.jl +++ b/ext/PDMatsSparseArraysExt/PDMatsSparseArraysExt.jl @@ -12,6 +12,11 @@ else const HAVE_CHOLMOD = isdefined(SuiteSparse, :CHOLMOD) end +# https://github.com/JuliaLang/julia/pull/29749 +if VERSION < v"1.1.0-DEV.792" + eachcol(A::AbstractVecOrMat) = (view(A, :, i) for i in axes(A, 2)) +end + if HAVE_CHOLMOD include("chol.jl") include("pdsparsemat.jl") diff --git a/src/utils.jl b/src/utils.jl index 98d2900..fecc3c8 100644 --- a/src/utils.jl +++ b/src/utils.jl @@ -125,7 +125,3 @@ else _ldiv!(Y::AbstractArray, s::Number, X::AbstractArray) = ldiv!(Y, s, X) end -# https://github.com/JuliaLang/julia/pull/29749 -if VERSION < v"1.1.0-DEV.792" - eachcol(A::AbstractVecOrMat) = (view(A, :, i) for i in axes(A, 2)) -end From 2dc005d4c46be9293a385f6f3549c3733ede8149 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Mon, 20 Nov 2023 20:08:25 +0100 Subject: [PATCH 06/14] Update README.md --- README.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/README.md b/README.md index 6235524..be14803 100644 --- a/README.md +++ b/README.md @@ -29,7 +29,7 @@ Elemenent types are in princple all Real types, but in practice this is limited ```julia struct PDMat{T<:Real,S<:AbstractMatrix,C<:Factorization} <: AbstractPDMat{T} mat::S # input matrix - chol::Cholesky{T,S} # Cholesky factorization of mat + chol::C # Cholesky factorization of mat end # Constructors From deb575d29dc54c3c9e01ac4deff3900af5e72288 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Tue, 21 Nov 2023 15:56:39 +0100 Subject: [PATCH 07/14] Improve type parameter Co-authored-by: Tim Holy --- src/pdmat.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pdmat.jl b/src/pdmat.jl index 6ba5eb3..1211e80 100644 --- a/src/pdmat.jl +++ b/src/pdmat.jl @@ -1,7 +1,7 @@ """ Full positive definite matrix together with a Cholesky factorization object. """ -struct PDMat{T<:Real,S<:AbstractMatrix,C<:Factorization} <: AbstractPDMat{T} +struct PDMat{T<:Real,S<:AbstractMatrix{T},C<:Factorization} <: AbstractPDMat{T} mat::S chol::C From c27aab2039d4044cb04a457a53fce22021c9e087 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Tue, 21 Nov 2023 23:15:58 +0100 Subject: [PATCH 08/14] Rename field `chol` to `fact` --- README.md | 4 +- ext/PDMatsSparseArraysExt/chol.jl | 1 - ext/PDMatsSparseArraysExt/pdsparsemat.jl | 36 +++++----- src/pdmat.jl | 84 ++++++++++++------------ test/pdmtypes.jl | 4 +- 5 files changed, 63 insertions(+), 66 deletions(-) diff --git a/README.md b/README.md index be14803..13b8ad8 100644 --- a/README.md +++ b/README.md @@ -27,9 +27,9 @@ Elemenent types are in princple all Real types, but in practice this is limited * `PDMat`: full covariance matrix, defined as ```julia -struct PDMat{T<:Real,S<:AbstractMatrix,C<:Factorization} <: AbstractPDMat{T} +struct PDMat{T<:Real,S<:AbstractMatrix{T},F<:Factorization} <: AbstractPDMat{T} mat::S # input matrix - chol::C # Cholesky factorization of mat + fact::F # factorization of mat end # Constructors diff --git a/ext/PDMatsSparseArraysExt/chol.jl b/ext/PDMatsSparseArraysExt/chol.jl index f04c084..195cf4d 100644 --- a/ext/PDMatsSparseArraysExt/chol.jl +++ b/ext/PDMatsSparseArraysExt/chol.jl @@ -1,4 +1,3 @@ - if isdefined(Base, :get_extension) const CholTypeSparse{T} = SparseArrays.CHOLMOD.Factor{T} else diff --git a/ext/PDMatsSparseArraysExt/pdsparsemat.jl b/ext/PDMatsSparseArraysExt/pdsparsemat.jl index b5c7d93..f6bef3f 100644 --- a/ext/PDMatsSparseArraysExt/pdsparsemat.jl +++ b/ext/PDMatsSparseArraysExt/pdsparsemat.jl @@ -1,12 +1,9 @@ """ Sparse positive definite matrix together with a Cholesky factorization object. """ -const PDSparseMat{T<:Real,S<:AbstractSparseMatrix,C<:CholTypeSparse} = PDMat{T,S,C} +const PDSparseMat{T<:Real,S<:AbstractSparseMatrix{T},C<:CholTypeSparse} = PDMat{T,S,C} function PDMats.PDMat(mat::AbstractSparseMatrix, chol::CholTypeSparse) - d = LinearAlgebra.checksquare(mat) - size(chol, 1) == d || - throw(DimensionMismatch("Dimensions of mat and chol are inconsistent.")) PDMat{eltype(mat),typeof(mat),typeof(chol)}(mat, chol) end Base.@deprecate PDMat{T,S}(d::Int, m::AbstractSparseMatrix{T}, c::CholTypeSparse) where {T,S} PDSparseMat{T,S,typeof(c)}(m, c) @@ -22,20 +19,18 @@ function Base.convert(::Type{PDMat{T}}, a::PDSparseMat) where {T<:Real} # So there is no point in recomputing `cholesky(mat)` and we just reuse # the existing Cholesky factorization mat = convert(AbstractMatrix{T}, a.mat) - return PDMat{T,typeof(mat),typeof(a.chol)}(mat, a.chol) + return PDMat{T,typeof(mat),typeof(a.fact)}(mat, a.fact) end ### Arithmetics -Base.:\(a::PDSparseMat{T}, x::AbstractVecOrMat{T}) where {T<:Real} = convert(Array{T},a.chol \ convert(Array{Float64},x)) #to avoid limitations in sparse factorization library CHOLMOD, see e.g., julia issue #14076 -Base.:/(x::AbstractVecOrMat{T}, a::PDSparseMat{T}) where {T<:Real} = convert(Array{T},convert(Array{Float64},x) / a.chol ) +Base.:\(a::PDSparseMat{T}, x::AbstractVecOrMat{T}) where {T<:Real} = convert(Array{T},a.fact \ convert(Array{Float64},x)) #to avoid limitations in sparse factorization library CHOLMOD, see e.g., julia issue #14076 +Base.:/(x::AbstractVecOrMat{T}, a::PDSparseMat{T}) where {T<:Real} = convert(Array{T},convert(Array{Float64},x) / a.fact ) ### Algebra -Base.inv(a::PDSparseMat{T}) where {T<:Real} = PDMat(inv(a.mat)) -LinearAlgebra.det(a::PDSparseMat) = det(a.chol) -LinearAlgebra.logdet(a::PDSparseMat) = logdet(a.chol) -LinearAlgebra.sqrt(A::PDSparseMat) = PDMat(sqrt(Hermitian(Matrix(A)))) +Base.inv(a::PDSparseMat) = PDMat(inv(a.mat)) +LinearAlgebra.cholesky(a::PDSparseMat) = a.fact ### whiten and unwhiten @@ -43,15 +38,15 @@ function PDMats.whiten!(r::AbstractVecOrMat, a::PDSparseMat, x::AbstractVecOrMat PDMats.@check_argdims axes(r) == axes(x) PDMats.@check_argdims a.dim == size(x, 1) # Can't use `ldiv!` due to missing support in SparseArrays - return copyto!(r, PDMats.chol_lower(a.chol) \ x) + return copyto!(r, PDMats.chol_lower(cholesky(a)) \ x) end function PDMats.unwhiten!(r::AbstractVecOrMat, a::PDSparseMat, x::AbstractVecOrMat) PDMats.@check_argdims axes(r) == axes(x) PDMats.@check_argdims a.dim == size(x, 1) # `*` is not defined for `PtL` factor components, - # so we can't use `chol_lower(a.chol) * x` - C = a.chol + # so we can't use `chol_lower(cholesky(a)) * x` + C = cholesky(a) PtL = sparse(C.L)[C.p, :] return copyto!(r, PtL * x) end @@ -64,8 +59,8 @@ end function PDMats.unwhiten(a::PDSparseMat, x::AbstractVecOrMat) PDMats.@check_argdims a.dim == size(x, 1) # `*` is not defined for `PtL` factor components, - # so we can't use `chol_lower(a.chol) * x` - C = a.chol + # so we can't use `chol_lower(cholesky(a)) * x` + C = cholesky(a) PtL = sparse(C.L)[C.p, :] return PtL * x end @@ -100,7 +95,7 @@ end function PDMats.invquad(a::PDSparseMat, x::AbstractVecOrMat) PDMats.@check_argdims a.dim == size(x, 1) - z = a.chol \ x + z = cholesky(a) \ x return x isa AbstractVector ? dot(x, z) : map(dot, eachcol(x), eachcol(z)) end @@ -108,9 +103,10 @@ function PDMats.invquad!(r::AbstractArray, a::PDSparseMat, x::AbstractMatrix) PDMats.@check_argdims eachindex(r) == axes(x, 2) PDMats.@check_argdims a.dim == size(x, 1) # Can't use `ldiv!` with buffer due to missing support in SparseArrays + C = cholesky(a) @inbounds for i in axes(x, 2) xi = view(x, :, i) - r[i] = dot(xi, a.chol \ xi) + r[i] = dot(xi, C \ xi) end return r end @@ -134,12 +130,12 @@ end function PDMats.X_invA_Xt(a::PDSparseMat, x::AbstractMatrix{<:Real}) PDMats.@check_argdims a.dim == size(x, 2) - z = a.chol \ collect(transpose(x)) + z = cholesky(a) \ collect(transpose(x)) return Symmetric(x * z) end function PDMats.Xt_invA_X(a::PDSparseMat, x::AbstractMatrix{<:Real}) PDMats.@check_argdims a.dim == size(x, 1) - z = a.chol \ x + z = cholesky(a) \ x return Symmetric(transpose(x) * z) end diff --git a/src/pdmat.jl b/src/pdmat.jl index 1211e80..ecd945a 100644 --- a/src/pdmat.jl +++ b/src/pdmat.jl @@ -1,18 +1,20 @@ """ -Full positive definite matrix together with a Cholesky factorization object. +Full positive definite matrix together with a factorization object. """ -struct PDMat{T<:Real,S<:AbstractMatrix{T},C<:Factorization} <: AbstractPDMat{T} +struct PDMat{T<:Real,S<:AbstractMatrix{T},F<:Factorization} <: AbstractPDMat{T} mat::S - chol::C - - PDMat{T,S,C}(m::AbstractMatrix{T}, c::Factorization) where {T,S,C} = new{T,S,C}(m,c) + fact::F + + function PDMat{T,S,C}(mat::S, fact::F) where {T,S<:AbstractMatrix{T},F<:Factorization} + d = LinearAlgebra.checksquare(mat) + if size(fact) != (d, d) + throw(DimensionMismatch("Dimensions of the matrix and the factorization are inconsistent.")) + end + new{T,S,F}(mat, fact) + end end function PDMat(mat::AbstractMatrix{T}, chol::Cholesky) where {T<:Real} - d = LinearAlgebra.checksquare(mat) - if size(chol, 1) != d - throw(DimensionMismatch("Dimensions of mat and chol are inconsistent.")) - end PDMat{T,typeof(mat),typeof(chol)}(mat, chol) end @@ -25,14 +27,17 @@ function Base.getproperty(a::PDMat, s::Symbol) end return getfield(a, s) end -Base.propertynames(::PDMat) = (:mat, :chol, :dim) +Base.propertynames(::PDMat) = (:mat, :fact, :dim) AbstractPDMat(A::Cholesky) = PDMat(A) +### Abstract type alias +const PDMatCholesky{T<:Real,S<:AbstractMatrix{T}} = PDMat{T,S,<:Cholesky} + ### Conversion Base.convert(::Type{PDMat{T}}, a::PDMat{T}) where {T<:Real} = a function Base.convert(::Type{PDMat{T}}, a::PDMat) where {T<:Real} - chol = convert(Factorization{float(T)}, a.chol) + chol = convert(Factorization{float(T)}, a.fact) mat = convert(AbstractMatrix{T}, a.mat) return PDMat{T,typeof(mat),typeof(chol)}(mat, chol) end @@ -43,7 +48,7 @@ Base.convert(::Type{AbstractPDMat{T}}, a::PDMat) where {T<:Real} = convert(PDMat Base.size(a::PDMat) = (a.dim, a.dim) Base.Matrix(a::PDMat) = Matrix(a.mat) LinearAlgebra.diag(a::PDMat) = diag(a.mat) -LinearAlgebra.cholesky(a::PDMat) = a.chol +LinearAlgebra.cholesky(a::PDMatCholesky) = a.chol ### Inheriting from AbstractMatrix @@ -60,60 +65,57 @@ end *(a::PDMat, c::Real) = PDMat(a.mat * c) *(a::PDMat, x::AbstractVector) = a.mat * x *(a::PDMat, x::AbstractMatrix) = a.mat * x -\(a::PDMat, x::AbstractVecOrMat) = a.chol \ x +\(a::PDMat, x::AbstractVecOrMat) = a.fact \ x function /(x::AbstractVecOrMat, a::PDMat) # /(::AbstractVector, ::Cholesky) is not defined if VERSION < v"1.9-" # return matrix for 1-element vectors `x`, consistent with LinearAlgebra - return reshape(x, Val(2)) / a.chol + return reshape(x, Val(2)) / a.fact else - return x isa AbstractVector ? vec(reshape(x, Val(2)) / a.chol) : x / a.chol + return x isa AbstractVector ? vec(reshape(x, Val(2)) / a.fact) : x / a.fact end end ### Algebra -Base.inv(a::PDMat) = PDMat(inv(a.chol)) -LinearAlgebra.det(a::PDMat) = det(a.chol) -LinearAlgebra.logdet(a::PDMat) = logdet(a.chol) +Base.inv(a::PDMat) = PDMat(inv(a.fact)) +LinearAlgebra.det(a::PDMat) = det(a.fact) +LinearAlgebra.logdet(a::PDMat) = logdet(a.fact) LinearAlgebra.eigmax(a::PDMat) = eigmax(a.mat) LinearAlgebra.eigmin(a::PDMat) = eigmin(a.mat) LinearAlgebra.sqrt(A::PDMat) = PDMat(sqrt(Hermitian(A.mat))) -function Base.kron( - A::PDMat{<:Real,<:AbstractMatrix,<:Cholesky}, - B::PDMat{<:Real,<:AbstractMatrix,<:Cholesky}, -) - return PDMat(kron(A.mat, B.mat), Cholesky(kron(A.chol.U, B.chol.U), 'U', A.chol.info)) +function Base.kron(A::PDMatCholesky, B::PDMatCholesky) + return PDMat(kron(A.mat, B.mat), Cholesky(kron(A.fact.U, B.fact.U), 'U', A.fact.info)) end ### (un)whitening -function whiten!(r::AbstractVecOrMat, a::PDMat, x::AbstractVecOrMat) +function whiten!(r::AbstractVecOrMat, a::PDMatCholesky, x::AbstractVecOrMat) @check_argdims axes(r) == axes(x) @check_argdims a.dim == size(x, 1) v = _rcopy!(r, x) return ldiv!(chol_lower(cholesky(a)), v) end -function unwhiten!(r::AbstractVecOrMat, a::PDMat, x::AbstractVecOrMat) +function unwhiten!(r::AbstractVecOrMat, a::PDMatCholesky, x::AbstractVecOrMat) @check_argdims axes(r) == axes(x) @check_argdims a.dim == size(x, 1) v = _rcopy!(r, x) return lmul!(chol_lower(cholesky(a)), v) end -function whiten(a::PDMat, x::AbstractVecOrMat) +function whiten(a::PDMatCholesky, x::AbstractVecOrMat) @check_argdims a.dim == size(x, 1) return chol_lower(cholesky(a)) \ x end -function unwhiten(a::PDMat, x::AbstractVecOrMat) +function unwhiten(a::PDMatCholesky, x::AbstractVecOrMat) @check_argdims a.dim == size(x, 1) return chol_lower(cholesky(a)) * x end ## quad/invquad -function quad(a::PDMat, x::AbstractVecOrMat) +function quad(a::PDMatCholesky, x::AbstractVecOrMat) @check_argdims a.dim == size(x, 1) aU_x = chol_upper(cholesky(a)) * x if x isa AbstractVector @@ -123,7 +125,7 @@ function quad(a::PDMat, x::AbstractVecOrMat) end end -function quad!(r::AbstractArray, a::PDMat, x::AbstractMatrix) +function quad!(r::AbstractArray, a::PDMatCholesky, x::AbstractMatrix) @check_argdims eachindex(r) == axes(x, 2) @check_argdims a.dim == size(x, 1) aU = chol_upper(cholesky(a)) @@ -136,7 +138,7 @@ function quad!(r::AbstractArray, a::PDMat, x::AbstractMatrix) return r end -function invquad(a::PDMat, x::AbstractVecOrMat) +function invquad(a::PDMatCholesky, x::AbstractVecOrMat) @check_argdims a.dim == size(x, 1) inv_aL_x = chol_lower(cholesky(a)) \ x if x isa AbstractVector @@ -146,7 +148,7 @@ function invquad(a::PDMat, x::AbstractVecOrMat) end end -function invquad!(r::AbstractArray, a::PDMat, x::AbstractMatrix) +function invquad!(r::AbstractArray, a::PDMatCholesky, x::AbstractMatrix) @check_argdims eachindex(r) == axes(x, 2) @check_argdims a.dim == size(x, 1) aL = chol_lower(cholesky(a)) @@ -161,39 +163,39 @@ end ### tri products -function X_A_Xt(a::PDMat, x::AbstractMatrix{<:Real}) +function X_A_Xt(a::PDMatCholesky, x::AbstractMatrix{<:Real}) @check_argdims a.dim == size(x, 2) - z = x * chol_lower(a.chol) + z = x * chol_lower(cholesky(a)) return Symmetric(z * transpose(z)) end -function Xt_A_X(a::PDMat, x::AbstractMatrix{<:Real}) +function Xt_A_X(a::PDMatCholesky, x::AbstractMatrix{<:Real}) @check_argdims a.dim == size(x, 1) - z = chol_upper(a.chol) * x + z = chol_upper(cholesky(a)) * x return Symmetric(transpose(z) * z) end -function X_invA_Xt(a::PDMat, x::AbstractMatrix{<:Real}) +function X_invA_Xt(a::PDMatCholesky, x::AbstractMatrix{<:Real}) @check_argdims a.dim == size(x, 2) - z = x / chol_upper(a.chol) + z = x / chol_upper(cholesky(a)) return Symmetric(z * transpose(z)) end -function Xt_invA_X(a::PDMat, x::AbstractMatrix{<:Real}) +function Xt_invA_X(a::PDMatCholesky, x::AbstractMatrix{<:Real}) @check_argdims a.dim == size(x, 1) - z = chol_lower(a.chol) \ x + z = chol_lower(cholesky(a)) \ x return Symmetric(transpose(z) * z) end ### Specializations for `Array` arguments with reduced allocations -function quad(a::PDMat{<:Real,<:Vector}, x::Matrix) +function quad(a::PDMatCholesky{T,Matrix{T}}, x::Matrix) where {T<:Real} @check_argdims a.dim == size(x, 1) T = typeof(zero(eltype(a)) * abs2(zero(eltype(x)))) return quad!(Vector{T}(undef, size(x, 2)), a, x) end -function invquad(a::PDMat{<:Real,<:Vector}, x::Matrix) +function invquad(a::PDMatCholesky{T,Matrix{T}}, x::Matrix) where {T<:Real} @check_argdims a.dim == size(x, 1) T = typeof(abs2(zero(eltype(x))) / zero(eltype(a))) return invquad!(Vector{T}(undef, size(x, 2)), a, x) diff --git a/test/pdmtypes.jl b/test/pdmtypes.jl index f20b672..aa719e3 100644 --- a/test/pdmtypes.jl +++ b/test/pdmtypes.jl @@ -69,7 +69,7 @@ using Test @test B == A @test (B === A) === (S === T) @test (B.mat === A.mat) === (S === T) - @test (B.chol === A.chol) === (S === T) + @test (B.fact === A.fact) === (S === T) end A = PDiagMat(ones(T, 2)) @@ -100,7 +100,7 @@ using Test @test (B.mat === A.mat) === (S === T) # CholMOD only supports Float64 and ComplexF64 type parameters! # Hence the Cholesky factorization is reused - @test B.chol === A.chol + @test B.fact === A.fact end end end From 6ce4b7795040da8a23de2a96e55b9ab479102b05 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Tue, 21 Nov 2023 23:20:02 +0100 Subject: [PATCH 09/14] Fix type parameter --- src/pdmat.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/pdmat.jl b/src/pdmat.jl index ecd945a..57774ae 100644 --- a/src/pdmat.jl +++ b/src/pdmat.jl @@ -5,7 +5,7 @@ struct PDMat{T<:Real,S<:AbstractMatrix{T},F<:Factorization} <: AbstractPDMat{T} mat::S fact::F - function PDMat{T,S,C}(mat::S, fact::F) where {T,S<:AbstractMatrix{T},F<:Factorization} + function PDMat{T,S,F}(mat::S, fact::F) where {T,S<:AbstractMatrix{T},F<:Factorization} d = LinearAlgebra.checksquare(mat) if size(fact) != (d, d) throw(DimensionMismatch("Dimensions of the matrix and the factorization are inconsistent.")) From bbedcb2fca14e47d02bf9ec30028fad16a41dfb0 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Tue, 21 Nov 2023 23:24:13 +0100 Subject: [PATCH 10/14] Fix variable conflicts with type parameters --- src/pdmat.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/pdmat.jl b/src/pdmat.jl index 57774ae..31d9d94 100644 --- a/src/pdmat.jl +++ b/src/pdmat.jl @@ -191,13 +191,13 @@ end function quad(a::PDMatCholesky{T,Matrix{T}}, x::Matrix) where {T<:Real} @check_argdims a.dim == size(x, 1) - T = typeof(zero(eltype(a)) * abs2(zero(eltype(x)))) - return quad!(Vector{T}(undef, size(x, 2)), a, x) + S = typeof(zero(T) * abs2(zero(eltype(x)))) + return quad!(Vector{S}(undef, size(x, 2)), a, x) end function invquad(a::PDMatCholesky{T,Matrix{T}}, x::Matrix) where {T<:Real} @check_argdims a.dim == size(x, 1) - T = typeof(abs2(zero(eltype(x))) / zero(eltype(a))) - return invquad!(Vector{T}(undef, size(x, 2)), a, x) + S = typeof(abs2(zero(eltype(x))) / zero(T)) + return invquad!(Vector{S}(undef, size(x, 2)), a, x) end From dca927383bd43b266ab76e44c9619a567fc9bee1 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Wed, 22 Nov 2023 00:35:48 +0100 Subject: [PATCH 11/14] Clean up switch to `fact` --- src/pdmat.jl | 2 +- test/pdmtypes.jl | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/src/pdmat.jl b/src/pdmat.jl index 31d9d94..484f40c 100644 --- a/src/pdmat.jl +++ b/src/pdmat.jl @@ -48,7 +48,7 @@ Base.convert(::Type{AbstractPDMat{T}}, a::PDMat) where {T<:Real} = convert(PDMat Base.size(a::PDMat) = (a.dim, a.dim) Base.Matrix(a::PDMat) = Matrix(a.mat) LinearAlgebra.diag(a::PDMat) = diag(a.mat) -LinearAlgebra.cholesky(a::PDMatCholesky) = a.chol +LinearAlgebra.cholesky(a::PDMatCholesky) = a.fact ### Inheriting from AbstractMatrix diff --git a/test/pdmtypes.jl b/test/pdmtypes.jl index aa719e3..f6decef 100644 --- a/test/pdmtypes.jl +++ b/test/pdmtypes.jl @@ -214,7 +214,7 @@ using Test for dim in (1, 5, 10) x = rand(dim, dim) M = PDMat(x' * x + I) - @test fieldnames(typeof(M)) == (:mat, :chol) + @test fieldnames(typeof(M)) == (:mat, :fact) @test propertynames(M) == (fieldnames(typeof(M))..., :dim) @test getproperty(M, :dim) === dim for p in fieldnames(typeof(M)) @@ -239,7 +239,7 @@ using Test if HAVE_CHOLMOD x = sprand(dim, dim, 0.2) M = PDMat(x' * x + I) - @test fieldnames(typeof(M)) == (:mat, :chol) + @test fieldnames(typeof(M)) == (:mat, :fact) @test propertynames(M) == (fieldnames(typeof(M))..., :dim) @test getproperty(M, :dim) === dim for p in fieldnames(typeof(M)) From 33a80f3b06eff4853f1ba8d01abe48789c4ae867 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Wed, 22 Nov 2023 01:38:17 +0100 Subject: [PATCH 12/14] Specialize `sqrt` for `PDSparseMat` --- ext/PDMatsSparseArraysExt/pdsparsemat.jl | 1 + src/pdmat.jl | 2 +- 2 files changed, 2 insertions(+), 1 deletion(-) diff --git a/ext/PDMatsSparseArraysExt/pdsparsemat.jl b/ext/PDMatsSparseArraysExt/pdsparsemat.jl index f6bef3f..129a33b 100644 --- a/ext/PDMatsSparseArraysExt/pdsparsemat.jl +++ b/ext/PDMatsSparseArraysExt/pdsparsemat.jl @@ -31,6 +31,7 @@ Base.:/(x::AbstractVecOrMat{T}, a::PDSparseMat{T}) where {T<:Real} = convert(Arr Base.inv(a::PDSparseMat) = PDMat(inv(a.mat)) LinearAlgebra.cholesky(a::PDSparseMat) = a.fact +Base.sqrt(A::PDSparseMat) = PDMat(sqrt(Hermitian(Matrix(A)))) ### whiten and unwhiten diff --git a/src/pdmat.jl b/src/pdmat.jl index 484f40c..82d63c2 100644 --- a/src/pdmat.jl +++ b/src/pdmat.jl @@ -83,7 +83,7 @@ LinearAlgebra.det(a::PDMat) = det(a.fact) LinearAlgebra.logdet(a::PDMat) = logdet(a.fact) LinearAlgebra.eigmax(a::PDMat) = eigmax(a.mat) LinearAlgebra.eigmin(a::PDMat) = eigmin(a.mat) -LinearAlgebra.sqrt(A::PDMat) = PDMat(sqrt(Hermitian(A.mat))) +Base.sqrt(A::PDMat) = PDMat(sqrt(Hermitian(A.mat))) function Base.kron(A::PDMatCholesky, B::PDMatCholesky) return PDMat(kron(A.mat, B.mat), Cholesky(kron(A.fact.U, B.fact.U), 'U', A.fact.info)) From a0acc4912c738295708dad62a8fcef68b46aebbc Mon Sep 17 00:00:00 2001 From: David Widmann Date: Mon, 5 Aug 2024 12:17:20 +0200 Subject: [PATCH 13/14] Add test from #207 --- test/pdmtypes.jl | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/test/pdmtypes.jl b/test/pdmtypes.jl index 9331671..ab4bac1 100644 --- a/test/pdmtypes.jl +++ b/test/pdmtypes.jl @@ -304,4 +304,12 @@ using Test @test_broken C - B isa Diagonal{Float64, Vector{Float64}} @test C - B ≈ Matrix(C) - Matrix(B) end + + # https://github.com/JuliaStats/PDMats.jl/pull/207 + @testset "PDMat from SymTridiagonal" begin + S = SymTridiagonal(fill(4, 4), fill(1, 3)) + M = @inferred(PDMat(S)) + @test M isa PDMat{Int,<:SymTridiagonal,<:Cholesky} + @test M == S + end end From 8c160df742fe71fed6f36be8a0b3bdf489ce6279 Mon Sep 17 00:00:00 2001 From: David Widmann Date: Mon, 5 Aug 2024 13:37:55 +0200 Subject: [PATCH 14/14] Run `SymTridiagonal` tests only on Julia >= 1.8 --- test/pdmtypes.jl | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/test/pdmtypes.jl b/test/pdmtypes.jl index ab4bac1..2d44b80 100644 --- a/test/pdmtypes.jl +++ b/test/pdmtypes.jl @@ -305,11 +305,14 @@ using Test @test C - B ≈ Matrix(C) - Matrix(B) end - # https://github.com/JuliaStats/PDMats.jl/pull/207 - @testset "PDMat from SymTridiagonal" begin - S = SymTridiagonal(fill(4, 4), fill(1, 3)) - M = @inferred(PDMat(S)) - @test M isa PDMat{Int,<:SymTridiagonal,<:Cholesky} - @test M == S + # Ref https://github.com/JuliaStats/PDMats.jl/pull/207 + # Note: `cholesky(::SymTridiagonal)` requires https://github.com/JuliaLang/julia/pull/44076 + if VERSION >= v"1.8.0-DEV.1526" + @testset "PDMat from SymTridiagonal" begin + S = SymTridiagonal(fill(4, 4), fill(1, 3)) + M = @inferred(PDMat(S)) + @test M isa PDMat{Int,<:SymTridiagonal,<:Cholesky} + @test M == S + end end end