Fixed

eigs() bug

5 years ago • updated by 5 years ago • 1

I noticed a pretty severe bug for the eigs() function. In particular, if I do:

A = sparse(mp(eye(4),34));
[evecs,evals] = eigs(A,1);

I expect an identity matrix to be returned for both evecs and evals, however, I get the following result after one call of the above code:

evecs = 

         0.6575264478049146996368762588811955    
         0.1350828341320901021606099191810989    
        0.09164029388649532440547542667314778    
         0.7355363042680420857014721043910981    


evals = 

    1

If I call it again I get an error:

Error using mp/subsref (line 1375)
Index exceeds matrix dimensions.

Error in mpeigs (line 523)
        D = ordeig(H(r,r));

Error in mp/eigs (line 3763)
            [varargout{1:nargout}] = mpeigs(varargin{:});

and in general it seems very sporadic, random, and most importantly, incorrect each time I call it. However, if I instead use the same function but without specifying the "1" as a parameter, e.g.:

A = sparse(mp(eye(4),34));
[evecs,evals] = eigs(A);

Everything seems to work as expected and I get identities returned for both evecs and evals. Obviously I don't need the eigenvalues and eigenvectors of the identity matrix but this is giving me little confidence in the eigenvectors and eigenvalues returned for the actual matrices I'm interested in. Am I missing something here or calling the function incorrectly? Also, I've noticed similar issues with the generalized eigenproblem equivalent of the above example.

Thank you.

Vote

Replies 1
Oldest first
- Newest first
- Oldest first

Fixed

5 years ago

The eigs produced wrong result / shows error because it didn't handle such small matrices (up to now).

Krylov iteration methods are robust when search subspace is large ( > max(2*evals,20)). That is why eigs must switch to full QR to compute spectrum of matrices smaller than some threshold. This safeguard was not included in the toolbox up to this moment.

I have just added this workaround to the mpeigs.m. I will send you new mpeigs.m by email in a few minutes. Just copy it to the toolbox directory (overwrite the old file). Next toolbox versions will include the new mpeigs.m in the bundle.

P.S.

This kind of special cases (small matrices, small K, special structure) - are not indicative of eigs performance on large/real-life matrices. Because all of these are corner cases, which are handled accordingly by separate codes.

The only way to be confident is to check math. properties of the results, e.g.:

mp.Digits(34);
A = sparse(eye(4,'mp'));
[V,D] = eigs(A);
norm(A*V - V*D,1)/norm(A,1)

Customer support service by UserEcho