Meilleure méthode itérative pour 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 la nombre d'étapes prévues avant l'absorption . De ce question Je sais que cela peut être calculé par l'équation suivante

(I - Q)t=1

qui peut être obtenu en utilisant le code 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 de l'itération de puissance. utilisé pour calculer le PageRank . Cette méthode me permettrait de calculer une approximation du nombre attendu d'étapes avant absorption dans un système de type mapreduce.

¿Est-ce que quelque chose de similaire existe ?

Answer 1

Si vous avez une matrice éparse, vérifiez si scipy.spare.linalg.spsolve travaux. Aucune garantie quant à la robustesse numérique, mais au moins pour des exemples triviaux, c'est significativement 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)

Timings pour une solution éparse :

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

Délais pour la solution dense :

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

Comme Tobias le mentionne dans les commentaires, je m'attendais à ce que d'autres solveurs soient plus performants que le générique, et c'est le cas pour les très grands systèmes. Pour ce petit exemple, la solution générique semble fonctionner suffisamment bien.

Answer 2

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

En utilisant la série de Neumann :

Si on multiplie chaque terme de la série par le vecteur 1 nous pourrions opérer séparément sur chaque ligne de la matrice Q et l'approximer successivement avec :

Voici le code python que j'utilise pour faire ce calcul :

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 voici le code utilisant Spark. Ce code s'attend à ce que Q soit un RDD où chaque ligne est un tuple (row_id, dict des poids pour cette ligne 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