6.4 The QR Algorithm for Eigenvalues

Lecture notes

1Reading¶

Material related to this page, as well as additional exercises, can be found in ALA 8.1.

2Learning Objectives¶

By the end of this page, you should know:

the definition of similar matrices,
the definition of the Schur decomposition of a matrix,
how to use the QR algorithm to find the eigenvalues of a matrix.

3Motivations¶

Before we introduce the QR algorithm, we introduce a few motivating definitions.

3.1Similar Matrices¶

First, we introduce the notion of similar matrices:

A key property of similar matrices is that they have the same eigenvalues:

To see why this is true, suppose that λ is an eigenvalue of $B$ , i.e., $B\vv v = \lambda \vv v$ for some nonzero $\vv v$ . Then,

\begin{align*} B \vv v = \lambda \vv v &\implies (P^{-1}AP) \vv v = \lambda \vv v\\ &\implies A(P \vv v) = P (\lambda \vv v) = \lambda (P\vv v) \end{align*}

(1)

meaning that $P\vv v$ is an eigenvector of $A$ , with the same eigenvalue of λ.

This motivates the following method: if we want to find the eigenvalues of $B$ , we could try to write it in the form $A = P^{-1}BP$ , where $B$ is a matrix for which it is easy to find the eigenvalues! This is exactly the motivation for the QR algorithm, which is an iterative algorithm for finding the Schur decomposition of $A$ .

We will revisit similar matrices in more detail a few sections down the line.

3.2The (Real) Schur Decomposition¶

Next, we define the Schur decomposition. The proof of this claim relies on material covered later in the course, so we won’t prove it here.

Note that, given a real Schur decomposition $A = Q^\top U Q$ , it’s easy to find the eigenvalues of $A$ ; we just find the eigenvalues of each diagonal block! To see why, recall Determinant Fact 4. Since $U$ is quasi-upper triangular, its diagonal blocks are all $1\times 1$ or $2\times 2$ matrices, which have easy to solve characteristic equations, so it’s easy to find their eigenvalues.

4The QR Algorithm¶

Equipped with these facts, the motivation behind the QR algorithm is as follows.

First, we find an approximation to a Schur decomposition $A = Q^\top U Q$ (where $Q$ is orthogonal, $U$ is quasi-upper triangular).
The eigenvalues of $A$ are thus the diagonal entries of $U$ .

4.1Approximating a Real Schur Decomposition¶

The question thus becomes, how do we find an approximate Schur decomposition for $A$ ? One idea is repeatedly apply the QR decomposition to iteratively find similar matrices to $A$ . In the case that all eigenvalues are distinct, this actually converges to the Schur decomposition.

Initialize:
$A_1 = A$
(3)
For $t = 1, 2, \dots$ , perform the following steps:
a. Compute the QR-factorization of $A_t$ :
$Q_t R_t = A_t$
(4)
where $Q_t$ is orthogonal and $R_t$ is upper triangular.
b. Define $A_{t+1}$ by reversing the factors:
$A_{t+1} = R_t Q_t = (Q_t^\top A_t) Q_t$
(5)

The mechanism of why this algorithm works can be described by two phenomena. First, each subsequent $A_{t+1}$ is similar to the previous $A_t$ so we know that eigenvalues are preserved over each iteration. Second, in each interation the elements below the diagonal get smaller and vice versa for those above (showing why this property holds is beyond the scope of this course). This means that over many iterations the elements below the diagonal effectively become 0, making $A_\infty$ a upper-triangular matrix.

5Pseudocode for the QR Algorithm¶

We are now ready to give the pseudocode for the QR Algorithm.

6Python Implementations of the QR Algorithm¶

We are now ready to use the QR Algorithm to find eigenvalues in Python.

6.1From Scratch¶

First, we’ll give an implementation of the Schur decomposition from scratch. We won’t reimplement a QR factorization algorithm, so if you’re interested in seeing the implementation for that, it can be found on this page (scroll down).

import numpy as np

def schur_form(A, iters=100): # A is a square matrix
    B = A
    for i in range(iters):
        Q, R = np.linalg.qr(B)
        B = R @ Q
    return B

print(np.round(schur_form(np.array([
    [2.0, 3.0, 1.0, 0.5, 4.0],
    [4.0, 5.0, 7.0, 0.1, 1.0],
    [5.0, 3.0, 6.0, 19.2, 9.0],
    [1.0, 4.0, 1.0, 4.0, 7.0],
    [3.0, 1.0, 6.0, 2.0, 6.0]
])), 2))

[[21.36 -6.11  9.82  4.82 -2.86]
 [ 0.   -2.14  8.72  3.1  -6.6 ]
 [-0.   -5.12 -1.82  1.99 -0.84]
 [ 0.   -0.    0.    4.34 -0.32]
 [ 0.   -0.    0.    0.    1.26]]

6.2Using `numpy.linalg.eigvals`¶

Next, we’ll use the numpy.linalg.eigvals function, which returns the eigenvalues of a square matrix. Under the hood, numpy.linalg.eigvals is based off of a more advanced version of the basic QR algorithm we outlined above.

import numpy as np

print(np.round(np.linalg.eigvals(np.array([
    [2.0, 3.0, 1.0, 0.5, 4.0],
    [4.0, 5.0, 7.0, 0.1, 1.0],
    [5.0, 3.0, 6.0, 19.2, 9.0],
    [1.0, 4.0, 1.0, 4.0, 7.0],
    [3.0, 1.0, 6.0, 2.0, 6.0]
])).T, 2))

[21.36+0.j   -1.98+6.68j -1.98-6.68j  1.26+0.j    4.34+0.j  ]

Verify for yourself that the real eigenvalues returned by numpy.linalg.eigvals are the diagonal $1\times 1$ blocks of the Schur form of the input matrix, which we found above. You can also verify that the complex eigenvalues returned by numpy.linalg.eigvals are the eigenvalues of the $2\times 2$ block diagonal.