inv(Σ) fails with StaticArrays

[3f6a]

Why this fails?

μ = @SVector randn(3)
Σ = @SMatrix randn(3,3)
Σ = Σ * Σ'
@assert Σ == Σ'
inv(MvNormal(μ, Σ).Σ)
# PosDefException: matrix is not Hermitian; Cholesky factorization failed.

If I do the same with normal arrays (non-static),

μ = randn(3)
Σ = randn(3,3)
Σ = Σ * Σ'
@assert Σ == Σ'
inv(MvNormal(μ, Σ).Σ) # works!

Therefore the issue is related to StaticArrays.

[andreasnoack]

I believe this issue belongs in StaticArrays. The inverse of the Cholesky here isn't symmetric

julia> issymmetric(MvNormal(μ, Σ).Σ.mat)
true

julia> issymmetric(inv(MvNormal(μ, Σ).Σ.chol))
false

This is probably because it's computed as two solves in https://github.com/JuliaArrays/StaticArrays.jl/blob/5913f330ac4834ef1521bc7faa626a0be1d6bb7a/src/cholesky.jl#L72C1-L74. Two solves is probably also more expensive than needed. A single solve followed by a symmetric rank-k update should be more efficient and ensure symmetry.

[devmotion]

AFAICT this issue (and the suggested fix) are already discussed in JuliaArrays/StaticArrays.jl#1218.

[astrozot]

I believe this issue is not relevant anymore and should be closed as the latest version of StaticArrays fixes this (see also JuliaArrays/StaticArrays.jl#1257).

[3f6a]

@astrozot But JuliaArrays/StaticArrays.jl#1218 seems to be still open.

[astrozot]

True @3f6a, I think this has just been overlooked. I added a comment JuliaArrays/StaticArrays.jl#1218 and hopefully some of the developers can verify and close that issue.

[devmotion]

I'm not sure if the issue should be considered solved either - as noted above, instead of averaging out the errors the implementation in StaticArrays should probably just be changed to a single solve.

[devmotion]

Ref JuliaArrays/StaticArrays.jl#1326

[3f6a]

@devmotion Could you please leave a quick comment clarifying why you are closing this ?
Is it because the actual issue should be tracked upstream in StaticArrays?
Or has this actually been fixed (e.g. I see JuliaArrays/StaticArrays.jl#1326 was merged a couple of days ago)? Thanks !

[devmotion]

The underlying issue (non-symmetric matrices returned by inv(::Cholesky{<:Real,<:StaticMatrix})) was fixed upstream in StaticArrays, initially in JuliaArrays/StaticArrays.jl#1258 (released in StaticArrays@1.9.15) and in an improved way in JuliaArrays/StaticArrays.jl#1326 (released in StaticArrays@1.9.16).