为什么我的 PCA 对旋转和轴交换不是不变的？
Why my PCA is not invariant to rotation and axis swap?
我有一个大小为 3x3x3 的体素 (np.array)，填充了一些值，这个设置对我来说很重要。我想要它的 rotation-invariant 表示。对于这种情况，我决定尝试 PCA 表示，它被认为是 invariant to orthogonal transformations. another
为简单起见，我进行了一些轴交换，但以防万一我弄错了，可以 np.rot90。
我将我的 3d 体素解释为一组加权 3d 立方体点向量，我错误地称之为“基础”，总共 27 个（所以这是 space 中的一组 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：Found the question similar to mine，但代码不可用
EDIT2：找到代码 InertiaRotate 并设法 monkey-do 以下内容：
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)
好像是分开用basis（坐标）和mass。结果与我上面的问题非常相似：转换后的数据的某些部分与符号匹配，我相信那些是立方体的边
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 的观点，即我将体素向量乘以权重并因此创建一些 non-cubic 点云的想法是不正确的。可能，我们确实需要采用旋转“基础”的解决方案（是的，这不是基础，而是一种确定点云的方法），当“基础”相同时，它将在稍后工作，但权重根据旋转到 3D 旋转。
根据 meTchaikovsky 提供的答案和参考资料，并找到 other answers 我们和我的朋友一起得出结论，rotate_func 来自上述分子包，试图发明一些约定来计算组件的标志。他们的解决方案尝试对前 2 个轴使用 3rd 力矩，对最后一个轴（？）使用行列式。我们尝试了另一种方法并成功地匹配了一半的表示：
# -*- 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
从上面的代码可以看出，我们尝试使用另一个基：范数最大的向量，范数最小的向量及其叉积。这里我们应该只有两个变体（叉积的方向）可以在以后修复，但我无法管理这个替代解决方案。
我希望有人能帮我完成这个并获得 rotation-invariant 体素表示。
编辑 3. 非常感谢我柴可夫斯基，但情况仍不明朗。我的问题最初在于处理 3d 体素，它们是 (3,3,3) numpy 数组。我们得出的结论是，为了找到不变的表示，我们只需要将 3d 体素固定为用于计算 cov 矩阵的权重，并在居中的“基础”（一些用于描述点云的向量）上应用旋转。
因此，当我们实现对“基础”旋转的不变性时，问题应该已经解决了：现在，当我们固定“基础”并使用旋转体素时，结果必须是不变的。令人惊讶的是，事实并非如此。在这里，我使用 basis2=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
此处 basis1=basis2（因此 pca=0 之前的基差），您可以看到 (0,0,0) 旋转为 0。但是旋转的体素会产生不同的结果。万一 scipy 做错了什么，我用 检查了方法，结果相同：
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（basis）而不管它，进一步居中x_f（features*basis）并用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 并分配权重（类似于刚性 body 的惯性）时，实际 分配权重到点 是一种外部信息，对于算法来说变成了 hard-defined。当我们通过应用旋转矩阵旋转基时，我们没有改变数组中向量的顺序，因此质量分配的顺序也没有改变。当我们开始旋转体素时，我们改变了质量的顺序，因此如果不对基础应用相同的变换，通常 PCA 算法将无法工作。因此，只有当我们有一些 3d 向量数组，通过一些旋转和相应的质量列表 re-arranged 进行变换时，我们才能使用 PCA 检测刚性 body 的旋转。否则，如果我们从点中分离出质量，那将是另一个 body。
那么它对我们有什么作用呢？它之所以有效，是因为我们的点在居中基础之后围绕中心完全对称。在这种情况下，质量的重新分配不会改变“body”，因为向量范数是相同的。在这种情况下，我们可以使用相同的（数值）basis2=basis1 来测试 24 次旋转和旋转的体素 2（旋转的点云立方体匹配，只是质量迁移）。这对应于质量点围绕立方体中心的点云的旋转。然后，PCA 将根据 body 的“惯性”以相同的方式对相同长度和不同质量的向量进行变换（在我们对分量的符号达成约定之后）。最后剩下的唯一事情就是对 pca 变换后的向量进行排序，因为它们在数组中的位置不同（因为我们的 body 被旋转了，质点改变了它们的位置）。这使我们丢失了一些与向量顺序相关的信息，但这看起来是不可避免的。
这是检查 24 次旋转的解决方案的代码。如果在理论上也应该在一般情况下工作，为更复杂的 objects 在更大的体素内旋转提供一些更接近的值：
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
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]]
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
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)
voxel2 = voxel1[ind2.reshape(27,)] #np.take(voxel10, ind2).reshape(27,1)
basis1 += 1e-6*np.random.normal(size=(27, 1)) # perturbation
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)
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 0.000, difference after pca 0.000 True
difference before pca 0.000, difference after pca 0.000 True
difference before pca 0.000, difference after pca 0.000 True
PS。我想提出更好的排序主题，以考虑体素中的零值，当 PCA 向量的整个第一列为零等时，这可能会混淆以前的方法。我建议按向量范数排序，乘以元素总和的符号.这是张量流 2 代码：
def infer_shape(x):
x = tf.convert_to_tensor(x)
# If unknown rank, return dynamic shape
if x.shape.dims is None:
return tf.shape(x)
static_shape = x.shape.as_list()
dynamic_shape = tf.shape(x)
ret = []
for i in range(len(static_shape)):
dim = static_shape[i]
if dim is None:
dim = dynamic_shape[i]
ret.append(dim)
return ret
def merge_last_two_dims(tensor):
shape = infer_shape(tensor)
shape[-2] *= shape[-1]
#shape.pop(1)
shape = shape[:-1]
return tf.reshape(tensor, shape)
def pca(inpt_voxel):
patches = tf.extract_volume_patches(inpt_voxel, ksizes=[1,3,3,3,1], strides=[1, 1,1,1, 1], padding="VALID")
features0 = patches[...,tf.newaxis]*basis
# centered basises
basis1_ = tf.ones(shape=tf.shape(patches[...,tf.newaxis]), dtype=tf.float32)*basis
basis1 = basis1_ - tf.math.divide_no_nan(tf.reduce_sum(features0, axis=-2), tf.reduce_sum(patches, axis=-1)[...,None])[:,:,:,:,None,:]
features = patches[...,tf.newaxis]*basis1
features_centered_basis = features - tf.reduce_mean(features, axis=-2)[:,:,:,:,None,:]
x = features_centered_basis
m = tf.cast(x.get_shape()[-2], tf.float32)
cov = tf.matmul(x,x,transpose_a=True)/(m - 1)
e,v = tf.linalg.eigh(cov,name="eigh")
proj = tf.matmul(x,v,transpose_b=False)
asign = tf.sign(proj)
max_ind = tf.argmax(tf.abs(proj),axis=-2)[:,:,:,:,None,:]
sign = tf.gather(asign,indices=max_ind, batch_dims=4, axis=-2)
sign = tf.linalg.diag_part(sign)
proj = proj * sign
# But we can have 1st coordinate zero. In this case,
# other coordinates become ambiguous
#s = tf.argsort(proj[...,0], axis=-1)
# sort by l2 vector norms, multiplied by signs of sums
sum_signs = tf.sign(tf.reduce_sum(proj, axis=-1))
norms = tf.norm(proj, axis=-1)
s = tf.argsort(sum_signs*norms, axis=-1)
proj = tf.gather(proj, s, batch_dims=4, axis=-2)
return merge_last_two_dims(proj)
首先，你的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对旋转和轴交换不不变？”的问题，来自this post，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 结果相同。
更新
我认为与其关注主成分的特征向量，不如关注投影。对于两个点云，即使它们本质上相同，特征向量也可能截然不同。因此，硬编码以某种方式让两组特征向量相同是一项非常困难的任务。
然而，基于this post，对于相同的点云，两组特征向量最多只能相差一个负号。因此，两组特征向量的投影也仅相差一个负号。这实际上为我们提供了一个优雅的解决方案，对于沿特征向量（主轴）的投影，我们需要做的就是切换投影的符号，使沿该主轴的绝对值最大的投影为正。
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) 进行比较是有意义的。
此外，scipy.ndimage 的方法 rotate 实际上改变了信息，你可以通过 运行ning 这个脚本看到
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，这个答案的第一部分很容易理解，因为点的索引不变，而我们只旋转基。为了理解这个答案的第二部分（回复EDIT 5），我们必须首先理解函数rotations24。此函数旋转的是索引而不是点的坐标，因此，如果我们在同一位置观察点，我们会感觉到点的位置发生了变化。
想到这里就不难理解了回复EDIT 5[=110=.
实际上，这个答案中的函数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))))
我有一个大小为 3x3x3 的体素 (np.array)，填充了一些值，这个设置对我来说很重要。我想要它的 rotation-invariant 表示。对于这种情况，我决定尝试 PCA 表示，它被认为是 invariant to orthogonal transformations. another
为简单起见，我进行了一些轴交换，但以防万一我弄错了，可以 np.rot90。
我将我的 3d 体素解释为一组加权 3d 立方体点向量，我错误地称之为“基础”，总共 27 个（所以这是 space 中的一组 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：Found the question similar to mine，但代码不可用
EDIT2：找到代码 InertiaRotate 并设法 monkey-do 以下内容：
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)
好像是分开用basis（坐标）和mass。结果与我上面的问题非常相似：转换后的数据的某些部分与符号匹配，我相信那些是立方体的边
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 的观点，即我将体素向量乘以权重并因此创建一些 non-cubic 点云的想法是不正确的。可能，我们确实需要采用旋转“基础”的解决方案（是的，这不是基础，而是一种确定点云的方法），当“基础”相同时，它将在稍后工作，但权重根据旋转到 3D 旋转。
根据 meTchaikovsky 提供的答案和参考资料，并找到 other answers 我们和我的朋友一起得出结论，rotate_func 来自上述分子包，试图发明一些约定来计算组件的标志。他们的解决方案尝试对前 2 个轴使用 3rd 力矩，对最后一个轴（？）使用行列式。我们尝试了另一种方法并成功地匹配了一半的表示：
# -*- 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
从上面的代码可以看出，我们尝试使用另一个基：范数最大的向量，范数最小的向量及其叉积。这里我们应该只有两个变体（叉积的方向）可以在以后修复，但我无法管理这个替代解决方案。
我希望有人能帮我完成这个并获得 rotation-invariant 体素表示。
编辑 3. 非常感谢我柴可夫斯基，但情况仍不明朗。我的问题最初在于处理 3d 体素，它们是 (3,3,3) numpy 数组。我们得出的结论是，为了找到不变的表示，我们只需要将 3d 体素固定为用于计算 cov 矩阵的权重，并在居中的“基础”（一些用于描述点云的向量）上应用旋转。
因此，当我们实现对“基础”旋转的不变性时，问题应该已经解决了：现在，当我们固定“基础”并使用旋转体素时，结果必须是不变的。令人惊讶的是，事实并非如此。在这里，我使用 basis2=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
此处 basis1=basis2（因此 pca=0 之前的基差），您可以看到 (0,0,0) 旋转为 0。但是旋转的体素会产生不同的结果。万一 scipy 做错了什么，我用
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（basis）而不管它，进一步居中x_f（features*basis）并用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 并分配权重（类似于刚性 body 的惯性）时，实际 分配权重到点 是一种外部信息，对于算法来说变成了 hard-defined。当我们通过应用旋转矩阵旋转基时，我们没有改变数组中向量的顺序，因此质量分配的顺序也没有改变。当我们开始旋转体素时，我们改变了质量的顺序，因此如果不对基础应用相同的变换，通常 PCA 算法将无法工作。因此，只有当我们有一些 3d 向量数组，通过一些旋转和相应的质量列表 re-arranged 进行变换时，我们才能使用 PCA 检测刚性 body 的旋转。否则，如果我们从点中分离出质量，那将是另一个 body。
那么它对我们有什么作用呢？它之所以有效，是因为我们的点在居中基础之后围绕中心完全对称。在这种情况下，质量的重新分配不会改变“body”，因为向量范数是相同的。在这种情况下，我们可以使用相同的（数值）basis2=basis1 来测试 24 次旋转和旋转的体素 2（旋转的点云立方体匹配，只是质量迁移）。这对应于质量点围绕立方体中心的点云的旋转。然后，PCA 将根据 body 的“惯性”以相同的方式对相同长度和不同质量的向量进行变换（在我们对分量的符号达成约定之后）。最后剩下的唯一事情就是对 pca 变换后的向量进行排序，因为它们在数组中的位置不同（因为我们的 body 被旋转了，质点改变了它们的位置）。这使我们丢失了一些与向量顺序相关的信息，但这看起来是不可避免的。
这是检查 24 次旋转的解决方案的代码。如果在理论上也应该在一般情况下工作，为更复杂的 objects 在更大的体素内旋转提供一些更接近的值：
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
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]]
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
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)
voxel2 = voxel1[ind2.reshape(27,)] #np.take(voxel10, ind2).reshape(27,1)
basis1 += 1e-6*np.random.normal(size=(27, 1)) # perturbation
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)
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 0.000, difference after pca 0.000 True
difference before pca 0.000, difference after pca 0.000 True
difference before pca 0.000, difference after pca 0.000 True
PS。我想提出更好的排序主题，以考虑体素中的零值，当 PCA 向量的整个第一列为零等时，这可能会混淆以前的方法。我建议按向量范数排序，乘以元素总和的符号.这是张量流 2 代码：
def infer_shape(x):
x = tf.convert_to_tensor(x)
# If unknown rank, return dynamic shape
if x.shape.dims is None:
return tf.shape(x)
static_shape = x.shape.as_list()
dynamic_shape = tf.shape(x)
ret = []
for i in range(len(static_shape)):
dim = static_shape[i]
if dim is None:
dim = dynamic_shape[i]
ret.append(dim)
return ret
def merge_last_two_dims(tensor):
shape = infer_shape(tensor)
shape[-2] *= shape[-1]
#shape.pop(1)
shape = shape[:-1]
return tf.reshape(tensor, shape)
def pca(inpt_voxel):
patches = tf.extract_volume_patches(inpt_voxel, ksizes=[1,3,3,3,1], strides=[1, 1,1,1, 1], padding="VALID")
features0 = patches[...,tf.newaxis]*basis
# centered basises
basis1_ = tf.ones(shape=tf.shape(patches[...,tf.newaxis]), dtype=tf.float32)*basis
basis1 = basis1_ - tf.math.divide_no_nan(tf.reduce_sum(features0, axis=-2), tf.reduce_sum(patches, axis=-1)[...,None])[:,:,:,:,None,:]
features = patches[...,tf.newaxis]*basis1
features_centered_basis = features - tf.reduce_mean(features, axis=-2)[:,:,:,:,None,:]
x = features_centered_basis
m = tf.cast(x.get_shape()[-2], tf.float32)
cov = tf.matmul(x,x,transpose_a=True)/(m - 1)
e,v = tf.linalg.eigh(cov,name="eigh")
proj = tf.matmul(x,v,transpose_b=False)
asign = tf.sign(proj)
max_ind = tf.argmax(tf.abs(proj),axis=-2)[:,:,:,:,None,:]
sign = tf.gather(asign,indices=max_ind, batch_dims=4, axis=-2)
sign = tf.linalg.diag_part(sign)
proj = proj * sign
# But we can have 1st coordinate zero. In this case,
# other coordinates become ambiguous
#s = tf.argsort(proj[...,0], axis=-1)
# sort by l2 vector norms, multiplied by signs of sums
sum_signs = tf.sign(tf.reduce_sum(proj, axis=-1))
norms = tf.norm(proj, axis=-1)
s = tf.argsort(sum_signs*norms, axis=-1)
proj = tf.gather(proj, s, batch_dims=4, axis=-2)
return merge_last_two_dims(proj)
首先，你的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对旋转和轴交换不不变？”的问题，来自this post，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 结果相同。
更新
我认为与其关注主成分的特征向量，不如关注投影。对于两个点云，即使它们本质上相同，特征向量也可能截然不同。因此，硬编码以某种方式让两组特征向量相同是一项非常困难的任务。
然而，基于this post，对于相同的点云，两组特征向量最多只能相差一个负号。因此，两组特征向量的投影也仅相差一个负号。这实际上为我们提供了一个优雅的解决方案，对于沿特征向量（主轴）的投影，我们需要做的就是切换投影的符号，使沿该主轴的绝对值最大的投影为正。
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) 进行比较是有意义的。
此外，scipy.ndimage 的方法 rotate 实际上改变了信息，你可以通过 运行ning 这个脚本看到
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，这个答案的第一部分很容易理解，因为点的索引不变，而我们只旋转基。为了理解这个答案的第二部分（回复EDIT 5），我们必须首先理解函数rotations24。此函数旋转的是索引而不是点的坐标，因此，如果我们在同一位置观察点，我们会感觉到点的位置发生了变化。
想到这里就不难理解了回复EDIT 5[=110=.
实际上，这个答案中的函数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))))