用左特征值求马尔可夫稳态（使用numpy或scipy）

Question

我需要使用一些python代码，利用Markov模型的转移矩阵的左特征向量来求出它们的稳态。在

在this question中已经确定西皮.利纳格.艾格无法提供所述的实际左特征向量，但在那里演示了修复。官方文件大多是无用的，像往常一样令人费解。在

比不正确的格式更大的问题是，产生的特征值不是以任何特定的顺序排列的（没有排序，每次都不同）。因此，如果你想找到对应于1个特征值的左特征向量，你必须去寻找它们，这就带来了它自己的问题（见下文）。数学是清楚的，但是如何让python计算这个值并返回正确的特征向量还不清楚。这个问题的其他答案，比如this one，似乎没有使用左特征向量，所以这些不可能是正确的解。在

This question提供了一个部分解，但它不能解释较大转移矩阵的无序特征值。所以，就用

leftEigenvector = scipy.linalg.eig(A,left=True,right=False)[1][:,0]
leftEigenvector = leftEigenvector / sum(leftEigenvector)

很接近，但一般不起作用，因为[:,0]位置的条目可能不是正确特征值的特征向量（在我的例子中通常不是）。在

好的，但是scipy.linalg.eig(A,left=True,right=False)的输出是一个数组，其中[0]元素是每个特征值的列表（不是按任何顺序排列），在{}位置后面是与这些特征值相对应的特征向量数组。在

我不知道一个很好的方法来通过特征值对整个事情进行排序或搜索，从而找出正确的特征向量（特征值为1的所有特征向量都被向量项的和规范化了）。我的想法是得到特征值的索引等于1，然后从特征向量数组中提取这些列。我的版本是缓慢和繁琐的。首先，我有一个函数（不太管用）来查找与值匹配的鞋楦位置：

^{pr2}$

然后我用它来得到匹配特征值=1的特征向量。在

M = transitionMatrix( G )
leftEigenvectors = scipy.linalg.eig(M,left=True,right=False)
unitEigenvaluePositions = findPositions(leftEigenvectors[0], 1.000)
steadyStateVectors = []
for i in unitEigenvaluePositions:
    thisEigenvector = leftEigenvectors[1][:,i]
    thisEigenvector / sum(thisEigenvector)
    steadyStateVectors.append(thisEigenvector)
print steadyStateVectors

但实际上这不管用。有一个特征值=1.00000000e+00 +0.00000000e+00j找不到，即使另外两个也找不到。在

我的期望是我不是第一个使用python来发现Markov模型的平稳分布的人。更精通/更有经验的人可能有一个可行的通用解决方案（无论是否使用numpy或scipy）。考虑到马尔可夫模型有多流行，我希望有一个库来完成这个任务，也许它确实存在，但我找不到。在

Answer 1

您链接到How do I find out eigenvectors corresponding to a particular eigenvalue of a matrix?并说它不计算左特征向量，但是您可以通过处理转置来修复它。在

例如

In [901]: import numpy as np

In [902]: import scipy.sparse.linalg as sla

In [903]: M = np.array([[0.5, 0.25, 0.25, 0], [0, 0.1, 0.9, 0], [0.2, 0.7, 0, 0.1], [0.2, 0.3, 0, 0.5]])

In [904]: M
Out[904]: 
array([[ 0.5 ,  0.25,  0.25,  0.  ],
       [ 0.  ,  0.1 ,  0.9 ,  0.  ],
       [ 0.2 ,  0.7 ,  0.  ,  0.1 ],
       [ 0.2 ,  0.3 ,  0.  ,  0.5 ]])

In [905]: eval, evec = sla.eigs(M.T, k=1, which='LM')

In [906]: eval
Out[906]: array([ 1.+0.j])

In [907]: evec
Out[907]: 
array([[-0.32168797+0.j],
       [-0.65529032+0.j],
       [-0.67018328+0.j],
       [-0.13403666+0.j]])

In [908]: np.dot(evec.T, M).T
Out[908]: 
array([[-0.32168797+0.j],
       [-0.65529032+0.j],
       [-0.67018328+0.j],
       [-0.13403666+0.j]])

要规范化特征向量（您知道应该是实数）：

^{pr2}$

如果您事先不知道特征值1的重数，请参阅@pv的注释，使用scipy.linalg.eig显示代码。下面是一个例子：

In [984]: M
Out[984]: 
array([[ 0.9 ,  0.1 ,  0.  ,  0.  ,  0.  ,  0.  ],
       [ 0.3 ,  0.7 ,  0.  ,  0.  ,  0.  ,  0.  ],
       [ 0.  ,  0.  ,  0.25,  0.75,  0.  ,  0.  ],
       [ 0.  ,  0.  ,  0.5 ,  0.5 ,  0.  ,  0.  ],
       [ 0.  ,  0.  ,  0.  ,  0.  ,  0.  ,  1.  ],
       [ 0.  ,  0.  ,  0.  ,  0.  ,  1.  ,  0.  ]])

In [985]: import scipy.linalg as la

In [986]: evals, lvecs = la.eig(M, right=False, left=True)

In [987]: tol = 1e-15

In [988]: mask = abs(evals - 1) < tol

In [989]: evals = evals[mask]

In [990]: evals
Out[990]: array([ 1.+0.j,  1.+0.j,  1.+0.j])

In [991]: lvecs = lvecs[:, mask]

In [992]: lvecs
Out[992]: 
array([[ 0.9486833 ,  0.        ,  0.        ],
       [ 0.31622777,  0.        ,  0.        ],
       [ 0.        , -0.5547002 ,  0.        ],
       [ 0.        , -0.83205029,  0.        ],
       [ 0.        ,  0.        ,  0.70710678],
       [ 0.        ,  0.        ,  0.70710678]])

In [993]: u = lvecs/lvecs.sum(axis=0, keepdims=True)

In [994]: u
Out[994]: 
array([[ 0.75, -0.  ,  0.  ],
       [ 0.25, -0.  ,  0.  ],
       [ 0.  ,  0.4 ,  0.  ],
       [ 0.  ,  0.6 ,  0.  ],
       [ 0.  , -0.  ,  0.5 ],
       [ 0.  , -0.  ,  0.5 ]])

In [995]: np.dot(u.T, M).T
Out[995]: 
array([[ 0.75,  0.  ,  0.  ],
       [ 0.25,  0.  ,  0.  ],
       [ 0.  ,  0.4 ,  0.  ],
       [ 0.  ,  0.6 ,  0.  ],
       [ 0.  ,  0.  ,  0.5 ],
       [ 0.  ,  0.  ,  0.5 ]])

Answer 2

好吧，在实施沃伦的解决方案时，我必须做出一些改变，我已经在下面列出了这些。基本上是一样的，所以他得到了所有的荣誉，但是numpy和scipy的数值近似的现实需要更多的按摩，我想这对其他将来尝试这样做的人会有帮助。我还将变量名改为超级noob友好型。在

如果我有什么问题或有进一步的改进建议（例如速度方面的改进），请告诉我。在

# in this case my Markov model is a weighted directed graph, so convert that nx.graph (G) into it's transition matrix
M = transitionMatrix( G )   

#create a list of the left eigenvalues and a separate array of the left eigenvectors
theEigenvalues, leftEigenvectors = scipy.linalg.eig(M, right=False, left=True)  

# for stationary distribution the eigenvalues and vectors are always real, and this speeds it up a bit
theEigenvalues = theEigenvalues.real                 
leftEigenvectors = leftEigenvectors.real

# set how close to zero is acceptable as being zero...1e-15 was too low to find one of the actual eigenvalues
tolerance = 1e-10 

# create a filter to collect the eigenvalues that are near enough to zero                               
mask = abs(theEigenvalues - 1) < tolerance           

# apply that filter
theEigenvalues = theEigenvalues[mask]                

# filter out the eigenvectors with non-zero eigenvalues
leftEigenvectors = leftEigenvectors[:, mask] 

# convert all the tiny and negative values to zero to isolate the actual stationary distributions    
leftEigenvectors[leftEigenvectors < tolerance] = 0  

# normalize each distribution by the sum of the eigenvector columns
attractorDistributions = leftEigenvectors / leftEigenvectors.sum(axis=0, keepdims=True)   

# this checks that the vectors are actually the left eigenvectors, but I guess it's not needed to usage 
#attractorDistributions = np.dot(attractorDistributions.T, M).T 

# convert the column vectors into row vectors (lists) for each attractor (the standard output for this kind of analysis)
attractorDistributions = attractorDistributions.T

# a list of the states in any attractor with the approximate stationary distribution within THAT attractor (e.g. for graph coloring)         
theSteadyStates = np.sum(attractorDistributions, axis=1)  

以一种简单的复制和粘贴格式将它们组合在一起：

^{pr2}$

通过对生成的Markov模型的分析，得到了一个吸引子（三个），稳态分布为0.19835218和0.80164782，而数学上的精确值为0.2和0.8。所以这是超过0.1%的折扣，对科学来说是一个很大的错误。这不是一个真正的问题，因为如果精度很重要的话，那么既然单个吸引子已经被识别，那么就可以使用矩阵子集对每个吸引子内部的行为进行更精确的分析。在