I tried to make my question clearer. The old digression can still be found at the end of this.
PROBLEM: I have an n by 3 matrix A (np.ndarray) of n points in three dimensions. We say those points are symmetric in relation to a transformation R (3 by 3 matrix) if they are stationary as an unordered set.
This means that A and A @ R.T differ by 1. at most a permutation and 2. after corrected for permutation, both matrices might differ by a numerical tolerance parameter: np.allclose(A, permuted(A @ R.T)) == True (I don't know permuted() beforehand and this surely depends on R).
QUESTION: how do I create the following function:
def is_symmetric(A, R, atol=1e-5) -> bool:
# checks symmetry as defined above, considering both numerical noise
# and permutation of vectors.
(Some discussion on possible ways and my doubts and attempts are to be found below.)
OLD DIGRESSION:
I want to check for symmetries in a collection of points, represented as vectors. That means checking whether those points are invariant regarding the application of a matrix transform in space, e.g., rotation or plane mirroring.
import numpy as np
R = np.array([[0, -1, 0], # 90° rotation around z as an example
[1, 0, 0],
[0, 0, 1]])
The point is that I accept a permutation of the vectors: as long as some initial position gets transformed to some other, preexisting, location, I'm fine. That means checking for transformed pairs of vectors from one to the other.
The naïve solution would be looping over the rows of A @ R.T (where A is a matrix whose rows are point positions) and try to match a transformed vector for each initial row of A, which seems to grow quadratically with the number of columns.
Another possibility would be preordering the vectors (say, by their coordinate values):
A = np.array([[1, 0, 0], # some points
[0, 1, 0],
[0, 0, 1],
[0.5, 0.5, 0]])
A = np.array(sorted(A, key=lambda v: (v[2], v[1], v[0]))) # sort by z, y, x values
# array([[1. , 0. , 0. ],
# [0.5, 0.5, 0. ],
# [0. , 1. , 0. ],
# [0. , 0. , 1. ]])
A_rotated = np.array(sorted(A @ R.T, key=lambda v: (v[2], v[1], v[0])))
# array([[-1. , 0. , 0. ], # no match
# [-0.5, 0.5, 0. ], # no match
# [ 0. , 1. , 0. ], # match
# [ 0. , 0. , 1. ]]) # match
(This approach does two sorts so O(n ln(n))?) The third idea would be to compare sets created from the original and rotated vectors. I have a gut feeling that this is as good as comparing sorted vectors.
But one thing remains: how to deal with approximate comparison? I want to accept two vectors v and w being equal if e.g. np.allclose(v, w) == True or equivalent (i.e., abs(v - w) < eps or similar):
np.allclose([1, 0, 0], [1, 0, 0])
# True
np.allclose([1, 0, 0], [1 + 1e-5, 0, 0], atol=1e-5)
# True
np.allclose([1, 0, 0], [1 + 1e-4, 0, 0], atol=1e-5)
# False
So here is the question: how do I (efficiently) compare two (unordered) sets of vectors for equality, taking numerical approximation (e.g. np.allclose) into account?
Here's a function that works using np.lexsort:
def is_symmetric(A, R, *args, **kwargs):
A = np.asanyarray(A)
A = A[np.lexsort(A.T)]
A_t = A @ np.asanyarray(R).T
A_t = A_t[np.lexsort(A_t.T)]
return np.allclose(A, A_t, *args, **kwargs)
Some results:
R = np.array([[0, -1, 0], # 90° rotation as an example
[1, 0, 0],
[0, 0, 1]])
is_symmetric([[0, 0, 0]], R)
# True
is_symmetric([[1, 0, 0],
[0, 1, 0],
[0, 0, 0]], R)
# False
is_symmetric([[1, 0, 0],
[0, 1, 0],
[0, 0, 0],
[-1, 0, 0]], R)
# False
is_symmetric([[1, 0, 0],
[0, 1, 0],
[0, 0, 0],
[-1, 0, 0],
[0, -1, 0]], R)
# True
Performance seems fine for 100000 random vectors:
A = np.random.rand(100000, 3)
%timeit is_symmetric(A, R)
# 82.2 ms ± 75.4 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)