使用scipy.sparse实现马尔可夫链稳态分布？

Question

使用scipy.sparse实现马尔可夫链稳态分布？

5

我有一个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)是单位矩阵），但(A-I)不会具有满秩，甚至整个系统可能是欠定的。此外，可能会生成具有负条目的特征向量，无法归一化为有效的概率分布。避免在p中出现负数将是不错的选择。

使用scipy.sparse.linalg.eigen.eigs不是一个解决方案：它不允许指定加性约束条件。（如果特征向量包含负条目，则规范化无法帮助。）此外，它与真实结果相差很大，并且有时无法收敛，性能比scipy.linalg.eig更差。（还使用了移位反演模式，它可以改善我想要的特征值类型的查找但不能改善它们的质量。如果不使用它，则会更加繁琐，因为我只感兴趣一个特定的特征值1。）
转换为密集矩阵并使用scipy.linalg.eig不是一个解决方案：除了负条目问题外，矩阵也太大了。
使用scipy.sparse.spsolve不是一个明显的解决方案：矩阵既可能不是正方形的（当将加性约束条件和特征向量条件组合在一起时），也可能不具有满秩（尝试以某种方式单独指定它们时），有时两者都不满足。

有没有一种很好的方法可以在Python中数值地获取表示为稀疏矩阵的Markov链的平稳状态？如果有一种获取详尽列表（并可能还包括接近平稳状态）的方法，那将是赞赏的，但这不是必需的。

- Anaphory

4个回答

1

这是解决可能存在不充分规定的矩阵方程，因此可以使用scipy.sparse.linalg.lsqr来完成。我不知道如何确保所有条目都为正数，但除此之外，它非常简单。

import scipy.sparse.linalg
states = A.shape[0]

# I assume that the rows of A sum to 1.
# Therefore, In order to use A as a left multiplication matrix,
# the transposition is necessary.
eigvalmat = (A - scipy.sparse.eye(states)).T
probability_distribution_constraint = scipy.ones((1, states))

lhs = scipy.sparse.vstack(
    (eigvalmat,
     probability_distribution_constraint))

B = numpy.zeros(states+1)
B[-1]=1

r = scipy.sparse.linalg.lsqr(lhs, B)
# r also contains metadata about the approximation process
p = r[0]

- Anaphory

1

针对静态解不唯一且可能不为非负数的观点进行说明。

这意味着您的马尔可夫链不是不可约的，您可以将问题分解为不可约的马尔可夫链。为此，您需要找到马尔可夫链的闭合通信类，这实际上是研究转移图的连通组件（维基百科建议使用一些线性算法来查找强连通组件）。此外，您可以放弃所有开放通信类，因为每个静态状态必须在这些状态上消失。

如果您有闭合通信类C_1，...，C_n，则您的问题有望分解为简单小块：每个闭合类C_i上的马尔可夫链现在是不可约的，因此受限制的转移矩阵M_i具有恰好一个特征值为0的特征向量，并且该特征向量仅具有正分量（请参见Perron-Frobenius定理）。因此，我们有恰好一个静态状态x_i。

现在，您完整的马尔可夫链的静态状态都是来自闭合类中x_i的线性组合。事实上，这些都是静态状态。

为了找到稳态 x_i，您可以依次应用 M_i，迭代将收敛于该状态（这也将保持您的归一化）。一般来说，很难确定收敛速度，但它为您提供了一种简单的方法来提高解决方案的准确性并验证解决方案。

- Helge Dietert

-1

使用幂迭代（例如）：http://www.google.com/search?q=power%20iteration%20markov%20chain

或者，您可以使用scipy.sparse.linalg.eig的shift-invert模式（即ARPACK）来查找接近1的特征值。 "指定"归一化不是必要的，因为您可以在之后规范化这些值。

- pv.

我使用了shift-invert模式，这可以提高我所需的特征值类型的查找效率，但并不能改善它们的质量。此外，仅仅对特征向量进行归一化是不起作用的，因为它们可能包含负数，需要通过找到特征向量之间的正确线性组合来解决——如果特征向量不精确，则会出现问题。 - Anaphory

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- pv. · Accepted Answer

可以在Google学术中找到几篇关于可能方法的文章摘要，以下是其中一篇：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()

编辑: 发现一个bug --- csgraph.connected_components 还会返回纯粹消退的组件，需要进行过滤。