找到具有左特征值的马尔可夫稳态（使用 numpy 或 scipy）
find Markov steady state with left eigenvalues (using numpy or scipy)
我需要使用一些 python 代码使用转移矩阵的左特征向量找到马尔可夫模型的稳态。
已在 this question 中确定 scipy.linalg.eig 无法提供所描述的实际左特征向量，但此处演示了修复方法。官方文档和往常一样大多无用且难以理解。
比不正确的格式更大的问题是生成的特征值没有按任何特定顺序排列（未排序且每次都不同）。因此，如果您想找到对应于 1 个特征值的左特征向量，您必须寻找它们，这会带来自己的问题（见下文）。数学很清楚，但是如何让 python 计算这个和 return 正确的特征向量还不清楚。这个问题的其他答案，比如 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 的特征向量都由向量条目的总和归一化。）我的想法是得到等于 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 来寻找马尔可夫模型平稳分布的人。更 proficient/experienced 的人可能有一个可行的通用解决方案（无论是否使用 numpy 或 scipy）。考虑到马尔可夫模型的流行程度，我希望有一个专门供他们执行此任务的库，也许它确实存在，但我找不到。
你链接到 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]])
归一化特征向量（你知道它应该是真实的）：
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 ]])
好吧，我在实施 Warren 的解决方案时不得不进行一些更改，我已将这些更改包括在下面。它基本上是一样的，所以他得到了所有的荣誉，但是 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% 的折扣，这对科学来说是一个很大的错误。这不是一个真正的问题，因为如果准确性很重要，那么既然已经识别了各个吸引子，就可以使用矩阵子集对每个吸引子内的行为进行更准确的分析。
我需要使用一些 python 代码使用转移矩阵的左特征向量找到马尔可夫模型的稳态。
已在 this question 中确定 scipy.linalg.eig 无法提供所描述的实际左特征向量，但此处演示了修复方法。官方文档和往常一样大多无用且难以理解。
比不正确的格式更大的问题是生成的特征值没有按任何特定顺序排列（未排序且每次都不同）。因此，如果您想找到对应于 1 个特征值的左特征向量，您必须寻找它们，这会带来自己的问题（见下文）。数学很清楚，但是如何让 python 计算这个和 return 正确的特征向量还不清楚。这个问题的其他答案，比如 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 的特征向量都由向量条目的总和归一化。）我的想法是得到等于 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 来寻找马尔可夫模型平稳分布的人。更 proficient/experienced 的人可能有一个可行的通用解决方案（无论是否使用 numpy 或 scipy）。考虑到马尔可夫模型的流行程度，我希望有一个专门供他们执行此任务的库，也许它确实存在，但我找不到。
你链接到 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]])
归一化特征向量（你知道它应该是真实的）：
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 ]])
好吧，我在实施 Warren 的解决方案时不得不进行一些更改，我已将这些更改包括在下面。它基本上是一样的，所以他得到了所有的荣誉，但是 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% 的折扣，这对科学来说是一个很大的错误。这不是一个真正的问题，因为如果准确性很重要，那么既然已经识别了各个吸引子，就可以使用矩阵子集对每个吸引子内的行为进行更准确的分析。