Discovering the QR Algorithm
The QR algorithm for finding eigenvalues is a deceptively simple piece of applied mathematics. Start with a square matrix \(\require{mathtools} A \coloneqq A_0\). Then, for the matrix \(A_k\), carry out a QR factorization \(A_k = Q_k R_k\). Set \(A_{k+1} = R_k Q_k\) and iterate to convergence. Magically, the result will be an upper-triangular matrix with the same eigenvalues as \(A\), and you can read the eigenvalues off the diagonal.
The goal of this post is to present an accessible proof of this remarkable fact. Rather than work backwards from the algorithm to the proof, we will discuss various methods of iteratively determining eigenvalues, until the QR algorithm falls out.
Step 1: Power Iteration
The standard way to get the dominant eigenvalue of \(A\) (that is, the eigenvalue with largest norm, assuming that one exists) is through power iteration. We start with a randomly chosen vector \(z_0\) and iteratively compute \(z_{k+1} = Az_k / \Vert Az_k \Vert\). In the case where \(A\) is diagonalizable with eigenvalues \(\vert \lambda_1 \vert > \vert \lambda_2 \vert > \cdots\), this converges geometrically to an eigenvector of \(A\) with eigenvalue \(\lambda_1\). To explain this, recall that, since \(A\) is diagonalizable, there’s a basis for the underlying space made up of eigenvectors. Writing \(z_0\) in terms of this basis and repeatedly applying \(A\) and normalizing will amplify the component parallel to the dominant eigenvector.
When \(A\) is not diagonalizable, this is slightly more complicated. However, by considering the Jordan canonical form, it can be shown that the power method will converge to the dominant eigenvector as long as one exists and there is some component of the starting vector that lies in the dominant generalized eigenspace. In short, in many cases, if you start with some vector and apply powers of \(A\), you’ll get an eigenvector of \(A\).
Step 2: Building Upper Triangular Matrices
So, if we can find one eigenvector, can we find them all? One potential way to do this, assuming that dominant eigenvectors exist throughout, would be to find the eigenvector \(v\), restrict \(A\) to \(\{v\}^\perp\) (the space of vectors perpendicular to \(v\)), find the dominant eigenvector on that, and repeat. The issue here is that \(\{v\}^\perp\) might not be invariant under \(A\). However, there is a variant of this idea that can be made rigorous. Our algorithm works like this.
- Find an eigenvector \(v_1\) of \(A\).
- For vectors \(v_1, \dots, v_k\), define \(U_k = \operatorname{span}(v_1, \dots, v_k)\) and let \(P_{U_k^\perp}\) be the projection onto \(U_k^\perp\).
- Since \(P_{U_k^\perp}A\) is an operator on \(U_k^\perp\), we can attempt to find a dominant eigenvector for it in \(U_k^\perp\). Call this eigenvector \(v_{k+1}\).
- Add \(v_{k+1}\) to our list of vectors, go back to step 2, and repeat until we have enough vectors to form a basis of the underlying vector space.
Crucially, since \(P_{U_k^\perp}Av_{k+1} = \lambda v_{k+1}\) for some \(\lambda\), \(Av_{k+1} \in U_{k+1}\). This is exactly the condition for \(A\) to be upper triangular in the basis \(v_1, \dots, v_n\). We can then recover the actual eigenvalues, and hence eigenvectors, of \(A\).
Although the power method does find eigenvectors, most of the eigenvectors we use it to find are not eigenvectors of \(A\) at all!
Step 3: Converging Multiple Vectors Simultaneously
Our method requires that we apply the power method to \(P_{U_k^\perp}A\), which requires applying \(P_{U_k^\perp}A\) to various vectors and then normalizing them. The most obvious way to do this is to compute \(P_{U_k^\perp}A\) as a matrix, but we don’t actually have to. Instead, we can apply \(A\), remove all components of the resulting vectors contained in \(U_k\), and normalize.
This is exactly what the Gram-Schmidt process does. Let \(z\) be some vector. Then we can apply \(P_{U_k^\perp}A\) to \(z\) by applying Gram-Schmidt to the list \(v_1, \dots, v_k, Az\) and taking the final vector of the result. Equivalently, we can consider the matrix with columns given by \(v_1, \dots, v_k, Az\), apply QR factorization, forget the R, and take the last column of the Q matrix.
The core to the QR algorithm, much as it might not look like it at first, is to converge all of these vectors simultaneously. Start with some matrix \(Q_0\). For any \(Q_k\), compute \(AQ_k\) and perform QR factorization to get \(Q_{k+1}R_{k+1}\). When we do this, the first column of \(Q_{k+1}\) is exactly the result of applying \(A\) to the first column of \(Q_k\) and normalizing. Therefore, under the usual power method assumptions, this column converges to \(v_1\). Once it has done so, or at least done so to within a negligible error, the second column experiences power method iteration and converges to an eigenvector of \(P_{U_1^\perp}A\). Once that has converged, the third column experiences power method iteration for \(P_{U_2^\perp}A\), and so on. In practice, rather than converging one at a time, all the vectors converge together very quickly. (There are analyses of the convergence rate of this method; I won’t go into them here.)
Step 4: An Algebraic Magic Trick
We now have a recurrence relation: \(AQ_{k-1} = Q_k R_k\). Furthermore, the actual matrix we want is \(A\) but written in the basis given by the columns of \(Q_k\). We call this \(A_k \coloneqq Q_k^* A Q_k\). Plugging in our recurrence relation,
\[A_{k-1} = Q_{k-1}^* A Q_{k-1} = Q_{k-1}^* Q_k R_k.\]Since \(Q_{k-1}^*\) and \(Q_k\) are both unitary, their product is as well, so this is a QR factorization of \(A_{k-1}\).
Now let’s try to compute \(A_k\). Using our recurrence relation once again, we can write \(A = Q_k R_k Q_{k-1}^*\). Therefore:
\[A_k = Q_k^* (Q_k R_k Q_{k-1}^*) Q_k = R_k (Q_{k-1}^* Q_k).\]This is exactly the iteration that we use in the QR algorithm: \(A_{k-1} = Q_k R_k\), \(A_k = R_k Q_k\). We’re done.