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

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.

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.