inverse of SMatrix sometimes fail with PDMat

I suspect that this bug is because the cholesky decomposition fills the lower-triangular elements with zeros (as mentioned in @194, and from here). I might be wrong.

If StaticArrays.cholesky is used to construct a PDMat object, then inv on this object can fail, when it needs to get the cholesky decomposition of the inverse, because the inverse can be slightly non-Hermitian.

Here is a small example.

julia> using LinearAlgebra, PDMats, StaticArrays

julia> J = [0.326392342118349 0.022376926902038786; 0.022376926902038786 5.002486325211338];

julia> inv(PDMat(MMatrix{2,2}(J)))
ERROR: PosDefException: matrix is not Hermitian; Cholesky factorization failed.
Stacktrace:
 [1] non_hermitian_error()
   @ StaticArrays ~/.julia/packages/StaticArrays/cZ1ET/src/cholesky.jl:2
 [2] #cholesky#532
   @ ~/.julia/packages/StaticArrays/cZ1ET/src/cholesky.jl:4 [inlined]
 [3] cholesky
   @ ~/.julia/packages/StaticArrays/cZ1ET/src/cholesky.jl:3 [inlined]
 [4] PDMat(mat::SMatrix{2, 2, Float64, 4})
   @ PDMats ~/.julia/packages/PDMats/bzppG/src/pdmat.jl:19
 [5] inv(a::PDMat{Float64, MMatrix{2, 2, Float64, 4}})
   @ PDMats ~/.julia/packages/PDMats/bzppG/src/pdmat.jl:81
 [6] top-level scope
   @ REPL[7]:1

yet J is very well behaved. There is no error if the values of J are slightly perturbed, or with a basic inverse:

julia> inv(PDMat(MMatrix{2,2}(J .+ [0 0; 0 -8e-15])))
2×2 PDMat{Float64, SMatrix{2, 2, Float64, 4}}:
  3.06474    -0.0137091
 -0.0137091   0.199962

julia> inv(J)
2×2 Matrix{Float64}:
  3.06474    -0.0137091
 -0.0137091   0.199962

Also, there is no error if we don't combine MMatrix, PDMat, and getting the inverse of that, as a PDMat:

julia> inv(MMatrix{2,2}(J))
2×2 MMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
  3.06474    -0.0137091
 -0.0137091   0.199962

julia> inv(cholesky(MMatrix{2,2}(J)))
2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
  3.06474    -0.0137091
 -0.0137091   0.199962

julia> inv(PDMat(J))
2×2 PDMat{Float64, Matrix{Float64}}:
 
 3.06474    -0.0137091
 -0.0137091   0.199962

The problem may come from the fact that the cholesky decomposition of a StaticArray stores a non-hermitian matrix:

julia> Jbad  = PDMat(MMatrix{2,2}(J)); Jbad.chol.factors # non-symmetric: 0 on lower off-diagonal
2×2 MMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 0.571308  0.0391679
 0.0       2.23628

julia> Jgood = PDMat(J); Jgood.chol.factors # not symmetric either actually!
2×2 Matrix{Float64}:
 0.571308   0.0391679
 0.0223769  2.23628

In any case, the problem above is because the inverse of J, computed using its cholesky decomposition, is seen as non-hermitian:

julia> Jbadinv = inv(Jbad.chol)
2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
  3.06474    -0.0137091
 -0.0137091   0.199962

julia> Jbadinv[2,1] == Jbadinv[1,2] # false: yet required by ishermitian
false

julia> ishermitian(Jbadinv)
false

julia> PDMat(Jbadinv)
ERROR: PosDefException: matrix is not Hermitian; Cholesky factorization failed.
...

If we slightly perturb the matrix, rounding errors don't appear, and the inverse is recognized as being hermitian:

julia> Jokay = PDMat(MMatrix{2,2}(J .+ [0 0; 0 -8e-15])); Jokay.chol.factors
2×2 MMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 0.571308  0.0391679
 0.0       2.23628

julia> Jokayinv = inv(Jokay.chol)
2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
  3.06474    -0.0137091
 -0.0137091   0.199962

julia> Jokayinv[2,1] == Jokayinv[1,2]
true

julia> ishermitian(Jokayinv)
true

julia> PDMat(Jokayinv); # okay: no error

This issue comes up frequently in what I do, actually. For now, I don't see any other workaround, other than giving up on using StaticArrays unfortunately. I would love to hear a work-around until this issue can be fixed.

Reproduce the PDMat inverse failure from the Julia examples, then inspect src/cholesky.jl at the reported factorization path and the PDMats pdmat.jl call sites in the stack trace. Compare the MMatrix and Matrix cases, including the reported cholesky factors and inverse symmetry. Done means the supplied example no longer raises a non-Hermitian PosDefException.