如何计算吸收马尔可夫链的基本矩阵？

Question

我有一个非常大的吸收马尔可夫链（从10个状态到数百万个状态）非常稀疏（大多数状态只能与4或5个其他状态相互作用）。

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

通常，我会通过计算(I-Q)^(-1)来完成这个任务，但我找不到一个实现稀疏矩阵求逆算法的好库！我看到了一些论文，其中大部分都是博士级别的工作。

我的大多数Google结果都指向了一些关于解线性（或非线性）方程组时为什么不应该使用矩阵求逆的帖子……我觉得那些帖子并没有什么帮助。计算基础矩阵是否类似于解方程组，我只是不知道如何用一个形式表达另一个？

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

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

或

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

如果可以提供Python解决方案，那将会非常棒（因为我的项目目前仍然是一个概念验证），但如果必须使用一些Fortran或C等编程语言，那也不是问题。

编辑：我刚刚意识到矩阵A的逆B可以定义为AB=I，其中I是单位矩阵。这可能允许我使用一些标准的稀疏矩阵求解器来计算逆矩阵……我得走了，所以请随意完成我的思路，我开始觉得这可能只需要一些非常基本的矩阵属性……

- Ryan Marcus

---

如果您需要Python解决方案，请标记为“python”。还有其他可能更或更少有用的Stack Exchange。 - Jared Farrish

我正在处理一些关于 PGM 的东西，想知道是否有一般计算的方法 - 不过对于稀疏矩阵还没有头绪，祝好运！ - argentage

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

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

可以使用类似于此处提出的方法来计算步数方差。

- Andrew Walker

Answer 2

你之所以听到不要使用矩阵求逆来解方程的建议，是因为数值稳定性。当你的矩阵具有特征值为零或接近零时，如果该特征值为零，则没有逆矩阵（如果接近零，则会出现数值稳定性问题）。因此，解决问题的方法是使用不需要存在逆矩阵的算法。解决方案是使用高斯消元法。这并不能提供一个完整的逆矩阵，而是使矩阵变换成行简化阶梯形式，它是上三角形式的推广。如果矩阵可逆，则结果矩阵的最后一行包含逆矩阵中的一行。因此，只需将要消元的最后一行排列成所需的行即可。

至于为什么I-Q总是可逆，我将留给你自己去理解。

- eh9