Non-symmetric inverse of Cholesky decomposition

[astrozot]

Because of rounding errors, in some cases explicitly symmetric matrices might produce non-symmetric inverses:

using StaticArrays, LinearAlgebra

c = cholesky(@SMatrix[π 3.1; 3.1 12.0])
d = inv(c)
d[1,2] - d[2,1] # returns -1.3877787807814457e-17

This can be a serious issue and it can break various codes: for example, cholesky(d) would fail. Note that the same code works without problems if one uses standard arrays instead of static ones:

c1 = cholesky([π 3.1; 3.1 12.0])
d1 = inv(c1)
d1[1,2] - d1[2,1] # returns 0.0

A fix is to force, in the definition of inv(::Cholesky)

StaticArrays.jl/src/cholesky.jl

Lines 72 to 74 in c4092a1

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

the symmetry of the result:

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