16 views (last 30 days)

Show older comments

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

Christine Tobler
on 9 Aug 2018

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.

Christine Tobler
on 13 Aug 2018

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

Christine Tobler
on 14 Aug 2018

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

Andrew Knyazev
on 21 Sep 2018

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!