`X_A_Xt` symmetric output breaks addition and subraction typing with static PDMats

[btmit]

using StaticArrays, PDMats

C = SMatrix{2,3}([0.366247  0.314705  0.0352938
        0.851331  0.495017  0.282046])
x = PDMat(SMatrix{2,2}([0.191844  0.183975
        0.183975  0.200666]))
y = PDMat(SMatrix{3,3}([0.0487753  0.0529655  0.0757205
        0.0529655  0.503523   0.478115
        0.0757205  0.478115   0.480115]))

x2 = x - X_A_Xt(y, C)

typeof(x2)  # Matrix{Float64}

Before the commit on 10/17, the output was type SMatrix{2,2,Float64} which I think is correct.

[andreasnoack]

I don't think the problem here is X_A_Xt but that (-)(::PDMat{<:SMatrix}, PDMat[<:SMatrix}) returns a Matrix. You can get the same problematic behavior just with x, i.e.

julia> x - x
2×2 Matrix{Float64}:
 0.0  0.0
 0.0  0.0

[btmit]

I'm not sure that's exactly this problem (but it may be another one). Only one of these two is a PDMat, not both. Though perhaps fixing the issue you highlight would fix the dispatch for this one too.

As a reference, both of the following cases "work".

typeof( Matrix(x) - X_A_Xt(y, C) )  # SMatrix - Symmetric{SMatrix} => SMatrix
typeof( x - X_A_Xt(y, C).data )  # PDMat{SMatrix} - SMatrix => SMatrix

[devmotion]

I think the first example

typeof( Matrix(x) - X_A_Xt(y, C) )  # SMatrix - Symmetric{SMatrix} => SMatrix

should return a Matrix. Seems wrong to get a SMatrix when a Matrix was requested, we should fix

PDMats.jl/src/pdmat.jl

Line 45 in ca14884

Base.Matrix(a::PDMat) = copy(a.mat)

Edit: This bug is already fixed by #188 but I opened a separate PR #196 since it seems this fix could be merged and released independent of the major changes in #188.

[btmit]

Thank you!

In general is there a function that returns whatever matrix backs an AbstractPDMat?

[devmotion]

No. Since AbstractPDMat <: AbstractMatrix, AbstractMatrix(A::AbstractPDMat) = A. AbstractPDMat are also not necessarily wrappers of other matrices, e.g., ScalMat only stores two scalar fields and PDiagMat a vector of the diagonal entries.

[btmit]

Understood.

Any suggestions or workarounds on the original issue?

[devmotion]

I think we can solve this issue, maybe not in all generality but at least for your example.

It would already be an improvement to define

Base.broadcastable(x::PDMat) = Base.broadcastable(x.mat)

Many generic methods (such as -(::AbstractArray, ::AbstractArray)) fall back to broadcasting. With this definition one gets

julia> x - x
2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 0.0  0.0
 0.0  0.0

However, it's not sufficient to fix the original example:

julia> x - X_A_Xt(y, C)
2×2 SizedMatrix{2, 2, Float64, 2, Matrix{Float64}} with indices SOneTo(2)×SOneTo(2):
 0.110047      0.000865491
 0.000865491  -0.210773

This is not a PDMats problem though, the same behaviour occurs without PDMats as well:

julia> @SMatrix(zeros(2, 2)) .- Symmetric(@SMatrix(zeros(2, 2)))
2×2 SizedMatrix{2, 2, Float64, 2, Matrix{Float64}} with indices SOneTo(2)×SOneTo(2):
 0.0  0.0
 0.0  0.0

AFAICT the problem is a missing definition

Base.BroadcastStyle(::Type{<:Symmetric{<:Any, <:StaticArray}}) = StaticArrays.StaticArrayStyle{2}()

in https://github.com/JuliaArrays/StaticArrays.jl/blob/72d2bd3538235c9162f630a5130112b83eaa0af7/src/broadcast.jl#L12-L13 (BTW should probably add Hermitian as well). With this definition your example gives

julia> x - X_A_Xt(y, C)
2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 0.110047      0.000865491
 0.000865491  -0.210773

[btmit]

I can confirm that #197 plus Base.BroadcastStyle(::Type{<:Symmetric{<:Any, <:StaticArray}}) = StaticArrays.StaticArrayStyle{2}() fixes this issue when x is either PDMat or PDiagMat.

Would you like me to open an issue in StaticArrays for Symmetric and Hermitian then link it here?

[devmotion]

I opened JuliaArrays/StaticArrays.jl#1220

[devmotion]

The issue is fixed on PDMats 0.11.30 and StaticArrays 1.7.0:

julia> using StaticArrays, PDMats

julia> C = SMatrix{2,3}([0.366247  0.314705  0.0352938
               0.851331  0.495017  0.282046])
2×3 SMatrix{2, 3, Float64, 6} with indices SOneTo(2)×SOneTo(3):
 0.366247  0.314705  0.0352938
 0.851331  0.495017  0.282046

julia> x = PDMat(SMatrix{2,2}([0.191844  0.183975
               0.183975  0.200666]))
2×2 PDMat{Float64, SMatrix{2, 2, Float64, 4}}:
 0.191844  0.183975
 0.183975  0.200666

julia> y = PDMat(SMatrix{3,3}([0.0487753  0.0529655  0.0757205
               0.0529655  0.503523   0.478115
               0.0757205  0.478115   0.480115]))
3×3 PDMat{Float64, SMatrix{3, 3, Float64, 9}}:
 0.0487753  0.0529655  0.0757205
 0.0529655  0.503523   0.478115
 0.0757205  0.478115   0.480115

julia> x2 = x - X_A_Xt(y, C)
2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 0.110047      0.000865491
 0.000865491  -0.210773

julia> typeof(x2)
SMatrix{2, 2, Float64, 4} (alias for SArray{Tuple{2, 2}, Float64, 2, 4})

[btmit]

I think the first example

typeof( Matrix(x) - X_A_Xt(y, C) )  # SMatrix - Symmetric{SMatrix} => SMatrix

should return a Matrix. Seems wrong to get a SMatrix when a Matrix was requested, we should fix

PDMats.jl/src/pdmat.jl

Line 45 in ca14884

Base.Matrix(a::PDMat) = copy(a.mat)

Edit: This bug is already fixed by #188 but I opened a separate PR #196 since it seems this fix could be merged and released independent of the major changes in #188.

I just noticed that PR #196 breaks the cov(::MvNormal) method in Distributions here. Is that an issue that should be raised in this package or Distributions?

smat = SMatrix{2,2}([0.191844  0.183975
               0.183975  0.200666])
x = PDMat(smat)
mv = MvNormal(x)

typeof(cov(mv)) == typeof(smat)  # false, it's a Matrix

[devmotion]

#196 fixed a bug, so IMO that's an issue with Distributions. IMO we should just merge JuliaStats/Distributions.jl#1373 - I was always hesitant to merge it since it changes the return type of cov but since AbstractPDMats are also AbstractArrays I think one should not expect any problems - and now it would even be a bugfix 😄