我有一个马尔可夫链作为一个大的稀疏矩阵scipy
矩阵{scipy.sparse.dok_matrix
格式的矩阵,但是转换成其他格式或者将其构造为csc_matrix
都可以。)
我想知道这个矩阵的任何平稳分布p
,它是特征值1
的特征向量。这个特征向量中的所有条目都应该是正的,加起来等于1,以表示概率分布。在
这意味着我需要系统的任何解决方案
(A-I) p = 0
,p.sum()=1
(其中I=scipy.sparse.eye(*A.shape)
是恒等矩阵),但是{p
中的负条目会很好。在
使用scipy.sparse.linalg.eigen.eigs
不是解决方案:
它不允许指定加法约束。(如果特征向量包含负项,则规范化没有帮助。)此外,它与实际结果有相当大的偏差,有时会出现收敛问题,表现得比scipy.linalg.eig
更糟糕。(另外,我使用了移位-反转模式,它可以改进我想要的特征值类型的查找,但不能改进其质量。如果我不使用它,那就更过分了,因为我只对一个特定的特征值感兴趣,1
。)
转换为稠密矩阵并使用scipy.linalg.eig
并不是一个解决方案:除了负输入问题外,矩阵太大。
使用scipy.sparse.spsolve
不是一个明显的解决方案:
矩阵要么不是正方形(当结合加性约束和特征向量条件时),要么不是满秩(当试图以某种方式分别指定它们时),有时两者都不是。
有没有一种好的方法可以用python数值求解一个以稀疏矩阵形式给出的Markov链的平稳状态? 如果有一种方法可以得到一个详尽的列表(也可能是几乎静止的状态),这是值得赞赏的,但不是必要的。在
解决了固定解不是唯一的,解可能不是非负的这一点。在
这意味着你的马尔可夫链不是不可约的,你可以把问题分解成不可约的马尔可夫链。为此,您需要找到Markov链的封闭通信类,这本质上是对转移图中的连接组件的研究(Wikipedia建议使用一些线性算法来查找强连接组件)。此外,您可以通过所有开放的通信类,因为每个静止状态必须消失在这些类上。在
如果你有你的封闭通讯类Cˉ1,…,你的问题有希望被分成几个小的简单部分:每个封闭类Cˉi上的马尔可夫链现在是不可约的,因此,受限转移矩阵M峈i只有一个特征值为0的特征向量,且该特征向量只有正分量(参见Perron-Frobenius定理)。因此我们只有一个稳态x
你的整个马尔可夫链的平稳态现在都是封闭类的x_i的线性组合。事实上,这些都是静止态。在
为了找到稳定态x,你可以连续地应用Méi,迭代会收敛到这个状态(这也会保持你的正规化)。一般来说,很难判断收敛速度,但它提供了一种简单的方法来提高解的精度和验证解。在
这是在求解一个可能未指定的矩阵方程,因此可以用
scipy.sparse.linalg.lsqr
来完成。我不知道如何确保所有条目都是正的,但除此之外,它非常简单。在谷歌学者(Google scholar)提供了几篇文章,总结了可能的方法,以下是其中一篇: http://www.ima.umn.edu/preprints/pp1992/932.pdf
下面所做的是上面@Helge Dietert关于先拆分为强连接组件的建议,以及上面链接的论文中的方法4。在
编辑:发现一个错误--csgraph.connected_组件还返回需要过滤掉的纯衰减组件。在
相关问题 更多 >
编程相关推荐