4.6. Powers of A

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

The brute force approach is to multiply \(\bf{A}\) by itself \(k-1\) times, which is slow if \(k\) and the size of the matrix (\(n\)) are large. If we are clever about the order of the multiplications, it could be reduced to \(\log_{2}(k)\) matrix multiplications. 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}}.\]

Diagonalization allows us to compute \(\mathbf{A}^k\) faster.

\[\begin{split}\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}\end{split}\]

\[\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).

\[\begin{split}\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}\end{split}\]

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]])

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 code).

• The eigenvalues and eigenvectors of a matrix are often complex. When a matrix has complex eigenvalues, the method for computing powers of \(\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.