Il modo migliore per calcolare la matrice fondamentale di una catena Markov assorbente?

Question

Ho una catena di Markov assorbente molto grande (scala alla dimensione del problema - da 10 stati a milioni) che è molto sparsa (la maggior parte degli stati può reagire solo a 4 o 5 altri stati).

Ho bisogno di calcolare una riga della matrice fondamentale di questa catena (la frequenza media di ogni stato dato uno stato iniziale).

Normalmente, lo farei calcolando (I - Q)^(-1), ma non sono stato in grado di trovare una buona libreria che implementa un algoritmo inverso di matrice sparsa! Ho visto alcuni documenti su di esso, la maggior parte di loro P.h.D. livello di lavoro.

La maggior parte dei risultati di Google mi indirizza a post che parlano di come non si dovrebbe utilizzare una matrice inversa quando si risolvono i sistemi di equazioni lineari (o non lineari) ... Non lo trovo particolarmente utile. Il calcolo della matrice fondamentale è simile alla soluzione di un sistema di equazioni e semplicemente non so come esprimerne uno nella forma dell'altro?

Quindi, mi pongo due domande specifiche:

Qual è il modo migliore per calcolare una riga (o tutte le righe) della inversa di una matrice sparsa?

O

Qual è il modo migliore per calcolare una riga della matrice fondamentale di una grande catena di Markov assorbimento?

Una soluzione Python sarebbe meravigliosa (come il mio progetto è ancora attualmente una dimostrazione di concetto), ma se devo sporcarmi le mani con un buon vecchio 'Fortran o C, non è un problema.

Modifica: ho appena realizzato che l'inversa B della matrice A può essere definita come AB = I, dove I è la matrice di identità. Questo potrebbe consentirmi di usare alcuni solutori di matrice sparse standard per calcolare l'inverso ... Devo scappare, quindi sentiti libero di completare il mio pensiero, che sto iniziando a pensare potrebbe richiedere solo una matrice davvero elementare proprietà ...

Answer 1

Partendo dal presupposto che ciò che si sta cercando di fare è lavorare fuori è la expected number of steps before absorbtion, l'equazione di "Catene di Markov Finite" (Kemeny e Snell), che è riprodotto su Wikipedia è:

o ampliare la matrice fondamentale

Riorganizzare:

che è nel formato standard per l'utilizzo di funzioni per i sistemi di equazioni lineari

Mettere in pratica la risoluzione per dimostrare la differenza di prestazioni (anche per i sistemi di molto inferiori a quelle che si' descrivendo).

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 

Presentando i due approcci alternativi per calcolare il numero di passi previsti.

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) 

Picking un esempio che è abbastanza grande per dimostrare i tipi di problemi che si verificano quando si calcola la matrice fondamentale:

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

produce i seguenti orari:

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

E:

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

Sono necessari ulteriori esperimenti per testare le implicazioni di prestazione per le matrici sparse, ma è chiaro che il calcolo dell'inverso è molto più lento di quanto ci si potrebbe aspettare.

Un approccio simile a quello qui presentato può essere utilizzato anche per la varianza del numero di passi

Answer 2

La ragione per cui stai ricevendo il consiglio di non utilizzare gli inversori di matrice per risolvere equazioni è a causa della stabilità numerica. Quando sei matrice ha autovalori che sono zero o quasi zero, hai problemi o dalla mancanza di una inversa (se zero) o stabilità numerica (se vicino allo zero). Il modo di affrontare il problema, quindi, consiste nell'utilizzare un algoritmo che non richiede che esista un'inversione. La soluzione è utilizzare Gaussian elimination. Ciò non fornisce un'inversione completa, ma piuttosto ti porta alla forma a fila di righe, una generalizzazione della forma triangolare superiore. Se la matrice è invertibile, l'ultima riga della matrice dei risultati contiene una riga dell'inverso. Quindi ordina che l'ultima riga che elimini sia la riga che desideri.

Lo lascerò a voi per capire perché I-Q è sempre invertibile.