如何获得三角形Toeplitz矩阵的真实特征值和特征向量

Question

我构造了100*100矩阵k，并希望使用numpy.linalg.eig对其化。

k=np.zeros((100,100))
np.fill_diagonal(k,-2)
np.fill_diagonal(k[1:,:-1],1.5)
np.fill_diagonal(k[:-1,1:],0.5)

当我尝试较小的矩阵时，例如

w,v=np.linalg.eig(k[:10,:10])

特征值w和特征向量v是真实的。但是，当我尝试更大的矩阵或整个完整矩阵

时

w,v=np.linalg.eig(k)

w和 v结果是复数，而虚构部分不可忽略。

我也尝试scipy.linalg.eig，它也有类似的问题。

我想拿出特征值和特征向量的自然对数。我的模型中没有复数的物理含义。

我只能拥有独立的实际数字特征值和特征向量？如果没有，如何将复杂的特征值和特征向量更改为python的真实值？

Answer 1

@daniel F和@ftp比我快，请参阅他们的评论，但是由于我将代码坐在这里，所以我最好分享：

import numpy as np
from scipy import sparse
def tri_toep_eig(a, b, c, n):
    evals = a + 2*np.sqrt(b*c) * np.cos(np.pi * np.arange(1, n+1) / (n+1))
    evecs = np.sin(np.outer(np.arange(1, n+1) * np.pi / (n+1),
                            np.arange(1, n+1))) 
        * np.sqrt(b/c)**np.arange(n)[:, None]
    return evals, evecs
def tri_toep(a, b, c, n):
    return sparse.dia_matrix((np.outer((b, a, c), np.ones((n,))),
                              (-1, 0, 1)), (n, n))
def check(a, b, c, n):
    evals, evecs = tri_toep_eig(a, b, c, n)
    tt = tri_toep(a, b, c, n)
    for eva, eve in zip(evals, evecs.T):
        assert np.allclose(tt @ eve, eva * eve)
check(-2, 0.5, 1.5, 100)

Answer 2

transpose修复特征值

显然，拉帕克讨厌非对称的三角形矩阵，其中较大的非对抗元素在对角线以下。使用较大元素在对角线上方的转置矩阵会导致实际特征值。(从理论上讲，矩阵及其转置具有相同的特征值。(除了真实之外，他们还同意分类和合理的错误的理论价值。

a, b, c, n = -2, 0.5, 1.5, 100
k = np.zeros((n, n))
np.fill_diagonal(k, a)
np.fill_diagonal(k[:-1, 1:], b)
np.fill_diagonal(k[1:, :-1], c)
theory = a - 2*np.sqrt(b*c) * np.cos(np.pi * np.arange(1, n+1) / (n+1))
computed = np.sort(np.linalg.eig(k.T)[0])
print(np.max(np.abs(theory - computed)))

这是打印6.183001044490766e-08，这是计算和理论特征值之间最大的差异。没有转座t，此错误最高为0.26。

特征向量

您还想要特征向量。np.eig返回的转置矩阵的特征向量是原始矩阵的左特征向量：也就是说，它们满足vec.dot(k) = lam*vec而不是k.dot(vec) = lam*vec。如果要获得原始矩阵的正确特征向量，请使用Scipy的eig：

from scipy import linalg as la
evals, right_evects, left_evects = (np.real(_) for _ in la.eig(k.T, left=True, right=True))

scipy的eigensolver与numpy的不同之处在于，它返回带有 +0j的特征向量和特征值；它认为它们复杂，但正确地评估了假想部分为0。我截断了上述假想部分。请注意，Scipy返回顺序是" EVALS，左，右"，但由于K被转移，我向左和右转。

让我们检查这些特征向量：

np.max([np.linalg.norm(k.dot(vec) - ev*vec) for ev, vec in zip(evals, right_evects.T)])

返回1.845213984555825e-14，还不错。带有特征向量的阵列上的转置是因为zip从矩阵中挑选行，我们需要列。

奖励内容

所以...问题解决了吗？好吧，我没有说我可以对角线矩阵。试图将左右特征向量形成的矩阵颠倒看起来像是一个失败的命题。倒数是可怕的。

另外，我们不应该过多地相信上述特征向量的测试。它给出了右特征向量的微小错误 ...让我们在错误的特征向量上尝试一下，，那些具有非平凡的想象中的部分。

wrong_evals, wrong_evects = np.linalg.eig(k)
np.max([np.linalg.norm(k.dot(vec) - ev*vec) for ev, vec in zip(wrong_evals, wrong_evects.T)])

这返回1.7136373586499598e-14。错误的特征向量甚至比真实的东西更好！