我试图确定一个100 × 100矩阵的(真实的)左特征向量和特征值。矩阵T的元素是从0到1的概率。矩阵T的每一行总和为1.0(100%概率)。列的总和不一定是任何特定的值。
scipy和numpy库返回复数-这显然是由于无法将特征多项式求解为真实的或算法被优化来做到这一点。我转向sympy库,它也返回复数作为解决方案。
使用sympy在这里描述是有意义的:Is it possible somehow to find complex eigenvalues using SymPy?
我的理解是,这将迫使答案为真实的。
下面是我的代码:
import numpy as np
from sympy import *
T = np.loadtxt('rep10_T_ij.dat', delimiter=' ')
T = np.asmatrix(T)
T = Matrix(T)
lam = symbols('\lambda')
char_poly = det(lam * eye(T.rows) - T)
roots_char_poly = solve(char_poly, lam)
display(char_poly)
display(roots(poly(char_poly, domain=CC)))
display([N(r) for r in roots_char_poly])
然而,这并没有在12小时内输出答案,任何数量的硬件都扔在它上面。
我在这里读到:https://github.com/sympy/sympy/issues/7200
100 × 100矩阵是否太贵?我可以在集群上向它扔更多的硬件,但还没能尝试这个。有解决办法吗?
1条答案
按热度按时间qgzx9mmu1#
在这种情况下,使用SymPy而不是NumPy或SciPy的优势在于,如果您想要精确或符号化地执行计算。确定特征值是精确真实的通常需要精确的算术。精确或符号计算比固定精度浮点要慢得多,所以与使用NumPy或SciPy相比,你应该期待事情会慢很多。如果你的初始矩阵是一个浮点数矩阵,那么它已经是近似的,在这种情况下,使用SymPy很有可能没有好处。
在任何情况下,这里有一些代码可以显示如何使用SymPy更有效地获得矩阵的确切真实的特征值。这通常只在你的原始矩阵有精确的有理数(不是浮点数)时才值得做,但它应该精确地获得真实的特征值:
这使用Berkowicz算法来计算特征多项式,这是操作中最慢的部分(
charpoly
方法)。这比使用行列式找到特征多项式要有效得多,但对于大型矩阵仍然很慢。对于N x N
稠密矩阵,如果底层系数字段具有O(1)
算术成本,则以这种方式找到特征多项式大致为O(N**4)
。不幸的是,精确有理数没有O(1)
的算术成本,所以实际上,当你扩展到更大的矩阵时,复杂度明显比O(N**4)
差。您可以了解不同大小所涉及的时间:
因此,已经找到40 x40矩阵的确切整数特征值需要30秒。扩展到100 x100会很慢,但可能是可能的。如果你想这样做,我建议安装
gmpy2
,因为它可以让SymPy中的精确有理数更快。最有可能的是,像这样计算精确的特征值并不是你想要做的,你应该使用NumPy/SciPy。如果它们给予的本征值虚部很小,那么你可以丢弃那些,如果你知道本征值应该是真实的。否则,只要初始矩阵仅是近似的(即,浮点数),那么上面所示的那种缓慢的精确计算可能是不值得的。