PDMatsExtras

This is a package for extra Positive (Semi-) Definated Matrix types. It is an extension to PDMats.jl.

It probably wouldn't exist, except Distributions.jl is currently very tied to the idea that the type of a covariance matrix should subtype AbstractPDMat. There is an issue open to change that. When that is resolve the matrix defined here may well move elsewhere, or cease to be required.

The Matrixes

PSDMat

A Positive Semi-Definite Matrix. It still subtypes AbstractPDMat. It's not quite as nice to work with as a truely positive definite matrix, since the math doesn't work out so well. But this is able to represent all covariences -- which must be positive semi-definate. You might not like the consequences,

julia> using LinearAlgebra, PDMatsExtras

julia> X = Float64[
               10   -9   2   4
               -9    8  -1  -4
                2   -1   1   1
                4   -4   1   6
           ];

julia> isposdef(X)
false

julia> PSDMat(X)
4×4 PSDMat{Float64, Matrix{Float64}}:
 10.0  -9.0   2.0   4.0
 -9.0   8.0  -1.0  -4.0
  2.0  -1.0   1.0   1.0
  4.0  -4.0   1.0   6.0


julia> # can also construct from a pivoted cholesky, even one of a rank deficient matrix (like this one)

julia> PSDMat(cholesky(X, Val(true); check=false))
4×4 PSDMat{Float64, Matrix{Float64}}:
 10.0  -9.0       2.0   4.0
 -9.0   9.26923  -1.0  -4.0
  2.0  -1.0       1.0   1.0
  4.0  -4.0       1.0   6.0

WoodburyPDMat

It is a positive definite Woodbury matrix. This is a special case of the Symmetric Woodbury Matrix (see WoodburyMatrices.jl's SymWoodbury type) which is given by A*D*A' + S for S and D being diagonal, which has the additional requirement that the diagonal matrices are also non-negative.

julia> using LinearAlgebra, PDMatsExtras

julia> A = Float64[
         2.0  2.0  -8.0   5.0  -1.0   2.0   6.0
         2.0  7.0  -1.0  -5.0  -4.0   8.0   7.0
        -2.0  9.0  -9.0  -5.0   9.0  -5.0  -3.0
         3.0  4.0  -6.0  -4.0   3.0  -3.0  -3.0
       ];

julia> D = Diagonal(Float64[1, 2, 3, 2, 2, 1, 5]);

julia> S = Diagonal(Float64[4, 2, 3, 6]);

julia> W = WoodburyPDMat(A, D, S)
4×4 WoodburyPDMat{Float64,Array{Float64,2},Diagonal{Float64,Array{Float64,1}},Diagonal{Float64,Array{Float64,1}}}:
 444.0  240.0   80.0   24.0
 240.0  498.0  -18.0  -33.0
  80.0  -18.0  694.0  382.0
  24.0  -33.0  382.0  259.0

julia> A*D*A' + S
4×4 Array{Float64,2}:
 444.0  240.0   80.0   24.0
 240.0  498.0  -18.0  -33.0
  80.0  -18.0  694.0  382.0
  24.0  -33.0  382.0  259.0