eigs function: incorrect eigenvalues

I am trying to use the function eigs to obtain several eigenvalues of the smallest magnitude of a singular sparse matrix. This funciton correctly obtaines the zero egienvalue and the corresponding eigenvector, but all non-zero eigenvalues are different from those obtained using the function eig. Moreover, calling eigs several times in a row produces different non-zero eigenvalues. I tried decreasing the tolerance and increasing the number of Lanczos vectors, but that didn't help.

I am attaching a mat file with the matrix. The code that I used for testing is below:

load('test.mat'); % this loads matrix A
[V,D] = eigs(A,5,0);
[V1,D1] = eig(full(A));

[Walter Roberson]

I seem to recall that eigs involves a random initialization, and that can matter when the eigenvalues are close together.

[Dmitry Kopelevich]

Yes, it does use a random inititialization, according to the documentation. Moreover, according to eig, some of the small-magnitude eigenvalues are very close together. Eigs does work well when I try to obtain the largest eigenvalues (which are further apart from each other).

Is there any way to tweak the eigs parameters so that it can handle eigenvalues that are close together?

[Walter Roberson]

Not a solution but informative https://www.mathworks.com/matlabcentral/answers/356200-eigs-runs-faster-for-more-eigenvalues-of-the-same-matrix

[Dmitry Kopelevich]

Thanks for pointing this out. Interestingly, eigs is able to correctly compute small eigenvalues if I start the search from the large-magntide eigenvalues, i.e. if I use something like

[V,D] = eigs(A,N-2);

Here, N is the size of matrix A. The large-magntiude eigenvalues of the test matrix are spaced further apart, which seems to help the iterations to converge even for the small-magnitude eigenvalues. Of course, at this point one might as well use eig instead of eigs.

[Walter Roberson]

Was Andrew's FEX contribution useful for you? https://www.mathworks.com/matlabcentral/fileexchange/48-lobpcg-m

[Dmitry Kopelevich]

Thanks for pointing it out. Unfortunately, this method is applicable only to Hermitian matrices. The matrices I need to analyze are not Hermitian.

[Dmitry Kopelevich]

I have just tried applying eigs to a large matrix that I need to analyze (not the test matrix I've uploaded earlier) and it does seem to work reliably. The eigenvalues are reproducible and the eigenvalue/eigenector pairs satisfy the eigenvalue equaiton. It is strange because my test case was essentially a scaled-down version of the matrix that I needed to analyze. But I am glad that it worked!

Answer 1

[Christine Tobler]

The problem is that the matrix A is badly conditioned:

>> cond(full(A))
ans =
     2.206837183464466e+16

This is just around the border of when EIGS (or even MLDIVIDE) would give a warning about the matrix being ill-conditioned. It doesn't happen to give that warning, but the results still come out wrong (I'll investigate exactly why the warning isn't given here).

The reason for the wrong results here is that for a sigma shift, EIGS operates on the inverse of (A - sigma*I): It computes the largest eigenvalues in absolute value of A\x instead of using A*x, and then inverts those eigenvalues to return the eigenvalues of A. But when A is ill-conditioned, A\x doesn't apply the inverse of a in exact arithmetic, because of round-off error.

Unfortunately, there isn't anything that will make this always come out right. For your particular matrix, I've found that moving the shift slightly to the right will make the matrix (A-sigma*I) be better conditioned, and give the correct eigenvalues:

     eigs(A, 5, 1e-4)

Another thing that works for this matrix is to request the largest eigenvalues by real value instead of the ones closest to zero (which happen to be the same ones for this matrix):

     eigs(A, 5, 'largestreal')

[Dmitry Kopelevich]

Thank you for the suggestion. For this particular type of matrices, 'largestreal' works well, since their eigenvales have to have non-positive real parts. (These matrices arise in a finite-difference scheme for a convection-diffusion operator which cannot have eigenvalues with positive real parts.)