8.5. Diagonalization and Powers of A

Here, we will establish a factoring of a matrix based on its eigenvalues and eigenvectors. This factoring is called both diagonalization and eigen-decomposition. Then, we will use the factoring to find a computationally efficient way to compute powers of the matrix (\mathbf{A}^k). Then in Change of Basis and Difference Equations, we will extend this to \mathbf{A}^k\,\bm{v}, where \bm{v} is a vector. This final result has direct application to difference equations, which are used to solve several types of problems.

8.5.1. Diagonalization

The eigenpairs of the n{\times}n matrix \bf{A} give us n equations.

\begin{array}{rl} \mathbf{A}\,\bm{x}_1 &= \lambda_1\,\bm{x}_1 \\ \mathbf{A}\,\bm{x}_2 &= \lambda_2\,\bm{x}_2 \\ &\vdots \\ \mathbf{A}\,\bm{x}_n &= \lambda_n\,\bm{x}_n \\ \end{array}

We wish to combine the n equations into one matrix equation.

Consider the product of \bf{A} and a matrix called \bf{X} containing the eigenvectors of \bf{A} as the columns.

\begin{array}{ll} \mathbf{A\,X} &= \mathbf{A} \spalignmat{\vertbar{} \vertbar{} {} \vertbar{}; \bm{x}_1 \bm{x}_2 \cdots{} \bm{x}_n; \vertbar{} \vertbar{} {} \vertbar{} } \\ &{} \\ &= \spalignmat{\vertbar{} \vertbar{} {} \vertbar{}; {\lambda_1\,\bm{x}_1} {\lambda_2\,\bm{x}_2} \cdots{} {\lambda_n\,\bm{x}_n}; \vertbar{} \vertbar{} {} \vertbar{} } \\ &= \spalignmat{\vertbar{} \vertbar{} {} \vertbar{}; \bm{x}_1 \bm{x}_2 \cdots{} \bm{x}_n; \vertbar{} \vertbar{} {} \vertbar{} } \spalignmat{\lambda_1 0 \cdots{} 0; 0 \lambda_2 \cdots{} 0; \vdots{} \vdots{} \ddots{} \vdots{}; 0 0 \cdots{} \lambda_n } \\ &= \mathbf{X\,\Lambda} \end{array}

The matrix \bf{\Lambda} is a diagonal eigenvalue matrix. If the matrix has linearly independent eigenvectors, which will be the case when all eigenvalues are different (no repeating \lambdas), then \bf{X} is invertible.

\begin{array}{c} \mathbf{A\,X} = \mathbf{X\,\Lambda} \\ \mathbf{X}^{-1}\mathbf{A\,X} = \mathbf{\Lambda} \\ {}\\ \boxed{\mathbf{A} = \mathbf{X}\,\mathbf{\Lambda}\,\mathbf{X}^{-1}} \end{array}

This factoring of \bf{A} is called the diagonalization factoring.

8.5.1.1. When does Diagonalization not work?

Notice that for the diagonalization factoring to work we need to find the inverse of the eigenvector matrix. So each eigenvector must be independent of the other eigenvectors.

One example of a matrix with repeating eigenvectors is a 2{\times}2 matrix with the same values on the forward diagonal and a zero on the backward diagonal. The determinant of the characteristic matrix has repeating eigenvalues, so it will also have repeating eigenvectors. Thus the eigenvector matrix is singular and may not be inverted.

\spaligndelims\vert\vert \spalignmat{{a-\lambda}, n; 0, {a-\lambda}} = (a - \lambda) (a - \lambda)

>> A = [5 7; 0 5]
A =
    5     7
    0     5
>> [X, L] = eig(A)
X =
    1.0000   -1.0000
         0    0.0000
L =
    5     0
    0     5
>> rank(X)
ans =
    1

8.5.1.2. Diagonalization of a Symmetric Matrix

Symmetric matrices have a simplified diagonalization because the matrix of eigenvectors of a symmetric matrix, \bf{S}, is orthogonal (The Schur Decomposition of Symmetric Matrices). Recall also from Matrix Transpose Properties that from the spectral theorem, orthogonal matrices have the property \mathbf{Q}^{T} = \mathbf{Q}^{-1}. Thus the diagonalization of a symmetric matrix is shown below and is also called the spectral decomposition of a symmetric matrix.

\boxed{\mathbf{S} = \mathbf{Q\,\Lambda\,Q}^T}.

Note

Recall that the columns of orthonormal matrices must be unit vectors (length of 1), which the eig function returns. Orthonormal columns are required to use a transposed matrix instead of the inverse matrix.

8.5.2. Powers of A

How does one compute matrix \bf{A} raised to the power k (\mathbf{A}^k)?

The brute force approach is to simply multiply \bf{A} by itself k - 1 times. This could be slow if k and the size of the matrix (n) are large.

If \bf{A} is a n{\times}n square matrix, the product \mathbf{A}\,\mathbf{A} requires n^3 multiply and addition operations (\mathcal{O}(n^3)). Thus \mathbf{A}^k has complexity of \mathcal{O}\left((k - 1)n^3\right); or if we are clever about the order of the multiplications, it could be reduced to \mathcal{O}\left(\log_{2}(k)n^3\right). That is:

\mathbf{A}^k = \underbrace{\underbrace{\underbrace{ \mathbf{A\,A}}_{\mathbf{A}^2}\, \mathbf{A}^2}_{\mathbf{A}^4} \, \mathbf{A}^4 \ldots \mathbf{A}^{k/2}}_{\mathbf{A}^{k}}.

Fortunately, diagonalization gives us a faster way to compute \mathbf{A}^k as a stand-alone matrix. The complexity of this method is \mathcal{O}(n^3).

\begin{array}{rcl} \mathbf{A}^2 &= \mathbf{X\,\Lambda\,X}^{-1}\: \mathbf{X\,\Lambda\,X}^{-1} &= \mathbf{X}\,\mathbf{\Lambda}^2\,\mathbf{X}^{-1} \\ \mathbf{A}^3 &= \mathbf{X}\,\mathbf{\Lambda}^2\,\mathbf{X}^{-1}\: \mathbf{X\,\Lambda\,X}^{-1} &= \mathbf{X}\,\mathbf{\Lambda}^3\,\mathbf{X}^{-1} \\ \end{array}

\boxed{\mathbf{A}^k = \mathbf{X\,\Lambda}^k\,\mathbf{X}^{-1}}

Because the \bf{\Lambda} matrix is diagonal, only the individual \lambda values need be raised to the k power (element-wise exponent).

\mathbf{\Lambda}^k = \spalignmat{\lambda_1^k 0 \cdots{} 0; 0 \lambda_2^k \cdots{} 0; \vdots{} \vdots{} \ddots{} \vdots{}; 0 0 \cdots{} \lambda_n^k}

Example Power of A

In this example, the matrix is symmetric, so the inverse simplifies to a transpose.

>> S = gallery('moler', 4)
S =
     1    -1    -1    -1
    -1     2     0     0
    -1     0     3     1
    -1     0     1     4
>> [Q, L] = eig(S)
Q =
    -0.8550    0.1910    0.3406   -0.3412
    -0.4348   -0.6421   -0.6219    0.1093
    -0.2356    0.6838   -0.4490    0.5248
    -0.1562   -0.2894    0.5437    0.7722
L =
    0.0334         0         0         0
         0    2.2974         0         0
         0         0    2.5477         0
         0         0         0    5.1215

>> Q * L.^5 * Q'
ans =
        425        -162        -639        -912
       -162         110         204         273
       -639         204        1022        1389
       -912         273        1389        2138
>> S^5
ans =
        425        -162        -639        -912
       -162         110         204         273
       -639         204        1022        1389
       -912         273        1389        2138

A few observations on computing powers of \bf{A}

  • Considerable computation is saved by using the element-wise exponent operator on \bf{\Lambda} (L in the MATLAB code). This gives the same result as a matrix exponent because \bf{\Lambda} is a diagonal matrix.
  • The eigenvalues and eigenvectors of a matrix are often complex. In the example, a symmetric matrix is used because such will always have only real eigenvalues and eigenvectors. When a matrix has complex eigenvalues, the method still works. After multiplying the eigenvalues by the eigenvectors, the imaginary values may not be exactly zero due to round-off errors, but they will be very small. Use the real function to get only the real component of the complex numbers.