Python 中的特征值:一个错误

Question

这里有两个关于特征向量和平方矩阵的特征值的假设。我相信两者都是真的：

如果矩阵是

1. 对称的并且只包含实值，那么它是一个埃尔米特矩阵，那么所有特征值都应该是实数，所有特征向量的所有分量也应该是实数。从埃尔米特矩阵计算特征向量和特征值时，结果中不应出现复数。

2. 从给定矩阵计算的给定特征值的特征向量应始终指向仅由矩阵和特征值确定的方向。用于计算它的算法对结果没有影响，只要算法正确实现即可。

但是，当您在 Python 中使用标准库来计算特征向量和特征值时，这两种假设都站不住脚。这些方法是否包含错误？

有四种不同的方法可以从埃尔米特矩阵计算特征值和特征向量：

1. numpy.linalg.eig
2. scipy.linalg.eig
3. numpy.linalg.eigh
4. scipy.linalg.eigh

#1 和 #2 可用于任何方阵(包括埃尔米特矩阵).
#3 和 #4 仅用于埃尔米特矩阵。据我所知，它们的目的只是它们运行得更快，但结果应该是相同的(只要输入真的是埃尔米特)。

但是这四种方法对于相同的输入提供了三种不同的结果。这是我用来测试所有四种方法的程序：

#!/usr/bin/env python3
import numpy as np
import scipy.linalg as la
A = [
[19, -1, -1, -1, -1, -1, -1, -1],
[-1, 19, -1, -1, -1, -1, -1, -1],
[-1, -1, 19, -1, -1, -1, -1, -1],
[-1, -1, -1, 19, -1, -1, -1, -1],
[-1, -1, -1, -1, 19, -1, -1, -1],
[-1, -1, -1, -1, -1, 19, -1, -1],
[-1, -1, -1, -1, -1, -1, 19, -1],
[-1, -1, -1, -1, -1, -1, -1, 19]
]
A = np.array(A, dtype=np.float64)
delta = 1e-12
A[5,7] += delta
A[7,5] += delta
if np.array_equal(A, A.T):
print('input is symmetric')
else:
print('input is NOT symmetric')
methods = {
'np.linalg.eig'  : np.linalg.eig,
'la.eig'         : la.eig,
'np.linalg.eigh' : np.linalg.eigh,
'la.eigh'        : la.eigh
}
for name, method in methods.items():
print('============================================================')
print(name)
print()
eigenValues, eigenVectors = method(A)

for i in range(len(eigenValues)):
print('{0:6.3f}{1:+6.3f}i '.format(eigenValues[i].real, eigenValues[i].imag), end=' |  ')
line = eigenVectors[i]
for item in line:
print('{0:6.3f}{1:+6.3f}i '.format(item.real, item.imag), end='')
print()
print('---------------------')
for i in range(len(eigenValues)):
if eigenValues[i].imag == 0:
print('real    ', end=' |  ')
else:
print('COMPLEX ', end=' |  ')
line = eigenVectors[i]
for item in line:
if item.imag == 0:
print('real    ', end='')
else:
print('COMPLEX ', end='')
print()
print()

这是它产生的输出：

input is symmetric
============================================================
np.linalg.eig
12.000+0.000i  |  -0.354+0.000i  0.913+0.000i  0.204+0.000i -0.013+0.016i -0.013-0.016i  0.160+0.000i -0.000+0.000i  0.130+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.183+0.000i  0.208+0.000i  0.379-0.171i  0.379+0.171i -0.607+0.000i  0.000+0.000i -0.138+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.182+0.000i  0.203+0.000i -0.468-0.048i -0.468+0.048i  0.153+0.000i  0.001+0.000i -0.271+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.182+0.000i  0.203+0.000i  0.657+0.000i  0.657-0.000i  0.672+0.000i -0.001+0.000i  0.617+0.000i 
20.000-0.000i  |  -0.354+0.000i -0.182+0.000i  0.203+0.000i -0.276+0.101i -0.276-0.101i -0.361+0.000i  0.001+0.000i -0.644+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.001+0.000i -0.612+0.000i -0.001+0.000i -0.001-0.000i  0.001+0.000i  0.706+0.000i -0.000+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.182+0.000i  0.203+0.000i -0.276+0.101i -0.276-0.101i -0.018+0.000i -0.000+0.000i  0.306+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.001+0.000i -0.612+0.000i -0.001+0.000i -0.001-0.000i  0.001+0.000i -0.708+0.000i  0.000+0.000i 
---------------------
real     |  real    real    real    COMPLEX COMPLEX real    real    real    
real     |  real    real    real    COMPLEX COMPLEX real    real    real    
real     |  real    real    real    COMPLEX COMPLEX real    real    real    
COMPLEX  |  real    real    real    real    real    real    real    real    
COMPLEX  |  real    real    real    COMPLEX COMPLEX real    real    real    
real     |  real    real    real    COMPLEX COMPLEX real    real    real    
real     |  real    real    real    COMPLEX COMPLEX real    real    real    
real     |  real    real    real    COMPLEX COMPLEX real    real    real    
============================================================
la.eig
12.000+0.000i  |  -0.354+0.000i  0.913+0.000i  0.204+0.000i -0.013+0.016i -0.013-0.016i  0.160+0.000i -0.000+0.000i  0.130+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.183+0.000i  0.208+0.000i  0.379-0.171i  0.379+0.171i -0.607+0.000i  0.000+0.000i -0.138+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.182+0.000i  0.203+0.000i -0.468-0.048i -0.468+0.048i  0.153+0.000i  0.001+0.000i -0.271+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.182+0.000i  0.203+0.000i  0.657+0.000i  0.657-0.000i  0.672+0.000i -0.001+0.000i  0.617+0.000i 
20.000-0.000i  |  -0.354+0.000i -0.182+0.000i  0.203+0.000i -0.276+0.101i -0.276-0.101i -0.361+0.000i  0.001+0.000i -0.644+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.001+0.000i -0.612+0.000i -0.001+0.000i -0.001-0.000i  0.001+0.000i  0.706+0.000i -0.000+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.182+0.000i  0.203+0.000i -0.276+0.101i -0.276-0.101i -0.018+0.000i -0.000+0.000i  0.306+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.001+0.000i -0.612+0.000i -0.001+0.000i -0.001-0.000i  0.001+0.000i -0.708+0.000i  0.000+0.000i 
---------------------
real     |  real    real    real    COMPLEX COMPLEX real    real    real    
real     |  real    real    real    COMPLEX COMPLEX real    real    real    
real     |  real    real    real    COMPLEX COMPLEX real    real    real    
COMPLEX  |  real    real    real    real    real    real    real    real    
COMPLEX  |  real    real    real    COMPLEX COMPLEX real    real    real    
real     |  real    real    real    COMPLEX COMPLEX real    real    real    
real     |  real    real    real    COMPLEX COMPLEX real    real    real    
real     |  real    real    real    COMPLEX COMPLEX real    real    real    
============================================================
np.linalg.eigh
12.000+0.000i  |  -0.354+0.000i  0.000+0.000i  0.000+0.000i -0.086+0.000i  0.905+0.000i -0.025+0.000i  0.073+0.000i  0.205+0.000i 
20.000+0.000i  |  -0.354+0.000i  0.000+0.000i -0.374+0.000i  0.149+0.000i -0.236+0.000i -0.388+0.000i  0.682+0.000i  0.206+0.000i 
20.000+0.000i  |  -0.354+0.000i  0.001+0.000i  0.551+0.000i  0.136+0.000i -0.180+0.000i  0.616+0.000i  0.317+0.000i  0.201+0.000i 
20.000+0.000i  |  -0.354+0.000i  0.001+0.000i -0.149+0.000i  0.719+0.000i -0.074+0.000i -0.042+0.000i -0.534+0.000i  0.207+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.005+0.000i  0.505+0.000i -0.386+0.000i -0.214+0.000i -0.556+0.000i -0.274+0.000i  0.203+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.707+0.000i -0.004+0.000i  0.002+0.000i  0.001+0.000i  0.002+0.000i -0.000+0.000i -0.612+0.000i 
20.000+0.000i  |  -0.354+0.000i  0.003+0.000i -0.529+0.000i -0.535+0.000i -0.203+0.000i  0.398+0.000i -0.262+0.000i  0.203+0.000i 
20.000+0.000i  |  -0.354+0.000i  0.707+0.000i  0.001+0.000i  0.001+0.000i  0.000+0.000i -0.005+0.000i -0.001+0.000i -0.612+0.000i 
---------------------
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
============================================================
la.eigh
12.000+0.000i  |  -0.354+0.000i  0.000+0.000i  0.000+0.000i -0.225+0.000i  0.882+0.000i  0.000+0.000i  0.065+0.000i -0.205+0.000i 
20.000+0.000i  |  -0.354+0.000i  0.000+0.000i -0.395+0.000i  0.332+0.000i -0.156+0.000i  0.227+0.000i  0.701+0.000i -0.205+0.000i 
20.000+0.000i  |  -0.354+0.000i  0.001+0.000i  0.612+0.000i  0.011+0.000i -0.204+0.000i -0.597+0.000i  0.250+0.000i -0.200+0.000i 
20.000+0.000i  |  -0.354+0.000i  0.001+0.000i -0.086+0.000i  0.689+0.000i  0.030+0.000i -0.054+0.000i -0.589+0.000i -0.205+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.005+0.000i  0.413+0.000i -0.264+0.000i -0.245+0.000i  0.711+0.000i -0.165+0.000i -0.205+0.000i 
20.000+0.000i  |  -0.354+0.000i -0.707+0.000i -0.004+0.000i -0.000+0.000i  0.001+0.000i -0.002+0.000i -0.001+0.000i  0.612+0.000i 
20.000+0.000i  |  -0.354+0.000i  0.003+0.000i -0.540+0.000i -0.542+0.000i -0.309+0.000i -0.290+0.000i -0.261+0.000i -0.205+0.000i 
20.000+0.000i  |  -0.354+0.000i  0.707+0.000i  0.001+0.000i -0.000+0.000i  0.001+0.000i  0.005+0.000i -0.001+0.000i  0.612+0.000i 
---------------------
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real    
real     |  real    real    real    real    real    real    real    real 

如您所见，numpy.linalg.eig和scipy.linalg.eig在其输出中生成复数，但它们不应该。这可以被接受为某种舍入误差，如果虚部的大小与实部的大小相比很小。但事实并非如此。产生的数字之一是-0.013+0.016i.在这里，虚部的量级甚至比实部高。

更糟糕的是：这四种方法会产生三种不同的结果。

所有四种方法仅计算一次特征值 12 和 7 乘以特征值 20。所有特征向量的长度始终为 1。这意味着，对于特征值 12，所有四种方法都应生成完全相同的特征向量。但只有numpy.linalg.eig和scipy.linalg.eig产生相同的输出。

以下是特征值 12 的特征向量的分量。仔细查看用箭头(<==)标记的线条。在这里，您可以找到三个不同的值，但这些值应该完全相等。如果你再看一眼，你会发现，只有在^第一行中，所有三个值都是相等的。在所有其他行中，您会发现 2 或 3 个不同的值。

numpy.linalg.eig  |                     |
scipy.linalg.eig  |  numpy.linalg.eigh  |  scipy.linalg.eigh
------------------+---------------------+-------------------
-0.354+0.000i  |      -0.354+0.000i  |      -0.354+0.000i
0.913+0.000i  |       0.000+0.000i  |       0.000+0.000i
0.204+0.000i  |       0.000+0.000i  |       0.000+0.000i
-0.013+0.016i  |      -0.086+0.000i  |      -0.225+0.000i   <===
-0.013-0.016i  |       0.905+0.000i  |       0.882+0.000i   <===
0.160+0.000i  |      -0.025+0.000i  |       0.000+0.000i   <===
-0.000+0.000i  |       0.073+0.000i  |       0.065+0.000i   <===
0.130+0.000i  |       0.205+0.000i  |      -0.205+0.000i

以下是我的问题：

1. 这怎么可能？
2. 这些是错误吗？
3. 其中一个结果是否正确？
4. 如果有一种方法可以提供正确的结果：是哪个？

---

p.s：以下是相关的版本信息：

• 我确实在iMac上运行了此代码(macOS Catalina版本10.15.7)
• python版本是3.8.5
• numpy 的版本是 1.19.5
• scipy 的版本是 1.6.0

---

这是numpy.show_config()的输出 (根据评论中的要求)：

blas_mkl_info:
NOT AVAILABLE
blis_info:
NOT AVAILABLE
openblas_info:
libraries = ['openblas', 'openblas']
library_dirs = ['/usr/local/lib']
language = c
define_macros = [('HAVE_CBLAS', None)]
blas_opt_info:
libraries = ['openblas', 'openblas']
library_dirs = ['/usr/local/lib']
language = c
define_macros = [('HAVE_CBLAS', None)]
lapack_mkl_info:
NOT AVAILABLE
openblas_lapack_info:
libraries = ['openblas', 'openblas']
library_dirs = ['/usr/local/lib']
language = c
define_macros = [('HAVE_CBLAS', None)]
lapack_opt_info:
libraries = ['openblas', 'openblas']
library_dirs = ['/usr/local/lib']
language = c
define_macros = [('HAVE_CBLAS', None)]
None

>附录(在提出问题后 1 天添加)对评论的反应：

• 实对称矩阵的复特征向量
  @Rethipher：谢谢！我确实阅读并理解了您链接到的问题(真正的对称矩阵可以有复杂的特征向量吗？)，我也确实阅读了答案，但我不理解它们。他们说"是"还是"否"？(反问，无需回答，见下一行)
  @Mark Dickinson & @bnaecker：谢谢你明确表示，我的假设是错误的。

• 实对称矩阵与埃尔米特矩阵
  @bnaecker：实数集是复数集的子集。那些等于它们自己的复共轭的复数称为实数。因此，实对称矩阵的集合是埃尔米特矩阵的子集。这很重要，因为numpy.linalg.eigh和scipy.linalg.eigh被设计用于处理埃尔米特矩阵。因为每个实对称矩阵都是埃尔米特矩阵，所以这些模块也可以用于我的目的。

• 混合行和列
  @Mark Dickinson & @bnaecker：谢谢，我认为你是对的。文档也这么说，我应该更仔细地阅读它。但是，即使您比较列而不是行，您仍然会发现这 4 种方法会产生 3 种不同的结果。但是，如果结果包含一个只能用 7 个真实基向量描述的 7 维子空间，我仍然觉得很奇怪，算法会产生复杂的基。

• "一个错误会令人惊讶">
  @bnaecker：这是真的，但令人惊讶的错误确实存在。(比如Heartbleed和其他一些。所以，这不是一个真正的论点。

• "I get reals" -"你的样本矩阵不包含浮点数">
  @Stef & @JohanC：对不起，你没有足够仔细地阅读我的程序。我在A[5,7]和A[7,5]中添加了1e-12值，以模拟在计算特征值和特征向量之前在我的真实应用程序中不可避免地出现的微小舍入误差。(我在这里发布的只是一个很小的测试程序，刚好足以证明这个问题。
  你是对的，Stef：在不添加这种微小噪音的情况下，我也得到了真正的结果。但只有百万分之一的微小变化会产生如此大的差异，我不明白为什么。

对DavidB2013的回答的反应：

我尝试了您建议的工具，得到了不同的结果。我想你也忘了在A[5,7]和A[7,5]中添加1e-12的小噪音.但是，所有结果仍然是真实的。我确实得到了这些特征值：

12.000000000000249
20
20.00000000000075
19.999999999999
20
20
20
20

以及这些特征向量：

0.3535533905932847   0.9128505045937204      0.20252576206455747  0.002673672081814904   -0.09302397289286794   -0.09302397289286794   -0.09302397289286794   -0.09302397289286794     
0.3535533905932848  -0.18259457246238117     0.20444330131542393 -0.00009386949436945406 -0.20415317121194954   -0.20415317121194954   -0.20415317121194954   -0.20415317121194954     
0.3535533905932848  -0.18259457246238117     0.20444330131542393 -0.00009386949436945406 -0.20415317121194954   -0.20415317121194954   -0.20415317121194954    0.9080920678356449     
0.3535533905932848  -0.18259457246238117     0.20444330131542393 -0.00009386949436945406 -0.20415317121194954    0.9080920678356449    -0.20415317121194954   -0.20415317121194954     
0.3535533905932848  -0.18259457246238117     0.20444330131542393 -0.00009386949436945406  0.9080920678356449    -0.20415317121194954   -0.20415317121194954   -0.20415317121194954     
0.35355339059324065 -0.00011103276380548543 -0.6116010247648269   0.7060012169461334      0.0005790869815273477  0.0005790869815273477  0.0005790869815273477  0.0005790869815273477     
0.3535533905932848  -0.18259457246238117     0.20444330131542393 -0.00009386949436945406 -0.20415317121194954   -0.20415317121194954    0.9080920678356449    -0.20415317121194954     
0.35355339059324054  0.0002333904819895115  -0.6131412438770024  -0.7082055415560993      0.0009655029234935232  0.0009655029234935232  0.0009655029234935232  0.0009655029234935232     

只有特征值 12 的向量与您的计算具有相同的值：(6 维中相差约 1.1e-14，其他两个维度相差约 3.3e-14，但我将其视为舍入误差。所有其他向量都存在显著差异(最小差异的大小为 0.02)。令我困惑的是，在输入矩阵的 2 个元素中，1e-12 的微小舍入误差会产生如此大的差异。

---

我用另一种方法(在 https://www.wolframalpha.com 的帮助下)计算了特征值，当我没有添加应该模拟舍入误差的微小增量值时，我只得到两个不同的特征值，即 12 和 20。

给定矩阵的特征多项式为：

(20 - λ)^7 * (12 - λ)

因此，它在 λ=12 处有一个根，在 λ=20 处有 7 个根，这 8 个根是 8 个特征值。所有这些都是实数。

当我添加微小的增量值时，我得到这个特征多项式：

(20 - λ)^5 * (19999999999999/1000000000000 - λ) * (1000000000000 λ^2 - 32000000000001 λ + 240000000000014)/1000000000000

它有这些根源：

λ=12.00000000000024999999999998 (rounded)
λ=19.999999999999 (exact value)
λ=20 (exact value)  
λ=20 (exact value)  
λ=20 (exact value)  
λ=20 (exact value)  
λ=20 (exact value)  
λ=20.00000000000075000000000002 (rounded)

同样，所有 8 个特征值都是实数。

然后我计算了特征向量。在不添加 1e-12 的情况下，我得到以下结果：

特征值 12 的向量：

v = (1,1,1,1,1,1,1,1)

这个向量的长度是 sqrt(8)，如果将向量乘以 1/sqrt(8)，你会得到其他计算的结果(每个维度为 0.35355339)。

但是特征值 20 的七个特征向量非常不同。它们是：

(-1,1,0,0,0,0,0,0)
(-1,0,1,0,0,0,0,0)
(-1,0,0,1,0,0,0,0)
(-1,0,0,0,1,0,0,0)
(-1,0,0,0,0,1,0,0)
(-1,0,0,0,0,0,1,0)
(-1,0,0,0,0,0,0,1)

即使将它们的长度设置为 1，它们也与所有其他结果不同，并且很容易看出它们是正确的。其他结果也是正确的，但我更喜欢这些简单的结果。

我还计算了具有微小噪声的版本的特征值。所有 8 个向量都非常接近无噪声结果，甚至 Wolfram Alpha 也将它们四舍五入到与以前完全相同的值。这正是我对计算特征值和特征向量的算法所期望的行为：

• 输入中的微小变化应该 - 如果可能的话 - 返回结果中的微小变化。

Answer 1

据我所知，假设 1 是正确的，但假设 2 不是。

实对称矩阵生成仅实数的特征值和特征向量。

但是，对于给定的特征值，关联的特征向量不一定是唯一的。

此外，对于实际上不是那么大或包含不是很小的数字的矩阵，舍入误差不应该那么重要。

为了进行比较，我通过 JavaScript 版本的 RG 运行了您的测试矩阵。F(实数通用，来自EISPACK库)：特征值和特征向量计算器

这是输出：

特征：

20
12
20
20
20
20
20
20

特征向量：

0.9354143466934854     0.35355339059327395     -0.021596710639534     -0.021596710639534     -0.021596710639534     -0.021596710639534     -0.021596710639533997     -0.021596710639533997
-0.1336306209562122     0.3535533905932738     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797
-0.1336306209562122     0.3535533905932738     0.9286585574999623     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797
-0.1336306209562122     0.3535533905932738     -0.15117697447673797     0.9286585574999623     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797
-0.1336306209562122     0.3535533905932738     -0.15117697447673797     -0.15117697447673797     0.9286585574999623     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797
-0.1336306209562122     0.3535533905932738     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     0.9286585574999623     -0.15117697447673797     -0.15117697447673797
-0.1336306209562122     0.3535533905932738     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     0.9286585574999622     -0.15117697447673797
-0.1336306209562122     0.3535533905932738     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     -0.15117697447673797     0.9286585574999622

没有虚构的成分。

为了确认或否认结果的有效性，您始终可以编写一个小程序，将结果代入原始方程。简单的矩阵和向量乘法。然后，您将确定输出是否正确。或者，如果他们错了，他们离正确答案有多远。