为什么我的 PCA 对旋转和轴交换不是不变的？
Slo*_*oke 5 python math numpy voxel pca
我有一个大小为 3x3x3 的体素（np.array），填充了一些值，这个设置对我来说是必不可少的。我想要它的旋转不变表示。对于这种情况，我决定尝试 PCA 表示，它被认为对正交变换是不变的。其他
为简单起见，我进行了一些轴交换，但万一我弄错了，可以有np.rot90.
我已经将我的 3d 体素解释为一组加权的 3d 立方体点向量，我错误地将其称为“基础”，总共 27 个（所以这是空间中的一些 3d 点，由向量表示，从立方体点获得，由体素缩放值）。
import numpy as np
voxel1 = np.random.normal(size=(3,3,3))
voxel2 = np.transpose(voxel1, (1,0,2)) #np.rot90(voxel1) #
basis = []
for i in range(3):
for j in range(3):
for k in range(3):
basis.append([i+1, j+1, k+1]) # avoid 0
basis = np.array(basis)
voxel1 = voxel1.reshape((27,1))
voxel2 = voxel2.reshape((27,1))
voxel1 = voxel1*basis # weighted basis vectors
voxel2 = voxel2*basis
print(voxel1.shape)
(27, 3)
然后我对那 27 个 3 维向量进行了 PCA：
def pca(x):
center = np.mean(x, 0)
x = x - center
cov = np.cov(x.T) / x.shape[0]
e_values, e_vectors = np.linalg.eig(cov)
order = np.argsort(e_values)
v = e_vectors[:, order].transpose()
return x.dot(v)
vp1 = pca(voxel1)
vp2 = pca(voxel2)
但结果在vp1和vp2是不同的。也许，我有一个错误（虽然我相信这是正确的公式），正确的代码必须是
x.dot(v.T)
但在这种情况下，结果很奇怪。转换后的数据的上下块相同，直到符号：
>>> np.abs(np.abs(vp1)-np.abs(vp2)) > 0.01
array([[False, False, False],
[False, False, False],
[False, False, False],
[ True, True, True],
[ True, True, True],
[ True, True, True],
[ True, True, True],
[ True, True, True],
[ True, False, True],
[ True, True, True],
[ True, True, True],
[ True, True, True],
[False, False, False],
[False, False, False],
[False, False, False],
[ True, True, True],
[ True, True, True],
[ True, True, True],
[ True, True, True],
[ True, True, True],
[ True, False, True],
[ True, True, True],
[ True, True, True],
[ True, True, True],
[False, False, False],
[False, False, False],
[False, False, False]])
我做错了什么？
我想要做的是找到我的加权体素的一些不变表示，比如根据惯性轴或主轴定位。如果有人帮助我，我将不胜感激。
UPD：发现与我的问题类似，但代码不可用
EDIT2：找到代码InertiaRotate并设法猴子执行以下操作：
import numpy as np
# https://github.com/smparker/orient-molecule/blob/master/orient.py
voxel1 = np.random.normal(size=(3,3,3))
voxel2 = np.transpose(voxel1, (1,0,2))
voxel1 = voxel1.reshape((27,))
voxel2 = voxel2.reshape((27,))
basis = []
for i in range(3):
for j in range(3):
for k in range(3):
basis.append([i+1, j+1, k+1]) # avoid 0
basis = np.array(basis)
basis = basis - np.mean(basis, axis=0)
def rotate_func(data, mass):
#mass = [ masses[n.lower()] for n in geom.names ]
inertial_tensor = -np.einsum("ax,a,ay->xy", data, mass, data)
# negate sign to reverse the sorting of the tensor
eig, axes = np.linalg.eigh(-inertial_tensor)
axes = axes.T
# adjust sign of axes so third moment moment is positive new in X, and Y axes
testcoords = np.dot(data, axes.T) # a little wasteful, but fine for now
thirdmoment = np.einsum("ax,a->x", testcoords**3, mass)
for i in range(2):
if thirdmoment[i] < 1.0e-6:
axes[i,:] *= -1.0
# rotation matrix must have determinant of 1
if np.linalg.det(axes) < 0.0:
axes[2,:] *= -1.0
return axes
axes1 = rotate_func(basis, voxel1)
v1 = np.dot(basis, axes1.T)
axes2 = rotate_func(basis, voxel2)
v2 = np.dot(basis, axes2.T)
print(v1)
print(v2)
它似乎分别使用基础（坐标）和质量。结果与我上面的问题非常相似：转换后的数据的某些部分与符号匹配，我相信这些是一些立方体的边
print(np.abs(np.abs(v1)-np.abs(v2)) > 0.01)
[[False False False]
[False False False]
[False False False]
[ True True True]
[ True True True]
[ True True True]
[ True True True]
[False False False]
[ True True True]
[ True True True]
[ True True True]
[ True True True]
[False False False]
[False False False]
[False False False]
[ True True True]
[ True True True]
[ True True True]
[ True True True]
[False False False]
[ True True True]
[ True True True]
[ True True True]
[ True True True]
[False False False]
[False False False]
[False False False]]
寻找一些解释。这段代码是为分子设计的，必须工作......
UPD：尝试从这 24 个向量中选择 3 个向量作为新的基础 - 具有最大范数的向量，具有最小范数的向量及其叉积。将它们组合到矩阵 V 中，然后使用公式 V^(-1)*X 来变换坐标，并得到同样的问题 - 旋转体素的结果向量集不相等。
UPD2：我同意 meTchaikovsky 的观点，即我将体素向量乘以权重从而创建一些非三次点云的想法是不正确的。可能，我们确实需要采用旋转“基础”的解决方案（是的，这不是基础，而是一种确定点云的方法），当“基础”相同时，稍后将起作用，但权重根据旋转到 3D 旋转。
根据 meTchaikovsky 提供的答案和参考资料，并找到其他答案，我们和我的朋友一起得出结论，rotate_func从上面提到的分子包中，试图发明一些计算分量符号的约定。他们的解决方案尝试对前 2 个轴使用第三个矩，对最后一个轴 (?) 使用行列式。我们尝试了另一种方法并成功匹配了一半的表示：
# -*- coding: utf-8 -*-
"""
Created on Fri Oct 16 11:40:30 2020
@author: Dima
"""
import numpy as np
from numpy.random import randn
from numpy import linalg as la
from scipy.spatial.transform import Rotation as R
np.random.seed(10)
rotate_mat = lambda theta: np.array([[np.cos(theta),-np.sin(theta),0.],[np.sin(theta),np.cos(theta),0.],[0.,0.,1.]])
def pca(feat, x):
# pca with attemt to create convention on sign changes
x_c =x- np.mean(x,axis=0)
x_f= feat*x
x_f-= np.mean(x_f, axis=0)
cov = np.cov(x_f.T)
e_values, e_vectors = np.linalg.eig(cov)
order = np.argsort(e_values)[::-1]
#print(order)
v = e_vectors[:,order]
v= v/np.sign(v[0,:])
if(la.det(v)<0):
v= -v
return x_c @ v
def standardize(x):
# take vector with biggest norm, with smallest and thir cross product as basis
x -= np.mean(x,axis=0)
nrms= la.norm(x, axis=1)
imin= argmin(nrms)
imax= argmax(nrms)
vec1= x[imin, :]
vec2= x[imax, :]
vec3= np.cross(vec1, vec2)
Smat= np.stack([vec1, vec2, vec3], axis=0)
if(la.det(Smat)<0):
Smat= -Smat
return(la.inv(Smat)@x.T)
angles = np.linspace(0.0,90.0,91)
voxel1 = np.random.normal(size=(3,3,3))
res = []
for angle in angles:
voxel2 = voxel1.copy()
voxel1 = voxel1.reshape(27,1)
voxel2 = voxel2.reshape(27,1)
basis1 = np.array([[i+1,j+1,k+1] for k in range(3) for j in range(3) for i in range(3)]).astype(np.double)
basis1 = basis1+1e-4*randn(27,3) # perturbation
basis2 = basis1 @rotate_mat(np.deg2rad(angle))
#voxel1 = voxel1*basis1
#voxel2 = voxel2*basis2
#print(angle,(np.abs(pca(voxel1) - pca(voxel2) )))
#gg= np.abs(standardize(basis1) - standardize(basis2) )
gg= np.abs(pca(voxel1, basis1) - pca(voxel1, basis2) )
ss= np.sum(np.ravel(gg))
bl= np.all(gg<1e-4)
print(angle,ss, bl)
#res.append(np.all(np.abs(pca(voxel1) - pca(voxel2) < 1e-6)))
del basis1, basis2
结果在 58 度角范围内都很好（但我们仍在试验 x、y 轴的旋转）。在那之后，我们有恒定的差异，这表明一些未计数的符号反转。这比不太一致的结果要好rotate_func：
0.0 0.0 True
1.0 1.1103280567106161e-13 True
2.0 5.150139890290964e-14 True
3.0 8.977126225544196e-14 True
4.0 5.57341699240722e-14 True
5.0 4.205149954378956e-14 True
6.0 3.7435437643664957e-14 True
7.0 1.2943967187158123e-13 True
8.0 5.400185371573149e-14 True
9.0 8.006410204958181e-14 True
10.0 7.777189536904011e-14 True
11.0 5.992073021576436e-14 True
12.0 6.3716122222085e-14 True
13.0 1.0120048110065158e-13 True
14.0 1.4193029076233626e-13 True
15.0 5.32774440341853e-14 True
16.0 4.056702432878251e-14 True
17.0 6.52062429116855e-14 True
18.0 1.3237663595853556e-13 True
19.0 8.950259695710006e-14 True
20.0 1.3795067925438317e-13 True
21.0 7.498727794307339e-14 True
22.0 8.570866862371226e-14 True
23.0 8.961510590826412e-14 True
24.0 1.1839169916779899e-13 True
25.0 1.422193407555868e-13 True
26.0 6.578778015788652e-14 True
27.0 1.0042963537887101e-13 True
28.0 8.438153062569065e-14 True
29.0 1.1299103064863272e-13 True
30.0 8.192453876745831e-14 True
31.0 1.2618492405483406e-13 True
32.0 4.9237819394886296e-14 True
33.0 1.0971028569666842e-13 True
34.0 1.332138304559801e-13 True
35.0 5.280024600049296e-14 True
从上面的代码可以看出，我们尝试使用另一个基：具有最大范数的向量，具有最小范数的向量及其叉积。在这里，我们应该只有两个变体（交叉产品的方向）可以稍后修复，但我无法管理此替代解决方案的工作。
我希望有人能帮我完成这个并获得体素的旋转不变表示。
编辑 3. 非常感谢 meTchaikovsky，但情况仍不清楚。我的问题最初在于处理 3d 体素，它们是 (3,3,3) numpy 数组。我们得出的结论是，为了找到不变的表示，我们只需要将 3d 体素固定为用于计算 cov 矩阵的权重，并在中心“基础”（一些用于描述点云的向量）上应用旋转。
因此，当我们实现“基”旋转的不变性时，问题应该已经解决了：现在，当我们固定“基”并使用旋转体素时，结果必须是不变的。令人惊讶的是，事实并非如此。在这里，我使用 base2=basis1 检查立方体的 24 次旋转（小扰动除外）：
import scipy.ndimage
def pca(feat, x):
# pca with attemt to create convention on sign changes
x_c = x - np.mean(x,axis=0)
x_f = feat * x
x_f -= np.mean(x_f,axis=0)
cov = np.cov(x_f.T)
e_values, e_vectors = np.linalg.eig(cov)
order = np.argsort(e_values)[::-1]
v = e_vectors[:,order]
# here is the solution, we switch the sign of the projections
# so that the projection with the largest absolute value along a principal axis is positive
proj = x_c @ v
asign = np.sign(proj)
max_ind = np.argmax(np.abs(proj),axis=0)[None,:]
sign = np.take_along_axis(asign,max_ind,axis=0)
proj = proj * sign
return proj
def rotate_3d(image1, alpha, beta, gamma):
# z
# The rotation angle in degrees.
image2 = scipy.ndimage.rotate(image1, alpha, mode='nearest', axes=(0, 1), reshape=False)
# rotate along y-axis
image3 = scipy.ndimage.rotate(image2, beta, mode='nearest', axes=(0, 2), reshape=False)
# rotate along x-axis
image4 = scipy.ndimage.rotate(image3, gamma, mode='nearest', axes=(1, 2), reshape=False)
return image4
voxel10 = np.random.normal(size=(3,3,3))
angles = [[x,y,z] for x in [-90,0,90] for y in [-90,0,90] for z in [-90,0,90]]
res = []
for angle in angles:
voxel2 = rotate_3d(voxel10, angle[0], angle[1], angle[2])
voxel1 = voxel10.reshape(27,1)
voxel2 = voxel2.reshape(27,1)
basis1 = np.array([[i+1,j+1,k+1] for k in range(3) for j in range(3) for i in range(3)]).astype(np.double)
basis1 += 1e-4*np.random.normal(size=(27, 1)) # perturbation
basis2 = basis1
original_diff = np.sum(np.abs(basis1-basis2))
gg= np.abs(pca(voxel1, basis1) - pca(voxel2, basis2))
ss= np.sum(np.ravel(gg))
bl= np.all(gg<1e-4)
print('difference before pca %.3f,' % original_diff, 'difference after pca %.3f' % ss,bl)
res.append(bl)
del basis1, basis2
print('correct for %.1f percent of time' % (100*(np.sum(res) / len(res))))
difference before pca 0.000, difference after pca 45.738 False
difference before pca 0.000, difference after pca 12.157 False
difference before pca 0.000, difference after pca 26.257 False
difference before pca 0.000, difference after pca 37.128 False
difference before pca 0.000, difference after pca 52.131 False
difference before pca 0.000, difference after pca 45.436 False
difference before pca 0.000, difference after pca 42.226 False
difference before pca 0.000, difference after pca 18.959 False
difference before pca 0.000, difference after pca 38.888 False
difference before pca 0.000, difference after pca 12.157 False
difference before pca 0.000, difference after pca 26.257 False
difference before pca 0.000, difference after pca 50.613 False
difference before pca 0.000, difference after pca 52.132 False
difference before pca 0.000, difference after pca 0.000 True
difference before pca 0.000, difference after pca 52.299 False
这里 base1=basis2（因此 pca=0 之前的基差），并且您可以看到 (0,0,0) 旋转为 0。但旋转体素给出不同的结果。如果 scipy 做错了什么，我已经用numpy.rot90检查了该方法，结果相同：
rot90 = np.rot90
def rotations24(polycube):
# imagine shape is pointing in axis 0 (up)
# 4 rotations about axis 0
yield from rotations4(polycube, 0)
# rotate 180 about axis 1, now shape is pointing down in axis 0
# 4 rotations about axis 0
yield from rotations4(rot90(polycube, 2, axis=1), 0)
# rotate 90 or 270 about axis 1, now shape is pointing in axis 2
# 8 rotations about axis 2
yield from rotations4(rot90(polycube, axis=1), 2)
yield from rotations4(rot90(polycube, -1, axis=1), 2)
# rotate about axis 2, now shape is pointing in axis 1
# 8 rotations about axis 1
yield from rotations4(rot90(polycube, axis=2), 1)
yield from rotations4(rot90(polycube, -1, axis=2), 1)
def rotations4(polycube, axis):
"""List the four rotations of the given cube about the given axis."""
for i in range(4):
yield rot90(polycube, i, axis)
def rot90(m, k=1, axis=2):
"""Rotate an array k*90 degrees in the counter-clockwise direction around the given axis"""
m = np.swapaxes(m, 2, axis)
m = np.rot90(m, k)
m = np.swapaxes(m, 2, axis)
return m
voxel10 = np.random.normal(size=(3,3,3))
gen = rotations24(voxel10)
res = []
for voxel2 in gen:
#voxel2 = rotate_3d(voxel10, angle[0], angle[1], angle[2])
voxel1 = voxel10.reshape(27,1)
voxel2 = voxel2.reshape(27,1)
basis1 = np.array([[i+1,j+1,k+1] for k in range(3) for j in range(3) for i in range(3)]).astype(np.double)
basis1 += 1e-4*np.random.normal(size=(27, 1)) # perturbation
basis2 = basis1
original_diff = np.sum(np.abs(basis1-basis2))
gg= np.abs(pca(voxel1, basis1) - pca(voxel2, basis2))
ss= np.sum(np.ravel(gg))
bl= np.all(gg<1e-4)
print('difference before pca %.3f,' % original_diff, 'difference after pca %.3f' % ss,bl)
res.append(bl)
del basis1, basis2
print('correct for %.1f percent of time' % (100*(np.sum(res) / len(res))))
我试图调查这个案例，唯一可能不相关的事情是我发现的：
voxel1 = np.ones((3,3,3))
voxel1[0,0,0] = 0 # if I change 0 to 0.5 it stops working at all
# mirrored around diagonal
voxel2 = np.ones((3,3,3))
voxel2[2,2,2] = 0
for angle in range(1):
voxel1 = voxel1.reshape(27,1)
voxel2 = voxel2.reshape(27,1)
basis1 = np.array([[i+1,j+1,k+1] for k in range(3) for j in range(3) for i in range(3)]).astype(np.double)
basis1 = basis1 + 1e-4 * randn(27,3) # perturbation
basis2 = basis1
# If perturbation is used we have
# difference before pca 0.000, difference after pca 0.000 True
# correct for 100.0 percent of time
# eigenvalues for both voxels
# [1.03417495 0.69231107 0.69235402]
# [0.99995368 0.69231107 0.69235402]
# If no perturbation applied for basis, difference is present
# difference before pca 0.000, difference after pca 55.218 False
# correct for 0.0 percent of time
# eignevalues for both voxels (always have 1.):
# [0.69230769 1.03418803 0.69230769]
# [1. 0.69230769 0.69230769]
目前不知道如何从那里进行。
编辑4：
我目前认为通过将体素旋转转换为基系数存在一些问题 voxel.reshape()
创建索引数组的简单实验
indices = np.arange(27)
indices3d = indices.reshape((3,3,3))
voxel10 = np.random.normal(size=(3,3,3))
assert voxel10[0,1,2] == voxel10.ravel()[indices3d[0,1,2]]
然后将其用于旋转
gen = rotations24(indices3d)
res = []
for ind2 in gen:
basis1 = np.array([[i+1,j+1,k+1] for k in range(3) for j in range(3) for i in range(3)]).astype(np.double)
voxel1 = voxel10.copy().reshape(27,1) #np.array([voxel10[i,j,k] for k in range(3) for j in range(3) for i in range(3)])[...,np.newaxis]
voxel2 = voxel1[ind2.reshape(27,)]
basis1 += 1e-4*np.random.normal(size=(27, 1)) # perturbation
basis2 = basis1[ind2.reshape(27,)]
original_diff = np.sum(np.abs(basis1-basis2))
gg= np.abs(pca(voxel1, basis1) - pca(voxel2, basis2))
ss= np.sum(np.ravel(gg))
bl= np.all(gg<1e-4)
print('difference before pca %.3f,' % original_diff, 'difference after pca %.3f' % ss,bl)
res.append(bl)
del basis1, basis2
print('correct for %.1f percent of time' % (100*(np.sum(res) / len(res))))
表明这些旋转是不正确的，因为在我看来旋转的体素和基础应该匹配：
difference before pca 0.000, difference after pca 0.000 True
difference before pca 48.006, difference after pca 87.459 False
difference before pca 72.004, difference after pca 70.644 False
difference before pca 48.003, difference after pca 71.930 False
difference before pca 72.004, difference after pca 79.409 False
difference before pca 84.005, difference after pca 36.177 False
编辑 5：好吧，所以我们至少要旋转 24 圈。起初，我们的 pca 函数中潜伏着轻微的逻辑变化。在这里，我们居中x_c（基础）并忘记它，进一步居中x_f（特征*基础）并使用 pca 对其进行转换。这可能是行不通的，因为我们的基础没有居中，并且乘以特征进一步增加了偏差。如果我们x_c先居中，然后乘以特征，一切都会好起来的。此外，之前我们proj = x_c @ v使用v计算的方式x_f在这种情况下是完全错误的，因为x_f和x_c以不同的中心为中心。
def pca(feat, x):
# pca with attemt to create convention on sign changes
x_c = x - np.mean(x,axis=0)
x_f = feat * x
x_f -= np.mean(x_f,axis=0)
cov = np.cov(x_f.T)
e_values, e_vectors = np.linalg.eig(cov)
order = np.argsort(e_values)[::-1]
v = e_vectors[:,order]
# here is the solution, we switch the sign of the projections
# so that the projection with the largest absolute value along a principal axis is positive
proj = x_f @ v
return proj
其次，正如我们已经发现的，我们需要对 pca 获得的向量进行排序，例如按第一列：
basis2 = basis1
original_diff = np.sum(np.abs(basis1-basis2))
a = pca(voxel1, basis1)
t1 = a[a[:,0].argsort()]
a = pca(voxel2, basis2)
t2 = a[a[:,0].argsort()]
gg= np.abs(t1-t2)
我们还发现的最后一件事是，reshape体素简单是错误的，它必须对应于旋转：
voxel2 = voxel1[ind2.reshape(27,)] #np.take(voxel10, ind2).reshape(27,1).
理解解决方案的另一个重要注释。当我们对 3d 向量（点云，由我们的基础定义）执行 PCA 并分配权重（类似于刚体的惯性）时，权重的实际分配给点是一种外部信息，对于算法来说是硬定义的。当我们通过应用旋转矩阵旋转基时，我们没有改变数组中向量的顺序，因此质量分配的顺序也没有改变。当我们开始旋转体素时，我们改变了质量的顺序，所以一般来说，如果没有对基础应用相同的变换，PCA 算法将无法工作。因此，只有当我们有一些 3d 向量数组，通过一些旋转进行转换并且质量列表相应地重新排列时，我们才能使用 PCA 检测刚体的旋转。否则，如果我们将质量从点中分离出来，那通常会是另一个物体。
那么它对我们有什么作用呢？它起作用是因为我们的点在居中基础上围绕中心完全对称。在这种情况下，质量的重新分配不会改变“身体”，因为向量范数是相同的。在这种情况下，我们可以使用相同的（数字）basis2=basis1用于测试 24 次旋转和旋转体素 2（旋转点云立方体匹配，只是质量迁移）。这对应于点云的旋转，质点围绕立方体的中心。PCA将根据身体的“惯性”以相同的方式转换具有相同长度和不同质量的向量（在我们对组件的符号达成约定之后）。唯一剩下的就是最后对 pca 变换向量进行排序，因为它们在数组中的位置不同（因为我们的身体旋转了，质点改变了它们的位置）。这让我们丢失了一些与向量顺序相关的信息，但这看起来是不可避免的。
这是检查 24 次旋转的解决方案的代码。如果理论上应该也适用于一般情况，为在更大体素内旋转的更复杂对象提供一些更接近的值：
import numpy as np
from numpy.random import randn
#np.random.seed(20)
def pca(feat, x):
# pca with attemt to create convention on sign changes
x_c = x - np.mean(x,axis=0)
x_f = feat * x_c
cov = np.cov(x_f.T)
e_values, e_vectors = np.linalg.eig(cov)
order = np.argsort(e_values)[::-1]
v = e_vectors[:,order]
# here is the solution, we switch the sign of the projections
# so that the projection with the largest absolute value along a principal axis is positive
proj = x_f @ v
asign = np.sign(proj)
max_ind = np.argmax(np.abs(proj),axis=0)[None,:]
sign = np.take_along_axis(asign,max_ind,axis=0)
proj = proj * sign
return proj
# must be correct /sf/ask/1066112561/
indices = np.arange(27)
indices3d = indices.reshape((3,3,3))
voxel10 = np.random.normal(size=(3,3,3))
assert voxel10[0,1,2] == voxel10.ravel()[indices3d[0,1,2]]
r
首先你的pca函数不正确，应该是
def pca(x):
x -= np.mean(x,axis=0)
cov = np.cov(x.T)
e_values, e_vectors = np.linalg.eig(cov)
order = np.argsort(e_values)[::-1]
v = e_vectors[:,order]
return x @ v
您不应该转置 ，e_vectors[:,order]因为我们希望v数组的每一列都是一个特征向量，因此，x @ v将是x这些特征向量上的投影。
其次，我认为你误解了轮换的含义。不是voxel1应该旋转，而是basis1. 如果进行旋转（通过转置）voxel1，您实际上所做的是重新排列网格点的索引，而点的坐标basis1不会改变。
为了旋转点（例如绕 z 轴），您可以首先定义一个函数来计算给定角度的旋转矩阵
rotate_mat = lambda theta: np.array([[np.cos(theta),-np.sin(theta),0.],[np.sin(theta),np.cos(theta),0.],[0.,0.,1.]])
使用此函数生成的旋转矩阵，您可以旋转数组basis1以创建另一个数组basis2
basis2 = basis1 @ rotate_mat(np.deg2rad(angle))
现在谈到你的问题的标题“为什么我的PCA对于旋转和轴交换不是不变的？”，从这篇文章来看，PCA结果不是唯一的，你实际上可以运行一个测试来看看这个
import numpy as np
np.random.seed(10)
rotate_mat = lambda theta: np.array([[np.cos(theta),-np.sin(theta),0.],[np.sin(theta),np.cos(theta),0.],[0.,0.,1.]])
def pca(x):
x -= np.mean(x,axis=0)
cov = np.cov(x.T)
e_values, e_vectors = np.linalg.eig(cov)
order = np.argsort(e_values)[::-1]
v = e_vectors[:,order]
return x @ v
angles = np.linspace(0,90,91)
res = []
for angle in angles:
voxel1 = np.random.normal(size=(3,3,3))
voxel2 = voxel1.copy()
voxel1 = voxel1.reshape(27,1)
voxel2 = voxel2.reshape(27,1)
basis1 = np.array([[i+1,j+1,k+1] for k in range(3) for j in range(3) for i in range(3)])
# basis2 = np.hstack((-basis1[:,1][:,None],basis1[:,0][:,None],-basis1[:,2][:,None]))
basis2 = basis1 @ rotate_mat(np.deg2rad(angle))
voxel1 = voxel1*basis1
voxel2 = voxel2*basis2
print(angle,np.all(np.abs(pca(voxel1) - pca(voxel2) < 1e-6)))
res.append(np.all(np.abs(pca(voxel1) - pca(voxel2) < 1e-6)))
print()
print(np.sum(res) / len(angles))
运行此脚本后，您将发现只有 21% 的情况下两个 PCA 结果相同。
更新
我认为您可以关注投影，而不是关注主成分的特征向量。对于两个点云，即使它们本质上相同，特征向量也可能截然不同。因此，硬编码以使两组特征向量相同是一项非常困难的任务。
然而，根据这篇文章，对于相同的点云，两组特征向量最多只能有一个负号不同。因此，两组特征向量的投影也仅相差一个负号。这实际上为我们提供了一个优雅的解决方案，对于沿特征向量（主轴）的投影，我们需要做的就是切换投影的符号，使得沿该主轴的绝对值最大的投影为正。
import numpy as np
from numpy.random import randn
#np.random.seed(20)
rotmat_z = lambda theta: np.array([[np.cos(theta),-np.sin(theta),0.],[np.sin(theta),np.cos(theta),0.],[0.,0.,1.]])
rotmat_y = lambda theta: np.array([[np.cos(theta),0.,np.sin(theta)],[0.,1.,0.],[-np.sin(theta),0.,np.cos(theta)]])
rotmat_x = lambda theta: np.array([[1.,0.,0.],[0.,np.cos(theta),-np.sin(theta)],[0.,np.sin(theta),np.cos(theta)]])
# based on https://en.wikipedia.org/wiki/Rotation_matrix
rot_mat = lambda alpha,beta,gamma: rotmat_z(alpha) @ rotmat_y(beta) @ rotmat_x(gamma)
deg2rad = lambda alpha,beta,gamma: [np.deg2rad(alpha),np.deg2rad(beta),np.deg2rad(gamma)]
def pca(feat, x):
# pca with attemt to create convention on sign changes
x_c = x - np.mean(x,axis=0)
x_f = feat * x
x_f -= np.mean(x_f,axis=0)
cov = np.cov(x_f.T)
e_values, e_vectors = np.linalg.eig(cov)
order = np.argsort(e_values)[::-1]
v = e_vectors[:,order]
# here is the solution, we switch the sign of the projections
# so that the projection with the largest absolute value along a principal axis is positive
proj = x_f @ v
asign = np.sign(proj)
max_ind = np.argmax(np.abs(proj),axis=0)[None,:]
sign = np.take_along_axis(asign,max_ind,axis=0)
proj = proj * sign
return proj
ref_angles = np.linspace(0.0,90.0,10)
angles = [[alpha,beta,gamma] for alpha in ref_angles for beta in ref_angles for gamma in ref_angles]
voxel1 = np.random.normal(size=(3,3,3))
res = []
for angle in angles:
voxel2 = voxel1.copy()
voxel1 = voxel1.reshape(27,1)
voxel2 = voxel2.reshape(27,1)
basis1 = np.array([[i+1,j+1,k+1] for k in range(3) for j in range(3) for i in range(3)]).astype(np.double)
basis1 = basis1 + 1e-4 * randn(27,3) # perturbation
basis2 = basis1 @ rot_mat(*deg2rad(*angle))
original_diff = np.sum(np.abs(basis1-basis2))
gg= np.abs(pca(voxel1, basis1) - pca(voxel1, basis2))
ss= np.sum(np.ravel(gg))
bl= np.all(gg<1e-4)
print('difference before pca %.3f,' % original_diff, 'difference after pca %.3f' % ss,bl)
res.append(bl)
del basis1, basis2
print('correct for %.1f percent of time' % (100*(np.sum(res) / len(res))))
通过运行此脚本可以看到，主轴上的投影是相同的，这意味着我们已经解决了 PCA 结果不唯一的问题。
回复编辑3
至于你提出的新问题，我认为你错过了一个重要的一点，那就是点云在主轴上的投影是不变的，而不是其他任何东西。因此，如果旋转voxel1并获得，它们是相同的，因为它们各自在点云主轴上的投影是相同的，与voxel2进行比较实际上没有太大意义。pca(voxel1,basis1)pca(voxel2,basis1)
此外，该方法rotate实际上scipy.ndimage更改了信息，正如您通过运行此脚本可以看到的那样
image1 = np.linspace(1,100,100).reshape(10,10)
image2 = scipy.ndimage.rotate(image1, 45, mode='nearest', axes=(0, 1), reshape=False)
image3 = scipy.ndimage.rotate(image2, -45, mode='nearest', axes=(0, 1), reshape=False)
fig,ax = plt.subplots(nrows=1,ncols=3,figsize=(12,4))
ax[0].imshow(image1)
ax[1].imshow(image2)
ax[2].imshow(image3)
输出图像是
可以看到旋转后的矩阵与原始矩阵不一样，原始矩阵的一些信息发生了变化。
回复编辑4
实际上，我们已经快到了，两个 PCA 结果不同，因为我们正在比较不同点的 PCA 分量。
indices = np.arange(27)
indices3d = indices.reshape((3,3,3))
# apply rotations to the indices, it is not weights yet
gen = rotations24(indices3d)
# construct the weights
voxel10 = np.random.normal(size=(3,3,3))
res = []
count = 0
for ind2 in gen:
count += 1
# ind2 is the array of indices after rotation
# reindex the weights with the indices after rotation
voxel1 = voxel10.copy().reshape(27,1)
voxel2 = voxel1[ind2.reshape(27,)]
# basis1 is the array of coordinates where the points are
basis1 = np.array([[i+1,j+1,k+1] for k in range(3) for j in range(3) for i in range(3)]).astype(np.double)
basis1 += 1e-4*np.random.normal(size=(27, 1))
# reindex the coordinates with the indices after rotation
basis2 = basis1[ind2.reshape(27,)]
# add a slight modification to pca, return the axes also
pca1,v1 = pca(voxel1,basis1)
pca2,v2 = pca(voxel2,basis2)
# sort the principal components before comparing them
pca1 = np.sort(pca1,axis=0)
pca2 = np.sort(pca2,axis=0)
gg= np.abs(pca1 - pca2)
ss= np.sum(np.ravel(gg))
bl= np.all(gg<1e-4)
print('difference after pca %.3f' % ss,bl)
res.append(bl)
del basis1, basis2
print('correct for %.1f percent of time' % (100*(np.sum(res) / len(res))))
运行此脚本，您会发现，对于每次旋转，两组主轴仅相差一个负号。两组 PCA 结果不同，因为旋转前后点云的索引不同（因为您对索引应用了旋转）。如果在比较之前对 pca 结果进行排序，您会发现两个 pca 结果完全相同。
概括
这个问题的答案可以分为两部分。在第一部分中，旋转应用于基础（点的坐标），而索引和相应的权重保持不变。在第二部分中，将旋转应用于指数，然后使用新指数重新排列权重和基础。对于这两部分，解pca函数是相同的
def pca(feat, x):
# pca with attemt to create convention on sign changes
x_c = x - np.mean(x,axis=0)
x_f = feat * x
x_f -= np.mean(x_f,axis=0)
cov = np.cov(x_f.T)
e_values, e_vectors = np.linalg.eig(cov)
order = np.argsort(e_values)[::-1]
v = e_vectors[:,order]
# here is the solution, we switch the sign of the projections
# so that the projection with the largest absolute value along a principal axis is positive
proj = x_f @ v
asign = np.sign(proj)
max_ind = np.argmax(np.abs(proj),axis=0)[None,:]
sign = np.take_along_axis(asign,max_ind,axis=0)
proj = proj * sign
return proj
这个函数的想法是，我们可以匹配主成分，而不是匹配主轴，因为毕竟主成分是旋转不变的。
基于这个函数pca，这个答案的第一部分很容易理解，因为点的索引没有改变，而我们只旋转基础。为了理解这个答案的第二部分（回复编辑5），我们必须首先理解这个函数rotations24。该函数旋转的是索引而不是点的坐标，因此，如果我们保持在同一位置观察点，我们会感觉到点的位置发生了变化。
考虑到这一点，就不难理解Reply to EDIT 5了。
实际上，这个答案中的函数pca可以应用于更一般的情况，例如（我们旋转索引）
num_of_points_per_dim = 10
num_of_points = num_of_points_per_dim ** 3
indices = np.arange(num_of_points)
indices3d = indices.reshape((num_of_points_per_dim,num_of_points_per_dim,num_of_points_per_dim))
voxel10 = 100*np.random.normal(size=(num_of_points_per_dim,num_of_points_per_dim,num_of_points_per_dim))
gen = rotations24(indices3d)
res = []
for ind2 in gen:
voxel1 = voxel10.copy().reshape(num_of_points,1)
voxel2 = voxel1[ind2.reshape(num_of_points,)]
basis1 = 100*np.random.rand(num_of_points,3)
basis2 = basis1[ind2.reshape(num_of_points,)]
pc1 = np.sort(pca(voxel1, basis1),axis=0)
pc2 = np.sort(pca(voxel2, basis2),axis=0)
gg= np.abs(pc1-pc2)
ss= np.sum(np.ravel(gg))
bl= np.all(gg<1e-4)
print('difference after pca %.3f' % ss,bl)
res.append(bl)
del basis1, basis2
print('correct for %.1f percent of time' % (100*(np.sum(res) / len(res))))
| 归档时间： |
|
| 查看次数： |
359 次 |
| 最近记录： |