Repository: JuliaLinearAlgebra/GenericLinearAlgebra.jl (julia)

Issue:
Convergence problems of some matrices
Hi,

I'm encountering some matrices that fail to converge when computing eigenvalues. Consider the following:

```julia
julia> t = [1 -2 1 -1 -1 0; 0 1 0 1 0 1; 1 -1 2 0 -1 0; 0 1 0 2 1 1; 1 0 1 0 0 0; 0 -1 1 -1 -2 0]
6×6 Array{Int64,2}:
 1  -2  1  -1  -1  0
 0   1  0   1   0  1
 1  -1  2   0  -1  0
 0   1  0   2   1  1
 1   0  1   0   0  0
 0  -1  1  -1  -2  0

julia> using LinearAlgebra

julia> using GenericLinearAlgebra

julia> eigvals(t)
6-element Array{Complex{Float64},1}:
 -7.294730590901527e-9 + 0.0im
 7.2947305672186945e-9 + 0.0im
    1.4999999999999996 - 0.8660254037844369im
    1.4999999999999996 + 0.8660254037844369im
    1.5000000000000018 - 0.8660254037844406im
    1.5000000000000018 + 0.8660254037844406im

julia> eigvals(big.(t))
ERROR: ArgumentError: iteration limit 600 reached
Stacktrace:
 [1] _schur!(::GenericLinearAlgebra.HessenbergFactorization{BigFloat,Array{BigFloat,2},GenericLinearAlgebra.Householder{BigFloat,S} where S<:(Union{DenseArray{T,1}, Base.ReinterpretArray{T,1,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, Base.ReshapedArray{T,1,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray}, SubArray{T,1,A,I,L} where L where I<:Tuple{Vararg{Union{Int64, AbstractRange{Int64}, Base.AbstractCartesianIndex},N} where N} where
A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, Base.ReshapedArray{T,N,A,MI} where MI<:Tuple{Vararg{Base.MultiplicativeInverses.SignedMultiplicativeInverse{Int64},N} where N} where A<:Union{Base.ReinterpretArray{T,N,S,A} where S where A<:Union{SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, SubArray{T,N,A,I,true} where I<:Union{Tuple{Vararg{Real,N} where N}, Tuple{AbstractUnitRange,Vararg{Any,N} where N}} where A<:DenseArray where N where T, DenseArray} where N where T, DenseArray}} where T)}; tol::BigFloat, debug::Bool, shiftmethod::Symbol, maxiter::Int64, kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/user/.julia/packages/GenericLinearAlgebra/sJLhE/src/eigenGeneral.jl:103
 [2] _schur! at /home/user/.julia/packages/GenericLinearAlgebra/sJLhE/src/eigenGeneral.jl:91 [inlined]
 [3] #_schur!#21 at /home/user/.julia/packages/GenericLinearAlgebra/sJLhE/src/eigenGeneral.jl:183 [inlined]
 [4] _schur! at /home/user/.julia/packages/GenericLinearAlgebra/sJLhE/src/eigenGeneral.jl:183 [inlined]
 [5] _eigvals!(::Array{BigFloat,2}; kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/user/.julia/packages/GenericLinearAlgebra/sJLhE/src/eigenGeneral.jl:260
 [6] _eigvals! at /home/user/.julia/packages/GenericLinearAlgebra/sJLhE/src/eigenGeneral.jl:260 [inlined]
 [7] eigvals!(::Array{BigFloat,2}; sortby::typeof(LinearAlgebra.eigsortby), kwargs::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}) at /home/user/.julia/packages/GenericLinearAlgebra/sJLhE/src/eigenGeneral.jl:274
 [8] eigvals! at /home/user/.julia/packages/GenericLinearAlgebra/sJLhE/src/eigenGeneral.jl:271 [inlined]
 [9] #eigvals#81 at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/eigen.jl:326 [inlined]
 [10] eigvals(::Array{BigInt,2}) at /buildworker/worker/package_linux64/build/usr/share/julia/stdlib/v1.4/LinearAlgebra/src/eigen.jl:326
 [11] top-level scope at REPL[33]:1


julia> GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(copy(BigFloat.(t)); maxiter=10000))
ERROR: ArgumentError: iteration limit 10000 reached
[stacktrace]

julia> GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(copy(BigFloat.(t)); debug=true))
Split! Subdiagonal element is:  3.144e-77 and istart now      4
block start is:      4, block end is:      6, d:  8.096e-02, t:  1.663e+00
Francis double shift! Subdiagonal is:  1.144e+00, last subdiagonal is:  8.664e-01
Split! Subdiagonal element is:  0.000e+00 and istart now      4
block start is:      4, block end is:      6, d:  1.835e-02, t:  1.380e+00
Francis double shift! Subdiagonal is: -7.721e-02, last subdiagonal is:  7.052e-01
Split! Subdiagonal element is:  0.000e+00 and istart now      4
block start is:      4, block end is:      6, d:  5.610e-06, t:  1.509e+00
Francis double shift! Subdiagonal is: -7.490e-04, last subdiagonal is:  1.118e+00
Split! Subdiagonal element is:  0.000e+00 and istart now      4
block start is:      4, block end is:      6, d: -7.790e-10, t:  1.665e+00
Francis double shift! Subdiagonal is: -3.435e-09, last subdiagonal is:  7.440e-01
Split! Subdiagonal element is:  0.000e+00 and istart now      4
block start is:      4, block end is:      6, d: -4.071e-19, t:  1.386e+00
Francis double shift! Subdiagonal is: -1.804e-18, last subdiagonal is:  7.020e-01
Split! Subdiagonal element is:  0.000e+00 and istart now      4
block start is:      4, block end is:      6, d: -9.812e-39, t:  1.469e+00
Francis double shift! Subdiagonal is:  3.850e-37, last subdiagonal is:  1.116e+00
Split! Subdiagonal element is:  0.000e+00 and istart now      4
block start is:      4, block end is:      6, d:  7.081e-76, t:  1.671e+00
Francis double shift! Subdiagonal is: -3.069e-75, last subdiagonal is:  7.506e-01
Split! Subdiagonal element is: 1.628e-150 and istart now      6
Bottom deflation! Block size is one. New iend is      5
Split! Subdiagonal element is:  0.000e+00 and istart now      4
Bottom deflation! Block size is two. New iend is      3
block start is:      1, block end is:      3, d:  1.000e+00, t:  2.000e+00
Wilkinson-like shift! Subdiagonal is: -1.000e+00, last subdiagonal is: -1.414e+00
block start is:      1, block end is:      3, d:  1.000e+00, t:  2.000e+00
Francis double shift! Subdiagonal is: -7.071e-01, last subdiagonal is: -1.000e+00
block start is:      1, block end is:      3, d:  1.000e+00, t:  2.000e+00
Francis double shift! Subdiagonal is: -1.000e+00, last subdiagonal is:  1.414e+00
block start is:      1, block end is:      3, d:  1.000e+00, t:  2.000e+00
Francis double shift! Subdiagonal is:  1.414e+00, last subdiagonal is: -7.071e-01
block start is:      1, block end is:      3, d:  1.000e+00, t:  2.000e+00
Francis double shift! Subdiagonal is: -7.071e-01, last subdiagonal is: -1.000e+00
block start is:      1, block end is:      3, d:  1.000e+00, t:  2.000e+00
Francis double shift! Subdiagonal is: -1.000e+00, last subdiagonal is:  1.414e+00
[more of the same type of debug messages]
Francis double shift! Subdiagonal is: -1.414e+00, last subdiagonal is: -7.071e-01
block start is:      1, block end is:      3, d:  1.000e+00, t:  2.000e+00
Francis double shift! Subdiagonal is: -7.071e-01, last subdiagonal is:  1.000e+00
block start is:      1, block end is:      3, d:  1.000e+00, t:  2.000e+00
Francis double shift! Subdiagonal is:  1.000e+00, last subdiagonal is: -1.414e+00
block start is:      1, block end is:      3, d:  1.000e+00, t:  2.000e+00
Wilkinson-like shift! Subdiagonal is: -1.414e+00, last subdiagonal is: -7.071e-01
ERROR: ArgumentError: iteration limit 600 reached
[stacktrace]
```

The above seems to show that indeed it is stuck in some infinite loop.

The analytical eigenvalues are `[0,0,1/2*(3-√3im),1/2*(3+√3im),1/2*(3-√3im),1/2*(3+√3im)]`, and this only seems to happen with non-diagonalizable matrices.

These are other matrices I have so far found that displays a similar behaviour:

```
[1 -2 1 -1 -1 0; 0 1 0 1 0 1; 1 -1 2 0 -1 0; 0 1 0 2 1 1; 1 0 1 0 0 0; 0 -1 1 -1 -2 0]
[1 0 -1 0 0 0; 0 1 1 1 -1 0; 1 0 -1 0 0 0; 1 -1 0 -1 1 1; 0 0 1 0 0 0; -1 0 1 0 0 0]
[1 -2 -1 1 -1 0; 0 1 0 -1 0 1; -1 1 2 0 1 0; 0 -1 -2 2 -1 -1; 1 0 -1 0 0 0; 0 -1 -1 1 0 0]
[2 0 -1 0 0 1; 0 2 0 -1 -1 0; -1 0 1 0 0 -1; 0 -1 0 1 1 0; 0 1 0 -1 0 0; -1 0 1 0 0 0]
[1 0 1 0 -1 0; -2 1 2 -1 1 1; -1 0 -1 0 1 0; 2 1 2 -1 -2 1; 1 0 1 0 -1 0; 1 -1 -2 1 0 -1]
```

Pull request:
Also reduce the default number of iterations. Add complicated
test matrices.

Fixes https://github.com/JuliaLinearAlgebra/GenericLinearAlgebra.jl/issues/67

Patch:
diff --git a/.travis.yml b/.travis.yml
index 67b441a..998992d 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -4,7 +4,6 @@ os:
   - linux
 
 julia:
-  - 1.0
   - 1
   - nightly
 
diff --git a/Project.toml b/Project.toml
index d283ee0..cbf80cb 100644
--- a/Project.toml
+++ b/Project.toml
@@ -15,4 +15,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 test = ["Quaternions", "Test"]
 
 [compat]
-julia = "1"
+julia = "1.3"
diff --git a/azure-pipelines.yml b/azure-pipelines.yml
deleted file mode 100644
index 8ca143a..0000000
--- a/azure-pipelines.yml
+++ /dev/null
@@ -1,70 +0,0 @@
-# Starter pipeline
-# Start with a minimal pipeline that you can customize to build and deploy your code.
-# Add steps that build, run tests, deploy, and more:
-# https://aka.ms/yaml
-
-jobs:
-- job: Linux
-
-  pool:
-    vmImage: 'Ubuntu 16.04'
-
-  strategy:
-    matrix:
-      Julia 1.0:
-        JULIA_VERSION: '1.0.1'
-
-  steps:
-  - script: |
-      wget -nv https://julialang-s3.julialang.org/bin/linux/x64/1.0/julia-$(JULIA_VERSION)-linux-x86_64.tar.gz
-      tar zxf julia-$(JULIA_VERSION)-linux-x86_64.tar.gz
-    displayName: 'Download and extract Julia 1.0'
-
-  - script: |
-      ./julia-$(JULIA_VERSION)/bin/julia --project=@. -e 'using Pkg; Pkg.instantiate()'
-      ./julia-$(JULIA_VERSION)/bin/julia --project=@. -e 'using Pkg; Pkg.test("GenericLinearAlgebra")'
-    displayName: 'Run the tests'
-
-- job: macOS
-
-  pool:
-    vmImage: 'macOS-10.13'
-
-  strategy:
-    matrix:
-      Julia 1.0:
-        JULIA_VERSION: '1.0.1'
-        JULIA_VERSION_SHORT: '1.0'
-
-  steps:
-  - script: |
-      wget -nv https://julialang-s3.julialang.org/bin/mac/x64/1.0/julia-$(JULIA_VERSION)-mac64.dmg
-      hdiutil attach julia-$(JULIA_VERSION)-mac64.dmg
-      ln -s /Volumes/Julia-$(JULIA_VERSION)/Julia-$(JULIA_VERSION_SHORT).app/Contents/Resources/julia/bin/julia
-    displayName: 'Download and extract Julia 1.0'
-
-  - script: |
-      ./julia --project=@. -e 'using Pkg; Pkg.instantiate()'
-      ./julia --project=@. -e 'using Pkg; Pkg.test("GenericLinearAlgebra")'
-    displayName: 'Run the tests'
-
-- job: Windows
-
-  pool:
-    vmImage: 'VS2017-Win2016'
-
-  strategy:
-    matrix:
-      Julia 1.0:
-        JULIA_VERSION: '1.0.1'
-
-  steps:
-  - powershell: |
-      wget https://julialang-s3.julialang.org/bin/winnt/x64/1.0/julia-$(JULIA_VERSION)-win64.exe -OutFile julia-$(JULIA_VERSION)-win64.exe
-      Start-Process -FilePath .\julia-$(JULIA_VERSION)-win64.exe -ArgumentList "/S /D=C:\julia" -NoNewWindow -Wait
-    displayName: 'Download and extract Julia 1.0'
-
-  - powershell: |
-      C:\julia\bin\julia --project=@. -e 'using Pkg; Pkg.instantiate()'
-      C:\julia\bin\julia --project=@. -e 'using Pkg; Pkg.test(\"GenericLinearAlgebra\")'
-    displayName: 'Run the tests'
diff --git a/src/eigenGeneral.jl b/src/eigenGeneral.jl
index 061b5b2..f1f7018 100644
--- a/src/eigenGeneral.jl
+++ b/src/eigenGeneral.jl
@@ -85,7 +85,7 @@ function _schur!(
     tol = eps(real(T)),
     debug = false,
     shiftmethod = :Francis,
-    maxiter = 100*size(H, 1),
+    maxiter = 30*size(H, 1),
     kwargs...) where T
 
     n = size(H, 1)
@@ -99,6 +99,7 @@ function _schur!(
 
     @inbounds while true
         i += 1
+        debug && println("iteration: $i")
         if i > maxiter
             throw(ArgumentError("iteration limit $maxiter reached"))
         end
@@ -137,35 +138,20 @@ function _schur!(
         else
             Hmm = HH[iend, iend]
             Hm1m1 = HH[iend - 1, iend - 1]
-            d = Hm1m1*Hmm - HH[iend, iend - 1]*HH[iend - 1, iend]
-            t = Hm1m1 + Hmm
-            t = iszero(t) ? eps(real(one(t))) : t # introduce a small pertubation for zero shifts
+            if iszero(i % 10)
+                # Use Eispack exceptional shift
+                β = abs(HH[iend, iend - 1]) + abs(HH[iend - 1, iend - 2])
+                d = (Hmm + β)^2 - Hmm*β/2
+                t = 2*Hmm + 3*β/2
+            else
+                d = Hm1m1*Hmm - HH[iend, iend - 1]*HH[iend - 1, iend]
+                t = Hm1m1 + Hmm
+            end
             debug && @printf("block start is: %6d, block end is: %6d, d: %10.3e, t: %10.3e\n", istart, iend, d, t)
 
             if shiftmethod == :Francis
-                # Run a bulge chase
-                if iszero(i % 10)
-                    # Vary the shift strategy to avoid dead locks
-                    # We use a Wilkinson-like shift as suggested in "Sandia technical report 96-0913J: How the QR algorithm fails to converge and how fix it".
-
-                    debug && @printf("Wilkinson-like shift! Subdiagonal is: %10.3e, last subdiagonal is: %10.3e\n", HH[iend, iend - 1], HH[iend - 1, iend - 2])
-                    _d = t*t - 4d
-                    if _d isa Real && _d >= 0
-                        # real eigenvalues
-                        a = t/2
-                        b = sqrt(_d)/2
-                        s = a > Hmm ? a - b : a + b
-                    else
-                        # complex case
-                        s = t/2
-                    end
-                    singleShiftQR!(HH, τ, s, istart, iend)
-                else
-                    # most of the time use Francis double shifts
-
-                    debug && @printf("Francis double shift! Subdiagonal is: %10.3e, last subdiagonal is: %10.3e\n", HH[iend, iend - 1], HH[iend - 1, iend - 2])
-                    doubleShiftQR!(HH, τ, t, d, istart, iend)
-                end
+                debug && @printf("Francis double shift! Subdiagonal is: %10.3e, last subdiagonal is: %10.3e\n", HH[iend, iend - 1], HH[iend - 1, iend - 2])
+                doubleShiftQR!(HH, τ, t, d, istart, iend)
             elseif shiftmethod == :Rayleigh
                 debug && @printf("Single shift with Rayleigh shift! Subdiagonal is: %10.3e\n", HH[iend, iend - 1])
 


Test patch:
diff --git a/test/eigengeneral.jl b/test/eigengeneral.jl
index 7df846d..91ae90b 100644
--- a/test/eigengeneral.jl
+++ b/test/eigengeneral.jl
@@ -93,4 +93,61 @@ end
     end
 end
 
+Demmel(η) = [0  1 0 0
+             1  0 η 0
+             0 -η 0 1
+             0  0 1 0]
+
+@testset "Demmel matrix" for t in (1e-10, 1e-9, 1e-8)
+    # See "Sandia technical report 96-0913J: How the QR algorithm fails to converge and how fix it"
+    A = Demmel(t)
+    vals = GenericLinearAlgebra._eigvals!(GenericLinearAlgebra._schur!(A, maxiter=35))
+    @test abs.(vals) ≈ ones(4)
 end
+
+function Hevil2(θ, κ, α, γ)
+    # Eq (13) and (14)
+    β = ω = 0.0
+    ν = cos(θ)*cos(2γ) + cos(α + β + ω)*sin(2γ)*κ/2
+    σ = 1 + κ*sin(2γ)*cos(α + β + ω - θ) + κ^2*sin(γ)^2
+    μ = -sin(θ)*cos(2γ) - sin(α + β + ω)*sin(2γ)*κ/2
+    ρ = sqrt(σ - ν^2)
+
+    return [ν    (cos(2θ) - ν^2)/ρ              μ/ρ*(cos(2θ) - ν^2 + ρ^2)/sqrt(ρ^2 - μ^2)      (-2*μ*ν - sin(2*θ))/sqrt(ρ^2 - μ^2)
+            ρ   -ν - sin(2*θ)*μ/ρ^2            -μ/ρ^2*(μ*sin(2*θ) + 2*ν*ρ^2)/sqrt(ρ^2 - μ^2)    -μ/ρ*(cos(2*θ) - ν^2 + ρ^2)/sqrt(ρ^2 - μ^2)
+            0    sin(2*θ)*sqrt(ρ^2 - μ^2)/ρ^2   ν + sin(2*θ)*μ/ρ^2                              (cos(2*θ) - ν^2)/ρ
+            0    0                              ρ                                               -ν]
+end
+
+@testset "Complicate matrix from Sandia technical report" begin
+    H = Hevil2(0.111866322512629152, 1.08867072154101741, 0.338146383137297168, -0.313987810419091240)
+
+    @test Float64.(abs.(eigvals(big.(H)))) ≈ ones(4)
+end
+
+@testset "Issue 67" for (A, λs) in (
+    ([1 -2 1 -1 -1 0; 0 1 0 1 0 1; 1 -1 2 0 -1 0; 0 1 0 2 1 1; 1 0 1 0 0 0; 0 -1 1 -1 -2 0]    ,
+        [0.0, 0.0, (3 - √big(3)*im)/2, (3 + √big(3)*im)/2, (3 - √big(3)im)/2, (3 + √big(3)im)/2]),
+    ([1 0 -1 0 0 0; 0 1 1 1 -1 0; 1 0 -1 0 0 0; 1 -1 0 -1 1 1; 0 0 1 0 0 0; -1 0 1 0 0 0]      ,
+        zeros(6)),
+    ([1 -2 -1 1 -1 0; 0 1 0 -1 0 1; -1 1 2 0 1 0; 0 -1 -2 2 -1 -1; 1 0 -1 0 0 0; 0 -1 -1 1 0 0],
+        [(1 - √big(3)*im)/2, (1 + √big(3)*im)/2, (1 - √big(3)im)/2, (1 + √big(3)im)/2, 2, 2]),
+    ([2 0 -1 0 0 1; 0 2 0 -1 -1 0; -1 0 1 0 0 -1; 0 -1 0 1 1 0; 0 1 0 -1 0 0; -1 0 1 0 0 0]    ,
+        ones(6)),
+    ([1 0 1 0 -1 0; -2 1 2 -1 1 1; -1 0 -1 0 1 0; 2 1 2 -1 -2 1; 1 0 1 0 -1 0; 1 -1 -2 1 0 -1] ,
+        [-1, -1, 0, 0, 0, 0]))
+
+    @testset "shouldn't need excessive iterations (30*n) in the Float64 case" begin
+        GenericLinearAlgebra._schur!(float(A))
+    end
+
+    # For BigFloats, many iterations are required for convergence
+    # Improving this is probably a hard problem
+    vals = eigvals(big.(A), maxiter=1500)
+
+    # It's hard to compare complex conjugate pairs so we compare the real and imaginary parts separately
+    @test sort(real(vals)) ≈ sort(real(λs)) atol=1e-25
+    @test sort(imag(vals)) ≈ sort(imag(λs)) atol=1e-25
+end
+
+end
\ No newline at end of file