用 numpy 的 eigh 和 svd 计算的特征向量不匹配
Eigenvectors computed with numpy's eigh and svd do not match
考虑奇异值分解M=USV*。然后M* M的特征值分解得到M* M= V (S* S) V*=VS* U* USV*。我希望通过显示 eigh
函数返回的特征向量与 svd
函数返回的特征向量相同来验证与 numpy 的相等性:
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))
最后一个断言失败。为什么?
只需使用小数字来调试您的问题。
从 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 在他下面的评论中指出的那样:
There might be other differences than just vectors multiplied by -1.
If any of the eigenvalues has a multi-dimensional eigenspace, you
might get an arbitrary orthonormal basis of that eigenspace, and to
such bases might be rotated against each other by an arbitraty
unitarian matrix
考虑奇异值分解M=USV*。然后M* M的特征值分解得到M* M= V (S* S) V*=VS* U* USV*。我希望通过显示 eigh
函数返回的特征向量与 svd
函数返回的特征向量相同来验证与 numpy 的相等性:
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))
最后一个断言失败。为什么?
只需使用小数字来调试您的问题。
从 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 在他下面的评论中指出的那样:
There might be other differences than just vectors multiplied by -1. If any of the eigenvalues has a multi-dimensional eigenspace, you might get an arbitrary orthonormal basis of that eigenspace, and to such bases might be rotated against each other by an arbitraty unitarian matrix