How to use eigenvalues and eigenvectors to compute $A^{1000}$

Question

$$A= \begin{bmatrix} 0.9 & 0.15 & 0.25 \\ 0.075 & 0.8 & 0.25 \\ 0.025 & 0.05 & 0.5 \\ \end{bmatrix} $$

I need to use a Python script to compute $A^{1000}$.

I can use numpy's linalg.eig function to find my eigenvalues and eigen vectors. Doing so should yield eigenvalues

$$\begin{bmatrix} 1& 0.741 & 0.459 \\ \end{bmatrix} $$

and eigenvectors $$ \begin{bmatrix} 0.891 & 0.737 & -0.276 \\ 0.445 & -0.673 & -0.528 \\ 0.089 & -0.063 & 0.803 \\ \end{bmatrix} $$

Given this information, I am unsure how to proceed with computing $A^{1000}$

I have tried to use $A^{1000} = PD^{1000}P^{-1}$, where

$$P= \begin{bmatrix} 0.891 & 0.737 & -0.276 \\ 0.445 & -0.673 & -0.528 \\ 0.089 & -0.063 & 0.803 \\ \end{bmatrix} $$

$$P^{-1}= \begin{bmatrix} 0.702 & 0.702 & 0.702 \\ 0.495 & -0.904 & -0.424 \\ -0.039 & -0.149 & 1.134 \\ \end{bmatrix} $$

and

$$D= \begin{bmatrix} 1 & 0 & 0 \\ 0.445 & 0.741 & 0 \\ 0 & & 1.134 \\ \end{bmatrix} $$

Answer 1

You seem to understand all the steps, so you just need the syntax. Note that this kind of question should be asked on Stack Overflow.

Solution with exponentiation by squaring (using the numpy matrix_power function) and solution with eigenvalues and eigenvectors.

import numpy as np
import numpy.linalg

a = np.matrix([[0.9, 0.15, 0.25], [0.075, 0.8, 0.25], [0.025, 0.05, 0.5]])

np.linalg.matrix_power(a, 1000)

d, p = np.linalg.eig(a)
p * np.diag(d**1000) * np.linalg.inv(p)

For both the output is the following matrix (the norm of the difference is close to $1.6\cdot10^{-13}$).

matrix([[0.625 , 0.625 , 0.625 ],
        [0.3125, 0.3125, 0.3125],
        [0.0625, 0.0625, 0.0625]])