Taking the power of a matrix is an important operation with applications in statistics, machine learning, and engineering. For example, solving linear ordinary differential equations, identifying the state of a Markov chain at time \(t\), or identifying the number of paths between nodes in a graph can all be solved using powers of matrices. In this quick post we’ll show how Matrix Diagonalization can be used to efficiently compute the power of a matrix.

If matrix \(M\) is an \(m \times m\) diagonalizable, then \(M^k\) can be calculated directly from the diagonalization \(M = P D P^{-1}\) as follows:

\[\begin{align} M^k &= M \times M \dots \times M \\ &= (P D P^{-1}) (P D P^{-1}) \dots (P D P^{-1}) \\ &= P D (P^{-1} P) D (P^{-1} P) \dots D P^{-1} \\ &= P D^k P^{-1} \end{align}\]

Therefore to calculate \(M^k\), we simply need to diagonalize \(M\) and re-matrix-multiply the diagonalization components after raising the diagonal matrix component \(D\) to the \(k\)-th power. Since \(D\) is a diagonal matrix, the \(k\)-th power is calculated by simply raising each element along the diagonal to the \(k\)-th power:

\[\begin{align} D^k &= \begin{bmatrix} d_{1} & & \\ & \ddots & \\ & & d_{m} \end{bmatrix}^k \\ &= \begin{bmatrix} d_{1}^k & & \\ & \ddots & \\ & & d_{m}^k \end{bmatrix} \end{align}\]

This trick allows us to calculate the matrix power by multiplying three matrices, rather than \(k\). Thus as \(k\) gets large, or the size of the matrix \(M\) grows, you get more and more gains in efficiency.

To demonstrate, let’s calculate the matrix power of a random matrix using brute force, the matrix diagonalization approach reviewed above, and we’ll also throw in results from numpy.linalg.matrix_power for completeness.

import numpy as np
np.random.seed(123)

# Generate a random 3 x 3 matrix
M = np.random.randn(3, 3)
k = 3  # power exponent

print('\nBrute Force:\n', eval("@".join([' M '] * k)))
# Brute Force:
#  [[-0.34077132 -0.70544947 -1.07778229]
#  [ 2.73462284 -0.71537115 -2.62514227]
#  [ 3.35955945  1.68986542 -4.1619396 ]]

# Diagonalize M via Eigenvalue Decomposition
D, P = np.linalg.eig(M)
D = np.diag(D)  # Put eigenvalues into a diagonal matrix

print('\nMatrix Diagonalization:\n', np.real(P @ D ** k @ np.linalg.inv(P)))
# Matrix Diagonalization:
#  [[-0.34077132 -0.70544947 -1.07778229]
#  [ 2.73462284 -0.71537115 -2.62514227]
#  [ 3.35955945  1.68986542 -4.1619396 ]]

print('\nnumpy.linalg.matrix_power:\n', np.linalg.matrix_power(M, k))
# numpy.linalg.matrix_power:
#  [[-0.34077132 -0.70544947 -1.07778229]
#  [ 2.73462284 -0.71537115 -2.62514227]
#  [ 3.35955945  1.68986542 -4.1619396 ]]

Works! 😁