When trying to compute the roots of a polynomial, I noticed that a common practice is to calculate the eigenvalues of the polynomial's companion matrix. For example, the "roots" function in matlab and numpy both adopt this method.
While these implementations are usually reliable, I found that the accuracy significantly deteriorates when the polynomials have multiple roots, especially with high multiplicity. On my desktop with numpy 1.21.2 (linked with openblas in the archlinux official repository):
>>> numpy.roots([1,4,6,4,1])
array([-1.00021255+0.j , -0.99999998+0.00021254j,
-0.99999998-0.00021254j, -0.99978748+0.j ])
>>> p = numpy.poly([2,3,4,5,5,5])
>>> numpy.roots(p)
array([4.99991037+0.00015525j, 4.99991037-0.00015525j,
5.00017925+0.j , 4. +0.j ,
3. +0.j , 2. +0.j ])
In the first example all the roots should have been -1. In the second example, the first three roots should have been 5. The errors in both examples are on the order of 1e-4, which is much larger than I expect for small matrices. My matlab (which uses intel MKL 2019.0.3) yields similar results. I tried generating companion matrices myself and call eig() to find their eigenvalues, and the above issue can be reproduced.
Interestingly I found that eig() works pretty well for general dense matrices with degenerate eigenvalues; companion matrices are the only type I found to have such problem. I did a few tests on several different machines and they all have the same issue.
I was wondering:
Is this an intrinsic problem of the eigenvalue algorithm, or is it just a machine/platform-dependent bug?
If it's intrinsic, then what would be a quick & robust way to compute polynomial roots? The companion matrix method is good enough in most scenarios but the above issue can be really annoying in some situations.