bugs in eig?

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

[John D'Errico]

Please don't start adding one answer after another just to make a comment.

Answer 1

[Christine Tobler]

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.

Answer 2

[Walter Roberson]

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.

[John D'Errico]

Moved an answer to a comment.

"Dear Walter, Thank you very much for your kindly reply. In my problem, both A and B are singular (det(A)=det(B)=0), but I think we can still calculate the eigenvalues and eigenvectors. Do you mean that the problem is rcond(X) being too small?

Thanks, Tang Laoya"

[Walter Roberson]

According to https://software.intel.com/en-us/forums/intel-integrated-performance-primitives/topic/293858

"[...] QR algorithm for solving eigenvalue problem. the algorithm are numerical unstable for badly conditioned matrix (det is close to 0). "

[tanglaoya]

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]

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.