How can I remove duplicates from a list of matrices based on some equivalence relation?

Question

Given some list of symmetric integer matrices, I want to remove all duplicates under the following equivalence relation:

Two k x k matrices M1, M2 are equivalent if there is some permutation s on {1,...,k} such that for all i and j in {1,...,k} we have M1_ij = M2_s(i)s(j), i.e. two matrices are "equal" if I can get one from the other by simultaneously permuting its rows and columns.

Unfortunately, my naive approach (while building the list, check if any of the permutations of the new matrix is already in the list) proves to be too slow.

Some probably faster alternative I could think of would be putting all matrices in the list, permute them to some "canonical permutation" and then remove the duplicates as described e.g. here. However, I am not sure how to achieve such a "canonical permutation" in code either.

To narrow this down some more: The matrices are relatively small (k <= 4), the list will contain some 5 or 6 digit number of matrices and the dtype of the matrices has to be some integer type (currently intc, but I could change that).

The order of the final list does not matter, neither does which representative of each equivalence class survives. The entire process may take some small number of hours if needed, but not days.

Is there some reasonably efficient way to achieve this? Did I (yet again) miss some cool NumPy or SciPy facility that would help me with this?

---

As requested, some small examples to demonstrate how that equivalence relation works:

The matrices {{1,1,1},{1,2,0},{1,0,3}} and {{1,1,1},{1,3,0},{1,0,2}} are equivalent, as the permutation {1,2,3}->{1,3,2} transforms one to the other.

The matrices {{1,1,1},{1,2,0},{1,0,3}} and {{1,1,0},{1,2,1},{0,1,3}} are not equivalent, you cannot change the position of the 1s without permuting the diagonal.

Answer 1

This is an algebraic answer. I suspect there should be a more satisfying combinatoric answer.

You say that two matrix M and M' are equivalent if there exists a permutation matrix P such that M' = P^{-1} M P.

Let's use the eigen decomposition of M and M' :

M = Q^{-1} D Q

M = Q'^{-1} D' Q'

Where D and D' are diagonal matrices containing the eigenvalues, and Q and Q' are orthogonal matrices.

We can rewrite the equality as:

D = D' up to permutation (i.e. the two matrices should have the same eigenvalues)

Q' = PQ

Testing for the second condition is easy. Given that Q is orthogonal, it amounts to checking whether the matrix dot(Q, Q'.T) is a permutation matrix, i.e. if it has only one "1" per row and column.

---

A draft for an algorithm is thus:

• Take M and M'
• Compute the eigen decomposition (Q, D) and (Q', D') of M and M' (using np.linalg.eigh)
• If they do not have the same eigenvalues (up to numerical precision of course), they are not equivalent
• Else, compute np.dot(Q, Q'.T) and test if it is a permutation matrix

I think that the bottleneck is the eigen decomposition, but you only have to do it once per matrix. Hopefully the first test will discard a lot of matrices quickly.

Hope this helps.

Answer 2

You can think of your matrices as representing the adjacency/weight matrix of a graph, and then test if the two graphs are isomorphic to each other. networkx has a convenient function (and can be installed via pip).

import numpy as np
import networkx as nx
from networkx.algorithms.isomorphism import numerical_edge_match

# create matrices
n = 4
a = np.random.randint(0, 10, size=(n,n))
a = a + a.T # i.e. symmetric
b = np.rot90(a, k=2) # i.e. a rotated by 180 degrees
c = np.ones((n,n), dtype=np.int) # counter-example

# create graphs
ga = nx.from_numpy_matrix(a)
gb = nx.from_numpy_matrix(b)
gc = nx.from_numpy_matrix(c)

# test if isomorphic
print "a isomorphic with b:", nx.is_isomorphic(ga, gb, edge_match=numerical_edge_match('weight', 1)) # True
print "a isomorphic with c:", nx.is_isomorphic(ga, gc, edge_match=numerical_edge_match('weight', 1)) # False

Answer 3

Just use your canonical approach. Search for the largest entry in the matrix, put it on the top right corner. Sort then the first colum and first row according it's entries.

A = np.array([[1,2,3,5],
     [3,6,2,6],
     [3,5,7,2],
     [1,3,6,3]])
a = np.where(A == np.amax(A))
sort_colums = np.argsort(A[a[0]].ravel())[::-1]
sort_rows = np.argsort(A[:,a[1]].ravel())[::-1]
Col_sorted = A[:,sort_colums]
Equiv_class = Col_sorted[sort_rows]
#returns [[7, 5, 3, 2],
          [6, 3, 1, 3],
          [3, 2, 1, 5],
          [2, 6, 3, 6]]

As pointed out in the comments, this only works, when the entries of the matrix occur only once. In case the do occur several times, but not that often, then one could adapt this method by generating several equivalence classes.