#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!!!

Leave a Reply

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