Apply a few fixes and improvements from PDMats.jl

Project.toml

@@ -14,7 +14,7 @@ SuiteSparse = "4607b0f0-06f3-5cda-b6b1-a6196a1729e9"
 ChainRulesCore = "0.9.17, 0.10, 1"
 Distributions = "0.23, 0.24"
 FiniteDifferences = "0.11, 0.12"
-PDMats = "0.9, 0.10, 0.11"
+PDMats = "0.11.19"
 Zygote = "0.5.5, 0.6"
 julia = "1"
 

src/PDMatsExtras.jl

@@ -13,4 +13,6 @@ include("psd_mat.jl")
 include("woodbury_pd_mat.jl")
 include("utils.jl")
 
+Base.@deprecate PSDMat{T,S}(d::Int, m::AbstractMatrix{T}, c::CholType{T,S}) where {T,S} PSDMat{T,S}(m, c)
+
 end

src/psd_mat.jl

@@ -17,27 +17,40 @@ References for discussion on supporting degenerate mvnormal distributions:
 Positive semi-definite matrix together with a CholeskyPivoted factorization object.
 """
 struct PSDMat{T<:Real,S<:AbstractMatrix} <: AbstractPDMat{T}
- dim::Int
  mat::S
  chol::CholType{T, S}
 
- PSDMat{T,S}(d::Int,m::AbstractMatrix{T},c::CholType{T,S}) where {T,S} = new{T, S}(d, m, c)
+ PSDMat{T,S}(m::AbstractMatrix{T},c::CholType{T,S}) where {T,S} = new{T, S}(m, c)
 end
 
 function PSDMat(mat::AbstractMatrix, chol::CholType)
- d = size(mat, 1)
+ d = LinearAlgebra.checksquare(mat)
  size(chol, 1) == d ||
  throw(DimensionMismatch("Dimensions of mat and chol are inconsistent."))
- PSDMat{eltype(mat),typeof(mat)}(d, mat, chol)
+ PSDMat{eltype(mat),typeof(mat)}(mat, chol)
 end
 
-PSDMat(mat::Matrix) = PSDMat(mat, cholesky(mat, Val(true); check=false))
+PSDMat(mat::Matrix) = PSDMat(mat, cholesky(mat, VERSION >= v"1.8.0-rc1" ? RowMaximum() : Val(true); check=false))
 PSDMat(mat::Symmetric) = PSDMat(Matrix(mat))
 PSDMat(fac::CholType) = PSDMat(Matrix(fac), fac)
 
+function Base.getproperty(a::PSDMat, s::Symbol)
+ if s === :dim
+ return size(getfield(a, :mat), 1)
+ end
+ return getfield(a, s)
+end
+Base.propertynames(::PSDMat) = (:mat, :chol, :dim)
+
 ### Conversion
-Base.convert(::Type{PSDMat{T}}, a::PSDMat) where {T<:Real} = PSDMat(convert(AbstractArray{T}, a.mat))
-Base.convert(::Type{AbstractArray{T}}, a::PSDMat) where {T<:Real} = convert(PSDMat{T}, a)
+Base.convert(::Type{PSDMat{T}}, a::PSDMat{T}) where {T<:Real} = a
+function Base.convert(::Type{PSDMat{T}}, a::PSDMat) where {T<:Real}
+ chol = convert(CholType{T}, a.chol)
+ S = typeof(chol.factors)
+ mat = convert(S, a.mat)
+ PSDMat{T,S}(mat, chol)
+end
+Base.convert(::Type{AbstractPDMat{T}}, a::PSDMat) where {T<:Real} = convert(PSDMat{T}, a)
 
 ### Basics
 
@@ -63,6 +76,7 @@ end
 ### Algebra
 
 Base.inv(a::PSDMat) = PSDMat(inv(a.chol))
+LinearAlgebra.cholesky(a::PSDMat) = a.chol
 LinearAlgebra.logdet(a::PSDMat) = logdet(a.chol)
 LinearAlgebra.eigmax(a::PSDMat) = eigmax(a.mat)
 LinearAlgebra.eigmin(a::PSDMat) = eigmin(a.mat)
@@ -92,8 +106,26 @@ end
 
 ### quadratic forms
 
-PDMats.quad(a::PSDMat, x::StridedVector) = dot(x, a * x)
-PDMats.invquad(a::PSDMat, x::StridedVector) = dot(x, a \ x)
+function PDMats.quad(a::PSDMat, x::AbstractVecOrMat)
+ if a.dim != size(x, 1)
+ throw(DimensionMismatch("Inconsistent argument dimensions."))
+ end
+ # https://github.com/JuliaLang/julia/commit/2425ae760fb5151c5c7dd0554e87c5fc9e24de73
+ if VERSION < v"1.4.0-DEV.92"
+ z = a.mat * x
+ return x isa AbstractVector ? dot(x, z) : map(dot, eachcol(x), eachcol(z))
+ else
+ return x isa AbstractVector ? dot(x, a.mat, x) : map(Base.Fix1(quad, a), eachcol(x))
+ end
+end
+
+function PDMats.invquad(a::PSDMat, x::AbstractVecOrMat)
+ if a.dim != size(x, 1)
+ throw(DimensionMismatch("Inconsistent argument dimensions."))
+ end
+ z = a.chol \ x
+ return x isa AbstractVector ? dot(x, z) : map(dot, eachcol(x), eachcol(z))
+end
 
 """
  quad!(r::AbstractArray, a::AbstractPDMat, x::StridedMatrix)

src/utils.jl

@@ -14,3 +14,8 @@ function submat(A::ScalMat, inds)
  checkbounds(Bool, A, inds) || throw(BoundsError(A, inds))
  return ScalMat(length(inds), A.value)
 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

test/psd_mat.jl

@@ -18,7 +18,7 @@ end
  verbose = 1
  @testset "Positive definite" begin
  M = TEST_MATRICES["Positive definite"]
- pivoted = cholesky(M, Val(true))
+ pivoted = cholesky(M, VERSION >= v"1.8.0-rc1" ? RowMaximum() : Val(true))
  C = PSDMat(M, pivoted)
  test_pdmat(
  C,
@@ -33,7 +33,7 @@ end
  @testset "Positive semi-definite" begin
  M = TEST_MATRICES["Positive semi-definite"]
  @test !isposdef(M)
- pivoted = cholesky(M, Val(true); check=false)
+ pivoted = cholesky(M, VERSION >= v"1.8.0-rc1" ? RowMaximum() : Val(true); check=false)
  C = PSDMat(M, pivoted)
  test_pdmat(
  C,

test/utils.jl

@@ -40,7 +40,7 @@
 
  @testset "PSDMat" begin
  SM = TEST_MATRICES["Positive semi-definite"]
- pivoted = cholesky(SM, Val(true); check=false)
+ pivoted = cholesky(SM, VERSION >= v"1.8.0-rc1" ? RowMaximum() : Val(true); check=false)
  M = PSDMat(SM, pivoted)
  M_dense = Matrix(M)