Improve accuracy of 2x2 symmetric/Hermitian eigendecomposition (#1158)

* Improve accuracy of 2x2 symmetric/Hermitian eigendecomposition * bump version * Update src/eigen.jl Co-authored-by: Yuto Horikawa <[email protected]> --------- Co-authored-by: Yuto Horikawa <[email protected]>
JuliaArrays · Apr 30, 2023 · dce8c1c · dce8c1c · mateuszbaran · Apr 30, 2023
1 parent b4e49ce
commit dce8c1c
Show file tree

Hide file tree

Showing 3 changed files with 8 additions and 10 deletions.
diff --git a/Project.toml b/Project.toml
@@ -1,6 +1,6 @@
 name = "StaticArrays"
 uuid = "90137ffa-7385-5640-81b9-e52037218182"
-version = "1.5.23"
+version = "1.5.24"
 
 [deps]
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"

diff --git a/src/eigen.jl b/src/eigen.jl
@@ -159,14 +159,10 @@ end
 @inline function _eig(::Size{(2,2)}, A::LinearAlgebra.RealHermSymComplexHerm{T}, permute, scale) where {T <: Real}
     a = A.data
     TA = eltype(A)
-
     @inbounds if A.uplo == 'U'
         if !iszero(a[3]) # A is not diagonal
             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.
+            tmp = norm(SVector((a[1] - a[4])/2, a[3]'))
             vals = SVector(t_half - tmp, t_half + tmp)
 
             v11 = vals[1] - a[4]
@@ -187,10 +183,7 @@ end
     else # A.uplo == 'L'
         if !iszero(a[2]) # A is not diagonal
             t_half = real(a[1] + a[4]) / 2
-            d = real(a[1] * a[4] - a[2]' * a[2]) # Should be real
-
-            tmp2 = t_half * t_half - d
-            tmp = tmp2 < 0 ? zero(tmp2) : sqrt(tmp2) # Numerically stable for identity matrices, etc.
+            tmp = norm(SVector((a[1] - a[4])/2, a[2]))
             vals = SVector(t_half - tmp, t_half + tmp)
 
             v11 = vals[1] - a[4]

diff --git a/test/eigen.jl b/test/eigen.jl
@@ -93,6 +93,11 @@ using StaticArrays, Test, LinearAlgebra
         @test (@inferred SVector{2,ComplexF64} eigvals(Symmetric(m1), Symmetric(m2))) ≈ eigvals(Symmetric(m1_a), Symmetric(m2_a))
     end
 
+    @testset "2×2, difficult case" begin
+        m1 = SA[1.6590891025248637 -2.7087777909606814e-7; -2.7087777909606814e-7 1.659089317183428]
+        @test norm(m1 - Array(eigen(m1))) < 1e-15
+    end
+
     @test_throws DimensionMismatch eigvals(SA[1 2 3; 4 5 6], SA[1 2 3; 4 5 5])
     @test_throws DimensionMismatch eigvals(SA[1 2; 4 5], SA[1 2 3; 4 5 5; 3 4 5])