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

Question

我有一个非常大的吸收马尔可夫链(扩展到问题大小 - 从10个状态到数百万个)非常稀疏(大多数状态只能对4或5个其他状态作出反应).

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

通常,我会通过计算来做到这一点(I - Q)^(-1),但我找不到一个实现稀疏矩阵逆算法的好库!我已经看过几篇论文,其中大部分都是博士级的工作.

我的大部分Google结果都指向了一些帖子,讨论在解决线性(或非线性)方程组时如何不应该使用矩阵逆...我觉得这没有特别有用.基本矩阵的计算是否类似于求解方程组,我根本不知道如何以另一个的形式表达一个？

所以,我提出两个具体问题:

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

要么

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

一个Python解决方案会很精彩(因为我的项目目前仍然是一个概念验证),但是如果我不得不用一些好的'Fortran或C弄脏,那不是问题.

编辑:我刚刚意识到矩阵A的逆B可以定义为AB = I,其中I是单位矩阵.这可能会允许我使用一些标准的稀疏矩阵解算器来计算逆...我得赶紧走了,可以随意完成我的思路,我刚开始想可能只需要一个非常基本的矩阵属性...

Answer 1

假设您想要做的是计算出吸收之前的预期步数，则维基百科上复制的“有限马尔可夫链”（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

和：

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

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

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