`eigen` not working properly for 2x2 Hermitian SMatrices

[guberger]

Hello! I noticed an incorrect behavior using eigen() on 2x2 Hermitian SMatrices:

using LinearAlgebra
using StaticArrays
A = @SMatrix [1.0 1e-10; 1e-10 1.0]
eigen(A)

outputs

Eigen{Float64, Float64, SMatrix{2, 2, Float64, 4}, SVector{2, Float64}}        
values:
2-element SVector{2, Float64} with indices SOneTo(2):
 1.0
 1.0
vectors:
2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 0.0  0.0
 1.0  1.0

The problem seems to be a numerical instability in

StaticArrays.jl/src/eigen.jl

Lines 161 to 176 in f59bceb

	t_half = real(a[1] + a[4]) / 2
	d = real(a[1] * a[4] - a[3]' * a[3]) # Should be real
	tmp2 = t_half * t_half - d
	tmp = tmp2 < 0 ? zero(tmp2) : sqrt(tmp2) # Numerically stable for identity matrices, etc.
	vals = SVector(t_half - tmp, t_half + tmp)
	v11 = vals[1] - a[4]
	n1 = sqrt(v11' * v11 + a[3]' * a[3])
	v11 = v11 / n1
	v12 = a[3]' / n1
	v21 = vals[2] - a[4]
	n2 = sqrt(v21' * v21 + a[3]' * a[3])
	v21 = v21 / n2
	v22 = a[3]' / n2

where vals should be [-1e-10, 1e-10] instead of [0.0, 0.0]. This has consequences on the rest of the evaluations, e.g., v11 is evaluated to 0.0 instead of -1e-10 and n1 is evaluated to 1e-10 instead of 2e-10 etc.

More info:

julia> versioninfo()
Julia Version 1.6.2
Commit 1b93d53fc4 (2021-07-14 15:36 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-7700 CPU @ 3.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-11.0.1 (ORCJIT, skylake)
(@v1.6) pkg> status StaticArrays
      Status `~\v1.6\Project.toml`
  [90137ffa] StaticArrays v1.2.12

[kylebeggs]

Just came across this bug myself... hopefully this fix gets merged soon! Thanks.

[KnutAM]

Fixed by #1158 ? On StaticArrays v1.6.0:

julia> A = @SMatrix [1.0 1e-10; 1e-10 1.0]
2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 1.0      1.0e-10
 1.0e-10  1.0

julia> eigen(A)
LinearAlgebra.Eigen{Float64, Float64, SMatrix{2, 2, Float64, 4}, SVector{2, Float64}}
values:
2-element SVector{2, Float64} with indices SOneTo(2):
 0.9999999999
 1.0000000001
vectors:
2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 -0.707107  0.707107
  0.707107  0.707107

[guberger]

Unfortunately not:

A = @SMatrix [1.0 1e-18; 1e-18 1.0]
2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 1.0      1.0e-18
 1.0e-18  1.0
julia> eigen(A)
Eigen{Float64, Float64, SMatrix{2, 2, Float64, 4}, SVector{2, Float64}}
values:
2-element SVector{2, Float64} with indices SOneTo(2):
 1.0
 1.0
vectors:
2×2 SMatrix{2, 2, Float64, 4} with indices SOneTo(2)×SOneTo(2):
 0.0  0.0
 1.0  1.0

The reason is the same as before: vals[1] = 1-1e-18.
Hence, because of rounding, v11 = 0 (instead of -1e-18) which causes the error.

ps: the values are correct (up to machine accuracy), but the vectors are not correct (even taking maching inacurracy into account) (they are not orthogonal or near orthogonal for a start).