yixuan / RSpectra

R Interface to the Spectra Library for Large Scale Eigenvalue and SVD Problems
http://cran.r-project.org/package=RSpectra
80 stars 12 forks source link

eigs_sym has negative eigenvalue #7

Closed gdkrmr closed 6 years ago

gdkrmr commented 6 years ago

I have a relatively large symmetric eigenvalue problem (ca. 12000x12000) and want the first couple of eigenvectors and values:

e <- eigs_sym(x, 6, opts = list(retvec = TRUE))
e$values

### [1]   68213.67   31349.46   24366.68   19671.44   17046.18 -102290.45

e <- eigs_sym(x, 5, opts = list(retvec = TRUE))
e$values

### [1]   68213.67   31349.46   24366.68   19671.44 -102290.45

Am I right, that the entire procedure is stable, but the last eigenvalue is garbage? e$nconv is 5 and 6 respectively.

Is the correct solution to compute one more eigenvector/value than necessary and throw the last one away?

privefl commented 6 years ago

Are you really sure your matrix is symmetric? I get huge negative eigenvalues when testing on non-symmetric matrices.

Negative eigenvalues should be dropped but they are usually very small (absolute value < 1e-8).

gdkrmr commented 6 years ago

shouldn't eigs_sym ignore one of the triangles and therefore it should be impossible to produce negative eigenvalues?

privefl commented 6 years ago

This is a very good remark. From the documentation, I see that the lower triangle is used by default.

Maybe then, it should be because there is no structure in the data?

See what I tested:

x <- matrix(0, 1000, 1000); x[] <- rnorm(length(x))
RSpectra::eigs_sym(x, k = 6)$values

x2 <- matrix(0, 1000, 4); x2[] <- rnorm(length(x2)); x2 <- tcrossprod(x2)
RSpectra::eigs_sym(x2, k = 6)$values
yixuan commented 6 years ago

It may be a convergence issue, especially with large matrices.

gdkrmr commented 6 years ago

indeed isSymmetric(x) was false, but checking, the asymmetry was tiny and I symmetrized it x <- (x + t(x)) / 2 and it did not help.

I was already suspicious that it converged so fast (I am asking for the first 21 EVs):

Browse[2]> e$nconv
[1] 21
Browse[2]> e$niter
[1] 4
Browse[2]> e$nops
[1] 75
gdkrmr commented 6 years ago

also setting tol to a smaller values does not change much:

Browse[2]>   e <- RSpectra::eigs_sym(dgx2, out_dim + 1, opts = list(tol = 1e-18, retvec = TRUE))
Browse[2]> e$nops
[1] 96
Browse[2]> e$niter
[1] 6
Browse[2]> e$values
 [1]   68213.667   31349.457   24366.681   19671.441   17046.183   13932.041
 [7]   11735.303   10842.215    9632.626    8062.106    7776.732    6621.357
[13]    6001.662    5667.498    5129.942    4880.293    4683.910   -5224.663
[19]   -6539.950   -8087.123 -102290.454

Browse[2]> isSymmetric(dgx2)
[1] TRUE
gdkrmr commented 6 years ago

Same session, asking for less eigenvalues, the negative eigenvalues are identical:

Browse[2]> e <- RSpectra::eigs_sym(dgx2, 15, opts = list(tol = 1e-18, retvec = TRUE))
Browse[2]> e$values
 [1]   68213.667   31349.457   24366.681   19671.441   17046.183   13932.041
 [7]   11735.303   10842.215    9632.626    8062.106    7776.732    6621.357
[13]   -6539.950   -8087.123 -102290.454
gdkrmr commented 6 years ago

and now asking for more:

Browse[2]> e <- RSpectra::eigs_sym(dgx2, 35, opts = list(tol = 1e-18, retvec = TRUE))
Browse[2]> e$values
 [1]   68213.667   31349.457   24366.681   19671.441   17046.183   13932.041
 [7]   11735.303   10842.215    9632.626    8062.106    7776.732    6621.357
[13]    6001.662    5667.498    5129.942    4880.293    4683.910    4392.650
[19]    4130.336    3974.102    3856.767    3637.668    3523.318    3335.982
[25]    3238.525    3085.208    3042.923   -3430.344   -3621.775   -4225.863
[31]   -4534.989   -5224.663   -6539.950   -8087.123 -102290.454

This strongly suggests that there is a problem on the last couple of eigenvalues!

yixuan commented 6 years ago

I have just updated Spectra and the accuracy of symmetric solvers have been improved. I'm going to update RSpectra accordingly, so do you have a publicly available matrix file to test on?

yixuan commented 6 years ago

OK, a second thought: from your output I am pretty sure that your matrix is not positive definite, so it contains negative eigenvalues. The default selection rule of eigs_sym is "LM", which means the largest absolute value. If you only need the largest positive eigenvalues, you need to choose which = "LA".

gdkrmr commented 6 years ago

Now it makes sense, it computed the eigenvalues with the largest absolute values, i.e. the most extreme positive and negative ones and the result was correct. Also: sorry for the confusion, somehow I had the wrong idea that symmetric matrices only had positive eigenvalues. Thanks for clearing that up!

gdkrmr commented 6 years ago

for me the issue is solved, I can give the data to @yixuan if it helps development.