I am working on a problem which involves the calculation of the eigenvalues (specifically the spectral radius) of matrices with a particular structure. These are specifically sparse, asymmetric matrices with binary entries (either 1.0 or 0.0) and only zeroes on the leading diagonal. I am able to calculate these eigenvalues using Arpack.eigs() for matrices of size 20x20 quite easily, and have been attempting to look at larger ones (128x128).
However, I have run into serious problems of accuracy/errors. By way of example (though I don’t think this is the only case in which the code generates inaccurate results), there is a subset of these which are nilpotent and I know analytically should have a spectrum of entirely zeroes, for which:
• Arpack.eigs() either inaccurately calculates small (but not very very small, i.e. of absolute value around 0.05) eigenvalues, or raises the error
Error: XYAUPD_Exception: Maximum number of iterations taken. All possible eigenvalues of OP has been found.
│ IPARAM(5) returns the number of wanted converged Ritz values.
Even when maxiter is raised to, e.g., 10000
• ArnoldiMethod and KrylovKit.eigsolve calculate similarly small but nontrivial eigenvalues
I have checked for different permutations of typing (e.g. trying the matrices with both float and int entries, and in both sparse and dense forms), and have seen no improvement.
I have managed to obtain correct - or at least sufficiently precise - results for these matrices using SciPy.jl, but it feels inelegant to use a python interface when Julia has native eigen-solvers, and it doesn’t explain what’s gone wrong. So why is Julia failing where python succeeds? Is there a technical aspect of calling these functions (e.g. similar to increasing maxiter) that I’ve missed that is necessary for precise handling of this kind of matrices, or is this a bug/limitation of the available Julia implementations?
These are tiny matrices. I would normally just use a dense eigensolver for matrices this small.
Are you trying to compute all the eigenvalues? For the spectral radius you only want the eigenvalue of the largest magnitude, no? How exactly are you calling eigs?
Conversely, if you are looking for the smallest magnitude eigenvalues, be aware that the Julia libraries require you to supply the matrix factorization yourself for inverse iteration, whereas Scipy may do this for you (which is convenient but a lot less flexible).
(For iterative eigensolvers I would tend to use KrylovKit.jl rather than Arpack these days.)
Eigenvalues can be notoriously ill-conditioned (i.e., sensitive to roundoff). A Jordan block is a well-known example, though defectiveness is not the essence of the matter.
When the answer can be so sensitive, it is usually not the right target for numerical computation (and in many applications, as well). The book Spectra and Pseudospectra is a detailed examination of the phenomenon.
I thought that defectiveness is precisely the essence of the matter — the condition number of an eigenvalue diverges if and only if it is becoming defective with another eigenvalue (i.e. two or more corresponding eigenvectors are becoming linearly dependent, i.e. it’s approaching a matrix with a ≥ 2x2 Jordan block for that eigenvalue).
In any case, it’s hard to say much about the OP’s problem without knowing more. (Having the extremal eigenvalue be nearly defective rarely happens by accident.)
If it’s a 128x128 sparse matrix, that’s tiny enough that the OP could just post it. Just do println(findnz(A))
Along these lines, similar intuition holds for several problems. There’s a nice paper by Demmel relating ill-conditioning to distance to an ill-posed problem for linear systems, eigenvalue problems, and polynomial zeros. The proofs are done for all these in a somewhat uniform way using differential inequaltities.
Yes, but in the case where an eigenvector basis is ill-conditioned, the same Bauer–Fike theorem implies that it only requires a small perturbation (relative to the norm of the matrix) to coalesce two eigenvalues and make the matrix defective, no?
Yes, when I said “large” in the title what I really meant was “problem emerges as the matrix gets larger” - I was if anything surprised things went wrong for such a comparatively small size. You’re right, my matrices came out of other functions in the code as sparse but I must have skipped over the dense solvers e.g. in LinearAlgebra, which do work as expected (which I could have sworn I’d already checked!! Just goes to show one should never be too dismissive of the trivial case). I think that would answer my main question insofar as how SciPy managed to outperform the methods I had tried. Definitely a case of asking for help and the answer immediately becoming obvious.
I ran some further checks and for Arpack it was only the nilpotent subset of matrices that had an error; I’m sufficiently convinced that this would be the property that caused issue. KrylovKit was inconsistent for several other cases (and even more inconsistent in the nilpotent case) but I think this is probably a case of trying to use a heavy tool for too simple of a problem (indeed I had only initially looked at this once Arpack started failing).
The ill-conditioned eigenvector basis is also associated with proximity to a defective matrix, as @stevengj was faster than me in noting. To quantify it for a randomly generated matrix, consider the following Julia code
using LinearAlgebra
n = 50
T = UpperTriangular(randn(n,n))
e, X = eigen(T)
Y = inv(X)
x1 = X[:, 1]
y1 = Y[1, :]
norm_r = norm(y1[2:n])
@show bound = 2*norm(T) / norm_r
This shows the bound on the distance to a matrix with repeated eigenvalues from Theorem 5 from the paper by Demmel that I linked to. Running it, I got