2

Consider a Markovian process with a diagonalizable transition matrix A such that A=PDP^-1, where D is a diagonal matrix with eigenvalues of A, and P is a matrix whose columns are eigenvectors of A.

To compute, for each state, the likelihood of ending up in each absorbing state, I'd like to raise the transition matrix to the power of n, with n approaching infinity:

A^n=PD^nP^-1

Which is the pythonic way of doing this in Numpy? I could naively compute the eigenvalues and eigenvectors of A, raising the eigenvalues to infinity. Due to my assumption that I have a transition matrix, we would have only eigenvalues equal to one (which will remain one), and eigenvalues between 0 and 1, which will become zero (inspired by this answer):

import numpy as np
from scipy import linalg    

# compute left eigenvalues and left eigenvectors
eigenvalues, leftEigenvectors = linalg.eig(transitionMatrix, right=False, left=True)

# for stationary distribution, eigenvalues and vectors are real (up to numerical precision)
eigenvalues = eigenvalues.real
leftEigenvectors = leftEigenvectors.real

# create a filter to collect the eigenvalues close to one
absoluteTolerance = 1e-10
mask = abs(eigenvalues - 1) < absoluteTolerance

# raise eigenvalues to the power of infinity
eigenvalues[mask] = 1
eigenvalues[~mask] = 0

D_raised = np.diag(eigenvalues)

A_raised = leftEigenvectors @ D_raised @ linalg.inv(leftEigenvectors)

Would this be the recommended approach, i.e., one that is both numerically stable and efficient?

Unis
  • 614
  • 1
  • 5
  • 17
  • Efficient exponentiation (matrix or what not) is a very common thing and can be done in log time – Julien Nov 27 '20 at 07:13
  • 1
    I think you'd want to ask [math.se] or [stats.se] instead of us. This isn't really an implementation question. It does seem to be sound with my limited knowledge of tensor theory though. – Daniel F Nov 27 '20 at 07:25
  • 1
    Only change I would make is `mask = np.isclose(eigenvalues, 1, atol = 1e-10)`, or even `eigenvalues = np.isclose(eigenvalues, 1, atol = 1e-10).astype(int)` – Daniel F Nov 27 '20 at 07:41
  • 1
    @DanielF Also adding [Computational Science](https://scicomp.stackexchange.com/) to the mix. – Lith Nov 27 '20 at 19:39
  • It depends on the size of the matrix. If you can afford the decomposition, use only the dominant left eigenvector and take the outer product with the constant vector of ones ( check row-column orientation here), i.e. the matrix has either the same rows or columns. Also, if you use right=true, you can get the inverse of the left matrix and it should be faster and more stable. Beware that the stationary distribution is not unique, if more than one EV is 1. If the matrix is very large, use the Power Method to get an arbitrarily good approximation of one EV with eigenvalue 1. – jhp Jan 31 '21 at 21:10

0 Answers0