¿Cuál es la mejor forma de calcular la matriz fundamental de una cadena de Markov absorbente?

Question

Tengo una cadena de Markov absorbente muy grande (escalas al tamaño del problema - de 10 estados a millones) que es muy escasa (la mayoría de los estados pueden reaccionar a solo 4 o 5 estados más).

Necesito calcular una fila de la matriz fundamental de esta cadena (la frecuencia promedio de cada estado dado un estado inicial).

Normalmente, haria esto calculando (I - Q)^(-1), pero no he podido encontrar una buena biblioteca que implemente un algoritmo de matriz escasamente inversa. He visto algunos artículos, la mayoría de ellos P.h.D. trabajo de nivel.

La mayoría de mis resultados de Google me indican que no se debe usar una matriz inversa cuando se resuelven sistemas de ecuaciones lineales (o no lineales) ... No me parece especialmente útil. ¿El cálculo de la matriz fundamental es similar a la resolución de un sistema de ecuaciones, y simplemente no sé cómo expresarlo en la forma del otro?

Por lo tanto, planteo dos preguntas específicas:

Cuál es la mejor manera de calcular una fila (o todas las filas) de la inversa de una matriz dispersa?

O

Cuál es la mejor manera de calcular una fila de la matriz fundamental de una gran cadena de Markov de absorción?

Una solución de Python sería maravillosa (ya que mi proyecto todavía es actualmente una prueba de concepto), pero si tengo que ensuciarme las manos con Fortran o C, eso no es un problema.

Editar: Me acabo de dar cuenta de que el inverso B de la matriz A se puede definir como AB = I, donde I es la matriz de identidad. Eso me puede permitir usar algunos solucionadores de matriz dispersos estándar para calcular el inverso ... Tengo que salir corriendo, así que siéntete libre de completar mi tren de pensamiento, que estoy empezando a pensar que solo podría requerir una matriz realmente elemental propiedad ...

Answer 1

Suponiendo que lo que estamos tratando de hacer es hacer ejercicio es la expected number of steps before absorbtion, la ecuación de "Cadenas de Markov finita" (Kemeny y Snell), que se reproduce en la Wikipedia es:

o la ampliación de la matriz fundamental

Reordenando:

que está en el formato estándar para el uso de funciones para la resolución de sistemas de ecuaciones lineales

llevarlo a la práctica para demostrar la diferencia en el rendimiento (incluso para sistemas mucho más pequeños que los que' re describiendo).

import networkx as nx 
import numpy 

def example(n): 
    """Generate a very simple transition matrix from a directed graph 
    """ 
    g = nx.DiGraph() 
    for i in xrange(n-1): 
     g.add_edge(i+1, i) 
     g.add_edge(i, i+1) 
    g.add_edge(n-1, n) 
    g.add_edge(n, n) 
    m = nx.to_numpy_matrix(g) 
    # normalize rows to ensure m is a valid right stochastic matrix 
    m = m/numpy.sum(m, axis=1) 
    return m 

Presentamos los dos enfoques alternativos para calcular el número de pasos esperados.

def expected_steps_fundamental(Q): 
    I = numpy.identity(Q.shape[0]) 
    N = numpy.linalg.inv(I - Q) 
    o = numpy.ones(Q.shape[0]) 
    numpy.dot(N,o) 

def expected_steps_fast(Q): 
    I = numpy.identity(Q.shape[0]) 
    o = numpy.ones(Q.shape[0]) 
    numpy.linalg.solve(I-Q, o) 

Recogiendo un ejemplo que es lo suficientemente grande como para demostrar los tipos de problemas que se producen en el cálculo de la matriz fundamental:

P = example(2000) 
# drop the absorbing state 
Q = P[:-1,:-1] 

produce los siguientes tiempos:

%timeit expected_steps_fundamental(Q) 
1 loops, best of 3: 7.27 s per loop 

Y:

%timeit expected_steps_fast(Q) 
10 loops, best of 3: 83.6 ms per loop 

Se requiere más experimentación para probar las implicaciones de rendimiento para las matrices dispersas, pero está claro que el cálculo de la inversa es mucho más lento de lo que cabría esperar.

un enfoque similar al que se presenta aquí también se puede utilizar para la varianza del número de pasos

Answer 2

La razón por la que recibe el consejo de no utilizar inversos de matriz para resolver ecuaciones es debido a la estabilidad numérica. Cuando su matriz tiene valores propios que son cero o cerca de cero, tiene problemas ya sea por falta de una inversa (si es cero) o de estabilidad numérica (si es cercana a cero). La forma de abordar el problema, entonces, es usar un algoritmo que no requiera que exista una inversa. La solución es usar Gaussian elimination. Esto no proporciona una inversión completa, sino que te lleva a la forma escalonada, una generalización de la forma triangular superior. Si la matriz es invertible, entonces la última fila de la matriz de resultados contiene una fila de la inversa. Así que simplemente arregle que la última fila que elimine sea la fila que desea.

Te dejo entender por qué I-Q siempre es invertible.