# Possible bug in KrylovKit.eigsolve

**URL:** <https://discourse.julialang.org/t/possible-bug-in-krylovkit-eigsolve/55962>\
**Category:** Performance\
**Created:** [February 24, 2021, 9:34pm UTC](https://discourse.julialang.org/t/possible-bug-in-krylovkit-eigsolve/55962 "2021-02-24T21:34:01Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![iamsuddhasattwa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/iamsuddhasattwa/32/7441_2.png) [@iamsuddhasattwa](https://discourse.julialang.org/u/iamsuddhasattwa)\
**Post date:** [February 24, 2021, 9:34pm UTC](https://discourse.julialang.org/t/possible-bug-in-krylovkit-eigsolve/55962/1 "2021-02-24T21:34:01Z")

</div>

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

```julia
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

```julia
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](https://discourse.julialang.org/t/bug-detected-in-arpack/55334/4) 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.

---

<div class="post-metadata">

**Author:** ![fgerick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fgerick/32/13228_2.png) [@fgerick](https://discourse.julialang.org/u/fgerick)\
**Post date:** [February 24, 2021, 9:47pm UTC](https://discourse.julialang.org/t/possible-bug-in-krylovkit-eigsolve/55962/2 "2021-02-24T21:47:58Z")

</div>

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

> [@Bug detected in Arpack](https://discourse.julialang.org/t/bug-detected-in-arpack/55334/1):
>
> A language meant to do high performance computing should be able to find 4 eigenvalues of a 16X16 symmetric matrix. I could do it myself on pen and paper.

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
julia> using LinearAlgebra

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

```

---

<div class="post-metadata">

**Author:** ![iamsuddhasattwa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/iamsuddhasattwa/32/7441_2.png) [@iamsuddhasattwa](https://discourse.julialang.org/u/iamsuddhasattwa)\
**Post date:** [February 24, 2021, 10:06pm UTC](https://discourse.julialang.org/t/possible-bug-in-krylovkit-eigsolve/55962/3 "2021-02-24T22:06:01Z")

</div>

> [@fgerick](#):
>
> 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.

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

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

```

Please note that the KrylovKit [documentation](https://jutho.github.io/KrylovKit.jl/v0.1/man/eig/) 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.

---

<div class="post-metadata">

**Author:** ![fgerick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fgerick/32/13228_2.png) [@fgerick](https://discourse.julialang.org/u/fgerick)\
**Post date:** [February 24, 2021, 10:12pm UTC](https://discourse.julialang.org/t/possible-bug-in-krylovkit-eigsolve/55962/4 "2021-02-24T22:12:18Z")

</div>

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
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

> [@iamsuddhasattwa](#):
>
> pretty well defined

---

<div class="post-metadata">

**Author:** ![fgerick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fgerick/32/13228_2.png) [@fgerick](https://discourse.julialang.org/u/fgerick)\
**Post date:** [February 24, 2021, 10:14pm UTC](https://discourse.julialang.org/t/possible-bug-in-krylovkit-eigsolve/55962/5 "2021-02-24T22:14:45Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![iamsuddhasattwa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/iamsuddhasattwa/32/7441_2.png) [@iamsuddhasattwa](https://discourse.julialang.org/u/iamsuddhasattwa)\
**Post date:** [February 24, 2021, 10:25pm UTC](https://discourse.julialang.org/t/possible-bug-in-krylovkit-eigsolve/55962/6 "2021-02-24T22:25:45Z")

</div>

Thank you very much for your time and the detailed work. I truly appreciate it 🙂  
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.

---

<div class="post-metadata">

**Author:** ![fgerick](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/fgerick/32/13228_2.png) [@fgerick](https://discourse.julialang.org/u/fgerick)\
**Post date:** [February 24, 2021, 10:38pm UTC](https://discourse.julialang.org/t/possible-bug-in-krylovkit-eigsolve/55962/7 "2021-02-24T22:38:00Z")

</div>

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
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!

---

<div class="post-metadata">

**Author:** ![iamsuddhasattwa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/iamsuddhasattwa/32/7441_2.png) [@iamsuddhasattwa](https://discourse.julialang.org/u/iamsuddhasattwa)\
**Post date:** [February 24, 2021, 10:49pm UTC](https://discourse.julialang.org/t/possible-bug-in-krylovkit-eigsolve/55962/8 "2021-02-24T22:49:12Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![juthohaegeman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/juthohaegeman/32/8620_2.png) [@juthohaegeman](https://discourse.julialang.org/u/juthohaegeman)\
**Post date:** [February 25, 2021, 12:01am UTC](https://discourse.julialang.org/t/possible-bug-in-krylovkit-eigsolve/55962/9 "2021-02-25T00:01:22Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![PetrKryslUCSD](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/petrkryslucsd/32/215825_2.png) [@PetrKryslUCSD](https://discourse.julialang.org/u/PetrKryslUCSD)\
**Post date:** [February 25, 2021, 12:52am UTC](https://discourse.julialang.org/t/possible-bug-in-krylovkit-eigsolve/55962/10 "2021-02-25T00:52:11Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![iamsuddhasattwa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/iamsuddhasattwa/32/7441_2.png) [@iamsuddhasattwa](https://discourse.julialang.org/u/iamsuddhasattwa)\
**Post date:** [March 10, 2021, 8:33pm UTC](https://discourse.julialang.org/t/possible-bug-in-krylovkit-eigsolve/55962/11 "2021-03-10T20:33:37Z")

</div>

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.

---

<div class="post-metadata">

**Author:** ![iamsuddhasattwa](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/iamsuddhasattwa/32/7441_2.png) [@iamsuddhasattwa](https://discourse.julialang.org/u/iamsuddhasattwa)\
**Post date:** [March 10, 2021, 8:36pm UTC](https://discourse.julialang.org/t/possible-bug-in-krylovkit-eigsolve/55962/12 "2021-03-10T20:36:06Z")

</div>

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