6.6. Calculating eigenvalues using QR decomposition#

Eigenvalues and eigenvectors feature prominently in the study of numerical methods for ODEs. Given a system of ODEs, the eigenvalues of the coefficient matrix provide information about the stability, divergence, oscillatory behavior, and constant solutions of the system.

Definition 6.5 (Eigenvalue)

Let \(A\) be an \(n \times n\) matrix then \(\lambda\) is an eigenvalue of \(A\) if there exists a non-zero vector \(\vec{v}\) such that

(6.4)#\[A \vec{v} = \lambda \vec{v}.\]

The vector \(\vec{v}\) is called the eigenvector associated with the eigenvalue \(\lambda\).

Rearranging equation (6.4) we have

\[\begin{split} \begin{align*} A \vec{v} &= \lambda \vec{v} \\ (A - \lambda I) \vec{v} &= \vec{0}, \end{align*} \end{split}\]

which has non-zero solutions for \(\vec{v}\) if and only if the determinant of the matrix \((A - \lambda I)\) is zero. Therefore we can calculate the eigenvaleus of \(A\) using

(6.5)#\[ \det(A - \lambda I) = 0.\]

For example, consider the eigenvalues of the matrix

\[\begin{split} A = \begin{pmatrix} 2 & 1 \\ 2 & 3 \end{pmatrix}.\end{split}\]

Using equation (6.5)

\[\begin{split} \begin{align*} \det \begin{pmatrix} 2 - \lambda & 1 \\ 2 & 3 - \lambda \end{pmatrix} &= 0 \\ \lambda^2 - 5 \lambda + 4 &= 0 \\ (\lambda - 1)(\lambda - 4) &= 0 \end{align*} \end{split}\]

so the eigenvalues are \(\lambda_1 = 4\) and \(\lambda_2 = 1\) [1]. The problem with using determinants to calculate eigenvalues is that it is too computationally expensive for larger matrices.

6.6.1. The QR algorithm#

The QR algorithm is a method of computing the eigenvalues of a square matrix. Let \(A_0\) be a matrix for which we wish to compute the eigenvalues and \(A_1, \ldots, A_k, A_{k+1}, \ldots\) be a sequence of matrices such that the \(k\)-th matrix in the sequence is calculated using \(A_{k+1} = R_kQ_k\) where \(Q_kR_k = A_k\) is the QR decomposition of \(A_k\). Since \(Q_k\) is orthogonal then \(Q^{-1} = Q^T\) and

\[ A_{k+1} = R_kQ_k = Q_k^{-1}Q_kR_kQ_k = Q^{-1}A_kQ_k = Q_k^\mathrm{T}A_kQ_k\]

The matrices \(Q_k^\mathrm{T}A_kQ_k\) and \(A_{k}\) are similar meaning that they have the same eigenvalues. As \(k\) gets larger the matrix \(A_k\) will converge to an upper triangular matrix where the diagonal elements contain eigenvalues of \(A_k\) (and therefore \(A_0\))

\[\begin{split} A_{k} = R_kQ_k = \begin{pmatrix} \lambda_1 & \star & \cdots & \star \\ 0 & \lambda_2 & \ddots & \vdots \\ \vdots & \ddots & \ddots & \star \\ 0 & \cdots & 0 & \lambda_n \end{pmatrix}. \end{split}\]

Algorithm 6.9 (The QR algorithm)

Inputs: An \(n \times n\) matrix \(A\) and an accuracy tolerance \(tol\).

Outputs: A vector \((\lambda_1, \lambda_2, \ldots, \lambda_n)\) containing the eigenvalues of \(A\).

  • For \(k = 1, 2, \ldots\) do

    • Calculate the QR decomposition of \(A\)

    • \(A_{old} \gets A\)

    • \(A \gets R Q\)

    • If \(\max(\operatorname{diag}(|A - A_{old}|)) < tol\)

      • Break

  • Return \(\operatorname{diag}(A)\)

Example 6.10

Use the QR algorithm to compute the eigenvalues of the matrix

\[\begin{split} A = \begin{pmatrix} 1 & 2 \\ 3 & 4 \end{pmatrix}, \end{split}\]

using an accuracy tolerance of \(tol = 10^{-4}\).

Solution (click to show)

Calculate the QR decomposition of \(A_0\)

\[\begin{split} \begin{align*} Q_0 &= \begin{pmatrix} -0.3162 & 0.9487 \\ -0.9487 & -0.3162 \end{pmatrix}, \\ R_0 &= \begin{pmatrix} -3.1623 & -4.4272 \\ 0 & 0.6325 \end{pmatrix} \end{align*} \end{split}\]

and calculate \(A_1 = R_0Q_0\)

\[\begin{split} \begin{align*} A_1 &= \begin{pmatrix} -3.1623 & -4.4272 \\ 0 & 0.6325 \end{pmatrix} \begin{pmatrix} -0.3162 & 0.9487 \\ -0.9487 & -0.3162 \end{pmatrix} \\ &= \begin{pmatrix} 5.2 & -1.6 \\ -0.6 & -0.2 \end{pmatrix} \end{align*} \end{split}\]

Calculate the QR decomposition of \(A_1\)

\[\begin{split} \begin{align*} Q_1 &= \begin{pmatrix} -0.9943 & -0.1146 \\ 0.1146 & -0.9934 \end{pmatrix}, \\ R_1 &= \begin{pmatrix} -5.2345 & 1.5665 \\ 0 & 0.3821 \end{pmatrix} \end{align*} \end{split}\]

and calculate \(A_2 = R_1Q_1\)

\[\begin{split} \begin{align*} A_2 &= \begin{pmatrix} -5.2345 & 1.5665 \\ 0 & 0.3821 \end{pmatrix} \begin{pmatrix} -0.9943 & -0.1146 \\ 0.1146 & -0.9934 \end{pmatrix} \\ &= \begin{pmatrix} 5.3796 & -0.9562 \\ 0.0438 & -0.3796 \end{pmatrix} . \end{align*} \end{split}\]

Calculate the maximum difference between the diagonal elements of \(A_1\) and \(A_2\)

\[\begin{split} \begin{align*} \max(|\operatorname{diag}(A_2 - A_1)|) &= \max \left ( \left| \begin{pmatrix} 5.3796 \\ -0.3796 \end{pmatrix} - \begin{pmatrix} 5.2 \\ -0.2 \end{pmatrix} \right| \right) \\ &= \max \begin{pmatrix} 0.1796 \\ 0.1796 \end{pmatrix} = 0.1796, \end{align*} \end{split}\]

since \(0.1796 > 10^{-4}\) we need to continue to iterate. The estimates of the eigenvalues obtained by iterating to an accuracy tolerance of \(10^{-4}\) are tabulated below.

\(k\)

\(\lambda_1\)

\(\lambda_2\)

Max difference

0

5.200000

-0.200000

4.20e+00

1

5.379562

-0.379562

1.80e-01

2

5.371753

-0.371753

7.81e-03

3

5.372318

-0.372318

5.65e-04

4

5.372279

-0.372279

3.90e-05

The exact values of the eigenvalues are \(\lambda_1 = (5 + \sqrt{33})/2 \approx 5.372281\) and \(\lambda_2 = (5 - \sqrt{33}) / 2 \approx -0.372281\) which shows the QR algorithm has calculated the eigenvalues correct to five decimal places.

6.6.2. Code#

The code below defines a function called eigvals() that uses the QR algorithm to compute the eigenvalues of a matrix.

def eigvals(A, tol=1e-6):
    for k in range(20):
        Q, R = qr_householder(A)
        A, Aold = np.matmul(R, Q), A
        if max(abs(np.diagonal(A - Aold))) < tol:
            break

    return np.diagonal(A)

function lambda = eigvals(A, tol)

for k = 1 : 20
    [Q, R] = qr_householder(A);
    Aold = A;
    A = R * Q;
    if max(abs(diag(A - Aold))) < tol
        break
    end
end

lambda = diag(A);

end