Python 中的特征值:一个错误?

问题描述

这里有两个关于方阵的特征向量和特征值的假设。我相信两者都是正确的:

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

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

但是当您在 Python 中使用标准库来计算特征向量和特征值时,这两个假设都不成立。这些方法是否包含错误

有四种不同的方法可以从 hermitian 矩阵计算特征值和特征向量:

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

#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,19,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.eigscipy.linalg.eig 在其输出生成复数,但它们不应该。如果虚部的大小与实部的大小相比微不足道,则可以将其视为某种舍入误差。但这种情况并非如此。产生的数字之一是-0.013+0.016i。这里虚部的量级甚至高于实部。

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

所有四种方法只计算一次 12 的特征值和 7 倍的特征值 20。并且所有特征向量的长度总是为 1。这意味着,所有四种方法都应该为特征值 12 产生完全相同的特征向量。但只有 {{ 1}} 和 numpy.linalg.eig 产生相同的输出

这里是特征值 12 的特征向量的组成部分。仔细查看标有箭头 (scipy.linalg.eig) 的线。在这里您可以找到三个不同的值,但这些值应该完全相等。如果你再看一眼,你会发现,只有在 1st 行中,三个值都是相等的。在所有其他行中,您会发现 2 或 3 个不同的值。

<==

这是我的问题:

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

ps:以下是相关版本信息:

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

这是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 输出 (根据评论中的要求):

numpy.show_config()

附录(在提出问题后 1 天添加

评论的反应:

  • 实对称矩阵的复特征向量
    @Rethipher:谢谢!我确实阅读并理解了the question you linked to实对称矩阵可以有复数特征向量吗?),我也阅读了答案,但我没有理解它们。他们说“是”还是“不是”? (反问,不用回答,见下一行)
    @Mark Dickinson 和 @bnaecker:感谢您明确指出我的假设是错误的。

  • 实对称矩阵与厄米矩阵
    @bnaecker:实数集是复数集的子集。那些等于它们自己的复共轭的复数称为实数。因此,实对称矩阵集是厄米矩阵的子集。这很重要,因为 numpy.linalg.eigh 和 scipy.linalg.eigh 旨在处理 hermitian 矩阵。因为每个实对称矩阵都是厄米矩阵,所以这些模块也可以用于我的目的。

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

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

  • “我得到实数” - “您的样本矩阵不包含浮点数”
    @Stef & @JohanC:抱歉,你没有仔细阅读我的程序。我将 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',None)] lapack_mkl_info: NOT AVAILABLE openblas_lapack_info: libraries = ['openblas',None)] lapack_opt_info: libraries = ['openblas',None)] None 的值添加1e-12A[5,7] 以模拟在计算特征值和特征向量之前不可避免地出现在我的真实应用程序中的微小舍入误差。 (我在这里发布的只是一个很小的测试程序,大到足以证明问题。)
    你是对的,Stef:不添加这个微小的噪音,我也得到了真实的结果。但只有百万分之一的微小变化会产生如此大的差异,我不明白为什么。

对 DavidB2013 回答的反应:

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

A[7,5]

和这些特征向量:

12.000000000000249
20
20.00000000000075
19.999999999999
20
20
20
20

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


我用另一种方法(在 https://www.wolframalpha.com 的帮助下)计算了特征值,当我没有添加应该模拟舍入误差的微小增量值时,我只得到两个不同的特征值,分别为 12 和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 处有一个根,在 λ=20 处有 7 个根,这 8 个根是 8 个特征值。都是实数。

当我添加微小的 delta 值时,我得到了这个特征多项式:

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

它有以下根源:

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

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

然后我计算了特征向量。不添加 1e-12 我得到这个结果:

特征值 12 的向量:

λ=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)

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

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

v = (1,1,1)

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

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

  • 输入的微小变化应该 - 只要有可能 - 返回结果的微小变化。

解决方法

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

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

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

此外,对于实际上不是那么大或包含不是非常小的数字的矩阵,舍入误差不应该如此显着。

为了进行比较,我通过 JavaScript 版本的 RG.F(Real General,来自 EISPACK 库)运行了您的测试矩阵:Eigenvalues and Eigenvectors Calculator

输出如下:

特征值:

   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

无虚部。

要确认或否认结果的有效性,您总是可以编写一个小程序,将结果插入到原始方程中。简单的矩阵和向量乘法。然后你肯定知道输出是否正确。或者,如果他们错了,他们离正确答案还有多远。