# #StackBounty: #python #numpy #math #pca #voxel Why my PCA is not invariant to rotation and axis swap?

### Bounty: 50

I have a voxel (np.array) with size 3x3x3, filled with some values, this setup is essential for me. I want to have rotation-invariant representation of it. For this case, I decided to try PCA representation which is believed to be invariant to orthogonal transformations. another

For simplicity, I took some axes swap, but in case I’m mistaken there can be `np.rot90`.

I have interpereted my 3d voxels as a set of weighted 3d cube point vectors which I incorrectly called "basis", total 27 (so that is some set of 3d point in space, represented by the vectors, obtained from cube points, scaled by voxel values).

``````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)
``````

Then I did PCA to those 27 3-dimensional vectors:

``````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)
``````

But the results in `vp1` and `vp2` are different. Perhaps, I have a mistake (though I beleive this is the right formula), and the proper code must be

`x.dot(v.T)`

But in this case the results are very strange. The upper and bottom blocks of the transofrmed data are the same up to the sign:

``````>>> 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]])
``````

What I’m doing wrong?

What I want to do is to find some invariant representation of my weighted voxel, something like positioning according to the axes of inertia or principal axes. I would really appreciate if someone helps me.

UPD: Found the question similar to mine, but code is unavailable

EDIT2: Found the code InertiaRotate and managed to monkey-do the following:

``````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)
``````

It seems to use basis (coordinates) and mass separately. The results are quite similar to my problem above: some parts of the transformed data match up to the sign, I believe those are some cube sides

`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]]

``````

Looking for some explanation. This code is designed for molecules, and must work…

UPD: Tried to choose 3 vectors as a new basis from those 24 – the one with biggest norm, the one with the smallest and their cross product. Combined them into the matrix V, then used the formula V^(-1)*X to transform coordinates, and got the same problem – the resulting sets of vectors are not equal for rotated voxels.

Get this bounty!!!

This site uses Akismet to reduce spam. Learn how your comment data is processed.