Difference between Mac and Linux for 2x2 non-symmetric matrices with real eigenvalues

[DavidMeadon]

There seems to be some difference between Mac and Linux for non-symmetric matrices with real eigenvalues:

# Mac
julia> A = [3 1; -1 5]
2×2 Matrix{Int64}:
  3  1
 -1  5

julia> eigen(A)
Eigen{ComplexF64, ComplexF64, Matrix{ComplexF64}, Vector{ComplexF64}}
values:
2-element Vector{ComplexF64}:
 4.0 - 9.009566630008139e-9im
 4.0 + 9.009566630008139e-9im
vectors:
2×2 Matrix{ComplexF64}:
 0.707107-0.0im         0.707107+0.0im
 0.707107-6.37073e-9im  0.707107+6.37073e-9im

julia> versioninfo()
Julia Version 1.11.2
Commit 5e9a32e7af2 (2024-12-01 20:02 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin24.0.0)
  CPU: 12 × Apple M2 Pro
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, apple-m2)
Threads: 1 default, 0 interactive, 1 GC (on 8 virtual cores)

# Linux
julia> A = [3 1; -1 5]
2×2 Matrix{Int64}:
  3  1
 -1  5

julia> eigen(A)
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
2-element Vector{Float64}:
 4.0
 4.0
vectors:
2×2 Matrix{Float64}:
 0.707107  -0.707107
 0.707107  -0.707107

julia> versioninfo()
Julia Version 1.11.2
Commit 5e9a32e7af2 (2024-12-01 20:02 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 64 × AMD EPYC 7502P 32-Core Processor
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, znver2)
Threads: 1 default, 0 interactive, 1 GC (on 64 virtual cores)

[stevengj]

That matrix is defective: in exact arithmetic, it has an eigenvalue $\lambda = 4$ with algebraic multiplicity 2 (i.e. a double root of the characteristic polynomial), but it only has a single linearly independent eigenvector [1, 1] (and multiples thereof), so it the eigenvalue has geometric multiplicity 1.

The problem is that repeated eigenvalues of defective matrices are extremely sensitive to roundoff errors — a random perturbation of order $\varepsilon$ causes the repeated eigenvalues to split and change by $O(\sqrt{\varepsilon})$.

So, slightly different rounding errors on different platforms/compilers/BLAS versions can cause this eigenvalue to split — either into two real eigenvalues or into two complex-conjugate eigenvalues. The platform-specific differences you are seeing are of precisely this nature, and are roughly comparable to the square root of the machine precision. It's pure luck that on one platform it happens to detect that the eigenvalue is exactly repeated.

This is inherent to the computation of repeated/ill-conditioned eigenvalues in finite-precision floating-point arithmetic, I'm afraid. There is a famous paper by Golub and Wilkinson (1976) on this issue. Nothing to do with Julia per se (these are just calls to the same LAPACK library that is used in Matlab, Scipy, R…) and probably nothing that can be practically addressed.

[DavidMeadon]

Yeah, I understand that the matrix is defective and that in general different platforms would give different solutions depending on rounding errors however I would have thought that for a 2x2 matrix, LAPACK would use the exact solution since there is a closed form solution. It seems that dlanv2 exists, so should Julia perhaps not call that if the size of the matrix is 2?

[stevengj]

If you want it to use different algorithms for 2x2 matrices, that's a different issue. But I suspect that specialized routines for 2x2 matrices belong in a different package, like StaticArrays (which currently doesn't support eigenvalues for non-Hermitian matrices: JuliaArrays/StaticArrays.jl#990).

(Of course, even with the closed-form quadratic-equation solution for 2x2 matrices you can have rounding errors if the terms aren't represented exactly.)