Julia: requiring type to signal Hermitian/Symmetric matrices

[tpapp]

(follow-up to a discussion on the forum)

This is a proposal to change the semantics of functions that require a symmetric/hermitian matrix as an argument.

Currently functions that require symmetric/hermitian matrices, like LinearAlgebra.cholesky, accept arbitrary <: AbstractMatrix arguments, then check for the desired symmetry. This leads to a lot of confusion among users since even if a calculation preserves symmetry in theory, in practice floating point error is very likely to break it. An example from that discussion is

julia> VERSION
v"1.2.0-DEV.219"

julia> using LinearAlgebra

julia> rot = [cos(π/4) -sin(π/4); sin(π/4) cos(π/4)]
2×2 Array{Float64,2}:
 0.707107  -0.707107
 0.707107   0.707107

julia> A = [25.0 0; 0 4.0]
2×2 Array{Float64,2}:
 25.0  0.0
  0.0  4.0

julia> cholesky(rot * A * rot')
ERROR: PosDefException: matrix is not Hermitian; Cholesky factorization failed.
Stacktrace:
 [1] checkpositivedefinite(::Int64) at /home/tamas/src/julia-git/usr/share/julia/stdlib/v1.2/LinearAlgebra/src/factorization.jl:11
 [2] #cholesky!#97(::Bool, ::Function, ::Array{Float64,2}, ::Val{false}) at /home/tamas/src/julia-git/usr/share/julia/stdlib/v1.2/LinearAlgebra/src/cholesky.jl:182
 [3] #cholesky#101 at ./none:0 [inlined]
 [4] cholesky at /home/tamas/src/julia-git/usr/share/julia/stdlib/v1.2/LinearAlgebra/src/cholesky.jl:275 [inlined] (repeats 2 times)
 [5] top-level scope at REPL[85]:1

but questions like this come up all the time.

The current workaround is to wrap a matrix in Symmetric or Hermitian. I am proposing that instead of working with general <: AbstractMatrix types, functions like cholesky only accept types which explicitly enforce symmetry. This would encourage a programming style that works uses wrapper types instead of (low-cost, but still significant) checks, and eliminate a surprises.

The change is breaking, even though a lot of existing code already wraps matrices so that would not be affected.

[andreasnoack]

Regardless, we could start by extending the error message to mention Hermitian.

[4aHxKzD]

Regardless, we could start by extending the error message to mention Hermitian.

Please do this. Just lost a bunch of hours trying to implement some Gaussian Processes code that just wouldn't work. I ended up re implementing in Python to find out that what I wanted was indeed possible and finally found this issue here. This could have been avoided if the documentation was a bit more clear :(