Convergence problems of some matrices

[johanbluecreek]

Hi,

I'm encountering some matrices that fail to converge when computing eigenvalues. Consider the following:

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]

[johanbluecreek]

Ah, I also wanted to mention that issue #29 is similar and may be related, but that one is fixed.

[andreasnoack]

Thanks for reporting this and sorry for the long wait. I had to revisit the Sandia report on convergence failures for the general eigenvalue problem and I realized that I'd misunderstood the conclusion. Hence, I'm changing the exceptional shift strategy in #69 and it seems to work well.

However, the matrices you are sharing seem very stressful in general when solving with BigFloats. The required number of iterations for Float64 seems to be reasonable based on my reading of the Sandia report. I guess the problem isn't that well covered when you go beyond Float64 so I don't know if there is a way to improve here. Bottom line is that the package can now handle all the matrices that you shared but it's necessary to increase the number of iterations.