Merge pull request #89 from haampie/fix-schur-nearly-repeated-eigenvalue-issue

Fix issue with nearly repeated eigenvalue in Schur factorization

Author: haampie

src/run.jl

@@ -333,6 +333,7 @@ function partition_schur_three_way!(R, Q, groups::AbstractVector{Int})
     mi = 1
     lo = 1
 
+    # Be very scared whenever eigenvalues are too close!
     while hi ≤ length(groups)
         group = groups[hi]
         blocksize = is_start_of_11_block(R, hi) ? 1 : 2

src/schurfact.jl

@@ -5,7 +5,7 @@ import LinearAlgebra: lmul!, rmul!
 import Base: Matrix
 
 @propagate_inbounds is_offdiagonal_small(H::AbstractMatrix{T}, i::Int, tol = eps(real(T))) where {T} = 
-    abs(H[i+1,i]) < tol*(abs(H[i,i]) + abs(H[i+1,i+1]))
+    abs(H[i+1,i]) ≤ tol*(abs(H[i,i]) + abs(H[i+1,i+1]))
 
 
 abstract type SmallRotation end
@@ -352,7 +352,7 @@ function local_schurfact!(H::AbstractMatrix{T}, start::Int, to::Int,
             H[from,from-1] = zero(T)
             to -= 1
         else
-            # Now we are sure we can work with a 2x2 block H[to-1:to,to-1:to]
+            # Now we are sure we can work with a 2×2 block H[to-1:to,to-1:to]
             # We check if this block has a conjugate eigenpair, which might mean we have
             # converged w.r.t. this block if from + 1 == to. 
             # Otherwise, if from + 1 < to, we do either a single or double shift, based on
@@ -361,22 +361,30 @@ function local_schurfact!(H::AbstractMatrix{T}, start::Int, to::Int,
             H₁₁, H₁₂ = H[to-1,to-1], H[to-1,to]
             H₂₁, H₂₂ = H[to  ,to-1], H[to  ,to]
 
-            # Matrix determinant and trace
-            d = H₁₁ * H₂₂ - H₂₁ * H₁₂
-            t = H₁₁ + H₂₂
-            discriminant = t * t - 4d
+            # Scaling to avoid losing precision in the case where we have nearly
+            # repeated eigenvalues.
+            scale = abs(H₁₁) + abs(H₁₂) + abs(H₂₁) + abs(H₂₂)
+            H₁₁ /= scale; H₁₂ /= scale; H₂₁ /= scale; H₂₂ /= scale
 
-            if discriminant ≥ zero(T)
+            # Trace and discriminant of small eigenvalue problem.
+            t = (H₁₁ + H₂₂) / 2
+            d = (H₁₁ - t) * (H₂₂ - t) - H₁₂ * H₂₁
+            sqrt_discr = sqrt(abs(d))
+
+            # Very important to have a strict comparison here!
+            if d < zero(T)
                 # Real eigenvalues.
                 # Note that if from + 1 == to in this case, then just one additional
                 # iteration is necessary, since the Wilkinson shift will do an exact shift.
 
                 # Determine the Wilkinson shift -- the closest eigenvalue of the 2x2 block
                 # near H[to,to]
-                sqr = sqrt(discriminant)
-                λ₁ = (t + sqr) / 2
-                λ₂ = (t - sqr) / 2
+                
+                λ₁ = t + sqrt_discr
+                λ₂ = t - sqrt_discr
                 λ = abs(H₂₂ - λ₁) < abs(H₂₂ - λ₂) ? λ₁ : λ₂
+                λ *= scale
+
                 # Run a bulge chase
                 single_shift_schur!(H, from, to, λ, Q)
             else
@@ -389,8 +397,8 @@ function local_schurfact!(H::AbstractMatrix{T}, start::Int, to::Int,
                     to -= 2
                 else
                     # Otherwise we do a double shift!
-                    sqr = sqrt(complex(discriminant))
-                    double_shift_schur!(H, from, to, (t + sqr) / 2, Q)
+                    complex_shift = scale * (t + sqrt_discr * im)
+                    double_shift_schur!(H, from, to, complex_shift, Q)
                 end
             end
         end

src/schursort.jl

@@ -62,6 +62,7 @@ struct CompletelyPivotedLU{T,N,TA<:SMatrix{N,N,T},TP}
     A::TA
     p::TP
     q::TP
+    singular::Bool
 end
 
 struct CompletePivoting end
@@ -76,6 +77,7 @@ function lu(A::SMatrix{N,N,T}, ::Type{CompletePivoting}) where {N,T}
     A = MMatrix(A)
     p = @MVector fill(N, N)
     q = @MVector fill(N, N)
+    singular = false
 
     # Maybe I should consider doing this recursively.
     for k = OneTo(N - 1)
@@ -102,6 +104,12 @@ function lu(A::SMatrix{N,N,T}, ::Type{CompletePivoting}) where {N,T}
 
         Akk = A[k,k]
 
+        # It has actually happened :(
+        if iszero(Akk)
+            singular = true
+            break
+        end
+
         for i = k+1:N
             A[i, k] /= Akk
         end
@@ -115,8 +123,12 @@ function lu(A::SMatrix{N,N,T}, ::Type{CompletePivoting}) where {N,T}
         end
     end
 
+    if iszero(A[N,N])
+        singular = true
+    end
+
     # Back to immutable land!
-    return CompletelyPivotedLU(SMatrix(A), SVector(p), SVector(q))
+    return CompletelyPivotedLU(SMatrix(A), SVector(p), SVector(q), singular)
 end
 
 function (\)(LU::CompletelyPivotedLU{T,N}, b::SVector{N}) where {T,N}
@@ -158,22 +170,21 @@ end
               T(0)          -B[1,2]       A[2,1]        A[2,2]-B[2,2]]
 
 """
-    sylv(A, B, C) → X
+    sylv(A, B, C) → X, singular
 
 Solve A * X - X * B = C for X, where A and B are 1×1 or 2×2 matrices.
 
 It works by recasting the Sylvester equation to a linear system 
 (I ⊗ A + Bᵀ ⊗ I) vec(X) = vec(C) of size 2 or 4, which is then solved by
 Gaussian elimination with complete pivoting.
-"""
-@inline sylv(A::SMatrix{1,1,T}, B::SMatrix{2,2,T}, C::SMatrix{1,2,T}) where {T} =
-    SMatrix{1,2,T}(lu(sylvsystem(A, B), CompletePivoting) \ SVector{2,T}(C))
-
-@inline sylv(A::SMatrix{2,2,T}, B::SMatrix{1,1,T}, C::SMatrix{2,1,T}) where {T} = 
-    SMatrix{2,1,T}(lu(sylvsystem(A, B), CompletePivoting) \ SVector{2,T}(C))
 
-@inline sylv(A::SMatrix{2,2,T}, B::SMatrix{2,2,T}, C::SMatrix{2,2,T}) where {T} =
-    SMatrix{2,2,T}(lu(sylvsystem(A, B), CompletePivoting) \ SVector{4,T}(C))
+If the eigenvalues of A and B are equal, then `singular = true`.
+"""
+@inline function sylv(A::SMatrix{N,N,T}, B::SMatrix{M,M,T}, C::SMatrix{N,M,T}) where {T,N,M}
+    fact = lu(sylvsystem(A, B), CompletePivoting)
+    rhs = SVector{N*M,T}(C)
+    SMatrix{N,M,T}(fact \ rhs), fact.singular
+end
 
 """
     swap22_rotations(X) → c₁, s₁, c₂, s₂, c₃, s₃, c₄, s₄
@@ -289,7 +300,10 @@ function swap22!(R::AbstractMatrix{T}, i::Integer, Q = NotWanted()) where {T}
         C = @SMatrix [R[i+0,i+2] R[i+0,i+3]; R[i+1,i+2] R[i+1,i+3]]
 
         # A * X - X * B = C
-        X = sylv(A, B, C)
+        X, singular = sylv(A, B, C)
+
+        # No need to swap if eigenvalues are indistinguishable
+        singular && return R
 
         # Rotations that upper triangularize X
         c₁,s₁, c₂,s₂, c₃,s₃, c₄,s₄ = swap22_rotations(X)
@@ -340,7 +354,10 @@ function swap21!(R::AbstractMatrix{T}, i::Integer, Q = NotWanted()) where {T}
         C = @SMatrix [R[i+0,i+2]; R[i+1,i+2]]
 
         # A * X - X * B = C
-        X = sylv(A, B, C)
+        X, singular = sylv(A, B, C)
+
+        # No need to swap if eigenvalues are indistinguishable
+        singular && return R
 
         # Rotations that upper triangularize X
         c₁,s₁, c₂,s₂ = swap21_rotations(X)
@@ -388,7 +405,10 @@ function swap12!(R::AbstractMatrix{T}, i::Integer, Q = NotWanted()) where {T}
         C = @SMatrix [R[i+0,i+1] R[i+0,i+2]]
 
         # A * X - X * B = C
-        X = sylv(A, B, C)
+        X, singular = sylv(A, B, C)
+
+        # No need to swap if eigenvalues are indistinguishable
+        singular && return R
 
         # Rotations that upper triangularize X
         c₁,s₁, c₂,s₂ = swap12_rotations(X)

test/schurfact.jl

@@ -2,7 +2,8 @@
 # Here we look at some edge cases.
 
 using Test, LinearAlgebra
-using ArnoldiMethod: eigenvalues, local_schurfact!, is_offdiagonal_small
+using ArnoldiMethod: eigenvalues, local_schurfact!, is_offdiagonal_small,
+                     NotWanted
 
 include("utils.jl")
 
@@ -108,6 +109,16 @@ include("utils.jl")
             # Test that the eigenvalues of H are the same before and after transformation
             @test sort!(eigvals(H), by=realimag) ≈ sort!(eigvals(H′), by=realimag)
         end
-
     end
 end
+
+@testset "Schur with nearly repeated eigenvalues" begin
+    # Matrix with nearly repeated eigenpair could converge very slowly or
+    # stagnate completely when there's a tiny perturbation in the computed
+    # shift.
+    mat(ε) = [2   0    0  ;
+              5ε  1-ε  2ε ;
+              0   3ε   1+ε]
+
+    @test local_schurfact!(mat(eps()))
+end

test/sort_schur.jl

@@ -251,4 +251,19 @@ end
     @test abs(λs_before[3]) ≈ abs(λs_after[1])
     @test opnorm(I - Q'Q, 1) < 10eps() # we should be able to get rid of this prefactor?
     @test opnorm(A * Q - Q * A′, 1) < opnorm(A, 1) * eps()
+end
+
+@testset "Identical eigenvalues should not blow up" begin
+    A = [1.0 2.0 3.0 4.0;
+         0.0 1.0 5.0 6.0;
+         0.0 0.0 1.0 7.0;
+         0.0 0.0 0.0 1.0]
+         
+    A′ = copy(A)
+    swap22!(A′, 1)
+    @test A == A′
+    swap12!(A′, 1)
+    @test A == A′
+    swap21!(A′, 1)
+    @test A == A′
 end

test/sylvester.jl

@@ -1,14 +1,40 @@
 using Test
 using ArnoldiMethod: sylv
 using StaticArrays
+using LinearAlgebra
 
 @testset "Tiny sylvester equation $T, $s" for T in (Float64, ComplexF64), s = ((2,2),(2,1),(1,2))
     p, q = s
     A = @SMatrix rand(T, p, p)
     B = @SMatrix rand(T, q, q)
     C = @SMatrix rand(T, p, q)
 
-    X = sylv(A, B, C)
+    X, singular = sylv(A, B, C)
 
     @test A * X - X * B ≈ C
+    @test !singular
+end
+
+@testset "Singular sylvester $T" for T in (Float64, ComplexF64)
+    let
+        A = @SMatrix T[1 2; 0 1]
+        B = @SMatrix T[1 3; 0 1]
+        C = @SMatrix rand(T, 2, 2)
+        X, singular = sylv(A, B, C)
+        @test singular
+    end
+    let
+        A = @SMatrix T[1]
+        B = @SMatrix T[1 3; 0 1]
+        C = @SMatrix rand(T, 1, 2)
+        X, singular = sylv(A, B, C)
+        @test singular
+    end
+    let
+        A = @SMatrix T[1 2; 0 1]
+        B = @SMatrix T[1]
+        C = @SMatrix rand(T, 2, 1)
+        X, singular = sylv(A, B, C)
+        @test singular
+    end
 end