The LinearAlgebra.jl package's `inv()` function does not preserve matrix symmetry

The LinearAlgebra.jl package’s inv() function does not preserve matrix symmetry.

while true
    R = randn((3,3))
    R[2,1] = R[1,2]; R[3,1] = R[1,3]; R[2,3]=R[3,2]
    @assert R == R'
    @assert inv(R) == inv(R)'
end

I’m about to open an issue for this, but in the meantime: is there a way to deal with this?

Note that the original poster on Slack cannot see your response here on Discourse. Consider transcribing the appropriate answer back to Slack, or pinging the poster here on Discourse so they can follow this thread.
(Original message :slack:) (More Info)

You can use Symmetric for this. It will make inv use the Bunch-Kaufman factorization which preserves symmetry (by only working on one of the triangles) instead of the default LU factorization where rounding errors will make the result slightly asymmetric. I.e.

julia> issymmetric(inv(R))
false

julia> issymmetric(inv(Symmetric(R)))
true

julia> typeof(inv(Symmetric(R)))
Symmetric{Float64, Matrix{Float64}}
5 Likes

Thanks :+1:

Also, floating point comparisons should use approximate equality.

You shouldn’t really be using matrix inversion anyway.