I have a matrix m that represents the probabilities transitioning from states to states.
E.g. for the sample below I will always get stuck in states 1,3,4, and state 2 I will randomly transition to one of the 4 states.
import numpy as np
m = np.eye(4)
m[1] = 0.25
print(m)
[[1. 0. 0. 0. ]
[0.25 0.25 0.25 0.25]
[0. 0. 1. 0. ]
[0. 0. 0. 1. ]]
How do I find a matrix representing the eventual end state following infinite transitions?
E.g. if I do this, I get intuitive result of states 1,3,4 --> 100% sticking in 1,3,4 but state 2 --> 1/3 chance ending up in all the others. Since all cases from state 2 eventually allocated evenly between 1,3,4 through multiple transitions.
t = m
for _ in range(100_000):
t = t @ t
print(t)
[[1. 0. 0. 0. ]
[0.33333333 0. 0.33333333 0.33333333]
[0. 0. 1. 0. ]
[0. 0. 0. 1. ]]
How can I calculate this without using repeated multiplications? I thought it corresponds to the eigenvector/eigenvalues of the matrix, but I get something very different when I calculate this.
np.linalg.eig(m)
[[0. , 0.9486833 , 0. , 0. ],
[1. , 0.31622777, 0.31622777, 0.31622777],
[0. , 0. , 0.9486833 , 0. ],
[0. , 0. , 0. , 0.9486833 ]]
Is there a methodology to calculate this using numpy? I need it to work for an arbitrary matrix, but there will be a known list of terminal states and positive probability of reaching these from all other states.
At the moment I am thinking of using the repeated multiplication method but it feels suboptimal and something there should be a function/trick that can calculate without looping.
I was reading this but didn't fully understand what the methodology is and how to implement it.
https://math.dartmouth.edu/archive/m20x06/public_html/Lecture14.pdf
I also looked in this question. People seemed to give some tips for hand-solving but not a general algorithm:
My friend pointed out the following trick.
Eigendecomposition means we can write the original matrix as
V x D x V^-1
Where D is a diagonal matrix with the eigenvalues, and V is the eigenvector.
If we multiply this by itself infinite times, it is
V x D^inf x V^-1
Which we can calculate in numpy using the below.
d, v = np.linalg.eig(m)
v @ np.diag(d >= 1).astype(int) @ np.linalg.inv(v)
Since if the diagonal values are < 1 they will tend to 0 as we multiply (assuming we have a matrix with valid probabilities, and all states can reach the terminal states).