scipy.sparse馬爾可夫鏈平穩分佈？

Question

我有一個馬爾可夫鏈作爲一個大的稀疏scipy矩陣A給出。 （I已經構造在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回報也完全腐爛的成分，這就需要被過濾掉。

Answer 2

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

或者，你可以使用scipy.sparse.linalg.eig的移反轉模式（這是ARPACK）尋找接近1特徵值「指定「規範化不是必要的，因爲你可以在之後對數值進行規範化。

Answer 3

這是解決一個可能未指定的矩陣方程，因此可以用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] 

Answer 4

解決固定解不是唯一的問題，並且解決方案可能不是非負的。

這意味着您的馬爾可夫鏈不是不可約的，您可以將問題分解爲不可約馬爾可夫鏈。爲此，您需要查找馬爾可夫鏈的閉合通信類，這本質上是對轉換圖的連通組件的研究（Wikipedia建議一些線性算法來查找強連接組件）。此外，你可以通過所有開放的通信類，因爲每一個靜止狀態都必須消失。

如果你有你的封閉通信類C_1，...，C_n，那麼你的問題有望被分解成小的簡單片斷：每個封閉類C_i上的馬爾可夫鏈現在是不可約的，所以受限制的轉換矩陣M_i恰好具有一個特徵值爲0的特徵向量，並且該特徵向量僅具有正分量（參見Perron-Frobenius定理）。因此我們有一個固定狀態x_i。

您的完整馬爾可夫鏈的平穩狀態現在都是來自您的封閉類的x_i的所有線性組合。事實上，這些都是靜止狀態。

爲了找到固定狀態x_i，您可以連續應用M_i並且迭代將收斂到該狀態（這也將保持您的標準化）。一般而言，很難說明收斂速度，但它可以讓您輕鬆地提高解決方案的準確性並驗證解決方案。