`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 https://github.com/JuliaStats/PDMats.jl/blob/ca14884f526e149b22d6749d80b22cb0485a8113/src/pdmat.jl#L45

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 https://github.com/JuliaArrays/StaticArrays.jl/pull/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

https://github.com/JuliaStats/PDMats.jl/blob/ca14884f526e149b22d6749d80b22cb0485a8113/src/pdmat.jl#L45

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 https://github.com/JuliaStats/Distributions.jl/pull/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 😄

[btmit]

Selfishly I would love to see that merged. It would definitely break some of our code, but would be a big improvement overall.

I'm guessing we might still need some function that would return the array that underlies the covariance (e.g. PDMat(::SMatrix) -> SMatrix, PDiagMat(::Vector) -> Diagonal, ScalMat -> Diagonal(Fill)? etc to help folks fix the issues that will arise. Alternatively, I could fix our issues if this function existed even without JuliaStats/Distributions.jl#1373

[btmit]

JuliaStats/Distributions.jl#572 had a nice discussion of whether AbstractPDMat should be an internal representation or something the user deals with. If there was some conversion method then cov() could preserve type without forcing AbstractPDMats on every application. This would probably be more palatable in general (and less breaking) than JuliaStats/Distributions.jl#1373, but I think either approach is better than cov returning a Matrix in all cases.

[devmotion]

Apart from the last comment, all the discussion in the issue happened before AbstractPDMat became a subtype of AbstractMatrix. I agree with https://github.com/JuliaStats/Distributions.jl/issues/572#issuecomment-277876750 that it should be fine to return structured AbstractArrays.