Notes - Numerical Analysis HT24, QR algorithm

• Set $A _ 1 = A$
• For $k = 1, 2, \cdots$
  • Calculate the $QR$ factorisation $A _ k = Q _ k R _ k$
  • Set $A _ {k+1} = R _ k Q _ k$

Then $A _ k$ converges to a diagonal matrix with the eigenvalues on the diagonal (if $A$ is not symmetric but still has eigenvalues of distinct modulus, then $A _ k$ converges to an upper triangular matrix. If the eigenvalues aren’t distinct, $A _ k$ might not converge at all).

In general, if you know the eigenvalues of a matrix, then it is straightforward to find eigenvectors by solving a linear system. (And likewise, if you know the eigenvectors, you can find the eigenvalues).

But in the case where $A$ is symmetric and converges to a diagonal matrix, if you let $Q^{(k)} = Q _ 1 \cdots Q _ k$, then

since $A _ k$ is diagonal (or almost diagonal), it follows that $Q^{(k)}$ consists of an orthonormal basis consisting of $A$’s eigenvectors.

Use induction on $k$. $A _ 1$ is similar to $A$ by definition. Then for $A _ {k+1}$:

and this is a similarity transformation.

• First reduce $A$ to tridiagonal form by orthogonal similarity transformations; this reduces the complexity of calculating the QR at each step. This can be done using Givens rotations or Householder reflectors.
• Use shifts, to instead calculate the eigenvalues of $A + sI$ for some $s$ chosen so that the convergence ratio is higher. A typical choice is the bottom right entry of the matrix, since this will approximately equal the smallest eigenvalue of the matrix, making the of the final entry to the smallest eigenvalue higher.
• “Deflate” the matrix $A$: Each time the bottom row and last column converge to $(0, 0, \cdots, \lambda)$ within some threshold, remove them from $A$. This reduces the size of the matrix, speeding up matrix operations, and it also means that when we pick a new shift using $a _ {(n-1), (n-1)}$, the speed of the convergence now improves for the second smallest eigenvalue, rather than for the smallest eigenvalue.

A full algorithm:

Reduce $A$ to tridiagonal form by orthogonal similarity transformations.

For $m = n, n -1, \cdots, 2$:

• While $a _ {m-1, m} >$ tolerance: (this detects when we should deflate)
  • $[Q, R] = \mathbf{QRFactorise}(A - a _ {m, m}I)$
  • $A = RQ + a _ {m, m} I$ (note this process preserves the fact the $A _ k$ are similar)
• Record eigenvalue $\lambda _ m = a _ {m, m}$
• Set $A$ to the leading $m-1$ by $m-1$ submatrix of $A$

Finally: record eigenvalue $\lambda _ 1 = a _ {1, 1}$.

We do this in several steps, each step building off the last:

1. Showing $A _ k$ are similar for all $k$
2. Showing $A _ {k+1} = (Q^{(k)})^\top A Q^{(k)}$, where $Q^{(k)} = Q _ 1 \cdots Q _ k$
3. Showing $A^k = Q^{(k)}R^{(k)}$, where $R^{(k)} = R _ k \cdots R _ 1$ (note $R^{(k)}$ is reversed)
4. Showing $(A _ k) _ {1,1}$ must converge to the largest eigenvalue (in modulus) by relating it to $A^k e _ 1$
5. Showing $(A _ k) _ {n,n}$ must converge to the smallest eigenvalue (in modulus) by relating it to $A^{-k}e _ n$
6. Finally showing that the rest of entries of $A$ must also converge

Showing $A _ k$ are similar for all $k$: Use induction on $k$. $A _ 1$ is similar to $A$ by definition. Then for $A _ {k+1}$:

Showing $A _ {k+1} = (Q^{(k)})^\top A Q^{(k)}$, where $Q^{(k)} = Q _ 1 \cdots Q _ k$: This follows from the above proposition by induction.

Showing $A^k = Q^{(k)}R^{(k)}$, where $R^{(k)} = R _ k \cdots R _ 1$: Again, use induction on $k$. $A^1 = Q^{(1)} R^{(1)}$ by definition. Then, for $A^{k+1}$.

Where does $(\star)$ come from? We have that $A _ {k+1} = (Q^{(k)})^\top A Q^{(k)}$ from the previos proposition. But $A _ {k+1}$ also equals $Q _ {k+1} R _ {k+1}$. Then

which, by premultiplying both sides by $Q^{(k)}$, we see

Showing $(A _ k) _ {1,1}$ must be the dominant eigenvalue by relating it to $A^k e _ 1$: Consider

since $R^{(k)}$ is upper-triangular, $R^{(k)} e _ 1$ will equal $[R^{(k)} _ {1,1}, 0, \cdots, 0] = R^{(k)} _ {1,1} e _ 1$ and so $Q^{(k)} R^{(k)} e _ 1$ is $R _ {1,1}^{(k)} q^{(k)} _ 1$ where $q _ 1^{(k)}$ is the first column of $Q^{(k)}$. But from the analysis of the power method, $A^k e _ 1$ must be parallel to the dominant eigenvector of $A$.

But now considering $A _ k$ (subscript instead of superscript), we know

so in fact $A _ k$ must have first entry approximately equal to $\lambda _ 1$, which can be seen by considering $Q^{(k-1)} A _ k = A Q^{(k-1)}$. We also know that the rate of convergence depends on $\left(\frac{\lambda _ 2}{\lambda _ 1}\right)^k$, where $\lambda _ 1$ is the largest eigenvalue in modulus, and $\lambda _ 2$ is the second largest in modulus.

(One concern here is that the analysis of the power method shows convergence to $\pm v _ 1$ where $v _ 1$ is the dominant eigenvector, but this doesn’t matter because $-v _ 1$ is still “the” dominant of eigenvector of $A$. Note also that here we are just using facts about $A^k$ to deduce that the first entry of $A _ k$ is $\lambda$, we never need to explicitly calculate $A _ k$).

Showing $(A _ k) _ {n,n}$ must converge to the smallest eigenvalue (in modulus) by relating it to $A^{-k}e _ n$: Now consider

where $c = (R^{(k)})^{-1} _ {n,n}$. But $A^{-k} e _ n$ is the inverse power method, so $q _ n^{(k)}$ must actually converge to the eigenvector corresponding to the smallest eigenvalue, and $(A _ k) _ {n,n}$ must converge to $\lambda _ n$. Furthermore, the rate of convergence depends on

Finally showing that the rest of entries of $A$ must also converge: From the above, we know that for large $k$, $A _ k$ looks approximately like

by the symmetry of $A _ k$. But because of all the zeroes, we can repeat the analysis on $A _ k’$ to see that the rest of the entries must also converge to a diagonal matrix (admittedly this isn’t very precise, but it is also not very precise in the course notes).