Python中的numpy/scipy特征向量似乎对马尔可夫链模型不正确

Question

我有一个大的（351,351）numpy过渡矩阵。我想使用numpy（我还尝试了具有完全相同功能的scipy）找到此矩阵的状态向量。

sstate = np.linalg.eig(T)[1][:,0]

我认为这应该给我主左特征值的特征向量。主左特征值是1+0j。这有些正确，主左特征值应该是1，但我对虚数不太了解。此外，sstate向量包含所有复数。现在，尝试检查是否正确，我进行以下矩阵乘法：

np.dot(sstate,T)

如果正确执行，这应该返回与'sstate'相同的向量。我不确定为什么这不起作用。虚数可能是问题吗？此外，这个转移矩阵可能不包含稳态向量。我的转移状态矩阵中的每一行和每一列都应该总和为1，然而，我发现舍入误差导致每一行和每一列的总和只约为1。

非常感谢您的帮助！

- mepstein1218

Answer 1

转移矩阵对称吗？如果不是，则考虑检查 T.T（转置），因为您需要确保查看正确的状态转换：您需要随机矩阵的左特征向量，但几乎所有开箱即用的科学软件包（包括numpy）默认计算右特征向量（这也是在处理此类问题时在教科书和其他材料中必须通过行向量预乘以而不是通常的矩阵列乘法的原因）。

也许还可以使用 sstate = sstate/sstate.sum() 确保概率总和为1，尽管存在舍入误差。

这里有一个使用numpy的示例an example with numpy. 从评论中添加关于右特征向量与左特征向量的详细信息： eig和类似它的函数可以计算右特征向量，如向量v，使得Av = (lambda)v，其中lambda是标量。但你需要的是矩阵A的左特征向量，满足v.T*A = (lambda)v.T，这不仅仅是右特征向量的转置或共轭。

因此，你需要基于A.T计算特征向量，但在检查状态向量是否真正稳定时，不要使用A.T进行计算。你需要查看np.dot(sstate, T)（确保sstate是行向量而非列向量），并对其进行评估（可能还需要考虑舍入误差问题并重新归一化）。

- ely

---

我刚刚检查了一下，矩阵并不对称，这很奇怪，因为它应该是对称的。我得去查一下这个问题。 - mepstein1218

转换并不一定要对称，但通常是这样的。例如，我链接的示例中的转换不对称，因为问题陈述不对称。 - ely

哦，好的，在检查了数学计算后，我的矩阵不应该是对称的。那么这会对numpy/scipy造成什么问题呢？ - mepstein1218

我尝试找到转置转移矩阵T.T的稳态向量，但即使如此，这个稳态向量乘以T.T也不等于它本身。 - mepstein1218

eig和类似的函数会计算出正确的特征向量，即满足Av = (lambda)v（其中lambda为标量）的向量v。但是你需要的是矩阵A的左特征向量，也就是满足v.T*A = (lambda)v.T的向量，这个向量不仅仅是右特征向量的转置或共轭。 - ely

你需要基于 A.T 计算特征向量，但在验证状态向量是否真的是稳态时，不要再使用 A.T 进行计算。你需要查看 np.dot(sstate, T)（验证 sstate 是行向量而不是列向量），并对其进行评估（可能还要考虑舍入误差的问题）。 - ely

Answer 2

因此，请查看从调用eig返回的1D数组（2元组中的第一个元素）； 这是特征值数组，正如您所看到的，它并没有按降序排列，您需要手动排序，然后将相同的顺序应用于特征向量数组。 您可能已经这样做了，但它不在您的代码片段中，也没有在OP中提到。

一旦您这样做了，就可以选择第一个特征向量：

>>> import numpy as NP
>>> from scipy import linalg as LA
>>> a = NP.random.rand(16).reshape(4, 4)

>>> E = LA.eig(a, left=True)

>>> evals, evecs = E

将特征值按降序排序

>>> idx = NP.argsort(evals)[::-1]

>>> eva = eva[idx]

将排序索引应用于特征向量矩阵

>>> eva[idx,]

选择第一个

>>> eva[0].real 

此外，使用来自scipy的linalg；如果机器上没有它，则NumPy安装程序包含一个通用BLAS可供构建

另外，如果您传递给eig的2D数组是稀疏的，则使用来自scipy.sparse的eig；它更快，特别是在您只需要一个特征向量的情况下。

- doug

---

我正在使用正确的特征值，这已经被检查过了。不过还是谢谢你的回答和提示 :) - mepstein1218