计算吸收马尔可夫链的基本矩阵的最佳方法？

我有一个非常大的吸收马尔可夫链（根据问题大小进行缩放——从 10 个状态到数百万个），它非常稀疏（大多数状态只能对 4 或 5 个其他状态做出反应）。

我需要计算该链的基本矩阵的一行（给定一个起始状态的每个状态的平均频率）。

通常，我会通过计算来做到这一点(I - Q)^(-1)，但是我一直没能找到一个好的库来实现稀疏矩阵逆算法！我看过几篇关于它的论文，其中大部分是博士论文。水平的工作。

我的大部分谷歌结果都指向我一些帖子，讨论在求解线性（或非线性）方程组时不应该使用矩阵逆......我觉得这不是特别有帮助。基本矩阵的计算是否类似于求解方程组，而我根本不知道如何用另一个形式来表达一个？

因此，我提出两个具体问题：

计算稀疏矩阵的逆矩阵的一行（或所有行）的最佳方法是什么？

OR

计算大型吸收马尔可夫链的基本矩阵的一行的最佳方法是什么？

Python 解决方案会很棒（因为我的项目目前仍处于概念验证阶段），但如果我必须亲自动手使用一些好的 Fortran 或 C，那不是问题。

编辑：我刚刚意识到矩阵 A 的逆 B 可以定义为 AB=I，其中 I 是单位矩阵。这可能允许我使用一些标准的稀疏矩阵求解器来计算逆矩阵......我必须逃跑，所以请随意完成我的思路，我开始认为这可能只需要一个真正的基本矩阵财产...

假设你想做的是锻炼吸收前的预期步数 http://en.wikipedia.org/wiki/Absorbing_Markov_chain#Expected_number_of_steps，维基百科上转载的“有限马尔可夫链”（Kemeny 和 Snell）的方程为：

或者扩展基本矩阵

重新排列：

这是使用函数求解线性方程组的标准格式

将其付诸实践以证明性能差异（即使对于比您所描述的系统小得多的系统）。

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

提出计算预期步数的两种替代方法。

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)

选择一个足够大的示例来演示计算基本矩阵时出现的问题类型：

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

产生以下时序：

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

And:

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

需要进一步的实验来测试稀疏矩阵的性能影响，但很明显，计算逆矩阵比您预期的要慢得多。

与此处介绍的方法类似的方法也可用于步骤数的方差