.. _powers: Diagonalization and Powers of **A** =================================== Here, we will establish a factoring of a matrix based on its eigenvalues and eigenvectors. The factoring is called *diagonalization* or the *eigendecomposition*. Then, we will use the factoring to find a computationally efficient way to compute powers of the matrix (:math:`\mathbf{A}^k`). .. index:: eigen-decomposition .. index:: diagonalization .. index:: powers of matrix .. _diagonalization: Diagonalization --------------- The eigenpairs of the :math:`n{\times}n` matrix :math:`\bf{A}` give us :math:`n` equations. .. math:: \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 :math:`n` equations into one matrix equation. Consider the product of :math:`\bf{A}` and the :math:`\bf{X}` matrix containing the eigenvectors. .. math:: \begin{array}{ll} \mathbf{A\,X} &= \mathbf{A} \begin{bmatrix} \vertbar{} & \vertbar{} & {} & \vertbar{} \\ \bm{x}_1 & \bm{x}_2 & \cdots{} & \bm{x}_n \\ \vertbar{} & \vertbar{} & {} & \vertbar{} \end{bmatrix} \\ &{} \\ &= \begin{bmatrix} \vertbar{} & \vertbar{} & {} & \vertbar{} \\ {\lambda_1\,\bm{x}_1} & {\lambda_2\,\bm{x}_2} & \cdots{} & {\lambda_n\,\bm{x}_n} \\ \vertbar{} & \vertbar{} & {} & \vertbar{} \end{bmatrix} \\[18pt] &= \begin{bmatrix} \vertbar{} & \vertbar{} & {} & \vertbar{} \\ \bm{x}_1 & \bm{x}_2 & \cdots{} & \bm{x}_n \\ \vertbar{} & \vertbar{} & {} & \vertbar{} \end{bmatrix} \begin{bmatrix} \lambda_1 & 0 & \cdots{} & 0 \\ 0 & \lambda_2 & \cdots{} & 0 \\ \vdots{} & \vdots{} & \ddots{} & \vdots{} \\ 0 & 0 & \cdots{} & \lambda_n \end{bmatrix} \\ &= \mathbf{X\,\Lambda} \end{array} The matrix :math:`\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 :math:`\lambda`\ s), then :math:`\bf{X}` is invertible. .. math:: \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} When Does Diagonalization Not Work? ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Notice that the inverse of the eigenvector matrix must exist for the eigendecomposition to be valid. So each eigenvector *must be independent* of the other eigenvectors. One example of a matrix with repeating eigenvectors is a :math:`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 roots, so it will also have repeating eigenvectors. Thus, the eigenvector matrix is singular and may not be inverted. .. math:: \begin{bmatrix} {a-\lambda} & n \\ 0 & {a-\lambda} \end{bmatrix} = (a - \lambda) (a - \lambda) :: In [36]: A = np.array([[5, 7], [0, 5]]) In [37]: L, X = eig(A); print(X) [[ 1. -1.] [ 0. 0.]] In [39]: np.linalg.matrix_rank(X) Out[39]: np.int64(1) .. _symdiag: Diagonalization of a Symmetric Matrix ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Symmetric matrices have a simplified diagonalization because the matrix of eigenvectors of a symmetric matrix, :math:`\bf{S}`, are orthogonal (:ref:`th-spect-decomp`). Recall also from :ref:`matTranspose` that from the spectral theorem, orthogonal matrices have the property :math:`\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. .. math:: \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. .. index:: spectral decomposition Powers of A ----------- How does one compute matrix :math:`\bf{A}` raised to the power :math:`k` (:math:`\mathbf{A}^k`)? The brute force approach is to multiply :math:`\bf{A}` by itself :math:`k-1` times, which is slow if :math:`k` and the size of the matrix (:math:`n`) are large. If we are clever about the order of the multiplications, it could be reduced to :math:`\log_{2}(k)` matrix multiplications. That is: .. math:: \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}}. Diagonalization allows us to compute :math:`\mathbf{A}^k` faster. .. math:: \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} .. math:: \boxed{\mathbf{A}^k = \mathbf{X\,\Lambda}^k\,\mathbf{X}^{-1}} Because the :math:`\bf{\Lambda}` matrix is diagonal, only the individual :math:`\lambda` values need be raised to the :math:`k` power (element-wise exponent). .. math:: \mathbf{\Lambda}^k = \begin{bmatrix} \lambda_1^k & 0 & \cdots{} & 0 \\ 0 & \lambda_2^k & \cdots{} & 0 \\ \vdots{} & \vdots{} & \ddots{} & \vdots{} \\ 0 & 0 & \cdots{} & \lambda_n^k \end{bmatrix} **Example Power of A** In this example, the matrix is symmetric, so the inverse simplifies to a transpose. :: In [54]: A = 2*np.random.randn(4, 4) In [55]: S = A.T @ A In [56]: L, Q = eig(S) In [57]: L = np.diag(L) In [58]: (Q @ L**5 @ Q.T) - S@S@S@S@S Out[58]: array([[-0.+0.j, 0.+0.j, -0.+0.j, -0.+0.j], [ 0.+0.j, 0.+0.j, -0.+0.j, -0.+0.j], [-0.+0.j, -0.+0.j, 0.+0.j, 0.+0.j], [-0.+0.j, -0.+0.j, 0.+0.j, 0.+0.j]]) .. topic:: A few observations on computing powers of :math:`\bf{A}` - Considerable computation is saved by using the element-wise exponent operator on :math:`\bf{\Lambda}` (``L`` in the code). - The eigenvalues and eigenvectors of a matrix are often complex. When a matrix has complex eigenvalues, the method for computing powers of :math:`\bf{A}` 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. .. index:: real