Shifted QR algorithm—why does the shift help?

Question

I read that a way to speed up the convergence rate of the QR algorithm is to shift the target matrix. It is not so clear to me why this helps. The convergence rate depends on the minimum gap between consecutive eigenvalues, but if we shift $A$ to $A + cI$, the gaps remain the same?

Answer 1

Convergence depends on the ratio between the eigenvalues, not on the difference. Oversimplifying: if $\lambda_1$, $\lambda_2$ are two eigenvalues and you shift by $\mu$, then the magic ratio is $\frac{\lambda_1-\mu}{\lambda_2-\mu}$. If $\mu$ is close to $\lambda_1$, convergence is fast.

Answer 2

(This is an extension of my comment on Federico Poloni's answer.)

For simplicity, I'll assume $A$ is a symmetric $n\times n$ matrix. Suppose $$A = \begin{bmatrix}A' & b \\ b^T & c\end{bmatrix},$$ where $c$ is a scalar, $b\in \mathbb{R}^{(n-1)\times 1}$, and $A'$ is whatever's left. Moreover, suppose $b$ is small, so the standard basis vector $e_n$ is nearly an eigenvector of $A$ with $c$ as the near-eigenvalue. Starting to do one step of QR iteration, we find an orthogonal $Q$ and upper triangular $R$ such that $A-cI=QR$. Since $A$ is symmetric, $A-cI=R^TQ^T$ as well. Rearrange this to get $Q = (A-cI)^{-1}R^T$. The important thing about this matrix equation is the last column: $R^T$ is lower-triangular, so its last column is a scalar multiple of $e_n$, so $q_n$, the last column of $Q$, is a scalar multiple of $(A-cI)^{-1}e_n$.

At this point it helps to make an interlude to look at things in a coordinate-independent way: note that $c=(e_n^TAe_n)/(e_n^Te_n)$. In general, for any approximate eigenvector $x$, $r=(x^TAx)/(x^Tx)$ (called the Rayleigh quotient) is the best approximation to the eigenvalue. Moreover, $x'=\alpha(A-rI)^{-1}x$ can be expected to be a better approximation to the eigenvector, where $\alpha$ is any convenient scalar. To see this, let $v_1,\ldots,v_n$ be orthonormal eigenvectors of $A$ with eigenvalues $\lambda_1,\ldots,\lambda_n$. $x$ is close to an eigenvector, so WLOG $x = v_1 + \sum_{i>1} a_iv_i$, where $\sum_{i>1} a_i^2 < \epsilon^2$ is small. In other words, $x = v_1 + O(\epsilon)$, where the $O(\epsilon)$ is orthogonal to $v_1$. The Rayleigh quotient is then $$r=\frac{x^TAx}{x^Tx} = \frac{\lambda_1 + O(\epsilon^2)}{1+O(\epsilon^2)} = \lambda_1 + O(\epsilon^2).$$ Then $$(A-rI)^{-1}x = \frac{v_1}{O(\epsilon^2)} + \sum_{i>1}\frac{a_i}{\lambda_i - \lambda_1 + O(\epsilon^2)}v_i = v_1/O(\epsilon^2) + O(\epsilon),$$ assuming $\lambda_1$ is a simple eigenvalue so the denominator in the sum stays bounded away from zero (and I don't see how to remove that assumption in this argument). We can renormalize this however we want, so making the coefficient of $v_1$ be 1, we get $$x' = \alpha(A-rI)^{-1}x = v_1 + O(\epsilon^3).$$ That justifies my claim that $\alpha(A-rI)^{-1}x$ is a better approximation to the eigenvector $v_1$ (in fact, cubically better!).

Returning to QR iteration, we have $q_n = (A-cI)^{-1}e_n$, where $q_n$ is the last column of $Q$ and $c$ is the Rayleigh quotient of $e_n$, so since $e_n$ is already nearly an eigenvector, $q_n$ is even closer. Finally, when we replace $A$ by $\hat{A} = RQ + cI$, all we're doing is changing our basis: $\hat{A} = Q^TAQ$. This makes the old $q_n$ the new $e_n$, so we're in the same situation we were in originally, except our approximation to the eigenvector is much better now.

Answer 3

This completes Polini's answer, which is perfectly right, but a bit 'elliptic'.

The QR algorithm is usually employed together with a preparation step: one puts the matrix $A$ in a unitarilly similar Hessenberg form $B$. This means that $b_{ij}=0$ unless $i\le j+1$. This preliminary step is cheap; its cost is $O(n^3)$ and does not exceed the cost of one step of QR. But it has a huge reward: the Hessenberg form is invariant under QR, and the complexity of each QR step is now reduced to $O(n^2)$.

Therefore we may assume wlog that $A$ is Hessenberg and also irreducible. If in addition its eigenvalues are of pairwise distinct moduli, the QR algorithm converges. See Theorem 13.3 of my book Matrices. Theory and applications. Springer GTM 216, 2nd edition 2010. Because the iterates are Hessenberg, thus close to triangular, it is not too much difficult to analyse how the convergence occurs. Actually, the tail (bottom right) converges faster, in the sense that $a_{n,n-1}\rightarrow0$ and $a_{nn}\rightarrow\lambda_n$ very fast, where $\lambda_n$ is the smallest eigenvalue. Thus $a_{nn}$ gives you an approximation of $\lambda_n$, which you can use to perform a shift. Then the ratio $\frac{\lambda_j-\mu}{\lambda_{j-1}-\mu}$ increases significantly and the convergence rate is enhanced. This may be used in two ways:

• to accelerate the calculation of $\lambda_n$,
• if you are happy with the approximation already obtained of $\lambda_n$, then you drop the last row and last colum and continue the QR algorithm with the remaining $(n-1)\times(n-1)$ block.

This is how black box software proceeds to compute the spectrum of a given matrix. This is also used to compute the roots of a polynomial, after having formed its companion matrix. Remark that the companion matrix is already of Hessenberg form.

Answer 4

As explained in Answer 3, the shift can be explained by the inverse power iteration. $x'=(A-rI)^{-1}x$ has larger component in $v_1$ if $r$ is closer to $\lambda_1$.

However in QR iteration, we don't use the inverse iteration, but the power iteration.

In the power iteration, $x'=(A-rI)x$ has smaller component in $v_1$ if $r$ is closer to $\lambda_1$. It seems the conclusion is opposite. I don't know what my misunderstanding is. I appreciated if anyone can help. Thanks.

Answer 5

A visualisation of the dynamics of the unshifted QR algorithm for PSD matrices is given here, with a remark as to how shifts help: https://www.youtube.com/watch?v=XGAaHoppcvE

The video is by me.

The argument is given for PSD matrices. For the set of PSD matrices similar to some PSD matrix $M$, there are multiple fixed points of the unshifted QR algorithm. Some of these fixed points are stable, while the others are unstable. The introduction of shifts eliminates the unstable fixed points, leaving behind only the stable fixed points. Unstable fixed points are undesirable because they compromise the speed of the algorithm. The existence of the unstable fixed point is an artifact of the continuousness of the function being iterated in the naive QR algorithm, as the video demonstrates by visualizing the naive LR algorithm as well. Therefore the shift must depend discontinuously on the matrix; but this does not compromise the practicality of employing shifts.

Feel free to leave a comment. Cheers!

[EDIT] Now part of a playlist on the QR algorithm's dynamics: https://www.youtube.com/playlist?list=PLvTRlE0qBBp35rv11Wp3WCWSixUisJqY4