inverse of SMatrix sometimes fail with PDMat

[cecileane]

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.

[bstkj]

I think this problem occurs because of how inv is defined for Cholesky decompositions of StaticArrays.jl types:

function inv(A::Cholesky{T,<:StaticMatrix{N,N,T}}) where {N,T}
    return A.U \ (A.U' \ SDiagonal{N}(I))
end

This sequence of operations to invert a positive-definite matrix does not guarantee that the inverse is also symmetric in the face of rounding errors. Using the same example as above:

julia> J = [0.326392342118349 0.022376926902038786; 0.022376926902038786 5.002486325211338]
2×2 Matrix{Float64}:
 0.326392   0.0223769
 0.0223769  5.00249

julia> Jchol = cholesky(J)
Cholesky{Float64, Matrix{Float64}}
U factor:
2×2 UpperTriangular{Float64, Matrix{Float64}}:
 0.571308  0.0391679
  ⋅        2.23628

julia> (Jchol.U \ (Jchol.U' \ I))[1,2]
-0.013709063214013707

julia> (Jchol.U \ (Jchol.U' \ I))[2,1]
-0.013709063214013708 # the off-diagonals should be the same!

Note that preserving symmetry after inversion is not an issue for Jchol::Cholesky{Float64, Matrix{Float64}} since this operation defaults to:

inv!(C::Cholesky{<:BlasFloat,<:StridedMatrix}) =
    copytri!(LAPACK.potri!(C.uplo, C.factors), C.uplo, true)

which is probably ensures that the return value is symmetric.

julia> inv(Jchol)[1,2]
-0.013709063214013708

julia> inv(Jchol)[2,1]
-0.013709063214013708

However, converting J to an MMatrix leads to the former inversion method being called. This then leads to downstream errors if the inverse is taken after MMatrix{2,2}(J) is wrapped in a PDMat instance, because PDMats.jl defaults to the inv method of the underlying matrix, then attempts to convert it back to PDMat.

inv(MMatrix{2,2}(J)) may not be symmetric, leading to the "not Hermitian" error for inv(PDMat(MMatrix{2,2}(J))).

I guess the inv method for Cholesky decomposition within StaticArrays.jl should be rewritten if we want to guarantee that the returned inverse is symmetric?