Why does PCA result change drastically with a slight change in the input?

Question

I am using PCA to reduce an Nx3 array to an Nx2 array. This is mainly because the PCA transformation (Nx2 matrix) is invariant to the rotations or translations performed on the original Nx3 array. Let's take the following as an example.

import numpy as np
from sklearn.decomposition import PCA
a = np.array([[0.5  , 0.5  , 0.5  ],
              [0.332, 0.456, 0.751],
              [0.224, 0.349, 0.349],
              [0.112, 0.314, 0.427]])
pca = PCA(n_components=2)
print(pca.fit_transform(a))

The following is the output. Note that we get the same output with print(pca.fit_transform(a-L)), L being any number, due to translational invariance. Same with rotations.

[[ 0.16752654  0.15593431]
 [ 0.20568992 -0.14688601]
 [-0.16899598  0.06364857]
 [-0.20422047 -0.07269687]]

Now, I give a very small perturbation (~1%) to the array a and perform PCA.

a_p = np.array([[0.51 , 0.53 , 0.52 ],
       [0.322, 0.452, 0.741],
       [0.217, 0.342, 0.339],
       [0.116, 0.31 , 0.417]])
pca = PCA(n_components=2)
print(pca.fit_transform(a_p))

The result is as follows. This is quite different from the PCA of the original array.

 [[-0.2056024 , -0.14346977]
 [-0.18563578  0.15627932]
 [ 0.17974942 -0.07001969]
 [ 0.21148876  0.05721014]]

I expected the PCA transformation of the perturbed array to be very similar to that of the original array, but the percentage change is huge. Why is this? Is there any way I could obtain a very similar PCA transformation for a slightly perturbed/shaked array?

I understand that I can get a similar PCA by performing only transform operation in the second case (e.g. pca.transform(a_p)), however, in this case, I lose rotational and translational invariance w.r.t. a_p.

The question is originally related to crystallography and my requirement is that the PCA (or whatever) transformation should not change significantly to a small change in the input as well as it should be invariant to the rotations and transformations of the input. Can anyone explain the above or suggest me an alternative method which serves my purpose?

EDIT

As mentioned in a comment I will try to explain the actual problem I am trying to solve and the properties of a good solution. Basically, I want to obtain a 2-dimensional representation of a 3D crystal structure (periodic unit cell shown below). This 2D representation should be invariant to the rotations or translations of the original 3D structure (because it's the same crystal anyway).

Let's say we have a unit cell with N number of atoms. The fractional coordinates of selected 4 atoms are indicated above. I first represent this crystal structure by an Nx3 matrix (A), 3 columns being x,y & z coordinates.

\begin{bmatrix}0&0&0\\0&1&0\\0&0&1\\0.5&0.5&0.5\\..&..&..\end{bmatrix}

Now, I want to find a rotationally and translationally independent representation of A. To do this, I perform PCA on A to obtain an Nx2 matrix (B). Due to the nature of PCA, I get rotational and translational invariance (i.e. B = PCA(A) = PCA(A-L) = PCA(A_rotated_by_x_degrees) where L and x are any numbers). However, here's my problem.

A small change in matrix A results in a completely different principle components (B). For example, if I change only the first element of first row from 0 to 0.5, the PCA result changes dramatically resulting in a completely new B. I want this change in B to reflect the slight change that I made in A (e.g. preserve the correlation between other unchanged atom coordinates/retain the other rows as they are, if possible). Can I expect this from PCA? Could you suggest a better solution? It doesn't necessarily have to be PCA.

EDIT2:

Following @whuber's comment I'd like to show how I verified rotational invariance of PCA ($\mathbb{R}^2$) for any rotation in the 3D structure ($\mathbb{R}^3$). Please correct me if I am doing anything wrong. The code below is partly adapted from this SO thread. Using the above array a as the input;

from numpy import cross, eye, dot
from scipy.linalg import expm, norm

def M(axis, theta):
    return expm(cross(eye(3), axis/norm(axis)*theta))

axis = [4,3,-1] # can be rotated around any axis
theta = 1.2 # can be any rotation
M0 = M(axis, theta)

a_rotated = []  # 'a' rotated by theta radians around the given axis
for vec in a:
    vec_r = dot(M0,vec)  # rotate each individual vector
    a_rotated.append(vec_r)

a_rotated = np.array(a_rotated)
print(a_rotated)

a_rotated is as follows;

[[ 0.8410511  -0.05506857  0.19899872]
 [ 0.80627738 -0.30144367  0.37577851]
 [ 0.49270955 -0.05237642  0.21970892]
 [ 0.42659993 -0.14804179  0.29927434]]

Now, I obtain the PCA of a_rotated.

print(pca.fit_transform(a_rotated))

This gives;

[[ 0.16752654  0.15593431]
 [ 0.20568992 -0.14688601]
 [-0.16899598  0.06364857]
 [-0.20422047 -0.07269687]]

which is exactly the same to the PCA of array a. Doesn't this mean PCA is invariant to the rotations of the input?