MATLAB: Eigs Function on Function Handle Not Converged

eigenvalueMATLAB

my little fun problem looks like:

F is some neat linear function handle. I need to find its largest real eigenvalues. (it basically looks like A*x where the construction of A involves unpleasant large matrix inversion and is huge (size of tens of thousands squared)

So I called eigs and it did not converge when specified to calculate 10 largestreal eigenvalues but converged when specified to calculate 10 smallestreal. I need both so any reason this might be it?

Or what is the threshold that determines convergence or not behind the curtains?

Best Answer

Hi Sam,

The most likely reason for the convergence problems is that the eigenvalues are close together (or even multiples). The convergence speed of the internal method (for the 'largestabs' case) depends on the ratio between the smallest chosen eigenvalue and the largest eigenvalue that was not chosen.

Because of this, increasing the number of eigenvalues you are asking for (or just the 'SubspaceDimension') can help with convergence.

Another factor is that the 'largestreal' and 'smallestreal' options have some issues compared to 'largestabs' and 'smallestabs'. The inner iteration tends to go for the largest eigenvalues by absolute value, and has to be called back on every outer iteration to go for largest or smallest real part instead. So if the spectrum of your matrix is larger in the imaginary axes than the real axes, that could also be a reason for the convergence problems. Unfortunately, there's not much I can think of to improve that case.

Related Solutions

MATLAB: Specify a tolerance for backslash operation or linsolve while A is non-square? (QR solve)

What do you mean by "tolerance"?

A QR solution has no need for a tolerance, nor does linsolve, as they are not iterative methods.

Anyway, you cannot specify the desired accuracy of a solution in a linear system like this, solved using QR. Just wanting a more accurate solution is not relevant. The solution exists, and is unique, or there are infinitely many solutions, all of which are equally good (or poor), or there is no exact solution, and that which qr/linsolve returns is entirely valid.

Perhaps you have a singular system, or one that is nearly singular. In that case, a tolerance is still not really relevant. Yes, you could decide which singular values to discard, in case of a pinv style solution. But that is not a "tolerance" in the way tolerances are traditionally used. And if your system is singular, then again, tolerances have no real meaning.

As far as your matrix being non-square, that just means (assuming you have an over-determined system) that you will never be able to achieve an exact solution. (Well, not entirely true. It IS possible that an over-determined system has an exact solution. But that is terribly rare.) So if no exact solution exists, then specifying a tolerance is again a meaningless goal.

If your problem is non-square in the sense that it is under-determined, then again, a tolerance is irrelevant. A solution may exist, or not. If one solution does exist, then there will be infinitely many solutions, all equally good. A tolerance is meaningless here.

I think that perhaps you are being confused, thinking that the solution to a linear system is like that of a nonlinear system, where you will employ convergence tolerances, etc. Or, perhaps you think that simply by cranking down on a tolerance, you can always gain an arbitrarily good solution. (Many people seem to think that. They are flat out wrong.) Or perhaps you are thinking that a solution like that provided by linsolve/qr is like an iterative solver, like that from lsqr.

If this has not cleared things up, then I will suggest you need to explain CLEARLY what you want to do, and why you think that a tolerance is a meaningful thing to ask for here. As well, you need to explain why it is that you think what I have said above does not apply.

MATLAB: In sparse matrix eigs(A,k,sigma), how to specify both shift and direction

With a numeric sigma, EIGS uses a linear system solve (A-sigma*I)\v instead of A*v to compute the largest eigenvalues of the inverse of A. You can do this directly by calling

[U, D] = eigs(@(x) A\x, size(A, 1), k, 'largestabs');
D = diag(1./diag(D));

which is equivalent to the call

[U, D] = eigs(A, k, 'smallestabs');

The first call can be changed to return only the largest real instead of the largest absolute value, by simply replacing 'largestabs' with 'largestreal' (note: If your MATLAB version is older than R2017b, use 'lm', 'sm', instead of 'largestabs' and 'smallestabs', and either 'la' or 'lr' instead of 'largestreal'). After inverting, this translate to find the eigenvalues with smallest (positive) real difference to the sigma.

Also, the call to A\x directly is very expensive, so you should compute a factorization and use that for the linear system solve instead.

In summary:

dA = decomposition(A);
[U, D] = eigs(@(x) dA\x, size(A, 1), k, 'largestreal');
D = diag(1./diag(D));

Caveats:

If sigma is to the right of the spectrum (that is, all eigenvalues have real part less than sigma), this method will find the leftmost eigenvalue, not the one closest to sigma.
The "largest real" algorithm is not as robust as the "largest abs" algorithm: There's a larger chance that some eigenvalues are not found even though they have the largest real part (particularly if another eigenvalue has a slightly smaller real part but much larger imaginary part).

Best Answer

Related Solutions

MATLAB: Specify a tolerance for backslash operation or linsolve while A is non-square? (QR solve)

MATLAB: In sparse matrix eigs(A,k,sigma), how to specify both shift and direction

Related Question