計算吸收馬爾可夫鏈基本矩陣的最佳迭代方法？

我有一個非常大的吸收馬爾可夫鏈。我想獲得這個鏈的基本矩陣來計算expected number of steps before absortion。從這個question我知道，這可通過以下公式來計算計算吸收馬爾可夫鏈基本矩陣的最佳迭代方法？

（I - Q）T = 1，它可以通過使用下面的Python代碼中獲得

：

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

然而，我想用計算PageRank時使用的某種類似於功率迭代方法的迭代方法來計算它。這個方法將允許我在mapreduce-like系統中計算在absortion之前的期望步數的近似值。

¿是否存在類似的情況？

來源

2017-07-18 Carlos Perez-Miguel

你看過像GMRES這樣的線性方程組的迭代求解器嗎？他們提供了一些錯誤估計來檢查你需要多少次迭代。 Python在['scipy.sparse.linalg']中提供了它們（https://docs.scipy.org/doc/scipy/reference/generated/scipy.sparse.linalg.gmres。html＃scipy.sparse.linalg.gmres） –

對不起，我誤解了你的問題！你看過基本矩陣的[Neumann series representation]（https://en.wikipedia.org/wiki/Absorbing_Markov_chain#Fundamental_matrix）嗎？可能僅使用有限數量的力量系列 –

中的加數來近似逆向感謝您的意見@TobiasRibizel。我認爲紐曼系列不是我正在尋找的，因爲它意味着將矩陣乘以幾次。但是，如果我正確理解了GMRES算法，則此解決方案可能以每個組件爲基礎並行化。 –

如果您有一個稀疏矩陣，請檢查scipy.spare.linalg.spsolve是否有效。沒有關於數值魯棒性的保證，但至少對於微不足道的例子來說，它比用密集矩陣求解要快得多。

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)

時機與稀疏解：

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

時機與密集的解決方案：

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

像托比亞斯提到的意見，我本來期望其它的求解器跑贏普通一個，而且它們可能適用於非常大的系統。對於這個玩具的例子，通用解決方案似乎工作得很好。

來源

2017-07-19 11:47:12

感謝您的回答安德魯，但是我正在尋找一種方法，可以並行計算解決方案的所有元素。遵循Tobias提出的GMRES方法，我發現了其他迭代方法，如[Jacobi方法]（https://en.wikipedia.org/wiki/Jacobi_method），可以滿足我的要求。 –

感謝@ tobias-ribizel的使用Neumann series的建議，我安排回答這個問題。如果我們的一部分從下面的公式：

$t=(I-Q)^-1 1$

使用紐曼系列：

$t=sum_0_inf(Q^k)1$

如果我們通過矢量我們可以單獨在操作乘以系列每學期矩陣的每一行Q並大致連續地與：

$t=sum_0_inf(Q*Q^k-1*1)$

這是Python代碼，我用它來計算這個

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

這是用放電的代碼。這段代碼預計Q是一個RDD，其中每一行都是一個元組（row_id，該矩陣行的權重字典）。

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

來源

2017-07-31 10:29:35

計算吸收馬爾可夫鏈基本矩陣的最佳迭代方法？

回答

相關問題