如何获得三对角 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

转置修正特征值

显然，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。

特征向量

您还需要特征向量。 np.eig 为转置矩阵返回的特征向量是原始矩阵的 left 特征向量：也就是说，它们满足 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 的特征解算器与 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 从矩阵中选择行，而我们需要列。

奖励内容

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

此外，我们不应该过于相信上述特征向量测试。它对正确的特征向量给出了一个微小的错误......让我们在错误的特征向量上尝试一下，那些具有非平凡虚部的。

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

Answer 2

@Daniel F 和 @FTP 比我快，看看他们的 cmets，但既然我有代码，我不妨分享一下：

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)