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

Posted 2021-04-26

我构建了一个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将复杂的特征值和特征向量更改为真实的？

答案

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

另一答案

Transpose fixes eigenvalues

显然，LAPACK讨厌非对称三对角矩阵，其中较大的非对角线元素位于对角线下方。使用转置矩阵，其中较大的元素在对角线上方，产生真实的特征值。 （从理论上讲，矩阵及其转置具有相同的特征值。）除了真实之外，它们与theoretical values一致，直到排序和合理的错误。

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。

Eigenvectors

你也想要特征向量。 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，left，right”，但由于k被换位，我左右切换。

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

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

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

Bonus content

那么...问题解决了吗？好吧，我没有说我可以对齐你的矩阵。试图反转由左或右特征向量形成的矩阵看起来像一个失败的命题;逆是可怕的。

另外，我们不应该过多地相信上述特征向量的测试。它为正确的特征向量提供了一个微小的误差...让我们在错误的特征向量上尝试它，那些具有非平凡的虚部。

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。错误的特征向量甚至比真实的更好！