eigs does not return the eigenvalues closest to shift sigma

Hi,

I'm using eigs with a shift SIGMA to compute a few interior eigenpairs of a (large and sparse) generalized eigenvalue problem with symmetric positive-definite matrices. Just for the record, I use a function handle that computes (A-SIGMA*I)\X , with A given by transforming the generalized eigenvalue problem into a standard eigenvalue problem (through Cholesky factorization of the mass matrix and LDL factorization of the shifted stiffness matrix). The value for SIGMA is (a + b)/2 with ]a,b[ the interval for which I want to compute all the eigenpairs. The number n of eigenpairs is determined by inertia count at a and inertia count at b (using LDL factorizations). This eigenvalue problem exhibits very clustered eigenvalues.

MATLAB successfully computes eigenvalues of this problem (I have a reference solution to verify the exactness). The issue is that some of them (let's say p of them) are out of the interval ]a,b[ (they are all to the left). And since MATLAB produces only n eigenpairs (as requested), p eigenpairs at the other end of the interval (to the right) are missing.

The default convergence tolerance of eigs is 1e-14 (which I'm using) and the above issue, equivalent to a shift in the eigenvalues, leads to a relative eigenvalue error in the order of up to 1e-2 . The shift aside, all the computed eigenvalues are exact; the p eigenvalues outside ]a,b[ are also exact. As a conclusion, the claim "If SIGMA is a real or complex scalar including 0, eigs finds the eigenvalues closest to SIGMA." is not satisfied here.

The function handle that returns (A-SIGMA*I)\X is correct since all the eigenvalues are correct. There are 87 eigenvalues in the interval ]2.395285793575430 , 2.402047359143023[ . There are 72 eigenvalues in the interval ]a,b[ = ]2.402047359143023 2.666048612636298[ . We can see that there is a cluster of eigenvalues at the left of ]a,b[ .

Has this issue been encountered before?

Many thanks. Olivier

Answer 1

[Christine Tobler]

I have not heard of this particular issue before, but I could see it happening when the eigenvalues are very clustered.

If you are able to share the matrices, I could take a look at what is going on here. But I expect the best way around this would be to choose a somewhat larger interval than the one you finally require, to have a bit of a buffer zone for these kinds of cases.

Answer 2

[Olivier Ezvan]

Thanks for your answer.

Yes, I ended up checking the number of eigenvalues inside the desired band and introducing a new shift in case some are missing.

In order to share the matrices with you, first, I assembled them (prior to that, I was not assembling them and I was using a function handle in eigs). They take a total of about 5 GB in memory.

Then, using the assembled matrices, I re-did the eigenvalue computation with the shift in the middle of ]a,b[ and eigs did return all the eigenvalues in the interval, unlike when I used the function handle. I should add that I had also computed all the eigenvalues with the function handle thanks to multiple shifts, and that the eigenvalues obtained with both methods coincide (1e-14 to 1e-11 relative error).

I should add that the issue encountered with my function handle was reproduced any time tested.

Olivier

[Christine Tobler]

Hi Olivier,

Thank you, this is interesting to know. I talked with a colleague, who suggested that when he did similar things, he would take the number of expected eigenvalues k, based on the inertia counts, of both shifts, and multiply it with a smallish factor (e.g. 1.2), before calling the iterative eigenvalue solver. He would then choose the eigenvalues inside the interval from the results, and be reasonably sure to have caught them all.

Now that I think about it, when eigenvalues are very tightly clustered around the interval's boundaries, it's quite likely for the round-off errors in the inertia computation, in the matrix factorization, and in the eigs function to result in slightly inconsisten results as to which eigenvalue is in the interval.

Christine

[Olivier Ezvan]

Hi Christine,

yes, it can be a solution to compute more eigenvalues than required, but we still cannot be sure in advance that we will obtain all the eigenvalues in the band. That is why I count the number of computed eigenvalues inside the band and compare with the count obtained through the 2 LDL inertias. If the count does not match, we have to either re-launch eigs for more eigenvalues or introduce a new shift. For my particular application, it is more efficient to just introduce a new shift whenever necessary. And, concerning the number of eigenvalues to request, it is more efficient to just consider "smallish factor = 0" (to request just as many eigenvalues as there are in the band) as the issue does not occur very often.

But what I would like to point out is the following. I believe that the round-off errors of the LDL factorizations for the 2 inertia counts are very small compared to the distance between any 2 clustered eigenvalues, for my application. The eigenvalues provided by eigs have an error that is negligible in comparison to the distance between any 2 clustered eigenvalues, for my application. When the issue happens, the eigenvalues provided by eigs correspond to an interval that is not centered around the shift, but shifted/deviated by a distance that is far greater than the 1e-14 eigs tolerance (I repeat, it can be up to 1e-2 for my application).

For my application, the eigenvalues computed by eigs are always consistent with any LDL inertia count. What is not consistent is which ones are computed, given the shift.

So, it seems to me that the issue is not round-off errors of eigenvalues and inertia counts compared to eigenvalue separation but rather, that it is related to the fact that the Krylov subspace somehow shifts/deviates from the shift. Which would mean that eigs solver is not 100% reliable in regards to the claim "If SIGMA is a real or complex scalar including 0, eigs finds the eigenvalues closest to SIGMA." That, or my application is somehow 'ill-posed'. I would like to have feedback to determine which one it is.

Thanks, Olivier

Answer 3

[Christine Tobler]

Hi Olivier,

Yes, whether it's preferable for performance to compute more eigenvalues than needed, or to factorize the matrix for several shifts will of course depend on your problem.

There isn't a 100% guarantee in eigs that it will find the eigenvalues closest to SIGMA - we check residuals of the computed eigenvalues, but we have to trust the Krylov subspace method to converge towards the largest eigenvalue of the shifted-and-inverted matrix. It's possible to construct cases where this will go wrong, but I don't think this is likely to be your case - I think the problem is more likely to be in the construction of this shifted-and-inverted matrix.

Since I think you said that eigs finds the correct eigenvalues when called with a matrix, but they are wrongly shifted when called with your function handle, I'm wondering if the problem could be with your function handle? Maybe it's using a different factorization than you were using, and avoid the problem that way? Can you tell me how this is constructed?

One way to check my last speculation would be to call eigs with your function handle, but using the 'lm' option instead of a shift. Does this have the same problem (that is, the eigenvalues that are returned are not the largest eigenvalues available)?

Thanks, Christine

[Olivier Ezvan]

Hi Christine,

yes, my function handle does not use the same factorization as eigs, because my application has a block structure with diagonal block that I exploit by performing a block LDL factorization of the shifted stiffness matrix, based on Schur complement.

I just did the following test: I constructed the function Y = (A-SIGMA*I)\X based on the direct LDL factorization of the shifted stiffness matrix by calling function 'ldl' (what eigs probably does behind the scene when given matrices as input instead of a function handle); the outputs of this function and of my function are identical (~1e-12 relative error). This is consistent with the fact that my function handle led to the correct eigenvalues.

The calculation I'm doing is a broadband eigenvalue calculation, by using 116 slices (shifts). The issue I'm talking about (wrong eigenvalues returned and consequently, missing eigenvalues) happened just for a few slices. I had also done the calculation with twice the number of slices and it reduced the number of missing eigenvalues, as, for a slice where the issue occurs, the window of computed eigenvalues would shift less. Thus, I don't know if it would be useful to try and compute the 'lm' eigenvalues with my function handle. Anyway, my function handle cannot do that, as it requires a real scalar SIGMA as input.

What I tested with eigs and the assembled matrices as inputs was the eigenvalue calculation for just one slice for which the issue had occured with my function handle. Maybe I should do the 116 slices with the assembled matrices as inputs and see if the issue occurs for any of them. I could also do that with the classic function handle that I mentioned above (based on direct ldl of the shifted stiffness matrix).

Olivier

Answer 4

[Christine Tobler]

Hi Olivier,

I'm not sure I know much more to try out. Note that EIGS is doing Y = (A-SIGMA*I)\X based on the LU factorization, not the LDL factorization. Maybe this would give you results that match what EIGS does when called with the matrix directly? I'd be really astonished if there was anything systematically different in which eigenvalues are returned based on using LDL or LU, though.

I don't understand how your function handle can take a real scalar SIGMA: don't you have to factorize the shifted matrix once, before you pass the function handle to EIGS? Certainly EIGS will not pass in any value SIGMA to the function handle it receives.

Christine

Answer 5

[Olivier Ezvan]

Christine,

I'm going to do such tests and let you know.

I do some preliminary calculations before passing the function handle to EIGS, yes, but the way I solve the linear system requires the knowledge of the shift SIGMA used. So I myself give SIGMA as an input to my function. In my function, SIGMA should be a finite real scalar. I don't know what SIGMA should be to correspond to 'LM', it seems to me that 'LM' would require a function that computes Y = A*X anyway?

Best,

Olivier

[Christine Tobler]

Oh, I see. In that case, you would give the same SIGMA as you normally give to EIGS, but tell EIGS that the SIGMA input is 'lm'. It will then treat your (A - SIGMA*I)^-1 function handle as if it's applying A*X, and return the eigenvalues of (A - SIGMA*I)^-1. You could then transform the eigenvalues returned back by doing 1./d + SIGMA, which is what EIGS does internally.

[Olivier Ezvan]

Ok, I haven't done that yet. What I just did is the eigenvalue computation of the faulty slice (one of the faulty slices, actually) with a function handle that uses the LDL factorization of the assembled shifted stiffness matrix instead of my function handle. This returned the same eigenvalues as with my function handle (relative error around 1e-15, and for some exactly 0 which surprises me quite a bit), I mean the shift issue happened exactly the same way.

You advised me to try LU instead of LDL. Which LU should I use? [L,U] = lu(A) or [L,U,P] = lu(A) or [L,U,P,Q] = lu(A) or [L,U,P,Q,R] = lu(A) ?

Thanks,

Olivier

Update1: I just did the computation with [V,D] = eigs(K,chol(M),n,SIGMA,'IsCholesky',1); and the spurious shift occured in the exact same way, same result, same eigenvalues.

Update2: I just did the computation with [V,D] = eigs(K,M,n,SIGMA); and the spurious shift occured in the exact same way, same result, same eigenvalues. This contradicts what I told you a few days ago and I found out why: a few days ago I used a slightly different shift. And the reason of the different shift is the following mistake:

eigenvalue = (2*pi*eigenfrequency)^2

The interval bounds can be defined either through the eigenvalues or the eigenfrequencies. The shift is the mean between the eigenvalues. I took the mean between the eigenfrequencies, then squared. Which is not the same.

Conclusion: a priori, eigs does return the expected eigenvalues, the problem was a mistake in calculating the shift SIGMA.

[Olivier Ezvan]

Ok, problem solved. I made a mistake in calculating the shift (see edited comment above). Now the calculation works perfectly.

Thanks for your time, Christine.

Answer 6

[Andrew Knyazev]

You may also want to try https://www.mathworks.com/matlabcentral/fileexchange/48-lobpcg-m passing your function handle to it. Compared to EIGS, which implements single-vector Lanczos and may thus miss some clustered eigenvalues, LOBPCG is a block method, i.e., very robust for clusters.