Python numpy/scipy 特征向量似乎不适合马尔可夫链模型

Python numpy/scipy eigenvectors seemingly not correct for markov chain model

我有一个很大的 (351,351) numpy 转换矩阵。我想使用 numpy 找到状态向量(我也试过 scipy 具有相同的功能)。

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

所以我相信这应该给我占主导地位的左特征值的特征向量。主左特征值为 1+0j。这在某种程度上是正确的,主要的左特征值应该是 1,我对此有点陌生,所以我不确定如何处理虚数。此外,sstate 向量包含所有复数。现在,尝试检查这是否正确,我执行以下矩阵乘法:

np.dot(sstate,T)

如果操作正确,这应该 return 与 'sstate.' 相同的矢量 我不确定为什么这不起作用。虚数可能是问题所在吗?另外,这个转移矩阵是否可能不包含稳态向量。我的过渡状态矩阵中的每一行和每一列的总和应该为 1,但是,我发现舍入误差导致每一行和列的总和仅为大约 1。

感谢任何帮助!

转移矩阵是对称的吗?如果不是,请考虑检查 T.T(转置),因为您需要确保您正在查看正确的状态转换:您需要随机的 left 特征向量矩阵,但几乎所有开箱即用的科学包(包括 numpy)默认计算 right 特征向量(这与教科书中的相同原因和你必须预乘行向量的东西而不是处理这些东西时通常的矩阵列乘法)。

也许还 sstate = sstate/sstate.sum() 以确保概率总和为 1,尽管有舍入。

这里是an example with numpy.

从评论中添加了关于右特征向量与左特征向量的详细信息:

eig 和类似的东西将计算 正确的 特征向量,如向量 v 这样 Av = (lambda)v 对于标量 lambda.不过,您需要的是 Aleft 特征向量,因此满足 v.T*A = (lambda)v.T 的东西不仅是右特征向量的转置或共轭。

所以你会想根据 A.T 计算特征向量,但你 不会 想在稍后检查状态向量时用 A.T 计算真的是静止的。您需要查看 np.dot(sstate, T)(验证 sstate 是行向量,而不是列),并对其进行评估(可能还有关于重新归一化以帮助四舍五入的另一点)。

所以看看调用 eig 返回的一维数组(二元组中的第一个元素);这是特征值数组,如您所见,它不是按降序排列的,您需要手动 sort 然后应用相同的顺序到特征向量数组。您可能已经这样做了,但它不在您的代码片段中,也没有在 OP 中提及。

一旦你这样做了,那么你就可以 select 第一个特征向量:

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

select第一个

>>> eva[0].real 

另外,使用 scipy 中的 linalg; NumPy 安装程序包含一个通用的 BLAS,如果它不在机器上,则可以用来构建它

此外,如果您传递给 eig 的二维数组是稀疏的,则使用 eig 来自 scipy.sparse;它要快得多,尤其是在你的情况下,因为你只需要一个特征向量。