Repository: JuliaLinearAlgebra/GenericLinearAlgebra.jl (julia)

Issue:
Eigen returns non-eigenvectors for certain matrices
I have [a test](https://github.com/moble/Quaternionic.jl/blob/4061d5b9f2aa84f4e3678f8c4c6e445578a31766/test/conversion.jl#L148-L152) that I'm trying to get working for BigFloats.  It somewhat randomly constructs 5,000 symmetric real 4x4 matrices of BigFloats and extracts the dominant eigenvector of each.  For 3 of these instances, I get results for which the answer is not an eigenvector at all (never mind with the given eigenvalue).  There doesn't seem to be anything special about these matrices.  I've extracted one to provide a concrete MWE.  The eigenvalues of this matrix are expected to be (and reported as) [-1, -1, -1, 3], so I don't know that there would be any numerical problem here.

Here's the MWE:
```julia
using LinearAlgebra
using GenericLinearAlgebra
using StaticArrays

M = Symmetric(@SMatrix[
    big"-0.4080898675832881399369478084264191594976530854542904557798567397269356887436951";
        big"-0.1032324294981949906363774065395184125237581835226155628209100984396171211818558";
        big"-1.0795157507124452910839896877334667387210301781514938067860918240771876343947";
        big"0.9172086645212876240254394768180975107502376572771647296150618931226550446699544";;
    big"-0.1032324294981949906363774065395184125237581835226155628209100984396171211818558";
        big"-0.9819956883377066621250198846550622559246996804965712336465013506629992739010227";
        big"0.1882735697944729855991976669864503854920622386133987141371224931350749728226066";
        big"-0.1599663084136352437739757607131301560774255778371317602542426234968564801904052";;
    big"-1.0795157507124452910839896877334667387210301781514938067860918240771876343947";
        big"0.1882735697944729855991976669864503854920622386133987141371224931350749728226066";
        big"0.9688026817149176598146701814747478080649943014810992426739997593840858865727305";
        big"-1.672789745967021000172452940954243617442140494364475046869527486458478435262502";;
    big"0.9172086645212876240254394768180975107502376572771647296150618931226550446699544";
        big"-0.1599663084136352437739757607131301560774255778371317602542426234968564801904052";
        big"-1.672789745967021000172452940954243617442140494364475046869527486458478435262502";
        big"0.4212828742060771422472975116067336073573584644697624467523583310058490760719874"
])

eigensystem = eigen(M)
```
Here, the supposed dominant eigenvector is
```julia
julia> eigensystem.vectors[:, 4]
4-element SVector{4, BigFloat} with indices SOneTo(4):
  0.5293186034906468513583044045578991187176767353991885285160593337162140639833333
  0.6076028636690888706511331534102457699010874314100936245656011216542549302322659
 -0.5375888975041810576235384352977728036722978161994470931619050808663662918353576
  0.2482715314732502633344123259979948815899264503590173204924838651280614259806134
```
and its reported eigenvalue (3) is nonzero, so
```julia
(M * eigensystem.vectors[:, 4]) ./ eigensystem.vectors[:, 4]
```
*should* be a vector with entries that are at least approximately equal.  Instead, I get
```julia
4-element SVector{4, BigFloat} with indices SOneTo(4):
  0.9999999999999999999999999999999999999999999999999999999999999999999999999999309
 -1.303869922808737255865818137316681872433951050681740949789081294567406034061292
  2.591452043646991636425616621765959618553871235916851252573571248344097579507887
  5.607430516593620629517859036040424284010328224975389025670154530756814862123736
```
It's interesting to note `eigensystem.vectors[:, 1]` is also not an eigenvector.

The dominant eigenvector (with eigenvalue 3) that I expect in this case is
```julia
x4 = [
    big"-0.3846784801677603137399569700744195708135504357934165621816525583943309062098062",
    big"0.06709007315224313025838456104648289261272108291672656267831069781563681483410404",
    big"0.7015701464776914343825844884840609580397271879944290132648084132383274894856047",
    big"-0.5960878446600964581372146249105424317644837852918997676276656289252753627271997",
]
```
and indeed I get the expected eigenvector behavior with the given matrix:
```julia
julia> (M * x4) ./ x4
4-element Vector{BigFloat}:
 3.0
 2.999999999999999999999999999999999999999999999999999999999999999999999999999862
 2.999999999999999999999999999999999999999999999999999999999999999999999999999965
 3.000000000000000000000000000000000000000000000000000000000000000000000000000035
```


---

```julia
julia> versioninfo()
Julia Version 1.9.0
Commit 8e630552924 (2023-05-07 11:25 UTC)
Platform Info:
  OS: macOS (x86_64-apple-darwin22.4.0)
  CPU: 8 × Intel(R) Core(TM) i7-4770HQ CPU @ 2.20GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-14.0.6 (ORCJIT, haswell)
  Threads: 1 on 8 virtual cores

(jl_6G5ebb) pkg> status
Status `/private/var/folders/76/v4_xw_mx1ml49pbk_yd5pnrh0000gn/T/jl_6G5ebb/Project.toml`
  [14197337] GenericLinearAlgebra v0.3.10
  [90137ffa] StaticArrays v1.5.25
  [37e2e46d] LinearAlgebra
```

Pull request description:
The old versions had some unstable corners.

Fixes #123

Patch:
diff --git a/src/eigenSelfAdjoint.jl b/src/eigenSelfAdjoint.jl
index a42e3e9..c0b6a2a 100644
--- a/src/eigenSelfAdjoint.jl
+++ b/src/eigenSelfAdjoint.jl
@@ -135,24 +135,26 @@ function eigvals2x2(d1::Number, d2::Number, e::Number)
     r2 = hypot(d1 - d2, 2 * e) / 2
     return r1 + r2, r1 - r2
 end
-function eigQL2x2!(d::StridedVector, e::StridedVector, j::Integer, vectors::Matrix)
+function eig2x2!(d::StridedVector, e::StridedVector, j::Integer, vectors::Matrix)
     dj = d[j]
-    dj1 = d[j+1]
+    dj1 = d[j + 1]
     ej = e[j]
     r1 = (dj + dj1) / 2
     r2 = hypot(dj - dj1, 2 * ej) / 2
-    λ = r1 + r2
-    d[j] = λ
-    d[j+1] = r1 - r2
+    λ⁺ = r1 + r2
+    λ⁻ = r1 - r2
+    d[j] = λ⁺
+    d[j + 1] = λ⁻
     e[j] = 0
-    c = ej / (dj - λ)
-    if isfinite(c) # FixMe! is this the right fix for overflow?
-        h = hypot(one(c), c)
-        c /= h
-        s = inv(h)
+    if iszero(ej)
+        c = one(λ⁺)
+        s = zero(λ⁺)
+    elseif abs(λ⁺ - dj) > abs(λ⁻ - dj)
+        c = -ej / hypot(ej, λ⁺ - dj)
+        s = sqrt(1 - c*c)
     else
-        c = one(c)
-        s = zero(c)
+        s = abs(ej) / hypot(ej, λ⁻ - dj)
+        c = copysign(sqrt(1 - s*s), -ej)
     end
 
     _updateVectors!(c, s, j, vectors)
@@ -173,14 +175,16 @@ function eigvalsPWK!(S::SymTridiagonal{T}; tol = eps(T)) where {T<:Real}
         end
         while true
             # Check for zero off diagonal elements
-            for be = blockstart+1:n
-                if abs(e[be-1]) <= tol * sqrt(abs(d[be-1])) * sqrt(abs(d[be]))
+            for be = (blockstart + 1):n
+                if abs(e[be - 1]) <= tol * sqrt(abs(d[be - 1])) * sqrt(abs(d[be]))
                     blockend = be - 1
                     break
                 end
                 blockend = n
             end
 
+            @debug "submatrix" blockstart blockend
+
             # Deflate?
             if blockstart == blockend
                 @debug "Deflate? Yes!"
@@ -200,11 +204,12 @@ function eigvalsPWK!(S::SymTridiagonal{T}; tol = eps(T)) where {T<:Real}
                 μ = d[blockstart] - sqrte / (μ + copysign(r, μ))
 
                 # QL bulk chase
+                @debug "Values before PWK QL bulge chase" e[blockstart] e[blockend-1] μ
                 singleShiftQLPWK!(S, μ, blockstart, blockend)
 
                 rotations = blockend - blockstart
                 iter += rotations
-                @debug "QL, blockstart, blockend, e[blockstart], e[blockend-1], μ, rotations" blockstart blockend e[blockstart] e[blockend-1] μ rotations
+                @debug "Values after PWK QL bulge chase" e[blockstart] e[blockend-1] rotations
             end
             if blockstart >= n
                 break
@@ -243,32 +248,37 @@ function eigQL!(
     @inbounds begin
         while true
             # Check for zero off diagonal elements
-            for be = blockstart+1:n
-                if abs(e[be-1]) <= tol * sqrt(abs(d[be-1])) * sqrt(abs(d[be]))
+            for be = (blockstart + 1):n
+                if abs(e[be - 1]) <= tol * sqrt(abs(d[be - 1])) * sqrt(abs(d[be]))
                     blockend = be - 1
                     break
                 end
                 blockend = n
             end
+
+            @debug "submatrix" blockstart blockend
+
             # Deflate?
             if blockstart == blockend
                 @debug "Deflate? Yes!"
                 blockstart += 1
             elseif blockstart + 1 == blockend
                 @debug "Deflate? Yes, but after solving 2x2 block"
-                eigQL2x2!(d, e, blockstart, vectors)
-                blockstart += 1
+                eig2x2!(d, e, blockstart, vectors)
+                blockstart += 2
             else
                 @debug "Deflate? No!"
                 # Calculate shift
-                μ = (d[blockstart+1] - d[blockstart]) / (2 * e[blockstart])
+                μ = (d[blockstart + 1] - d[blockstart]) / (2 * e[blockstart])
                 r = hypot(μ, one(T))
                 μ = d[blockstart] - (e[blockstart] / (μ + copysign(r, μ)))
 
                 # QL bulk chase
+                @debug "Values before bulge chase" e[blockstart] e[blockend - 1] d[blockstart] μ
                 singleShiftQL!(S, μ, blockstart, blockend, vectors)
-                @debug "QL, blockstart, blockend, e[blockstart], e[blockend-1] μ" blockstart blockend e[blockstart] e[blockend-1] μ
+                @debug "Values after bulge chase" e[blockstart] e[blockend - 1] d[blockstart]
             end
+
             if blockstart >= n
                 break
             end
@@ -291,32 +301,35 @@ function eigQR!(
     @inbounds begin
         while true
             # Check for zero off diagonal elements
-            for be = (blockstart+1):n
-                if abs(e[be-1]) <= tol * sqrt(abs(d[be-1])) * sqrt(abs(d[be]))
+            for be = (blockstart + 1):n
+                if abs(e[be - 1]) <= tol * sqrt(abs(d[be - 1])) * sqrt(abs(d[be]))
                     blockend = be - 1
                     break
                 end
                 blockend = n
             end
+
+            @debug "submatrix" blockstart blockend
+
             # Deflate?
             if blockstart == blockend
                 @debug "Deflate? Yes!"
                 blockstart += 1
             elseif blockstart + 1 == blockend
                 @debug "Deflate? Yes, but after solving 2x2 block!"
-                eigQL2x2!(d, e, blockstart, vectors)
-                blockstart += 1
+                eig2x2!(d, e, blockstart, vectors)
+                blockstart += 2
             else
                 @debug "Deflate? No!"
                 # Calculate shift
-                μ = (d[blockend-1] - d[blockend]) / (2 * e[blockend-1])
+                μ = (d[blockend - 1] - d[blockend]) / (2 * e[blockend - 1])
                 r = hypot(μ, one(T))
-                μ = d[blockend] - (e[blockend-1] / (μ + copysign(r, μ)))
+                μ = d[blockend] - (e[blockend - 1] / (μ + copysign(r, μ)))
 
                 # QR bulk chase
+                @debug "Values before QR bulge chase" e[blockstart] e[blockend - 1] d[blockend] μ
                 singleShiftQR!(S, μ, blockstart, blockend, vectors)
-
-                @debug "QR, blockstart, blockend, e[blockstart], e[blockend-1], d[blockend], μ" blockstart blockend e[blockstart] e[blockend-1] d[blockend] μ
+                @debug "Values after QR bulge chase" e[blockstart] e[blockend - 1] d[blockend]
             end
             if blockstart >= n
                 break


Test patch:
diff --git a/test/eigenselfadjoint.jl b/test/eigenselfadjoint.jl
index 4599db2..a5a7508 100644
--- a/test/eigenselfadjoint.jl
+++ b/test/eigenselfadjoint.jl
@@ -118,4 +118,43 @@ Base.isreal(q::Quaternion) = q.v1 == q.v2 == q.v3 == 0
         T = T + T'
         @test issorted(eigvals(T))
     end
+
+    @testset "issue 123" begin
+        M = Hermitian([
+            big"-0.4080898675832881399369478084264191594976530854542904557798567397269356887436951";
+            big"-0.1032324294981949906363774065395184125237581835226155628209100984396171211818558";
+            big"-1.0795157507124452910839896877334667387210301781514938067860918240771876343947";
+            big"0.9172086645212876240254394768180975107502376572771647296150618931226550446699544";;
+            big"-0.1032324294981949906363774065395184125237581835226155628209100984396171211818558";
+            big"-0.9819956883377066621250198846550622559246996804965712336465013506629992739010227";
+            big"0.1882735697944729855991976669864503854920622386133987141371224931350749728226066";
+            big"-0.1599663084136352437739757607131301560774255778371317602542426234968564801904052";;
+            big"-1.0795157507124452910839896877334667387210301781514938067860918240771876343947";
+            big"0.1882735697944729855991976669864503854920622386133987141371224931350749728226066";
+            big"0.9688026817149176598146701814747478080649943014810992426739997593840858865727305";
+            big"-1.672789745967021000172452940954243617442140494364475046869527486458478435262502";;
+            big"0.9172086645212876240254394768180975107502376572771647296150618931226550446699544";
+            big"-0.1599663084136352437739757607131301560774255778371317602542426234968564801904052";
+            big"-1.672789745967021000172452940954243617442140494364475046869527486458478435262502";
+            big"0.4212828742060771422472975116067336073573584644697624467523583310058490760719874"
+        ])
+        F = eigen(M)
+        @test M * F.vectors ≈ F.vectors * Diagonal(F.values)
+
+        @testset "2x2 cases, d: $d, e: $e" for (d, e) in (
+            ([1.0, 1.0], [0.0]),
+            ([0.0, 0.0], [1.0]),
+            ([1.0, 1.0], [1e-20]),
+            ([1e-20, 1e-30], [0.5]),
+            ([1e-20, 1e-30], [1e-10]),
+            ([1e-20, 1e-30], [1e-40]),
+            ([1e20, 1e30], [0.5]),
+            ([1e20, 1e30], [1e10]),
+            ([1e20, 1e30], [1e40])
+        )
+            V = Matrix{Float64}(I, 2, 2)
+            c, s = GenericLinearAlgebra.eig2x2!(d, e, 1, V)
+            @test hypot(c, s) ≈ 1
+        end
+    end
 end