`eigen` returns `ComplexF64` eigenvalues for a Hermitian matrix

[yakovbraver]

Consider calling eigen for a ComplexF64 matrix which happens to be real-diagonal:

julia> D = ComplexF64[1 0; 0 2];
julia> eigen(D)
Eigen{ComplexF64, ComplexF64, Matrix{ComplexF64}, Vector{ComplexF64}}
values:
2-element Vector{ComplexF64}:
 1.0 + 0.0im
 2.0 + 0.0im
vectors:
2×2 Matrix{ComplexF64}:
 1.0+0.0im  0.0+0.0im
 0.0+0.0im  1.0+0.0im

The eigenvalues are returned in a Vector{ComplexF64} while I would expect a Vector{Float64}, as is the case for any Hermitian matrix:

julia> M = ComplexF64[1 1; 1 2];
julia> eigen(M)
Eigen{ComplexF64, Float64, Matrix{ComplexF64}, Vector{Float64}}
values:
2-element Vector{Float64}:
 0.38196601125010515
 2.618033988749895
vectors:
2×2 Matrix{ComplexF64}:
 -0.850651+0.0im  0.525731+0.0im
  0.525731+0.0im  0.850651+0.0im

Moreover, eigvals does return a Vector{Float64} for the above diagonal matrix D:

julia> eigvals(D)
2-element Vector{Float64}:
 1.0
 2.0

Seems inconsistent?

versioninfo

```julia Julia Version 1.8.3 Commit 0434deb161 (2022-11-14 20:14 UTC) Platform Info: OS: Windows (x86_64-w64-mingw32) CPU: 8 × Intel(R) Core(TM) i7-6700K CPU @ 4.00GHz WORD_SIZE: 64 LIBM: libopenlibm LLVM: libLLVM-13.0.1 (ORCJIT, skylake) Threads: 8 on 8 virtual cores Environment: JULIA_NUM_THREADS = 8 ```

[jishnub]

For consistent results, you may use the Hermitian wrapper

julia> eigen(Hermitian(D))
Eigen{ComplexF64, Float64, Matrix{ComplexF64}, Vector{Float64}}
values:
2-element Vector{Float64}:
 1.0
 2.0
vectors:
2×2 Matrix{ComplexF64}:
 -1.0+0.0im  0.0+0.0im
  0.0+0.0im  1.0+0.0im

[mikmoore]

Indeed, Hermitian should be used when you know that your matrix has that structure. However, the point stands that we get different output types from these two inputs whose only difference is that one is diagonal and one is hermitian (but neither have a type indicating them as such).

julia> eigen(ComplexF64[1 0; 0 2])
Eigen{ComplexF64, ComplexF64, Matrix{ComplexF64}, Vector{ComplexF64}}
values:
2-element Vector{ComplexF64}:
 1.0 + 0.0im
 2.0 + 0.0im
vectors:
2×2 Matrix{ComplexF64}:
 1.0+0.0im  0.0+0.0im
 0.0+0.0im  1.0+0.0im

julia> eigen(ComplexF64[1 1; 1 2])
Eigen{ComplexF64, Float64, Matrix{ComplexF64}, Vector{Float64}}
values:
2-element Vector{Float64}:
 0.38196601125010515
 2.618033988749895
vectors:
2×2 Matrix{ComplexF64}:
 -0.850651+0.0im  0.525731+0.0im
  0.525731+0.0im  0.850651+0.0im

It seems that first we check for diagonal-ness and then for hermitian-ness and have different codepaths and return types for each. The diagonal case is trivial but can't generally be assumed to have real eigenvalues, hence the complex return type. The only reasonable way I see to "fix" this and achieve type stability is that hermitian inputs (but not Hermitian-typed inputs, obviously) should return Complex eigenvalues. We can't fix this other direction -- by making complex diagonal matrices with real-valued entries return real eigenvalues -- because that would be a type instability.

[jishnub]

It would certainly be good to have eigen and eigvals be consistent:

julia> D = ComplexF64[1 0; 0 2];

julia> using LinearAlgebra

julia> eigvals(D)
2-element Vector{Float64}:
 1.0
 2.0

julia> eigen(D).values
2-element Vector{ComplexF64}:
 1.0 + 0.0im
 2.0 + 0.0im

[yakovbraver]

Thanks for the tip!

We can't fix this other direction -- by making complex diagonal matrices with real-valued entries return real eigenvalues -- because that would be a type instability.

But this is how eigvals operates currently:

julia> eigvals(ComplexF64[1 0; 0 2])
2-element Vector{Float64}:
 1.0
 2.0

[mikmoore]

If eigvals operates that way currently, then it is also unstable. I'd suggest it be changed as well.

In summary, I propose that all Complex-valued inputs except Hermitian (but not excepting Symmetric) should have Complex eigenvalues with either eigen and eigvals.

[yakovbraver]

I tend to agree.