Eigenvalues of rank deficient matrices are not sorted correctly.

[jpfairbanks]

If L is rank deficient (it is a combinatorial laplacian of a graph) then the 0 eigenvalue appears with the LR() eigenvalues.

For example:

julia> L = sparse(lapl)
10×10 SparseMatrixCSC{Float64,Int64} with 40 stored entries:
  [1 ,  1]  =  2.0
  [5 ,  1]  =  -1.0
  [6 ,  1]  =  -1.0
  [2 ,  2]  =  3.0
  [3 ,  2]  =  -1.0
  [4 ,  2]  =  -1.0
  [7 ,  2]  =  -1.0
  [2 ,  3]  =  -1.0
  [3 ,  3]  =  3.0
  [6 ,  3]  =  -1.0
  [8 ,  3]  =  -1.0
  [2 ,  4]  =  -1.0
  [4 ,  4]  =  3.0
  [9 ,  4]  =  -1.0
  [10,  4]  =  -1.0
  [1 ,  5]  =  -1.0
  [5 ,  5]  =  3.0
  [8 ,  5]  =  -1.0
  [10,  5]  =  -1.0
  [1 ,  6]  =  -1.0
  [3 ,  6]  =  -1.0
  [6 ,  6]  =  4.0
  [9 ,  6]  =  -1.0
  [10,  6]  =  -1.0
  [2 ,  7]  =  -1.0
  [7 ,  7]  =  1.0
  [3 ,  8]  =  -1.0
  [5 ,  8]  =  -1.0
  [8 ,  8]  =  3.0
  [9 ,  8]  =  -1.0
  [4 ,  9]  =  -1.0
  [6 ,  9]  =  -1.0
  [8 ,  9]  =  -1.0
  [9 ,  9]  =  4.0
  [10,  9]  =  -1.0
  [4 , 10]  =  -1.0
  [5 , 10]  =  -1.0
  [6 , 10]  =  -1.0
  [9 , 10]  =  -1.0
  [10, 10]  =  4.0

julia> partialeigen(partialschur(L, which=LR())[1])
([1.05398e-15, 0.574619, 5.89222, 5.72277, 4.85151, 3.81628, 3.04859, 1.54485, 2.08427, 2.46491], [-0.316228 -0.294971 … 0.140312 0.371051; -0.316228 0.339792 … -0.0655438 0.372631; … ; -0.316228 -0.135193 … 0.169911 -0.163748; -0.316228 -0.156412 … 0.330283 -0.162148])

[jpfairbanks]

Actually there is a slightly bigger problem:

`julia julia> partialeigen(partialschur(L, which=SR())[1]) ([-1.72814e-16, 5.89222, 5.72277, 0.574619, 1.54485, 3.81628, 3.04859, 2.08427, 2.46491, 4.85151], [0.316228 -0.188541 … 0.371051 0.0566662; 0.316228 0.0299987 … 0.372631 -0.601047; … ; 0.316228 -0.46855 … -0.163748 -0.348467; 0.316228 -0.269598 … -0.162148 -0.00893568])

julia> partialeigen(partialschur(L, which=LR())[1]) ([-4.42121e-16, 5.89222, 5.72277, 0.574619, 4.85151, 3.81628, 1.54485, 3.04859, 2.46491, 2.08427], [-0.316228 -0.188541 … -0.371051 -0.140312; -0.316228 0.0299987 … -0.372631 0.0655438; … ; -0.316228 -0.46855 … 0.163748 -0.169911; -0.316228 -0.269598 … 0.162148 -0.330283])

I'm getting the LR eigenvalues when I pass in the SR value for which.

[haampie]

Whenever the maxdim param is of the size of the matrix you get all eigenvalues. In this case: nev = 6 is the default parameter, and then automatically maxdim = max(2nev, order of matrix) = 10. So you get a full schur / eigen decomposition.

Currently I'm working on sorting the eigenvalues by the which target already during the partialschur call in #81, and also returning just nev eigenvalues even if more have converged.

[haampie]

Upon reading this issue again: do you mean the method should not converge to the 0 eigenvalue?

Current master orders the eigenvalues the right way! :)

[jarlebring]

I cannot reproduce this error. I don't see what output was not as expected in the original post. Whatever it was, was maybe solved in #81. A MWE would have been nice from the start. Here it is:

julia> using LinearAlgebra, SparseArrays;
julia> L_I=vec([1  5  6  2  3  4  7  2  3  6  8  2  4  9  10  1  5  8  10  1  3  6  9  10  2  7  3  5  8  9  4  6  8  9  10  4  5  6  9  10]);
julia> L_J=vec([1  1  1  2  2  2  2  3  3  3  3  4  4  4  4  5  5  5  5  6  6  6  6  6  7  7  8  8  8  8  9  9  9  9  9  10  10  10  10  10]);
julia> L_F=vec([2.0  -1.0  -1.0  3.0  -1.0  -1.0  -1.0  -1.0  3.0  -1.0  -1.0  -1.0  3.0  -1.0  -1.0  -1.0  3.0  -1.0  -1.0  -1.0  -1.0  4.0  -1.0  -1.0  -1.0  1.0  -1.0  -1.0  3.0  -1.0  -1.0  -1.0  -1.0  4.0  -1.0  -1.0  -1.0  -1.0  -1.0  4.0]);
julia> L=sparse(L_I,L_J,L_F);
julia> partialeigen(partialschur(L, which=LR())[1])
([5.89222, 5.72277, 4.85151, 3.81628, 3.04859, 2.46491], [-0.188541 0.14886 … 0.086437 -0.371051; 0.0299987 0.0885386 … 0.484552 -0.372631; … ; -0.46855 -0.449346 … -0.349383 0.163748; -0.269598 0.64257 … -0.14963 0.162148])

This seems to be the correct solution

julia> -sort(-eigvals(Matrix(L)))[1:6]
6-element Array{Float64,1}:
 5.892222827045327 
 5.722769111699081 
 4.851506514189556 
 3.816275449388274 
 3.048585741536626 
 2.4649075018504543

Same for SR:

julia> partialeigen(partialschur(L, which=SR())[1])
([3.11071e-17, 0.574619, 1.54485, 2.08427, 2.46491, 3.04859], [0.316228 -0.294971 … 0.371051 0.086437; 0.316228 0.339792 … 0.372631 0.484552; … ; 0.316228 -0.135193 … -0.163748 -0.349383; 0.316228 -0.156412 … -0.162148 -0.14963])
julia> sort(eigvals(Matrix(L)))[1:6]
6-element Array{Float64,1}:
 -1.2382770441425396e-15
  0.5746192521010607    
  1.5448476819141115    
  2.0842659202755205    
  2.4649075018504543    
  3.048585741536626     

Everything as expected as far as I can see. I would propose to close this issue unless more details are provided.