LAPACKException(22) for 32-bit storage and large matrices

[smataigne]

For sufficiently large matrices and 32-bit storage, LAPACKException(22) arises when eigen is called. eigen calls stegr! which throws the exception.

[smataigne]

[stevengj]

According to the dstegr documentation, 22 corresponds to error code -2 from dlarrv:

Problem in DLARRF when computing the RRR of a child. When a child is inside a tight cluster, it can be difficult to find an RRR. A partial remedy from the user's point of view is to make the parameter MINRGP smaller and recompile. However, as the orthogonality of the computed vectors is proportional to 1/MINRGP, the user should be aware that he might be trading in precision when he decreases MINRGP.

[stevengj]

Can you convert the matrix to Float64 and find the eigenvalues, and find the minimum separation?

[stevengj]

In particular, if you generate a random real skew-Hermitian matrix, the eigenvalues always come in ± pairs. If, on top of that, you have a near degeneracy, this will lead to a near fourfold degeneracy of the eigenvalues, which LAPACK might not be built to deal with?

[stevengj]

(If the problem is the near degeneracy, a simple trick would be to shift the tridiagonal matrix by e.g. the Frobenius norm to break up the ± pairs).

[smataigne]

I tried to shift the SymTridiagonal matrix as you proposed but it didn't work, the error is still thrown by stegr!. The following discussion also reports the problem https://bytemeta.vip/repo/Jutho/KrylovKit.jl/issues/57

[smataigne]

However, scaling the matrix by its norm(thus scaling its eigenvalues) did work instead of shifting its eigenvalues. Remains to know if scaling a matrix by its norm is a good idea. If the norm is to big, it could overflow, if it is less than one, it squeezes the eigenvalues instead of getting these far from each other.

(If the problem is the near degeneracy, a simple trick would be to shift the tridiagonal matrix by e.g. the Frobenius norm to break up the ± pairs).

[stevengj]

Scaling the matrix doesn't make any sense to me. The only thing that I can think of is that LAPACK has an erroneous absolute tolerance somewhere, which it really shouldn't — the algorithms to should be invariant under overall scaling.

I wouldn't do a multiplicative scaling, even if it happens to work for some mysterious reason, because it seems like just piling bugs on top of bugs.

[stevengj]

(It seems like, since we know all of the eigenvalues come in pairs, there should be a way to reduce it to an eigenproblem of half the size somehow.)

I would maybe look again through some of the papers on skew-symmetric eigensolvers in case they have any suggestions about the degenerate magnitudes.

[smataigne]

Scaling the matrix doesn't make any sense to me. The only thing that I can think of is that LAPACK has an erroneous absolute tolerance somewhere, which it really shouldn't — the algorithms to should be invariant under overall scaling.

I wouldn't do a multiplicative scaling, even if it happens to work for some mysterious reason, because it seems like just piling bugs on top of bugs.

Scaling didn't make any sense to me either and I don't have any idea why it allowed stegr! to compile. I will have a further look at possibilities to get rid of this LAPACKException

[smataigne]

The error only occurs for random SkewHermTridiagonal 32bit-storage matrices. Non-random matrices don't get the exception.
Moreover, the problem is created by the numbers inside the SkewHermTridiagonal(rand(Float32, n)) matrix. I mean by that that casting the matrix in Float64 when one calls stegr! doesn't fix the problem while the SkewHermTridiagonal(rand(Float64, n)) never creates errors.
This paper: https://netlib.org/lapack/lawnspdf/lawn163.pdf addresses the problem, but the solutions they bring ask to modify what stegr! does.

[stevengj]

It seems like this will have to be fixed in LAPACK itself.

[stevengj]

(Alternatively, you'd have to implement your own eigensolver, e.g. with a version of QR iterations, I guess. In principle this might be able to outperform the generic symtridiagonal eigensolver because you can take advantage of the diagonals being zero and the eigenvalues coming in ±λ pairs, maybe ala algorithm 530.)

[smataigne]

I think that trying to implement a fast eigensolver for tridiagonal skew-symmetric matrices could be a very nice challenge for my last week in Boston. As the package becomes quite complete, it is a very interesting problem to address, which is more on the "math side" of this package.

[smataigne]

#95 new skew-symmetric eigensolver was added. It is faster for tridiagonal matrices than the symmetric solver. However it cannot always work where LAPACK fails. For dense matrices, it is faster for matrices of size less than 200(on my machine), then it becomes slightly slower due to the hessenberg reduction that is slower. The skew-symmetric mv product is thus the last stone to beat LAPACK in single threaded.

[smataigne]

#95

[stevengj]

Can this issue be closed if you are no longer using the LAPACK solver?