马尔可夫链固定分布与scipy.sparse？

Question

我有一个Markov链给出一个大的稀疏scipy矩阵A.(我已经以scipy.sparse.dok_matrix格式构建了矩阵,但转换为其他矩阵或构造它csc_matrix很好.)

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

这意味着我想要系统的任何解决方案 (A-I) p = 0,p.sum()=1(在哪里I=scipy.sparse.eye(*A.shape)是idententy矩阵),但(A-I)不是完全排名,甚至整个系统可能都是不确定的.另外,可以生成具有否定条目的特征向量,其不能被归一化为有效概率分布.防止负面输入p会很好.

• 使用scipy.sparse.linalg.eigen.eigs不是解决方案:它不允许指定附加约束.(如果特征向量包含负数条目,则归一化没有帮助.)此外,它与真实结果有很大差异,有时会出现收敛问题,表现比较差scipy.linalg.eig.(另外,我使用了shift-invert模式,它改进了我想要的特征值的类型,但不是它们的质量.如果我不使用它,它甚至更具有杀伤力,因为我只对一个特定的特征值感兴趣,1.)

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

• 使用scipy.sparse.spsolve不是一个明显的解决方案:矩阵要么不是正方形(当组合加性约束和特征向量条件时),要么不是满秩(当试图以某种方式单独指定它们时),有时也不是.

有没有一种很好的方法可以使用python以数字形式获得马尔可夫链的静态状态作为稀疏矩阵？如果有办法获得详尽的清单(也可能是几乎静止的状态),那是值得赞赏的,但并非必要.

Answer 1

谷歌学者可以找到几篇关于可能方法摘要的文章,其中一篇:http: //www.ima.umn.edu/preprints/pp1992/932.pdf

下面的内容是@Helge Dietert上面提出的关于首先拆分为强连接组件的建议,以及上面链接的文章中的#4方法的组合.

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也返回纯粹腐烂的组件,需要过滤掉.