Necessity of the hard symmetry/hermiticity checks in `eigh(!)` decompositions

[leburgel]

When using the Hermitian eigenvalue decompositions eigh_full and eigh_trunc, I've run into the practical issue that I almost always have to manuall 'resymmetrize' my input matrix in order to be able to use these decompositions. If I don't do this, numerical inaccuracies always trigger the input checks in the implementation of the LAPACK wrappers. There's currently no way of using the Hermitian decompositions 'at my own risk' for not entirely Hermitian objects. This is quite different behavior from the LinearAlgebra.LAPACK routines, where there's no such strict input checking.

For example, I think it would not be unreasonable to use the Hermitian decompositions on objects of the kind:

using LinearAlgebra
using MatrixAlgebraKit

a = randn(ComplexF64, 10, 10)
a = (a + a') / 2
ap = a + 1e-14 * norm(a, Inf) * randn(eltype(a), size(a)...) # perturb slightly
@show ishermitian(ap) # false
@show norm(ap - ap') # very small

However, eigh_full using a divide and conquer (or any other) algorithm just errors out inYALAPACK.heevd!

alg = LAPACK_DivideAndConquer()
D, V = eigh_full(ap, alg)

Stacktrace:
 [1] heevd!(A::Matrix{ComplexF64}, W::Vector{Float64}, V::Matrix{ComplexF64}; uplo::Char)
   @ MatrixAlgebraKit.YALAPACK ~/.julia/packages/MatrixAlgebraKit/xv5dv/src/yalapack.jl:1300
 [2] heevd!
   @ ~/.julia/packages/MatrixAlgebraKit/xv5dv/src/yalapack.jl:1288 [inlined]
 [3] eigh_full!
   @ ~/.julia/packages/MatrixAlgebraKit/xv5dv/src/implementations/eigh.jl:86 [inlined]
 [4] eigh_full!(A::Matrix{ComplexF64}, alg::LAPACK_DivideAndConquer{@NamedTuple{}})
   @ MatrixAlgebraKit ~/.julia/packages/MatrixAlgebraKit/xv5dv/src/algorithms.jl:244
 [5] eigh_full(A::Matrix{ComplexF64}, alg::LAPACK_DivideAndConquer{@NamedTuple{}})
   @ MatrixAlgebraKit ~/.julia/packages/MatrixAlgebraKit/xv5dv/src/algorithms.jl:234
 [6] top-level scope
   @ REPL[84]:1

At the same time, I can use the routines from LinearAlgebra.LAPACK just fine on this kind of object,

D, V = LAPACK.syevd!('V', 'U', copy(ap)))
D

10-element Vector{Float64}:
 -3.407857755214439
 -2.387292210584882
 -1.45936482040303
 -0.68936395419045
 -0.21846056222030638
  0.6404302497673092
  1.5579324860696977
  2.3320360765034693
  3.2679939158673332
  3.650456704365599

Would it be reasonable to either remove these input checks, make them less strict by checking with respect to some tolerance, or at the very least replace the error with a warning so people know that what they're doing might not be technically correct and they acting at their own risk?

[lkdvos]

I dove a bit deeper into the way LinearAlgebra handles this, let me summarize what is going on there:

• eigen simply checks for exact hermiticity, and dispatches to the hermitian/non-hermitian implementation based on that.
• eigen(::Hermitian) has to copy the input before it can call eigen!(::Hermitian), and in that copy they only copy the upper/lower triangular part
• eigen!(::Hermitian), calls the corresponding LAPACK routine, which takes in which triangular part it has to consider, and ignores the other.

Given that our eigh calls out to LAPACK as well, effectively we are ignoring the lower triangular part of the input.

---

So, I can think of a couple solutions to this problem:

• eigh is "use-at-your-own-risk", so we disable the checks (possibly changing them for debug messages) and just use the LAPACK implementations as is. This is very convenient if you know what you are doing, but a bit nontransparent for figuring out that we are silently ignoring the lower triangular part if you select a lapack backend.
• eigh(A) projects onto the hermitian part instead of simply copying, eigh!(A) keeps its check
• We provide a higher-level wrapper like LinearAlgebra that automatically dispatches to the correct eig or eigh.
• The LAPACK algorithms add a keyword argument for the ishermitian tolerance.

I'm a bit on the fence about what we should do, especially since the projection done by project_hermitian really isn't the same as the projection done by ignoring the lower triangular part, and it does seem like silently ignoring large differences might be dangerous.
I think that because we don't have an interface for automatically selecting the right kind of eigensolver, I'm less inclined of making the eigh implementations fully "use at your own risk", so maybe we should just invest into that?

[mtfishman]

What about algorithm wrappers where you can either select a tolerance for how far off of Hermitian the input is allowed to be and/or select some way to project the input to be Hermitian (either project_hermitian or ignoring the upper or lower triangular part)? It seems better to be strict by default and require the user to be explicit about what they want, I find the LinearAlgebra.jl interface to be opaque and hard to remember, which is convenient until you have an algorithm that is sensitive to these choices...

Maybe a compromise could be to choose a tolerance that is comparable to the one LAPACK uses (say something close to eps) but then allow users to change that tolerance, but again the issue is that it is hard to remember which kind of projection is used and ideally we would allow users to choose the projection.

[lkdvos]

LAPACK doesn't really have a tolerance, you basically just give it a matrix and tell it to only look at the upper or lower triangular part, and it will ignore the other part. It therefore also doesn't error, because you really need to explicitly tell it which half to look at, so the interface doesn't allow inputting a nonhermitian matrix.

We can copy that behavior by simply not checking anything, in which case eigh with the LAPACK routines basically becomes equivalent to looking only at the upper triangular part.

I'm not sure we really need to handle the kind of projection used from within eigh, mostly because:

• We can't actually control what LAPACK does, so we don't have that much control over the kind of projection used.
• eigh_full!(project_hermitian!(A)) seems reasonable to me as an interface, if you want control over the kind of projection it might be better to be explicit about it. (in principle if you are using a LAPACK algorithm you could consider only projecting the upper triangular part, but that might be taking it a bit far)

I think the compromise is probably the easiest to implement and maintain, where indeed we just set a tolerance that is reasonable, document what LAPACK does, and allow to customize this with the additional keyword argument

[mtfishman]

I think having a tolerance keyword argument with a reasonable default along with having users explicitly call project_hermitian is a good plan. If the input falls within the tolerance would it just use the upper triangular part?