计算吸收马尔可夫链基本矩阵的最佳方法? [英] Best way to calculate the fundamental matrix of an absorbing Markov Chain?
问题描述
我有一个非常大的吸收马尔可夫链(规模到问题大小-从10个州到数百万个)，该链非常稀疏(大多数州只能对其他4个或5个州做出反应).
I have a very large absorbing Markov chain (scales to problem size -- from 10 states to millions) that is very sparse (most states can react to only 4 or 5 other states).
我需要计算该链的基本矩阵的一行(给定一个起始状态的每个状态的平均频率).
I need to calculate one row of the fundamental matrix of this chain (the average frequency of each state given one starting state).
通常，我会通过计算(I - Q)^(-1)来做到这一点，但是我一直找不到能够实现稀疏矩阵逆算法的好的库！我看过几篇论文，其中大多数是P.h.D.水平的工作.
Normally, I'd do this by calculating (I - Q)^(-1), but I haven't been able to find a good library that implements a sparse matrix inverse algorithm! I've seen a few papers on it, most of them P.h.D. level work.
我在Google上的大部分搜索结果都指向一些帖子，内容是关于在求解线性(或非线性)方程组时不应该使用矩阵逆的问题……我认为这没有什么特别的帮助.基本矩阵的计算是否类似于求解方程组，而我根本不知道如何以另一种形式来表达?
Most of my Google results point me to posts talking about how one shouldn't use a matrix inverse when solving linear (or non-linear) systems of equations... I don't find that particularly helpful. Is the calculation of the fundamental matrix similar to solving a system of equations, and I simply don't know how to express one in the form of the other?
所以，我提出两个具体问题:
So, I pose two specific questions:
计算稀疏矩阵逆的一行(或所有行)的最佳方法是什么?
OR
计算大吸收性马尔可夫链的基本矩阵行的最佳方法是什么?
Python解决方案将是很棒的(因为我的项目目前仍是概念验证)，但是如果我不得不动手使用Fortran或C语言，那不是问题.
A Python solution would be wonderful (as my project is still currently a proof-of-concept), but if I have to get my hands dirty with some good ol' Fortran or C, that's not a problem.
我刚刚意识到矩阵A的逆B可以定义为AB = I，其中I是单位矩阵.这可能使我可以使用一些标准的稀疏矩阵求解器来计算逆数...我必须要跑了，所以请随意完成我的思路，我开始认为这可能只需要一个真正的基本矩阵财产...
I just realized that the inverse B of matrix A can be defined as AB=I, where I is the identity matrix. That may allow me to use some standard sparse matrix solvers to calculate the inverse... I've got to run off, so feel free to complete my train of thought, which I'm starting to think might only require a really elementary matrix property...
推荐答案
Assuming that what you're trying to do is work out is the expected number of steps before absorbtion, the equation from "Finite Markov Chains" (Kemeny and Snell), which is reproduced on Wikipedia is:
或扩展基本矩阵
重新排列:
使用函数求解线性方程组的标准格式
Which is in the standard format for using functions for solving systems of linear equations
将其付诸实践以证明性能上的差异(即使是比您所描述的系统小得多的系统).
Putting this into practice to demonstrate the difference in performance (even for much smaller systems than those you're describing).
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
介绍了两种用于计算预期步骤数的替代方法.
Presenting the two alternative approaches for calculating the number of expected steps.
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)
选择一个足以说明计算基本矩阵时出现的问题类型的示例:
Picking an example that's big enough to demonstrate the types of problems that occur when calculating the fundamental matrix:
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
需要进一步的实验来测试稀疏矩阵的性能影响，但是很明显，计算逆函数比您期望的要慢得多.
Further experimentation is required to test the performance implications for sparse matrices, but it's clear that calculating the inverse is much much slower than what you might expect.
与此处介绍的方法类似的方法也可以用于步数的变化
A similar approach to the one presented here can also be used for the variance of the number of steps
这篇关于计算吸收马尔可夫链基本矩阵的最佳方法?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！
