找到具有左特征值的马尔可夫稳态（使用 numpy 或 scipy）

Question

我需要使用一些 python 代码使用转换矩阵的左特征向量找到马尔可夫模型的稳态。

在这个问题中已经确定scipy.linalg.eig 无法提供所描述的实际左特征向量，但在那里演示了修复。官方文档和往常一样，大多无用且难以理解。

比错误格式更大的问题是生成的特征值没有任何特定的顺序（没有排序并且每次都不同）。所以如果你想找到对应于 1 个特征值的左特征向量，你必须寻找它们，这会带来它自己的问题（见下文）。数学很清楚，但如何让 python 计算它并返回正确的特征向量尚不清楚。这个问题的其他答案，比如这个，似乎没有使用左特征向量，所以那些不是正确的解决方案。

这个问题提供了部分解决方案，但它没有考虑较大转换矩阵的无序特征值。所以，只需使用

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 来寻找马尔可夫模型平稳分布的人。更精通/经验丰富的人可能有一个可行的通用解决方案（无论是否使用 numpy 或 scipy）。考虑到马尔可夫模型的流行程度，我希望有一个专门供它们执行此任务的库，也许它确实存在，但我找不到。

Answer 1

您链接到如何找出与矩阵的特定特征值相对应的特征向量？并说它不计算左特征向量，但您可以通过使用转置来解决这个问题。

例如，

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