Fix 2x2 symmetric eigenvalue computation (#124)

Author: andreasnoack

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/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