Is the QR algorithm stable?

Question

I always thought the QR algorithm was backwards stable. Then I tried the following:

Let $$A = \begin{pmatrix} 1& 0\\ 0& 2 \end{pmatrix},\ E= \begin{pmatrix} 0 & 1 \\ 1 & 0 \end{pmatrix}$$ and $\varepsilon = 10^{-15}$

For the matrix $ B = A +\varepsilon E$ I computed the eigenvalues via QR algorithm. I plotted $|(B_k)_{21}|$ where $k$ means the iteration number and I also plotted $|(B_k)_{11}-1|$ which is just the difference between the matrix iterate and the first eigenvalue $(\lambda_1 = 1)$ of $A$.

But the behaviour actually indicates that the QR algorithm is not stable. I only perturbed the Matrix in the off diagonal elements by $10^{-15}$ and yet it took $100$ iterates to get back to "normal".

How can this behaviour be explained?

Answer 1

Without perturbation, the QR decomposition of $A$ is trivial, $Q=I$ and $R=A$ (or depending on the method used, $[Q,R]=[-I,-A]$). Thus it is not surprising that the QR iteration will at first move slowly when the off-diagonal entries are only minimally perturbed. That it moves at all rests on the tendency of the QR algorithm without shifts to arrange the eigenvalues in decreasing order of magnitude.

The change of the off-diagonal elements is approximately by a factor of $\frac{λ_2}{λ_1}$. In the initial order, this gives a factor $2$ for each iteration, approximately doubling the element. As $10^{-15}\approx 2^{-50}$, it needs about 50 doublings to get the off-diagonal elements in the same magnitude as the diagonal elements. After that the order of the eigenvalues switches, so that the factor rapidly becomes $1/2$. It needs another 50 iterations to get the off-diagonal elements to numerically zero.

A equivalence transform that switches rows and columns with a rotation is the rotation by $90°$. The half-way point is the rotation by $45°$. The transformed matrix $A_k$, $k$ around $50$, at that point is $$ \frac12\pmatrix{1&1\\-1&1}\pmatrix{1&0\\0&2}\pmatrix{1&-1\\1&1} =\frac12\pmatrix{1&2\\-1&2}\pmatrix{1&-1\\1&1} =\frac12\pmatrix{3&1\\1&3} $$ The off-diagonal entry $\frac12$ corresponds to the maximum in the second graph.