bugs in eig?

1 次查看（过去 30 天）
tanglaoya
tanglaoya 2018-5-3
评论： Christine Tobler 2018-5-3
Dear all,
I am using matlab's eig function to solve eigenvalue problem. I found that for two group of matrices that are nearly identify (max(abs(A1-A2))<1.0e-16), but the eigenvectors have great difference. Could anyone help me to take a look at it?
You can reproduce the problem by the following steps:
1) download the attached files;
2) run the following matlab commands:
load A1;load B1;[c1,d1]=eig(A1,B1);load A2;load B2;[c2,d2]=eig(A2,B2);max(abs(A1(:)-A2(:))),max(abs(B1(:)-B2(:))),max(abs(c1(:)-c2(:)))
Thanks,
Thang Laoya
  1 个评论
John D'Errico
John D'Errico 2018-5-3
Please don't start adding one answer after another just to make a comment.

请先登录，再进行评论。

请先登录，再回答此问题。

回答（2 个）

Christine Tobler
Christine Tobler 2018-5-3
Both answers are correct. The eigenvectors of a matrix are not uniquely determined: Every eigenvector can be multiplied by an arbitrary number, and is still a valid eigenvector of this matrix. But when two matrices are slightly different, this can completely change the scaling that is chosen for each eigenvector.
You can see that both answers are completely correct by verifying that they satisfy the defintion: A*U = B*U*D
>> norm(A1*c1 - B1*c1*d1)
ans =
9.3235e-11
>> norm(A2*c2 - B2*c2*d2)
ans =
9.5119e-11
This shows that both results are correct (up to a reasonable precision).
The reason that the two matrices c1 and c2 are very different is:
  • The scaling of two matching eigenvectors is not the same. Some eigenvectors must be multiplied with -1, or with a complex number of absolute value 1, so that they can match their counterpart.
  • Order of eigenvalues: EIG does not guarantee any order of the eigenvalues. In this case, they seem to be in the same order in the beginning, but towards the end, many eigenvalues are very close together: a small change can mean that they are sorted in a different order.
  • Multiple eigenvalues: If several eigenvalues are very close together, they are treated as one eigenvalue with multiple eigenvectors. Any linear combination of these eigenvectors is a valid eigenvector, and any valid basis of the space spanned by these eigenvectors may be returned.
It's best to write algorithms that do not depend on the scaling or order of the eigenvalues that EIG returns, as these are completely arbitrary and can change for different releases or on different machines.
Walter Roberson
Walter Roberson 2018-5-3
The matrices are not well conditioned.
>> rcond(A1)
ans =
4.67833932391801e-08
>> rcond(A2)
ans =
4.67833932391802e-08
>> rcond(B1)
ans =
2.4799219889766e-07
>> rcond(B2)
ans =
2.4799219889766e-07
>> help rcond
rcond LAPACK reciprocal condition estimator.
rcond(X) is an estimate for the reciprocal of the
condition of X in the 1-norm obtained by the LAPACK
condition estimator. If X is well conditioned, rcond(X)
is near 1.0. If X is badly conditioned, rcond(X) is
near EPS.
  4 个评论
显示 2更早的评论
tanglaoya
tanglaoya 2018-5-3
Dear Walter, Thank you very much for your kindly reply. Do you have any suggestion on how to get accurate eigenvalues and eigenvectors of this kind of matrices?
Thanks, Tang Laoya
Christine Tobler
Christine Tobler 2018-5-3
The matrices here are ill-conditioned, but this does not significantly affect the results of eig.
I agree with Michael Chuvelev's reply on the Intel forum thread linked by Walter: The algorithm used in eig is backwards-stable, meaning that a small disturbance in the inputs has a small effect on the error norm(A*U - B*U*D) (proportional to the condition of the eigenvalue problem). Of course, this small disturbance can still have a large effect on the eigenvectors in U, or on the order of the eigenvalues in D.
On the side: The determinant should not be used to determine the condition of a matrix, as it is not a good indication of the difficulty of a problem ( doc det ). Use cond or rcond instead.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Linear Algebra 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by