具有scipy.sparse的马尔可夫链平稳分布? [英] Markov chain stationary distributions with scipy.sparse?

问题描述

我有一个马尔可夫链，它是一个大的稀疏scipy矩阵A. (我已经以scipy.sparse.dok_matrix格式构造了矩阵，但是转换为其他矩阵或将其构造为csc_matrix很好.)

I have a Markov chain given as a large sparse scipy matrix A. (I've constructed the matrix in scipy.sparse.dok_matrix format, but converting to other ones or constructing it as csc_matrix are fine.)

我想知道此矩阵的任何平稳分布p，它是特征值1的特征向量.该特征向量中的所有条目都应为正并加1，以表示概率分布.

I'd like to know any stationary distribution p of this matrix, which is an eigenvector to the eigenvalue 1. All entries in this eigenvector should be positive and add up to 1, in order to represent a probability distribution.

这意味着我想要该系统的任何解决方案 (A-I) p = 0，p.sum()=1(其中I=scipy.sparse.eye(*A.shape)是恒等矩阵)，但是(A-I)不能完全排名，甚至整个系统都可能无法确定.另外，可以生成带有负项的特征向量，不能将其归一化为有效的概率分布.防止在p中出现否定条目会很不错.

This means I want any solution for the system (A-I) p = 0, p.sum()=1 (where I=scipy.sparse.eye(*A.shape) is the idententy matrix), but (A-I) will not be of full rank, and even the whole system may be under-determined. In addition, eigenvectors with negative entries can be generated, which cannot be normalized to be valid probability distributions. Preventing negative entries in p would be nice.

• 使用scipy.sparse.linalg.eigen.eigs不是解决方案: 它不允许指定加性约束. (如果特征向量包含否定项，归一化将无济于事.)此外，它与真实结果有很大出入，有时会出现收敛问题，其表现要比scipy.linalg.eig差. (此外，我使用了移位反转模式，该模式可以改善查找所需特征值的类型的能力，但不能提高其质量.如果不使用它，它的作用就更大了，因为我只对一个特定的特征值感兴趣，.)

• Using scipy.sparse.linalg.eigen.eigs is not a solution: It does not permit specifying the additive constraint. (Normalizing does not help if the eigenvectors contain negative entries.) Furthermore, it deviates quite a bit from the true results, and sometimes has problems converging, behaving worse than scipy.linalg.eig. (Also, I used shift-invert mode, which improves finding the types of eigenvalues I want, but not their quality. If I don't use it, it's even more overkill, because I am only interested in one particular eigenvalue, 1.)

转换为密集矩阵并使用scipy.linalg.eig并不是解决方案:除了负输入问题外，矩阵太大.

Converting to dense matrices and using scipy.linalg.eig is not a solution: In addition to the negative entry problem, the matrix is too large.

使用scipy.sparse.spsolve并不是一个显而易见的解决方案: 矩阵不是正方形(当将加性约束和特征向量条件组合在一起时)，或者不是满秩的(当试图以某种方式分别指定它们时)，有时也不是.

Using scipy.sparse.spsolve is not an obvious solution: The matrix is either not square (when combining the additive constraint and the eigenvector condition) or not of full rank (when trying to specify them separately in some way), sometimes neither.

是否存在一种很好的方法来使用python数字地获得以稀疏矩阵形式给出的马尔可夫链的平稳状态? 如果有一种方法可以获取详尽的列表(可能还包括几乎静止的状态)，则表示赞赏，但这不是必需的.

Is there a good way to numerically get a stationary state of a Markov chain given as sparse matrix using python? If there is a way to get an exhaustive list (and possibly also nearly stationary states), that's appreciated, but not necessary.

推荐答案

Google学者可以找到几篇文章，其中概要介绍了可能的方法，以下是其中之一: http://www.ima.umn.edu/preprints/pp1992/932.pdf

Several articles with summaries of possible approaches can be found with Google scholar, here's one: http://www.ima.umn.edu/preprints/pp1992/932.pdf

下面要做的是上述@Helge Dietert提出的关于首先拆分为强连接的组件的建议以及上面链接的论文中的方法4的组合.

What's done below is a combination of the suggestion by @Helge Dietert above on splitting to strongly connected components first, and approach #4 in the paper linked above.

import numpy as np
import time

# NB. Scipy >= 0.14.0 probably required
import scipy
from scipy.sparse.linalg import gmres, spsolve
from scipy.sparse import csgraph
from scipy import sparse 


def markov_stationary_components(P, tol=1e-12):
    """
    Split the chain first to connected components, and solve the
    stationary state for the smallest one
    """
    n = P.shape[0]

    # 0. Drop zero edges
    P = P.tocsr()
    P.eliminate_zeros()

    # 1. Separate to connected components
    n_components, labels = csgraph.connected_components(P, directed=True, connection='strong')

    # The labels also contain decaying components that need to be skipped
    index_sets = []
    for j in range(n_components):
        indices = np.flatnonzero(labels == j)
        other_indices = np.flatnonzero(labels != j)

        Px = P[indices,:][:,other_indices]
        if Px.max() == 0:
            index_sets.append(indices)
    n_components = len(index_sets)

    # 2. Pick the smallest one
    sizes = [indices.size for indices in index_sets]
    min_j = np.argmin(sizes)
    indices = index_sets[min_j]

    print("Solving for component {0}/{1} of size {2}".format(min_j, n_components, indices.size))

    # 3. Solve stationary state for it
    p = np.zeros(n)
    if indices.size == 1:
        # Simple case
        p[indices] = 1
    else:
        p[indices] = markov_stationary_one(P[indices,:][:,indices], tol=tol)

    return p


def markov_stationary_one(P, tol=1e-12, direct=False):
    """
    Solve stationary state of Markov chain by replacing the first
    equation by the normalization condition.
    """
    if P.shape == (1, 1):
        return np.array([1.0])

    n = P.shape[0]
    dP = P - sparse.eye(n)
    A = sparse.vstack([np.ones(n), dP.T[1:,:]])
    rhs = np.zeros((n,))
    rhs[0] = 1

    if direct:
        # Requires that the solution is unique
        return spsolve(A, rhs)
    else:
        # GMRES does not care whether the solution is unique or not, it
        # will pick the first one it finds in the Krylov subspace
        p, info = gmres(A, rhs, tol=tol)
        if info != 0:
            raise RuntimeError("gmres didn't converge")
        return p


def main():
    # Random transition matrix (connected)
    n = 100000
    np.random.seed(1234)
    P = sparse.rand(n, n, 1e-3) + sparse.eye(n)
    P = P + sparse.diags([1, 1], [-1, 1], shape=P.shape)

    # Disconnect several components
    P = P.tolil()
    P[:1000,1000:] = 0
    P[1000:,:1000] = 0

    P[10000:11000,:10000] = 0
    P[10000:11000,11000:] = 0
    P[:10000,10000:11000] = 0
    P[11000:,10000:11000] = 0

    # Normalize
    P = P.tocsr()
    P = P.multiply(sparse.csr_matrix(1/P.sum(1).A))

    print("*** Case 1")
    doit(P)

    print("*** Case 2")
    P = sparse.csr_matrix(np.array([[1.0, 0.0, 0.0, 0.0],
                                    [0.5, 0.5, 0.0, 0.0],
                                    [0.0, 0.0, 0.5, 0.5],
                                    [0.0, 0.0, 0.5, 0.5]]))
    doit(P)

def doit(P):
    assert isinstance(P, sparse.csr_matrix)
    assert np.isfinite(P.data).all()

    print("Construction finished!")

    def check_solution(method):
        print("\n\n-- {0}".format(method.__name__))
        start = time.time()
        p = method(P)
        print("time: {0}".format(time.time() - start))
        print("error: {0}".format(np.linalg.norm(P.T.dot(p) - p)))
        print("min(p)/max(p): {0}, {1}".format(p.min(), p.max()))
        print("sum(p): {0}".format(p.sum()))

    check_solution(markov_stationary_components)


if __name__ == "__main__":
    main()

编辑:发现了一个错误--- csgraph.connected_components还返回纯衰减的分量，需要将其滤除.

EDIT: spotted a bug --- csgraph.connected_components returns also purely decaying components, which need to be filtered out.

这篇关于具有scipy.sparse的马尔可夫链平稳分布?的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！