如何解决为什么我的PCA不变于旋转和轴交换？
我有一个尺寸为3x3x3的体素（np.array），其中填充了一些值，此设置对我来说至关重要。我想要它的旋转不变表示。对于这种情况，我决定尝试PCA表示，该表示被认为是invariant to orthogonal transformations。 another
为简单起见，我进行了一些轴交换，但是如果我弄错了，可能会有np.rot90。
我将3d体素作为一组加权的3d立方体点向量进行了穿插处理，我错误地将其称为“基础”，总共27个（因此，这是一组空间中的3d点，由从立方体点获得的向量表示，通过体素值缩放）。
import numpy as np
voxel1 = np.random.normal(size=(3,3,3))
voxel2 = np.transpose(voxel1,(1,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,[ True,True,True],False]])
我做错了什么？
我想要做的是找到我的加权体素的一些不变表示，例如根据惯性轴或主轴进行定位。如果有人帮助我，我将不胜感激。
UPD：Found the question similar to mine，但代码不可用
EDIT2：找到了代码InertiaRotate，并成功地完成了以下操作：
import numpy as np
# https://github.com/smparker/orient-molecule/blob/master/orient.py
voxel1 = np.random.normal(size=(3,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,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：我同意柴可夫斯基的观点，即将体素矢量乘以权重从而创建一些非立方点云的想法是错误的。也许，我们确实需要采取旋转“基础”的解决方案（是的，这不是基础，而是确定点云的一种方法），当“基础”相同时，权重将根据3D旋转。
根据meTchaikovsky提供的答案和参考，发现other answers，我们和我的朋友一起得出的结论是，上述分子包装中的rotate_func试图发明一些计算符号的惯例。组件。他们的解决方案尝试将第三个轴的第三个力矩用作最后一个轴的行列式（？）。我们尝试了另一种方法，并成功使一半的表示形式匹配：
# -*- 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.,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))
res = []
for angle in angles:
voxel2 = voxel1.copy()
voxel1 = voxel1.reshape(27,1)
voxel2 = voxel2.reshape(27,1)
basis1 = np.array([[i+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）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,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,mode='nearest',axes=(0,1),reshape=False)
# rotate along y-axis
image3 = scipy.ndimage.rotate(image2,2),reshape=False)
# rotate along x-axis
image4 = scipy.ndimage.rotate(image3,gamma,axes=(1,reshape=False)
return image4
voxel10 = np.random.normal(size=(3,3))
angles = [[x,y,z] for x in [-90,90] for y in [-90,90] for z in [-90,90]]
res = []
for angle in angles:
voxel2 = rotate_3d(voxel10,angle[0],angle[1],angle[2])
voxel1 = voxel10.reshape(27,1)
basis1 = np.array([[i+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 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。但是旋转的体素会产生不同的结果。万一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,2)
yield from rotations4(rot90(polycube,-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)
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,axis)
m = np.rot90(m,k)
m = np.swapaxes(m,axis)
return m
voxel10 = np.random.normal(size=(3,3))
gen = rotations24(voxel10)
res = []
for voxel2 in gen:
#voxel2 = rotate_3d(voxel10,1)) # perturbation
basis2 = basis1
original_diff = np.sum(np.abs(basis1-basis2))
gg= np.abs(pca(voxel1,basis2
print('correct for %.1f percent of time' % (100*(np.sum(res) / len(res))))
我试图调查此案，而唯一可能无关的发现是：
voxel1 = np.ones((3,3))
voxel1[0,0] = 0 # if I change 0 to 0.5 it stops working at all
# mirrored around diagonal
voxel2 = np.ones((3,3))
voxel2[2,2] = 0
for angle in range(1):
voxel1 = voxel1.reshape(27,1)
basis1 = np.array([[i+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]
当前不知道如何从那里继续。
EDIT4：
我目前正在考虑将体素旋转通过voxel.reshape()转换为基本系数存在一些问题
创建索引数组的简单实验
indices = np.arange(27)
indices3d = indices.reshape((3,3))
voxel10 = np.random.normal(size=(3,3))
assert voxel10[0,1,2] == voxel10.ravel()[indices3d[0,2]]
然后将其用于旋转
gen = rotations24(indices3d)
res = []
for ind2 in gen:
basis1 = np.array([[i+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,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,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个旋转和旋转的voxel2（旋转的点云立方体匹配，只是质量会迁移）。这对应于具有围绕立方体中心的质量点的点云的旋转。然后，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,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)
proj = proj * sign
return proj
# must be correct https://stackoverflow.com/questions/15230179/how-to-get-the-linear-index-for-a-numpy-array-sub2ind
indices = np.arange(27)
indices3d = indices.reshape((3,2]]
rot90 = np.rot90
def rotations24(polycube):
# imagine shape is pointing in axis 0 (up)
# 4 rotations about axis 0
yield from rotations4(polycube,axis)
return m
gen = rotations24(indices3d)
res = []
for ind2 in gen:
basis1 = np.array([[i+1,1)
voxel2 = voxel1[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,0].argsort()]
a = pca(voxel2,0].argsort()]
gg= np.abs(t1-t2)
ss= np.sum(np.ravel(gg))
bl= np.all(gg<1e-4)
print('difference before pca %.3f,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
PS。我想提出一个更好的排序主题，以考虑体素中的零值，当PCA向量的整个第一列为零时，这可能会混淆先前的方法，依此类推。我提议按向量范数排序，再乘以元素总和的符号。这是tensorflow 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,1],strides=[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)[:,:]
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),:]
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,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.,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),1.]])
def pca(x):
x -= np.mean(x,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][:,-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))
运行此脚本后，您将看到两个PCA结果相同的可能性只有21％。
更新
我认为，您可以将注意力集中在投影上，而不是专注于主成分的特征向量。对于两点云，即使它们本质上相同，特征向量也可能截然不同。因此，以某种方式使两组特征向量相同的硬编码是一项非常困难的任务。
但是，基于this post，对于相同的点云，两组特征向量仅在负号之前可以不同。因此，两组特征向量的投影也仅在负号之前不同。实际上，这为我们提供了一种优雅的解决方案，对于沿特征向量（主轴）的投影，我们要做的就是切换投影的符号，以使沿该主轴的绝对值最大的投影为正。 >
import numpy as np
from numpy.random import randn
#np.random.seed(20)
rotmat_z = lambda theta: np.array([[np.cos(theta),1.]])
rotmat_y = lambda theta: np.array([[np.cos(theta),np.sin(theta)],1.,[-np.sin(theta),np.cos(theta)]])
rotmat_x = lambda theta: np.array([[1.,-np.sin(theta)],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,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,gamma] for alpha in ref_angles for beta in ref_angles for gamma in ref_angles]
voxel1 = np.random.normal(size=(3,3))
res = []
for angle in angles:
voxel2 = voxel1.copy()
voxel1 = voxel1.reshape(27,1)
basis1 = np.array([[i+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结果不唯一的问题。
回复EDIT 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,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)
输出图像是
如您所见，旋转后的矩阵与原始矩阵不同，因此原始矩阵的某些信息已更改。
回复EDIT 4
实际上，我们差不多到了，两个pca结果是不同的，因为我们正在比较不同点的pca成分。
indices = np.arange(27)
indices3d = indices.reshape((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))
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,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,order]
# here is the solution,axis=0)
proj = proj * sign
return proj
此功能的想法是，我们不必匹配主轴，而是可以匹配主分量，因为毕竟它们是旋转不变的主分量。
基于此功能pca，此答案的第一部分很容易理解，因为点的索引不变，而我们仅旋转基础。为了理解此答案的第二部分（答复EDIT 5 ），我们必须首先了解函数rotations24。此功能旋转索引而不是点的坐标，因此，如果我们在观察点的位置相同，则会感觉到点的位置已更改。
考虑到这一点，不难理解 回复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))
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,basis2
print('correct for %.1f percent of time' % (100*(np.sum(res) / len(res))))
版权声明：本文内容由互联网用户自发贡献，该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务，不拥有所有权，不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容， 请发送邮件至 dio@foxmail.com 举报，一经查实，本站将立刻删除。
