【发布时间】:2021-04-11 14:26:07
【问题描述】:
这里有两个关于方阵特征向量和特征值的假设。我相信两者都是正确的:
-
如果一个矩阵是对称的并且只包含实数,那么它是一个 Hermitian 矩阵,那么所有特征值都应该是实数并且所有特征向量的所有分量也应该是实数。从 Hermitian 矩阵计算特征向量和特征值时,结果中不应出现复数。
-
从给定矩阵计算的给定特征值的特征向量应始终指向仅由矩阵和特征值确定的方向。用于计算它的算法对结果没有影响,只要算法实现正确即可。
但是,当您在 Python 中使用标准库计算特征向量和特征值时,这两个假设都不成立。这些方法是否包含错误?
有四种不同的方法可以从 Hermitian 矩阵计算特征值和特征向量:
#1 和#2 可用于任何方阵(包括 Hermitian 矩阵)。
#3 和 #4 仅适用于 Hermitian 矩阵。据我所知,他们的目的只是为了让他们跑得更快,但结果应该是一样的(只要输入真的是 Hermitian)。
但是,这四种方法会为相同的输入提供三种不同的结果。这是我用来测试所有四种方法的程序:
#!/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 和特征值 20 的 7 倍。所有特征向量的长度始终为 1。这意味着,所有四种方法都应该为特征值 12 产生完全相同的特征向量。但只有 @987654334 @ 和 scipy.linalg.eig 产生相同的输出。
这里是特征值 12 的特征向量的分量。仔细查看标有箭头 (<==) 的线。在这里您可以找到三个不同的值,但这些值应该完全相同。如果您再看一眼,您会发现,只有在 1st 行中,所有三个值都是相等的。在所有其他行中,您会发现 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
这是我的问题:
- 这怎么可能?
- 这些是错误吗?
- 其中一项结果是否正确?
- 如果有一种方法可以提供正确的结果:它是哪一种?
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 天后添加)
对 cmets 的反应:
-
实对称矩阵的复特征向量
@Rethither:谢谢!我确实阅读并理解了the question you linked to(一个实数对称矩阵可以有复杂的特征向量吗?),我也阅读了答案,但我不明白。他们说“是”还是“不是”? (修辞问题,无需回答,见下一行)
@Mark Dickinson 和@bnaecker:感谢您明确指出我的假设是错误的。 -
实对称矩阵与 Hermitian 矩阵
@bnaecker:实数集是复数集的子集。那些等于它们自己的复共轭的复数称为实数。因此,实对称矩阵的集合是 Hermitian 矩阵的子集。这很重要,因为 numpy.linalg.eigh 和 scipy.linalg.eigh 旨在处理 Hermitian 矩阵。而且因为每个实对称矩阵都是 Hermitian 矩阵,所以这些模块也可以用于我的目的。 -
混合行和列
@Mark Dickinson 和@bnaecker:谢谢,我认为你是对的。文件也这么说,我应该更仔细地阅读它。但即使您比较列而不是行,您仍然会发现这 4 种方法会产生 3 种不同的结果。但如果结果包含一个 7 维子空间,只能用 7 个实基向量来描述,我还是觉得奇怪,算法会产生一个复基。 -
“一个错误会令人惊讶”
@bnaecker:这是真的,但确实存在令人惊讶的错误。 (比如Heartbleed 和其他一些人。)所以,这并不是一个真正的论点。 -
“我得到实数” - “你的样本矩阵不包含浮点数”
@Stef & @JohanC:抱歉,你没有仔细阅读我的程序。我在A[5,7]和A[7,5]中添加了1e-12的值,以模拟在计算特征值和特征向量之前在我的实际应用程序中不可避免地出现的微小舍入误差。 (我在这里发布的只是一个很小的测试程序,足以证明问题。)
你是对的,Stef:不添加这种微小的噪音,我也得到了真正的结果。但是只有百万分之一的微小变化会产生如此大的差异,我不明白为什么。
对 DavidB2013 回答的反应:
我尝试了您建议的工具,但得到了不同的结果。我想你也忘了将1e-12 的小噪音添加到A[5,7] 和A[7,5]。但是,所有结果仍然是真实的。我确实得到了这些特征值:
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 的帮助下)计算了特征值,当我没有添加应该模拟舍入误差的微小 delta 值时,我只得到两个不同的特征值,即 12 和 20。
给定矩阵的特征多项式为:
(20 - λ)^7 * (12 - λ)
因此,它在 λ=12 处有 1 个根,在 λ=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 也将它们四舍五入到与以前完全相同的值。这正是我对计算特征值和特征向量的算法所期望的行为:
- 输入中的微小变化应该(尽可能)返回结果中的微小变化。
【问题讨论】:
-
找到了一个类似于你所问的答案,但它更数学here。我会小心特征值/特征向量计算。浮点数和数值问题通常会潜入并给出与手工计算的特征向量/特征值对不同的答案。另请注意,您的答案显示为 0.0i。尝试格式化数据以显示更多有效数字,您可能会看到一个表示舍入错误的小数字。
-
只是在你的第一个假设中关于术语的注释。 Hermitian 矩阵根据定义包含复杂元素。如果所有值都是实数,则称为对称。
-
你的假设(2)是错误的。虽然由矩阵表示的真正数学对象具有唯一的特征分解,但实际的机器表示的二维数组却没有。这些算法是迭代的,并使用有限精度浮点数。这些引入的小错误可以减少或纠正,但永远不会完全消除。
-
不要堆积,但假设(2)在数学上也是错误的:如果特定特征值的特征空间的维数大于 1,则特征向量不再确定为标量倍数。 (这已经发生在原始示例中:从数学上讲,矩阵
A的特征值恰好为20,特征空间为五维,以及两个不同的特征值,它们接近但不完全等于20,并且接近(但又不完全等于)12的特征值。) -
@JohanC 他将
1e-12添加到A[5,7]和A[7,5]
标签: python numpy scipy linear-algebra eigenvector