C语言和Python中特征向量计算程序的结果不同

Question

我注意到对于一个全1的4x4矩阵，使用numpy.linalg.eig在Python中得出了不同的特征分解结果。

matrix = numpy.ones((M,M), dtype=float);
values, vectors = numpy.linalg.eig(matrix);

Python 结果：

V1: [-0.866025 +0.288675 +0.288675 +0.288675]
V2: [+0.500000 +0.500000 +0.500000 +0.500000]
V3: [+0.391955 +0.597433 -0.494694 -0.494694]
V4: [+0.866025 -0.288675 -0.288675 -0.288675]

使用LAPACK DSYEV在C语言中：

#define NN 4
#define LDA NN
void main(){
    int n = NN, lda = LDA, lwork=NN*NN*NN*NN*NN, info;
    char both = 'V';
    char uplo = 'U';
    double w[NN*NN];
    double work[NN*NN*NN*NN*NN];
    double a[LDA*NN] = {
       1,  1,  1,  1,
       1,  1,  1,  1,
       1,  1,  1,  1,
       1,  1,  1,  1 
    };

    dsyev_(&both, &uplo, &n, a, &lda, w, work, &lwork, &info);
    return;
}

C DSYEV Result:

V1: +0.000596 +0.000596 -0.707702 +0.706510 
V2: +0.500000 +0.500000 -0.499157 -0.500842 
V3: +0.707107 -0.707107 -0.000000 +0.000000 
V4: +0.500000 +0.500000 +0.500000 +0.500000  

在C语言中使用LAPACK DGEEV:

#define NN 4
#define LDA NN
#define LDVL NN
#define LDVR NN
void main() {
    char compute_left = 'V';
    char compute_right = 'V';
    int n = NN, lda = LDA, ldvl = LDVL, ldvr = LDVR, info, lwork=2*NN*NN;
    double work[2*NN*NN];
    double wr[NN], wi[NN], vl[LDVL*NN], vr[LDVR*NN];
    double a[LDA*NN] = {
       1,  1,  1,  1,
       1,  1,  1,  1,
       1,  1,  1,  1,
       1,  1,  1,  1 
    };
    dgeev_( &compute_left, &compute_right, &n, a, &lda, wr, wi, vl, &ldvl, vr, &ldvr, work, &lwork, &info );
    return;
}

C DGEEV 结果：

V1: -0.866025 +0.288675 +0.288675 +0.288675 
V2: -0.500000 -0.500000 -0.500000 -0.500000 
V3: -0.000000 -0.816497 +0.408248 +0.408248 
V4: -0.000000 -0.000000 -0.707107 +0.707107 

结果都不一样！

所以我有两个主要问题:

1. 为什么？这是由于1矩阵中的退化吗？
2. 如何复制Python在C中的结果？

任何见解都将不胜感激。

- The Dude

---

结果似乎偏差太大，这可能是由于库的数值不稳定性造成的。也许它期望某些矩阵已经被初始化为零或其他值？ - Paul Ogilvie

@PaulOgilvie 我刚试过事先将所有相关数组初始化为零，但我仍然得到相同的答案。 - The Dude

Answer 1

所有都是正确的。您的矩阵有两个特征值，4和0。对于4的特征空间是由[1,1,1,1]张成的线性空间，所有列表中都会出现它的倍数。对于0的特征空间是三维空间x_1 + x_2 + x_3 + x_4 = 0。每种方法都为这个子空间提供了不同的基础，除了numpy，它只给出了张成二维子空间的向量，原因不明。

在我看来，DGEEV的结果是您报告的最好的结果，因为它以合理的阶梯形式给出了0特征空间的正交基。

- Nick Matteo

---

我想这大概就是这样。有没有办法让我在C中得到与numpy相同的基础呢？ - The Dude

@TheDude：据我所知没有。有无限多种可能的基础，实现选择任何一个的特定原因。 - Nick Matteo

@TheDude：实际上，注意到numpy给出的结果并不是一个基础。在你给出的结果中，V4只是负的V1。在我的系统上，它重复了一个向量两次，甚至没有符号变化。 - Nick Matteo

Answer 2

四个特征值中有三个是0（在你的Python脚本中尝试打印values）。由于矩阵的所有元素都相同（为1），任何元素之和为零的向量都将是对应于零特征值的有效特征向量。选择这个向量的方式并不重要，因此不同的软件找到零特征值的不同特征向量并不重要。您应该确认这些特征向量确实具有元素之和为0。

最后一个特征值是4（非零），这意味着相应的特征向量必须具有相同的（非零）元素。这些元素的确切值取决于特征向量的归一化，这在您的示例中似乎也有所不同。

总的来说，一切都还好，只是你的矩阵特征向量非常不唯一。另一个解决方案，我认为更令人满意的是在Wolfram Alpha中找到的。

更新

通常，一个M乘以M的矩阵有M个特征值。如果其中两个（或更多）相同，则存在无限数量的可能实现相应特征向量，因为从这两个（称为v1和v2）中，我们可以通过v3 = a*v1 + b*v2构造出一个新的特征向量，其中a和b是任意常数。

- jmd_dk

---

确实，这就是答案。但是我如何确保在不同平台上获得相同的基础？或者除了深入实现的细节之外，没有其他方法吗？ - The Dude

你可以在获得特征向量后按照自己的需求进行归一化。这样，无论实现方式如何，你都可以获得相同特征值4的特征向量。我不知道有什么方法可以使另外三个特征向量无论实现方式如何都保持一致。 - jmd_dk