Le meilleur moyen itératif de calculer la matrice fondamentale d'une chaîne de Markov absorbante?

Question

J'ai une très grande chaîne de Markov absorbante. Je veux obtenir la matrice fondamentale de cette chaîne pour calculer le expected number of steps before absortion. A partir de ce question je sais que cela peut être calculé par l'équation

(I - Q) t = 1

qui peut être obtenu en utilisant le code de python suivant:

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

Cependant, Je voudrais le calculer en utilisant une sorte de méthode itérative similaire à la méthode d'itération de puissance utilisé pour calculer le PageRank. Cette méthode me permettrait de calculer une approximation du nombre attendu d'étapes avant l'absorption dans un système de type mapreduce.

¿Est-ce que quelque chose de similaire existe?

Answer 1

Si vous avez une matrice clairsemée, vérifiez si scipy.spare.linalg.spsolve fonctionne. Aucune garantie sur la robustesse numérique, mais au moins pour les exemples triviaux, il est nettement plus rapide que de résoudre avec des matrices denses.

import networkx as nx 
import numpy as np 
import scipy.sparse as sp 
import scipy.sparse.linalg as spla 

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/np.sum(m, axis=1) 
    return m 

A = sp.csr_matrix(example(2000)[:-1,:-1]) 
Ad = np.array(A.todense()) 

def sp_solve(Q): 
    I = sp.identity(Q.shape[0], format='csr') 
    o = np.ones(Q.shape[0]) 
    return spla.spsolve(I-Q, o) 

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

synchronisations pour la solution clairsemée:

%timeit sparse_solve(A) 
1000 loops, best of 3: 1.08 ms per loop 

Timings pour la solution dense:

%timeit dense_solve(Ad) 
1 loops, best of 3: 216 ms per loop 

Comme Tobias mentionne dans les commentaires, je me attendais d'autres solveurs de surperformer le générique, et ils peuvent pour de très grands systèmes. Pour cet exemple de jouet, la résolution générique semble assez bien fonctionner.

Answer 2

Je suis arrivé à cette réponse grâce à la suggestion de @ tobias-ribizel d'utiliser le Neumann series. Si nous nous séparons de l'équation suivante:

Utilisation de la série Neumann:

Si l'on multiplie chaque terme de la série par le vecteur nous pourrions fonctionner séparément sur chaque rangée de la matrice Q et d'approcher successivement avec:

Ce code python J'utilise pour calculer ceci:

def expected_steps_iterative(Q, n=10): 
    N = Q.shape[0] 
    acc = np.ones(N) 
    r_k_1 = np.ones(N) 
    for k in range(1, n): 
     r_k = np.zeros(N) 
     for i in range(N): 
      for j in range(N): 
       r_k[i] += r_k_1[j] * Q[i, j] 
     if np.allclose(acc, acc+r_k, rtol=1e-8): 
      acc += r_k 
      break 
     acc += r_k 
     r_k_1 = r_k 
    return acc 

Et ceci est le code en utilisant Spark. Ce code s'attend à ce que Q soit un RDD où chaque ligne est un tuple (row_id, dict de poids pour cette rangée de la matrice).

def expected_steps_spark(sc, Q, n=10): 
    def dict2np(d, sz): 
     vec = np.zeros(sz) 
     for k, v in d.iteritems(): 
      vec[k] = v 
     return vec 
    sz = Q.count() 
    acc = np.ones(sz) 
    x = {i:1.0 for i in range(sz)} 
    for k in range(1, n): 
     bc_x = sc.broadcast(x) 
     x_old = x 
     x = Q.map(lambda (u, ol): (u, reduce(lambda s, j: s + bc_x.value[j]*ol[j], ol, 0.0))) 
     x = x.collectAsMap() 
     v_old = dict2np(x_old, sz) 
     v = dict2np(x, sz) 
     acc += v 
     if np.allclose(v, v_old, rtol=1e-8): 
      break 
    return acc