I'm writing some C++ matrix library for hobby. For computing eigenvalues and eigenvectors, I referred the following "Francis double step QR algorithm":
In particular, page 82 of https://people.inf.ethz.ch/arbenz/ewp/Lnotes/chapter4.pdf
This is my code: https://github.com/frozenca/Ndim-Matrix/blob/main/LinalgOps.h#L1045
My implementation correctly computes eigenvalues, but fails to compute eigenvectors, except the first eigenvector corresponding to the largest eigenvalue.
For $\begin{pmatrix} -2 & 2 & 3 & 4 \\ -9 & 7 & 5 & 5 \\ -5 & 2 & 6 & 6 \\ -7 & 2 & 8 & 9 \end{pmatrix}$, Wolframalpha gives the following result:
$\lambda_{1} = 12.6009, \lambda_{2} = 4.2635, \lambda_{3} = 2.5983, \lambda_{4} = 0.5373$ with corresponding eigenvectors
$v_1 = (0.520545, 0.706921, 0.728858, 1)$
$v_2 = (-2.02909, -7.68532, -0.446183, 1)$
$v_3 = (8.65707, 12.3871, 3.67796, 1) $
$v_4 = (-0.748529, -0.608655, -1.56064, 1)$
My code outputs the following:
(12.6009,0) {(0.343092,0), (0.465934,0), (0.480393,0), (0.659103,0)}
Av(normalized) : {(0.343092,0), (0.465934,0), (0.480393,0), (0.659103,0)}
(4.26349,0) {(0.101326,0), (0.836911,0), (-0.197859,0), (-0.500164,0)}
Av(normalized) : {(-0.180481,0), (0.234039,0), (-0.48548,0), (-0.822777,0)}
(2.59833,0) {(0.924778,0), (-0.243128,0), (-0.270181,0), (-0.112592,0)}
Av(normalized) : {(-0.203286,0), (-0.674777,0), (-0.41863,0), (-0.572799,0)}
(0.537284,0) {(0.12964,0), (-0.152886,0), (0.810603,0), (-0.55022,0)}
Av(normalized) : {(-0.276698,0), (-0.77434,0), (0.503774,0), (0.264664,0)}
,0 part is the imaginary part.
I collected the column vectors of the $Q$ matrix upon convergence of the algorithm.
The eigenvalues are correct, and the first eigenvector $\mathbf{v}$ is also correct - being equal to the normalized $\mathbf{A}\mathbf{v}$. But other eigenvectors seem just broken. Really weird that only one eigenvector is correct.
Is the paper wrong? Should I use another algorithm for computing eigenvectors? Or just is my implementation broken?
