Why my PCA is not invariant to rotation and axis swap?

Question

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.

---

UPD2: I agree with meTchaikovsky that my idea of multiplying voxel vectors by weights and thus creating some non-cubic point cloud was incorrect. Probably, we indeed need to take the solution for rotated "basis"(yes, this is not a basis, but rather a way to determine point cloud) which will work later when "basis" is the same, but the weights are rotated according to the 3D rotation.

Based on the answer and the reference provided by meTchaikovsky, and finding other answers we together with my friend came to conclusion that rotate_func from molecular package mentioned above tries to invent some convention for computing the signs of the components. Their solution tries to use 3rd moment for the first 2 axes and determinant for the last axis (?). We tried a bit another approach and succeeded to have half of the representations matching:

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

The results are good up to 58 degree angle (yet we're still experimenting with rotation of x, y axes). After that we have constant difference which indicates some uncounted sign reverse. This is better than the less consistent result of 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

From the code above, you can see that we tried to use another basis: vector with the biggest norm, vector with the smallest and their cross product. Here we should have only two variants (direction of the cross product) which could be later fixed, but I couldn't manage this alternative solution to work.

I hope that someone can help me finish this and obtain rotation-invariant representation for voxels.

---

EDIT 3. Thank you very much meTchaikovsky, but the situation is still unclear. My problem initially lies in processing 3d voxels which are (3,3,3) numpy arrays. We reached the conclusion that for finding invariant representation, we just need to fix 3d voxel as weights used for calculating cov matrix, and apply rotations on the centered "basis" (some vectors used for describing point cloud).

Therefore, when we achieved invariance to "basis" rotations, the problem should have been solved: now, when we fix "basis" and use rotated voxel, the result must be invariant. Surprisingly, this is not so. Here I check 24 rotations of the cube with basis2=basis1 (except small perturbation):

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

Here basis1=basis2 (hence basis difference before pca=0), and you can see 0 for (0,0,0) rotation. But rotated voxels give different result. In case scipy does something wrong, I've checked the approach with numpy.rot90 with the same result:

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

I tried to investigate this case, and the only perhaps irrelevant thing I found the following:

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]

Currently don't know how to proceed from there.

---

EDIT4:

I'm currently thinking that there is some problem with voxel rotations transformed into basis coefficients via voxel.reshape()

Simple experiment with creating array of indices

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

And then using it for rotations

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

Shows that those rotations are not correct, because on my opinion rotated voxel and basis should match:

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

---

EDIT 5: Okaaay, so here we go at least for 24 rotations. At first, we had a slight change of logic lurked into our pca function. Here we center x_c (basis) and forget about it, further centering x_f (features*basis) and transforming it with pca. This does not work perhaps because our basis is not centered and multiplication by features further increased the bias. If we center x_c first, and multiply it by features, everything will be Ok. Also, previously we had proj = x_c @ v with v computed from x_f which was totally wrong in this case, as x_f and x_c were centered around different centers.

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

Secondly, as we already found, we need to sort vectors obtained by pca, for example by the first column:

    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)

And the last thing we also discovered already, is that simple reshape is wrong for voxel, it must correspond to rotation:

voxel2 = voxel1[ind2.reshape(27,)] #np.take(voxel10, ind2).reshape(27,1).

One more important comment to understand the solution. When we perform PCA on the 3d vectors (point cloud, defined by our basis) with weights assigned (analogously to the inertia of the rigid body), the actual assignment of the weights to the points is sort of external information, which becomes hard-defined for the algorithm. When we rotated basis by applying rotation matrices, we did not change the order of the vectors in the array, hence the order of the mass assignments wasn't changed too. When we start to rotate voxel, we change the order of the masses, so in general PCA algorithm will not work without the same transformation applied to the basis. So, only if we have some array of 3d vectors, transformed by some rotation AND the list of masses re-arranged accordingly, we can detect the rotation of the rigid body using PCA. Otherwise, if we detach masses from points, that would be another body in general.

So how does it work for us then? It works because our points are fully symmetric around the center after centering basis. In this case reassignment of the masses does not change "the body" because vector norms are the same. In this case we can use the same (numerically) basis2=basis1 for testing 24 rotations and rotated voxel2 (rotated point cloud cubes match, just masses migrate). This correspond to the rotation of the point cloud with mass points around the center of the cube. PCA will transform vectors with the same lengths and different masses in the same way according to the body's "inertia" then (after we reached convention on the signs of the components). The only thing left is to sort the pca transformed vectors in the end, because they have different position in the array (because our body was rotated, mass points changed their positions). This makes us lose some information related to the order of the vectors but it looks inevitable.

Here is the code which checks the solution for 24 rotations. If should theoretically work in the general case as well, giving some closer values for more complicated objects rotated inside a bigger voxel:

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 https://mcmap.net/q/674442/-how-to-get-the-linear-index-for-a-numpy-array-sub2ind
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. I want to propose better ordering theme to take into account zero values in the voxel which might confuse previous approach when entire first column of PCA vectors is zero, etc. I propose to sort by vector norms, multiplied by the sign of the sum of elements. Here is tensorflow 2 code:

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)

Rundell

Answer 1

Firstly, your pca function is not correct, it should be

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

You shouldn't transpose the e_vectors[:,order] because we want each column of the v array is an eigenvector, therefore, x @ v will be projections of x on those eigenvectors.

Secondly, I think you misunderstand the meaning of rotation. It is not voxel1 that should be rotated, but the basis1. If you rotate (by taking transposition) voxel1, what you really do is to rearrange the indices of grid points, while the coordinates of the points basis1 are not changed.

In order to rotate the points (around the z axis for example), you can first define a function to calculate the rotation matrix given an angle

rotate_mat = lambda theta: np.array([[np.cos(theta),-np.sin(theta),0.],[np.sin(theta),np.cos(theta),0.],[0.,0.,1.]])

with the rotation matrix generated by this function, you can rotate the array basis1 to create another array basis2

basis2 = basis1 @ rotate_mat(np.deg2rad(angle))

Now it comes to the title of your question "Why my PCA is not invariant to rotation and axis swap?", from this post, the PCA result is not unique, you can actually run a test to see this

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

After you run this script, you will see that in only 21% of times the two PCA results are the same.

---

UPDATE

I think instead of focusing on the eigenvectors of the principal components, you can instead focus on the projections. For two clouds of points, even though they are essentially the same, the eigenvectors can be drastically different. Therefore, hardcoding in order to somehow let the two sets of eigenvectors to be the same is a very difficult task.

However, based on this post, for the same cloud of points, two sets of eigenvectors can be different only up to a minus sign. Therefore, the projections upon the two sets of eigenvectors are also different only up to a minus sign. This actually offers us an elegant solution, for the projections along an eigenvector (principal axis), all we need to do is to switch the sign of the projections so that the projection with the largest absolute value along that principal axis is positive.

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

As you can see by running this script, the projections on the principal axis are the same, this means we have resolved the issue of PCA results being not unique.

---

Reply to EDIT 3

As for the new issue you raised, I think you missed an important point, it is the projections of the cloud of points onto the principal axes that are invariant, not anything else. Therefore, if you rotate voxel1 and obtain voxel2, they are the same in the sense that their own respective projections onto the principal axes of the cloud of points are the same, it actually does not make too much sense to compare pca(voxel1,basis1) with pca(voxel2,basis1).

Furthermore, the method rotate of scipy.ndimage actually changes information, as you can see by running this script

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)

The output image is

As you can see the matrix after rotation is not the same as the original one, some information of the original matrix is changed.

---

Reply to EDIT 4

Actually, we are almost there, the two pca results are different because we are comparing pca components for different points.

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

Running this script, you will find, for each rotation, the two sets of principal axes are different only up to a minus sign. The two sets of pca results are different because the indices of the cloud of points before and after rotation are different (since you apply rotation to the indices). If you sort the pca results before comparing them, you will find the two pca results are exactly the same.

---

Summary

The answer to this question can be divided into two parts. In the first part, the rotation is applied to the basis (the coordinates of points), while the indices and the corresponding weights are unchanged. In the second part, the rotation is applied to the indices, then the weights and the basis are rearranged with the new indices. For both of the two parts, the solution pca function is the same

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

The idea of this function is, instead of matching the principal axes, we can match the principal components since it is the principal components that are rotationally invariant after all.

Based on this function pca, the first part of this answer is easy to understand, since the indices of the points are unchanged while we only rotate the basis. In order to understand the second part of this answer (Reply to EDIT 5), we must first understand the function rotations24. This function rotates the indices rather than the coordinates of the points, therefore, if we stay at the same position observing the points, we will feel that the positions of the points are changed.

With this in mind, it is not hard to understand Reply to EDIT 5.

Actually, the function pca in this answer can be applied to more general cases, for example (we rotate the indices)

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

Automatism

Thank you very much! I have updated my post. I tried to invent some convention on the signs of the components, which I believe rotation_func does. Also I tried to use three manually created axes, but with no success. I believe this task can be solved – Rundell

@Rundell You are welcome! :) I have updated my post, the task has been solved :) – Automatism

Well, either I'm missing something or there's still a problem when we go back from rotating basis to the initial problem with (3,3,3) voxels. When I fix basis2=basis1 and use rotated voxel, I lose invariance. Or perhaps I have a mistake which I cannot find. Could you please help me a bit more to understand this? – Rundell

@Rundell Hi Slowpoke, I have updated my post and addressed your EDIT 3 in Reply to EDIT 3. – Automatism

Thank you very much again! I still hesitate to close this question, because I haven't solved my problem yet, though I found another possible source of mistake and described it in EDIT4 - .reshape() operation might be not consistent with basis rotations, though it worked somehow when voxel was fixed and basis was rotated. I am trying to find a way to connect basis and voxel rotations and will try to update my results on this later when I'm done – Rundell

@Rundell You are welcome Slowpoke! I also lean a lot of things in solving this problem :) I must apologize first for the mistake in the function pca, the line proj = x_c @ v should be proj = x_f @ v, I have corrected it in the updated post. Actually, we are almost there, if you sort the two pca results, you will find they are exactly the same :) – Automatism

Phew, looks like I've managed to make it working. You can check EDIT 5. Not sure this this the last point in this question, but it seems to pass validation with 24 rotations when cube is matching and the difference is very small. Could you please update your answer perhaps, so I can close the question and others may find the solution easily? I tried to summarize my understanding so far of how it works, not sure if it's clear enough to be helpful – Rundell

@Rundell I have updated my post, in which I showed that the approach also works for general cases. From the last part of code in your post, I guess you might still misunderstand something, rotations24 is a function to rotate the indices, therefore if you rearrange the weights with the new indices, you also need to rearrange the basis with the same indices, otherwise a[a[:,0].argsort()] will not help. – Automatism

Yes, everything you wrote is correct. Indeed, we need to apply the same indices rotation for the basis in general case. But, as I wrote in those two vast paragraphs in EDIT5, there is the only case when we can use non-rotated basis. This is because the centered basis is sort of spherically-symmetric around the center of coordinates and reassigning indices becomes not very necessary for 24 rotations when rearranged cube matches previous one. In this case we have "computational" similarity between those two bodies and we can determine the same invariant features by sorting the output vectors. – Rundell