I've implemented a version of the double-shift QR algorithm featured in this report from ETH Zurich (Begins on page 77). The algorithm takes advantage of the Implicit Q theorem by applying an orthogonal-similarity transformation to our original matrix (in Hessenberg) form. This introduces a 'bulge' which we then chase out of the matrix by applying successive householder reflectors. The end goal is to covert the original matrix to real schur form (Although not truly, since complex eigenvalues form as blocks along the diagonal). The pseudocode for the algorithm is as follows:
(I've implemented this particular algorithm in C)
However, when trying to solve companion matrices to retrieve the eigenvalues, I find it often gets things wrong. Even their example results with a $6 \times 6$ couldn't solve for all the eigenvalues:
Compare this to my results, which have much the same mistake:
5.000000 6.000000 0.390635 -1.005053 -7.905002 7.702966
-6.000000 5.000000 -10.791799 -2.710116 -14.579765 -19.663045
0.000000 0.000000 3.491643 -2.656106 -4.440683 9.373136
0.000000 0.000000 3.843327 -1.491643 -10.829108 -5.067750
0.000000 0.000000 -0.000000 0.000000 4.000000 4.922432
-0.000000 -0.000000 0.000000 0.000000 0.000000 3.000000
Even smaller matrices pose problems (though they are close) with very low tolerances:
Input: (Eigenvalues are: 4, -3, 2, -1)
2.000000 13.000000 -14.000000 -24.000000
1.000000 0.000000 0.000000 0.000000
0.000000 1.000000 0.000000 0.000000
0.000000 0.000000 1.000000 0.000000
Output:
2.968869 11.491388 21.498124 17.801229
0.535591 -1.968869 -3.241184 -2.692433
-0.000000 0.000000 2.000000 1.078408
0.000000 -0.000000 -0.000000 -1.000000
Can someone suggest what sort of improvements can be made to obtain better results? I've tried implementing another algorithm I found, but it provided worse results. Even when I ran it as-is through Matlab instead of my own C-version of the software.