From 6670fd275aa5676ef9a25731fd704b2774ed5ff6 Mon Sep 17 00:00:00 2001
From: Guillaume Berger <47389731+guberger@users.noreply.github.com>
Date: Sat, 4 Sep 2021 14:10:03 +0200
Subject: [PATCH 1/2] Update eigen.jl

Addressing issue https://github.com/JuliaArrays/StaticArrays.jl/issues/956 on the numerical instability of the `eigen` function for Hermitian 2x2 SMatrix. The previous implementation also failed when the `A` is not Hermitian but `Hermitian(A)` is passed as argument to `eigen`.
---
 src/eigen.jl | 38 +++++++++++++++++++-------------------
 1 file changed, 19 insertions(+), 19 deletions(-)

diff --git a/src/eigen.jl b/src/eigen.jl
index e66aea1d..38a362f7 100644
--- a/src/eigen.jl
+++ b/src/eigen.jl
@@ -158,20 +158,20 @@ end
 
     @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.
+            t_half = (real(a[1]) + real(a[4])) / 2
+            s_half = (real(a[1]) - real(a[4])) / 2
+            tmp2 = s_half^2 + abs2(a[3]) # expansion of (tr(A)^2 - 4*det(A)) / 4
+            # use abs2 to ensure that it is real and > 0 (not sure about x'x being always > 0 if a[3] ≠ 0)
+            tmp = sqrt(tmp2) # normally, tmp > 0 if tmp2 > 0
             vals = SVector(t_half - tmp, t_half + tmp)
 
-            v11 = vals[1] - a[4]
-            n1 = sqrt(v11' * v11 + a[3]' * a[3])
+            v11 = s_half - tmp
+            n1 = sqrt(abs2(v11) + abs2(a[3]))
             v11 = v11 / n1
             v12 = a[3]' / n1
 
-            v21 = vals[2] - a[4]
-            n2 = sqrt(v21' * v21 + a[3]' * a[3])
+            v21 = s_half + tmp
+            n2 = sqrt(abs2(v21) + abs2(a[3]))
             v21 = v21 / n2
             v22 = a[3]' / n2
 
@@ -182,27 +182,27 @@ end
         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.
+            t_half = (real(a[1]) + real(a[4])) / 2
+            s_half = (real(a[1]) - real(a[4])) / 2 
+            tmp2 = s_half^2 + abs2(a[2]) # expansion of (tr(A)^2 - 4*det(A)) / 4
+            # use abs2 to ensure that it is real and > 0 (not sure about x'x being always > 0 if a[3] ≠ 0)
+            tmp = sqrt(tmp2) # normally, tmp > 0 if tmp2 > 0
             vals = SVector(t_half - tmp, t_half + tmp)
 
-            v11 = vals[1] - a[4]
-            n1 = sqrt(v11' * v11 + a[2]' * a[2])
+            v11 = s_half - tmp
+            n1 = sqrt(abs2(v11) + abs2(a[3]))
             v11 = v11 / n1
             v12 = a[2] / n1
 
-            v21 = vals[2] - a[4]
-            n2 = sqrt(v21' * v21 + a[2]' * a[2])
+            v21 = s_half + tmp
+            n2 = sqrt(abs2(v21) + abs2(a[3]))
             v21 = v21 / n2
             v22 = a[2] / n2
 
             vecs = @SMatrix [ v11  v21 ;
                               v12  v22 ]
 
-            return Eigen(vals,vecs)
+            return Eigen(vals, vecs)
         end
     end
 

From f5d94d8a8c2d9d24d814f88c326f291120dc8303 Mon Sep 17 00:00:00 2001
From: "guberger@PC-1r0-013" <guillaume.berger@uclouvain.be>
Date: Sat, 4 Sep 2021 17:28:13 +0200
Subject: [PATCH 2/2] Changes _eigvals for 2x2 Hermitian matrices

---
 src/eigen.jl | 21 +++++----------------
 1 file changed, 5 insertions(+), 16 deletions(-)

diff --git a/src/eigen.jl b/src/eigen.jl
index 38a362f7..5e1efc4a 100644
--- a/src/eigen.jl
+++ b/src/eigen.jl
@@ -31,22 +31,11 @@ end
 
 @inline function _eigvals(::Size{(2,2)}, A::LinearAlgebra.RealHermSymComplexHerm{T}, permute, scale) where {T <: Real}
     a = A.data
-
-    if A.uplo == 'U'
-        @inbounds t_half = real(a[1] + a[4])/2
-        @inbounds d = real(a[1]*a[4] - a[3]'*a[3]) # Should be real
-
-        tmp2 = t_half*t_half - d
-        tmp2 < 0 ? tmp = zero(tmp2) : tmp = sqrt(tmp2) # Numerically stable for identity matrices, etc.
-        return SVector(t_half - tmp, t_half + tmp)
-    else
-        @inbounds t_half = real(a[1] + a[4])/2
-        @inbounds d = real(a[1]*a[4] - a[2]'*a[2]) # Should be real
-
-        tmp2 = t_half*t_half - d
-        tmp2 < 0 ? tmp = zero(tmp2) : tmp = sqrt(tmp2) # Numerically stable for identity matrices, etc.
-        return SVector(t_half - tmp, t_half + tmp)
-    end
+    @inbounds t_half = (real(a[1]) + real(a[4])) / 2
+    @inbounds s_half = (real(a[1]) - real(a[4])) / 2
+    @inbounds tmp2 = A.uplo == 'U' ? s_half^2 + abs2(a[3]) : s_half^2 + abs2(a[2]) # expansion of (tr(A)^2 - 4*det(A)) / 4
+    tmp = sqrt(tmp2) # normally, tmp > 0 if tmp2 > 0
+    return SVector(t_half - tmp, t_half + tmp)
 end
 
 @inline function _eigvals(::Size{(3,3)}, A::LinearAlgebra.RealHermSymComplexHerm{T}, permute, scale) where {T <: Real}