# shifted QR algorithm--why does the shift help ?

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 ?

-
I am not specialist in the subject. By chance I read today "Francis's Algorithm" by David S. Watkins in the Monthly of May 2011. Seems related to your question... –  Luis H Gallardo May 16 '11 at 18:23
Thanks. I will check this article. –  vanvu May 16 '11 at 23:56
Watkins actually wrote a lot of pedagogical articles on QR and other eigenvalue methods; as I recall most of those were in SIAM Reviews... –  J. M. May 17 '11 at 0:37

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.

-
Well,that is exactly the point I would like to raise. How does one find $\mu$ close to $\lambda_1$ but (relatively) far from $\lambda_2$ while $\lambda_1$ and $\lambda_2$ are originally close to each other. If one can do this, that means one can more or less compute $\lambda_1$ itself, which is sort of the purpose of the algorithm. So the problem seems a bit circular, doesn't it ? –  vanvu May 16 '11 at 21:49
@vanvu: A common strategy (that gives local quadratic convergence (cubic for symmetric matrices)) is to use the last entry in the matrix (negated) as the shift. A better strategy is to use one of the eigenvalues of the last 2x2 block of the matrix. (IIRC, this leads to global convergence on symmetric matrices.) It makes sense if you think of the matrix as being already close to triangular, in which case we'd expect the diagonal to approximate the eigenvalues. –  Darsh Ranjan May 16 '11 at 22:16
@ranjan I understand that this is sort of strategy people choose to do in practice, but fail to see (mathematically) why ? There is no apparent reason that the last entry is a good $\mu$ in the above sense. Is there a good theorem here or mainly heuristic ? –  vanvu May 16 '11 at 23:39
@vanvu: Well, the Wilkinson shift is rather ad hoc (Jim Wilkinson says so himself); but it certainly works for most matrices "in the wild". In the symmetric case, Wilkinson showed that the convergence rate is at most cubic, so his heuristic works well here. In the unsymmetric case... it may well be that there are better shifts than Wilkinson's, but I've not heard about them and LAPACK still uses the Wilkinson shift along with a bunch of elaborate trickery (e.g. "aggressive deflation"). –  J. M. May 17 '11 at 0:35
I think you mean "at least cubic", J.M.? –  Federico Poloni May 17 '11 at 6:48

(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.

-
dx.doi.org/10.1016/0024-3795(68)90017-7 is a paper by Wilkinson where he demonstrates the cubical convergence of the QR algorithm in the symmetric tridiagonal case; both the Rayleigh shift and his shift lead to cubic convergence, but his shift has a slight edge. –  J. M. May 17 '11 at 3:19

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

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 an $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 proceed 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.

-
Thanks. I see it a bit clearer now. Btw, what happens if there are multiple eigenvalues ? –  vanvu May 17 '11 at 16:02
The QR algorithm does not treat efficiently eigenvalues that have the same modulus. Example 1: if $A$ is unitary (then all ev have modulus $1$), the iterates coincide with $A$. Example 2: if $A$ has real entries, the iterates remain real and cannot converge to a triangular form if there is at least one pair of complex conjugate ev. Both cases can be cured with a shift! I do not know what really happens when one ev has multiplicity, but this case cannot be avoided just by a shift. –  Denis Serre May 18 '11 at 10:50

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.

-
Why the shift helps can be explained by considering the QR of $H-rI=QR$. Since H is unreduced Hessenberg, it's first $n-1$ columns are linearly independent for any $r$. But if $r=\lambda_1$, $H-rI$ is singular. Then $r_{nn}=0$ and the last row of $R$ is zero. Then the last row of $H'+rI=RQ$ is zero except $nn$ component is $r$. So $r=\lambda_1$ or close to $\lambda_1$ will make $H'+rI$ $n-1,n$ component zero or closer to zero. This explains why the shift help. But this explanation doesn't fit well with the power iteration, where the component of $x'=(A-rI)x$ in $v_1$ actually decreases. –  Tongru Huo Jun 7 '11 at 17:29