使用numpy的eigh和svd计算的特征向量不匹配。

16

考虑奇异值分解M=USV*,则M* M的特征值分解为M* M = V(S* S)V*=VS*U*USV*。我希望通过numpy验证这个等式,即通过展示eigh函数返回的特征向量与svd函数返回的相同来验证。

import numpy as np
np.random.seed(42)
# create mean centered data
A=np.random.randn(50,20)
M= A-np.array(A.mean(0),ndmin=2)

# svd
U1,S1,V1=np.linalg.svd(M) 
S1=np.square(S1)
V1=V1.T  

# eig
S2,V2=np.linalg.eigh(np.dot(M.T,M))
indx=np.argsort(S2)[::-1]
S2=S2[indx]
V2=V2[:,indx]

# both Vs are in orthonormal form
assert np.all(np.isclose(np.linalg.norm(V1,axis=1), np.ones(V1.shape[0])))
assert np.all(np.isclose(np.linalg.norm(V1,axis=0), np.ones(V1.shape[1])))
assert np.all(np.isclose(np.linalg.norm(V2,axis=1), np.ones(V2.shape[0])))
assert np.all(np.isclose(np.linalg.norm(V2,axis=0), np.ones(V2.shape[1])))

assert np.all(np.isclose(S1,S2))
assert np.all(np.isclose(V1,V2))

最后一个断言失败了。为什么?


你可以将一个正数加到所有对角线元素上,即令M2=M+a*I,其中a足够大使得M2是半正定的。然后进行奇异值分解和特征值分解应该会更一致。 - Shaohua Li
1个回答

14

只需使用小的数值来调试您的问题。

A=np.random.randn(3,2)开始,而不是您的大小为(50,20)的大矩阵。

在我的随机情况中,我发现

v1 = array([[-0.33872745,  0.94088454],
   [-0.94088454, -0.33872745]])

并且对于v2

v2 = array([[ 0.33872745, -0.94088454],
   [ 0.94088454,  0.33872745]])

它们只相差一个符号,即使将向量归一化为单位模,向量也可能是相差一个符号的。

现在,如果您尝试这个技巧

assert np.all(np.isclose(V1,-1*V2))

对于您原来的大矩阵,它失败了…… 再次失败也没关系。发生的情况是一些向量已经被乘以-1,而有些向量则没有。

检查向量之间相等的正确方法是:

assert allclose(abs((V1*V2).sum(0)),1.)

确实,为了感受这个过程的工作原理,您可以打印出这个量:

(V1*V2).sum(0)

这确实取决于向量是+1还是-1

array([ 1., -1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,
    1., -1.,  1.,  1.,  1., -1., -1.])

编辑:大多数情况下会发生这种情况,特别是从随机矩阵开始。但要注意,如果一个或多个特征值具有维数大于1的特征空间,则此测试可能会失败,如@Sven Marnach在下面的评论中指出:

除了向量乘以-1之外,可能还会有其他差异。 如果任何一个特征值具有多维特征空间,则您可能会获得该特征空间的任意正交基础,并且这样的基础可以通过任意酉矩阵相互旋转。


@matus 好的,我有点迷惑了 :) 但我相信你的判断,所以我会删除我的评论,以免混淆未来的读者。干杯! - BartoszKP
可能存在除向量乘以-1之外的其他差异。 如果任何特征值具有多维特征空间,则可以获得该特征空间的任意正交基,且这些基可以通过任意酉矩阵相互旋转。 - Sven Marnach
@SvenMarnach,这是一个非常有价值的观点。我会编辑帖子以指出这个警告。 - gg349

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接