eigs bug for 0 as lowest eigenvalue in parts of the matlab versions

7 vues (au cours des 30 derniers jours)
pages larry
pages larry le 18 Mai 2019
Modifié(e) : Matt J le 20 Mai 2019
I'm eigs a sparse diagonal matrix that has 0 as the lowest eigenvalue but in some new versions of matlab it gets 1 as the wrong answer.
This is a code that get wrong results in matlab 2017b and 2018b2017b and 2018b, while gets the right answer in matlab 2016a and the online matlab. I'm comparing eig and eigs for identical matrix. The right answer will be always 0 0, and the wrong answer will always be 1 and 0. Sparse eigs has missed 0 as my lowest eigenvalue.
So is there any way to gurantee it right?
%%% my code
clearvars
clc
matrix = diag(0:20);
% matrix(1,1)=1;
% matrix(2,2)=0;
sp_matrix = sparse(matrix);
Eg_sp = eigs(sp_matrix,1,'sa');
Es = eig(matrix);
Es = sort(Es);
Eg_full = Es(1);
fprintf('Eg by sparse matrix is %f\n',Eg_sp)
fprintf('Eg by full matrix is %f\n',Eg_full)
  1 commentaire
Walter Roberson
Walter Roberson le 18 Mai 2019
R2019a for Mac gives -9.04753446879501e-17 for the Eg_sp, but R2018b for Mac gives 1 - 1.4432899320127e-15

Connectez-vous pour commenter.

Réponse acceptée

Matt J
Matt J le 20 Mai 2019
Modifié(e) : Matt J le 20 Mai 2019
From Tech Support:
The issue that you are seeing is because of a bug which was introduced in R2017b, while implementing the refactoring algorithm for 'eigs'.
Instead of building the Krylov subspace as {x, A*x, A^2*x, ...}, it was being built as {A*x, A^2*X, A^3*x, ...}. This means that the search space never contains any components in the direction of exact zero eigenvalues. This will in practice only be noticeable for cases where a matrix has an all-zero row, since otherwise numerical error will always keep a small component in the direction of the zero eigenvector.
The bug has subsequently been fixed in R2019a, and hence updating to R2019a should resolve the issue.
Otherwise, a workaround for the same is as shown below:
sigma = 1; % Any value which is not an exact eigenvalue of A
[U, D] = eigs(sp_matrix - sigma*speye(size(sp_matrix)), 1, 'smallestreal');
D = D + sigma*speye(size(D)); % Revert the shift
I hope this helps in resolving the issue. Meanwhile, I am closing this case for now. Please write back to me if you have any follow up questions/ comments regarding the same. I shall be happy to reopen the case and assist you further.
Please preserve the Reference ID in further correspondence on this query. This allows our systems to automatically associate your reply to the appropriate Case.
If you have a new technical support question, please submit a new request here:
Sincerely,
Abhi Baruah
Mathworks Technical Support

Plus de réponses (2)

Matt J
Matt J le 18 Mai 2019
That does seem like an interesting bug. This seems to work though
Eg_sp = eigs(sp_matrix,1,1e-12);
  1 commentaire
pages larry
pages larry le 18 Mai 2019
Thank you for your replying. But this bug was found in my exact diagonalization program, where the ground state energy in the program is not neccessarily near 0. So I cannot do that.

Connectez-vous pour commenter.


David Goodmanson
David Goodmanson le 18 Mai 2019
Modifié(e) : David Goodmanson le 18 Mai 2019
Hi pl,
2018b does not list 'sa' as an option for eigs, so I am speculating that that 'sa' is legacy. However, compariing 'sa' to 'smallestabs',
matrix = diag(0:20);
sp_matrix = sparse(matrix);
Eg_sp = eigs(sp_matrix,1,'sa')
Eg_sp = 1.0000 % incorrect
Eg_sp = eigs(sp_matrix,1,'smallestabs') % this errors out,
% error message suggests that the smallest eigenvalue is 0
% checking for small nonzero smallest eigenvalue
sp_matrix(1,1)=1e-10;
Eg_sp = eigs(sp_matrix,1,'smallestabs')
Eg_sp = 1.0000e-10
sp_matrix(1,1)=1e-100;
Eg_sp = eigs(sp_matrix,1,'smallestabs')
Eg_sp = 1.0000e-100 % but with a warning that the matrix is almost singular
so it appears that you could trap for zero eigenvalue, and then the other cases are covered.
  2 commentaires
pages larry
pages larry le 18 Mai 2019
In matlab 2018b both 'sa' and 'smallestreal' will encounter this bug and give the same answer.
David Goodmanson
David Goodmanson le 18 Mai 2019
well, that sucks.

Connectez-vous pour commenter.

Catégories

En savoir plus sur Logical dans Help Center et File Exchange

Produits


Version

R2017b

Community Treasure Hunt

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

Start Hunting!

Translated by