Possible bug in KrylovKit.eigsolve

I am trying to find the first 4 eigenvalues of a symmetric 8X8 matrix, arranged in increasing order of real part. Symmetry means that the eigenvalues are real. Moreover, as you will see from the output, it is a block diagonal matrix with two identical blocks. the eigenvalues are pretty easy to calculate actuall.y

using Serialization;
using KrylovKit;
K = deserialize("K_problematic.txt");
display(K);
A =K-K'; print("\nnorm of K-K' is ", norm(A), ", hence K is symmetric\n");
(vals, vecs, info) = KrylovKit.eigsolve(K, 4, :SR);
display(vals');

However, I get the following error

8×8 view(::Array{Float64,2}, 1:8, 1:8) with eltype Float64:
  3.0  -1.0  -1.0  -1.0  -0.0  -0.0  -0.0  -0.0
 -1.0   3.0  -1.0  -1.0  -0.0  -0.0  -0.0  -0.0
 -1.0  -1.0   3.0  -1.0  -0.0  -0.0  -0.0  -0.0
 -1.0  -1.0  -1.0   3.0  -0.0  -0.0  -0.0  -0.0
 -0.0  -0.0  -0.0  -0.0   3.0  -1.0  -1.0  -1.0
 -0.0  -0.0  -0.0  -0.0  -1.0   3.0  -1.0  -1.0
 -0.0  -0.0  -0.0  -0.0  -1.0  -1.0   3.0  -1.0
 -0.0  -0.0  -0.0  -0.0  -1.0  -1.0  -1.0   3.0

norm of K-K' is 0.0, hence K is symmetric
┌ Warning: Invariant subspace of dimension 2 (up to requested tolerance `tol = 1.0e-12`), which is smaller than the number of requested eigenvalues (i.e. `howmany == 4`); setting `howmany = 2`.
└ @ KrylovKit ~/.julia/packages/KrylovKit/OLgKs/src/eigsolve/lanczos.jl:30
1×2 Adjoint{Float64,Array{Float64,1}}:
 0.0  4.0

A few days ago I pointed out a bug in Arpack and people got pretty upset about it. The conversation started to steer towards whether it counts as a bug in the package or the language. For me, that labelling is a technicality and not of any importance to an user. A user such as me is currently confused about how to compute eigenvalues of small matrices in Julia. I tried two packages and they failed. If I am making a mistake in the above commands, please let me know. If there is no mistake from my part, then it is a bug. If Julia has no packages that can handle eigenvalue computation of 8X8 symmetric matrices, then it adds up to a deficiency in Julia. Fighting over whether it is a shortcoming of Julia or a package wont help an user.

Okay to be fair, the first post you made was not very delicately formulated also. I agree that sometimes discussions here tend towards nitpicking the way questions are formulated, but you can’t go ahead and say things like

without people getting upset.

I’m sure for the matrices you are considering there is a good algebraic reason, why trying to find the eigenvalues of smallest real part is causing issues. But I’m not qualified to judge that. KrylovKit seems to give a reason already for the Matrix in this thread here.

If you’re looking at such small matrices, you can always use the standard library package LinearAlgebra and use

julia> using LinearAlgebra

julia> sort(real(eigvals(K)))[1:4]
4-element Vector{Float64}:
 -1.1102230246251565e-16
 -1.1102230246251565e-16
  3.9999999999999982
  3.9999999999999982
2 Likes

The matrix that you see is a Laplacian matrix, D-A associated with a graph. D is a diagonal matrix containing degrees, and A is the symmetric adjacency matrix containing 0s and 1s. Graph Laplacians arise all the time in data-analysis.

KrylovKit seems to give a reason already for the Matrix in this thread here.

Did you mean the warning message

Warning: Invariant subspace of dimension 2...not found

Please note that the KrylovKit documentation does not place any constraint on howmany. The problem is pretty well defined and simple from a mathematical point of view. When working with a general symmetric matrix, there is no way a user can detect these invariant spaces on his/her own. It is precisely the job of the algorithm.

So I disagree with you that KrylovKit provides a handy explanation. It merely tells the user of its own inadequacy. Lastly,

If you’re looking at such small matrices, you can always use the standard library package LinearAlgebra and use …

Thank you for the suggestion. But unfortunately, I was running algorithms on small matrices as a test run. I am finally interested in 1E5 or 1E6 dimensional matrices.

Okay so let’s look at your matrix here with pen and paper, i.e. SymPy.

The first block of your Matrix is the same as the second one, so both of these submatrices have the same eigenvalues already.

If you look at one submatrix only:

julia> using SymPy

julia> @syms λ
(λ,)

julia> A = Int.(K[1:4,1:4])
4×4 Matrix{Int64}:
  3  -1  -1  -1
 -1   3  -1  -1
 -1  -1   3  -1
 -1  -1  -1   3

julia> #characteristic polynomial:
       factor(det(A-one(A)*λ))
         3
λ⋅(λ - 4)

You see that it has only two eigenvalues λ = 0,4. Also, you already see a problem with the matrix as it is singular (det(K)=0).

So your problem is not

If you have trouble with matrices that are beyond what you present here, then you have to provide this and show that there should be no such problems for an algorithm.

Thank you very much for your time and the detailed work. I truly appreciate it :slight_smile:
However, according to my understanding, the singularity of the matrix does NOT affect the question of its eigenvalues being well defined.
Obviously, the eigenvalues are 0 6-times and 4 2-times. So I was expecting the first 4 smallest eigenvalues being reported as 0, and any 4 linearly independent vectors of the null space being reported as the eigenvectors.
In fact, in any eigenvalue problem where a λ has mutiplicity higher than 1, there is no unique choice of corresponding eigenvectors, the routine needs to converge to a spanning set of eigenvectors in the null space of K-λ. Note that singularity of K just means that λ=0 is an eigenvalue.

Suppose a is a number not in the spectrum of K. Then K-aI is definitely nonsingular. But it wouldnt change the nature or complexity of the eigenvalue problem. That is why I think that the singularity is not an issue.

Perhaps the choice :SR is itself problematic, as then one would need Krylov iterations for the inverse operator.

Okay so you have to be a bit clearer next time on what you are expecting to see and why. It was not obvious from your initial post.

For the issue: you’re expecting 2 eigenvalues of 0 and 6 eigenvalues of 4 for the example here.

when using eigen you already find that. Using KrylovKit.jl does result in just finding one of degenerate eigenpairs, i do not know why. Arpack.jl has a bug in finding the lowest real eigenvalues, as you have seen before. ArnoldiMethod.jl however, gives (I believe) exactly what you want:

julia> using ArnoldiMethod

julia> partialeigen(partialschur(K,nev=4,which=SR())[1])[1]
4-element Vector{Float64}:
 0.0
 2.0220868699519785e-16
 3.999999999999999
 4.0

although it sometimes returns only 3…
I guess the problem lies within the fact that K is singular for Arnoldi iterations.
If your matrices are not singular you can use ArnoldiMethod as it matches the performance of Arpack you’re good to go for your large problems!

The ArnoldiMethod ? I will try to learn about it, thanks.
didnt mention the expected result as my issue was with the warning message itself.

Author of KrylovKit here. For a symmetric matrix, the method that KrylovKit uses is known as the (implicitly restarted) Lanczos method, which is a special case of the Arnoldi method that explicitly uses the symmetric structure.

To answer your first questions. Clearly Julia or any other scientific language provides generic and stable methods to compute eigenvalues of symmetric 8 x 8 matrices, in fact it is provided via LAPACK. These methods compute all eigenvalues (although there are also specialized methods to compute a range of eigenvalues or specific eigenvalues in the case of symmetric matrices). All of these are accessible via Julia’s LinearAlgebra.eigen.

Once you have large problems of which you only want a few eigenvalues, you should know at least a bit about the kind of methods that you will use to solve this problem. Not to the extent that you know how to implement them, but at least such that you know the basic philosophy behind a method and when you can use it and when you can not. KrylovKit in particular uses methods based on Krylov subspaces, which means that you start from a given initial vector v, and generate a subspace spanned by v, A*v, A^2*v, A^3*v, and so forth. Within this subspace, an approximation to the eigenvalues and eigenvectors is sought. Generically, this works well for extremal eigenvalues, i.e. if you were to plot all the eigenvalues in the complex plane, Krylov methods are generically good (without shift-and-invert strategies) only for eigenvalues on the periphery of that set of points. For a hermitian / real symmetric matrix, these eigenvalues lie on the real line, and the extremal eigenvalues are the most positive and most negative ones.

What can go wrong is of course that v, if you were to decompose it into the basis of eigenvectors, has no component along a certain eigenvector. In that case, no component will be generated in the Krylov subspace, and so you can never recover that eigenvector in that Krylov subspace, even if it is for example the one of largest magnitude eigenvalue. This is mostly a theoretical concern, as numerical floating point errors will typically generate small (order machine precision) components which then grow in the Krylov subspace. However, for your example, if you start with a vector v which only has nonzero components in its first four rows, you will never generate components in the last four rows.

A related problem is degenerate eigenvalues. Suppose a certain eigenspace is two-dimensional. The starting vector has a certain component in that eigenspace, but that selects one direction in that eigenspace. A second linearly independent direction in that eigenspace cannot in theory not be generated in the Krylov subspace, and so in theory Krylov subspace methods cannot find degenerate eigenvalues. In practice, once the first eigenvector is well approximated and deflated, a second linearly independent component can again be generated from numerical noise, but it’s kind of bad practice to rely on this behaviour.

There is one aspect that sets KrylovKit apart from other packages such as Arpack. Packages like Arpack do not require that A is an explict matrix (that would often defeat the purpose of using them), but they do assume that all the vectors take the standard format of a Vector. Once the Krylov subspace v, A*v, A^2*v, … , A^k * v becomes an invariant subspace, this means that A^(k+1)* v is linearly dependent on the previous vectors, and so after orthogonalizing it with respect to the previous ones, the zero vector is obtained. Arpack would (could? I am not entirely sure it does this) then just continue with a new random vector, again orthogonalized against the already found invariant subspace, but this time yielding a new nonzero linearly independent direction.

However, KrylovKit does not assume even that vectors in the problem are standard Vector types. It does not even assume that they have a finite length (or thus, dimension), nor does it rely on being able to fill those vectors with random numbers or anything. This is useful in many cases. But as a consequence, once an invariant subspace is obtained, and A^(k+1) * v after orthogonalization has norm zero, there is no way for KrylovKit to automatically generate a new random direction. And thus, it just stops and prints a nice warning.

14 Likes

Krylov methods must work with blocks of vectors to resolve geometric multiplicity.

I discovered yet another problem with the algorithm. I took a 16X16 block diagonal, real symmetric matrix in which each block is 4X4 and has 1 as an eigenvalu. The top eigenvalue is 1, with multiplicity 4. However krylovkit only detects the multiplicity of 1 to be two. So it means that krylovkit may not be able to detect all the top eigenvalues too. If you wish, I can provide you with the code and output.

Thank you. Although I am sure I fully explained and displayed the problem with the package in my initial post.