12

我有一个非常大的吸收马尔可夫链（扩展到问题规模——从 10 个状态到数百万个）非常稀疏（大多数状态只能对 4 或 5 个其他状态做出反应）。

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

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

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

所以，我提出两个具体问题：

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

或者

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

Python 解决方案会很棒（因为我的项目目前仍然是一个概念验证），但是如果我不得不用一些好的 Fortran 或 C 来弄脏我的手，那不是问题。

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

4

2 回答 2

4

假设您要做的是计算出吸收之前的预期步数，那么维基百科上转载的“有限马尔可夫链”（Kemeny 和 Snell）的方程式是：

t=N1

或扩展基本矩阵

t=(智商)^-1 1

重新排列：

(智商) t = 1

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

A x = b

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

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

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

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

于 2014-07-15T13:38:02.380 回答
3

您得到建议不要使用矩阵求逆来求解方程的原因是因为数值稳定性。当您的矩阵具有零或接近零的特征值时，您会因缺乏逆（如果为零）或数值稳定性（如果接近零）而遇到问题。那么，解决这个问题的方法是使用一种不需要逆存在的算法。解决方法是使用高斯消元法。这并没有提供完整的逆，而是让您进入行梯形形式，即上三角形形式的推广。如果矩阵是可逆的，则结果矩阵的最后一行包含逆矩阵。因此，只需安排您消除的最后一行是您想要的行。

我会让你明白为什么I-Q总是可逆的。

于 2012-11-08T15:34:01.150 回答