10
$\begingroup$

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?

$\endgroup$
2
  • 3
    $\begingroup$ 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... $\endgroup$ May 16, 2011 at 18:23
  • 1
    $\begingroup$ Watkins actually wrote a lot of pedagogical articles on QR and other eigenvalue methods; as I recall most of those were in SIAM Reviews... $\endgroup$ May 17, 2011 at 0:37

5 Answers 5

10
$\begingroup$

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.

$\endgroup$
6
  • $\begingroup$ 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 ? $\endgroup$
    – vanvu
    May 16, 2011 at 21:49
  • $\begingroup$ @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. $\endgroup$ May 16, 2011 at 22:16
  • $\begingroup$ @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 ? $\endgroup$
    – vanvu
    May 16, 2011 at 23:39
  • $\begingroup$ @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"). $\endgroup$ May 17, 2011 at 0:35
  • $\begingroup$ I think you mean "at least cubic", J.M.? $\endgroup$ May 17, 2011 at 6:48
6
$\begingroup$

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

$\endgroup$
2
6
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ Thanks. I see it a bit clearer now. Btw, what happens if there are multiple eigenvalues ? $\endgroup$
    – vanvu
    May 17, 2011 at 16:02
  • 1
    $\begingroup$ 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. $\endgroup$ May 18, 2011 at 10:50
0
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ 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. $\endgroup$
    – Tongru Huo
    Jun 7, 2011 at 17:29
  • 1
    $\begingroup$ You refer to "Answer 3", but the order of answers can change over time (depending on which ordering you use). I guessed which answer you meant, based on its use of inverse power iteration—namely, @DarshRanjan's. I hope that was correct. $\endgroup$
    – LSpice
    Mar 16, 2022 at 18:43
0
$\begingroup$

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

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.