Is this a bug in eig?

This issue came up with the roots function (see roots_ on File Exchange), where eig was failing on the companion matrix. Following are two test cases; the first one works and the second one doesn't. (I reported this to tech support, but it's not listed in the official Matlab bug reports.)

a = [-1,-1,0;1,0,0;0,1,0];

a(1,3) = 10e-32;
[v,d] = eig(a);
disp(abs(a*v-v*d))
%    1.0e-15 *
%    0.055511151231258   0.055511151231258   0.000000000000000
%    0.166533453693773   0.166533453693773   0.000000000000000
%    0.138777878078145   0.138777878078145                   0
disp(num2str(diag(d)))
%  -0.5+0.86603i
%  -0.5-0.86603i
% 1e-31+0i 

a(1,3) = 9e-32;
[v,d] = eig(a);
disp(abs(a*v-v*d))
%    0.000000000000000   0.000000000000000   0.000000000000000
%    0.000000000000000   0.000000000000000   0.000000000000000
%    0.707106781186548   0.707106781186548   0.000000000000000
disp(num2str(diag(d)))
% -0.5+0.86603i
% -0.5-0.86603i
%    0+0i 

Answer 1

[Matt J]

It's fine. You just need to disable balancing:

a = [-1,-1,0;1,0,0;0,1,0];


for e=[10e-32, 9e-32]
    
 a(1,3) = e;
 [v,d] = eig(a,'nobalance');

 Discrepancy = abs(a*v-v*d)

end

Discrepancy = 3×3

1.0e+-15 * 0.3608 0.3608 0.0342 0.3554 0.3554 0.0220 0.2355 0.2355 0.2282

Discrepancy = 3×3

1.0e+-15 * 0.3608 0.3608 0.0342 0.3554 0.3554 0.0220 0.2355 0.2355 0.2282

[Kenneth Johnson]

The eig documentation says "[V,D] = eig(A) produces a diagonal matrix D of eigenvalues and a full matrix V whose columns are the corresponding eigenvectors so that A*V = V*D." Also, "eig(A,'balance') is the same as eig(A)." The behavior fails to conform to the documentation.

The "nobalance" option "sometimes gives more accurate results" but this is not a minor accuracy discrepancy; it is a complete failure. (Discrepancy of order 1.) Shouldn't the eig algorithm be able to detect that it has failed?

[Matt J]

Shouldn't the eig algorithm be able to detect that it has failed?

How? It has no way of knowing whether you consider discrepancies of order 1 large.

[Matt J]

Kenneth Johnson's comment moved here:

The algorithm can use any reasonable default tolerance; this is an obvious failure.

[Matt J]

If I had to guess, I would say there were 2 considerations in their decision not to issue a warning.

(1) The choice of the criterion. Error on the order of 1 wouldn't be considered as significant if the eigenvalues were on the order of 1e8. But then what error criterion would you use? Error relative to the maximum eigenvalue? To the median eigenvalue? Everyone will have a different prefered criterion.

(2) Computational cost. For large matrices, the evaluation of a*v-v*d is significant extra computation. Not all users wish to compromise on speed in the interest of playing it safe. Therefore, they leave it up to the user to do their own post-checking.

[Kenneth Johnson]

(1) The eig argument and the eigenvectors are of order 1 and the error is of order 1, obviously a failure by any error criterion. Also, the behavior is discontinuous across a(1,3) = 9e-32 or 10e-32, suggesting that there is some errant branch logic in the code. This test case is derived from roots([1,1,1,-10e-32]), which correctly give a root of 10e-32, whereas roots([1,1,1,-9e-32]) incorrectly gives a root of zero. roots(c) should never give a zero root when c(end)~=0. If 'nobalance' is used in roots, then the answer is wrong, so 'nobalance' does not actually improve accuracy of the eigenvalues in this case.

(2) "The answer is wrong, but it computes fast!" Efficiency is no excuse for a bug.

[Matt J]

The eig argument and the eigenvectors are of order 1

I assume you mean eigenvalues..

and the error is of order 1, obviously a failure by any error criterion.

In this case maybe, but the software has to choose a policy applicable to a range of cases.

[Kenneth Johnson]

I meant eigenvectors (v). The error, a*v-v*d, is proportional to v. (But in the roots problem it is the eigenvalue that is of concern.)

"... the software has to choose a policy applicable to a range of cases." No excuse for a complete failure in this case.

[Matt J]

But the documentation warns you that it may occur precisely in the situations that your matrix represents. Is your beef that the code didn't issue a warning, or that it didn't autoselect the correct settings for you?

[Catalytic]

I for one would expect a warning, similar to what mldivide() used to do if it deemed that it was doing a singular matrix inversion. Or maybe eig should have an exitflag output similar to what Optimization Toolbox functions have

[V,D, exitflag] = eig(___)

That way, the overhead of post-checks might be incurred only if 3 output arguments were requested.

[Kenneth Johnson]

That explains it. In this case I think it should issue a warning when the result is so obviously wrong.

It seems that the balance function isn't really suitable for use in roots. It forces the scaling factors to be powers of 2, but if I use a scaling matrix t = diag([1,1,1/sqrt(a(1,3))]), then it works fine for my particular test case. I suspect that roots accuracy could be significantly improved by using a better adapted balancing function.

[Bruno Luong]

Interesting. I would say the 'balance' is a damn dangeruous option that can produce such obviously wrong result, especially it affects the two eigen space that is correctly balanced at the first place.

Since when this option is introduced?

[Bruno Luong]

Google I found this page where there is a paper discuss issue with banancing. Someone claims (last post) he has developped a correct balancing method in Lapack 3.5.0.

[Kenneth Johnson]

For my test case it doesn't work without balancing, but Matlab's balancing doesn't work. There is a discussion of matrix balancing in Numerical Recipes in C++.