3x3 matrix takes too many QR iterations

[haampie]

This 3x3 matrix with (nearly) repeated eigenvalues takes an astonishing 2700 iterations of the QR algorithm:

julia> H = [-9.000000022477964 -4.51417416984849e-5 -0.03599511897656601; 2.555189460215141e-12 -9.00000552392542 -0.004386727237710469; 0.0 6.956006220233481e-9 -8.999994476095713]

julia> wrap(H::AbstractMatrix{T}) where T = GenericLinearAlgebra.HessenbergFactorization{T,typeof(H),Vector{GenericLinearAlgebra.Householder{T}}}(H, [])

julia> HH = wrap(copy(H)); GenericLinearAlgebra._schur!(HH, maxiter=100000)

Haven't checked yet how LAPACK handles it.

[haampie]

I can bring it down to 19 iterations by doing a single shift Wilkinson shift.

LAPACK seems to also always do a double shift in case of real arithmetic? I've never fully understood it.

In ArnoldiMethod.jl my implementation was to do double shift if complex conjugate, and single shift with eigenvalues closest to bottom right corner matrix value if real. That seems to work alright -- maybe by accident.

[andreasnoack]

Are you using the latest version? If you look at line 139 in EigenGeneral.jl then you can see I've introduced an exceptional shift for every tenth iteration. I thought it would be sufficient to avoid cycles but maybe it's not sufficient.

[haampie]

Tried commit 8255d4e

Now getting 577 iterations. I believe I haven't changed a thing, maybe I didn't properly repr(H) when making the issue.

But 577 is still a lot for a 3x3 matrix :p.

IIUC there are two differences compared to reference LAPACK that may contribute to this:

• They use a double shift with the same eigenvalue repeated (god knows why), instead of shifting both eigenvalues of the 2x2 block separately
• They use scaling tricks when computing the column (H-λ₁I)(H-λ₂I)e₁ -- could check if there's catastrophic cancellation here, cause I hit that earlier with similar matrices.

However, I still have to check if LAPACK struggles with this matrix too or not.

[andreasnoack]

However, I still have to check if LAPACK struggles with this matrix too or not.

@haampie did you reach a conclusion here?