logo标签列表

使用numpy或scipy查找马尔可夫稳态的左特征值。

pythonnumpyscipyeigenvectormarkov-chains
4
我需要使用转移矩阵的左特征向量来找到马尔可夫模型的稳态，需要使用一些Python代码。在这个问题中已经确定scipy.linalg.eig不能提供实际的左特征向量，但是在那里演示了一个修复方法。官方文档通常是无用和难以理解的。比格式不正确更大的问题是产生的特征值没有任何特定顺序（未排序且每次都不同）。因此，如果您想找到与1特征值相对应的左特征向量，则必须寻找它们，这会带来自己的问题（请参见下文）。数学很清楚，但如何让Python计算并返回正确的特征向量并不清楚。其他回答这个问题的答案，如这个答案，似乎没有使用左特征向量，因此那些不能成为正确的解决方案。在这个问题中提供了部分解决方案，但它没有考虑更大的转移矩阵的无序特征值。因此，只需使用
leftEigenvector = scipy.linalg.eig(A,left=True,right=False)[1][:,0]
leftEigenvector = leftEigenvector / sum(leftEigenvector)

[:,0]位置上的元素可能不是正确特征值对应的特征向量（在我的情况中通常不是），因此 is close，但通常并不适用。

好的，但是 scipy.linalg.eig(A,left=True,right=False) 的输出是一个数组，其中 [0] 元素是每个特征值（无序）的列表，并且根据相应的顺序在位置 [1] 后跟随一个特征向量数组。

我不知道如何从整个数组中按特征值排序或搜索以提取正确的特征向量（所有特征值为1的特征向量都除以向量元素之和标准化）。 我的想法是获取等于1的特征值的索引，然后从特征向量数组中提取这些列。 我的版本比较慢且繁琐。 首先我有一个函数（不完全有效），可以找到与给定值匹配的位置。

# Find the positions of the element a in theList
def findPositions(theList, a):
  return [i for i, x in enumerate(theList) if x == a]

然后我用它这样的方式获取与特征值 = 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查找马尔可夫模型的稳态分布的人。或许有更熟练/经验丰富的人已经找到了通用解决方案（无论是使用NumPy还是SciPy）。考虑到马尔可夫模型有多受欢迎，我期望有一个专门用于执行此任务的库，也许它存在，但我找不到。

我对马尔可夫链分析并不是特别熟悉，但听起来你已经掌握了大部分。在特定的特征值中搜索所需的值似乎是一个可行的解决方案。如果性能是一个问题，也许可以发布你正在使用的代码，这样其他人就可以看一下了？ - Marco Tompitak
性能不是一个严重的问题，只要它实际上工作就可以了。目前在搜索与1相等的特征值列表时无法正常工作：它会错过其中一个。另外两个也不对应实际稳态，所以可能特征向量条目的顺序也以某种方式混乱，而文档没有指示。或者我的矩阵有多个单独的吸引子状态。 - Aaron Bramson
2
可能是一个精度问题？你可以尝试使用x < a + epsilon && x > a - epsilon这样的表达式，其中epsilon是一个适当的小数。 - Marco Tompitak
2
w，vl = eig（P，right = False，left = True）; tol = 1e-15; v1 = vl [:，abs（w-1）<tol] .T.conj（）; assert numpy.allclose（v1.dot（P），v1） - pv.
2个回答
6
你提供了一个链接：如何找到矩阵特定特征值对应的特征向量？，并且说它没有计算左特征向量，但是你可以通过使用转置矩阵来解决这个问题。
例如：
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]])

为了将特征向量标准化（你应该知道它应该是实数）：
In [913]: u = (evec/evec.sum()).real

In [914]: u
Out[914]: 
array([[ 0.18060201],
       [ 0.36789298],
       [ 0.37625418],
       [ 0.07525084]])

In [915]: np.dot(u.T, M).T
Out[915]: 
array([[ 0.18060201],
       [ 0.36789298],
       [ 0.37625418],
       [ 0.07525084]])

如果您事先不知道特征值 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 ]])
1

好的，当我在实施沃伦的解决方案时，我不得不做一些更改，我已经在下面包含了这些内容。基本上它是相同的，所以他获得了所有的荣誉，但是使用numpy和scipy进行数值近似的现实需要更多的调整，我认为对于其他试图在将来做到这一点的人来说，这些调整会很有帮助。我还将变量名称更改为非常适合新手的名称。

如果我做错了什么或者有进一步推荐的改进（例如加速），请告诉我。

# 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)

将所有内容整合成易于复制粘贴的格式：

将所有内容整合成易于复制粘贴的格式：

M = transitionMatrix( G ) 
theEigenvalues, leftEigenvectors = scipy.linalg.eig(M, right=False, left=True)  
theEigenvalues = theEigenvalues.real                 
leftEigenvectors = leftEigenvectors.real
tolerance = 1e-10            
mask = abs(theEigenvalues - 1) < tolerance 
theEigenvalues = theEigenvalues[mask]    
leftEigenvectors = leftEigenvectors[:, mask] 
leftEigenvectors[leftEigenvectors < tolerance] = 0  
attractorDistributions = leftEigenvectors / leftEigenvectors.sum(axis=0, keepdims=True)   
attractorDistributions = attractorDistributions.T
theSteadyStates = np.sum(attractorDistributions, axis=0)

使用这种分析方法在生成的马尔可夫模型中产生了一个吸引子（总共三个），其稳态分布为0.19835218和0.80164782，与数学上准确值0.2和0.8相比有超过0.1％的误差，对于科学来说是一个很大的误差。尽管如此，这并不是真正的问题，因为如果精度很重要，那么现在已经确定了各个吸引子，可以使用矩阵子集对每个吸引子内的行为进行更准确的分析。
看起来我需要更多的帮助。这个代码在去年12月底完美地运行，每次都为手动生成的马尔可夫矩阵提供相同且正确的特征向量。现在在相同的矩阵上重新运行相同的代码，会产生几组不同和不正确的左特征向量...而且每次运行结果都不同。它之前一直很好用，没有任何改变，现在却出了问题并且不一致。这怎么可能？有人有类似的经验吗？有什么修复建议吗？ - Aaron Bramson
最后，如果这还不够奇怪的话，当我从实验室电脑运行代码时，每次都能正常工作，但是从我的笔记本电脑运行时会出现上述问题。相同的代码，相同的WinPython x64 2.7.9.4安装文件和相同的操作系统。有人知道可能是什么原因导致这种差异吗？ - Aaron Bramson
你能否更新这个答案，包括一个完整的重现案例。了解可能出错的第一步是自己运行代码，但你所发布的不完整或不自包含。 - talonmies
好的，代码就是上面那些加上导入numpy和scipy，我会在这个问题中添加一个具有多个吸引子的特定矩阵，以生成一个最小运行示例。 - Aaron Bramson
“transitionMatrix” 是我认为最重要的缺失之一。 - talonmies
我已经解决了问题，或者说是某种程度上解决了，但唯一的解决方案是不使用scipy。上述代码是正确的，但scipy在寻找马尔可夫模型的稳态分布方面并不可靠，应该避免使用（至少对于这个目的来说）。 - Aaron Bramson
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接