# Listing matrices up to symmetry

I am interested in the equivalence relation on N x N binary matrices, in which two matrices are equivalent if one can be obtained by rotating/reflecting the other. I would like to obtain a list

binclasses = { {B1, m1},...{Bk, mk} },

where Bi is a representative element of class i, and mi is the size of that class. (Of course the requirement of binary matrices is somewhat artificial here, so the question could be rephrased for matrices taking values in a finite set). For N =2, I am looking to obtain the list

binclasses = {{{{0, 0}, {0, 0}}, 1}, {{{0, 0}, {0, 1}}, 4}, {{{0, 0}, {1, 1}}, 4}, {{{0, 1}, {1, 0}}, 2}, {{{0, 1}, {1, 1}}, 4}, {{{1, 1}, {1, 1}}, 1}}

My thought process so far has been as follows:

1. Define a list binlist of all N x N binary matrices.
2. Choose an element currmat from binlist, and generate a list currsym of all rotations and reflections of currmat.
3. Define a new list binminsym which is binlist with all elements of currsym removed.
4. Add the pair {currmat, Length[currsym]} to a new list binclasses.
5. Set binlist = binminsym, return to step 1.

I am guessing that this is not a particularly efficient approach, and am curious to know how one could better go about the problem.

EDIT A second update, I have improved on the original code, however, it still runs slowly (original edit left below this for information). In this variation I first group binary matrices by the number of 1s, and the number of adjacent 1s; this means that fewer list elements need be checked against each other. At the moment for the 4 x 4 case, this version takes a little over 4 minutes to identify the 8548 equivalence classes, compared to the 17 minutes of the previous method.

(* sTotal function gives total of all matrix entries *)
sTotal[matrix_] :=
Sum[matrix[[i, j]], {i, 1, Length[matrix]}, {j, 1, Length[matrix]}];

(* eTotal function gives total number of adjacent non-zero entries *)

halo[matrix_] :=
Module[{cmat},
cmat = ConstantArray[0, {Length[matrix] + 1, Length[matrix] + 1}];
cmat[[1 ;; Length[matrix], 1 ;; Length[matrix] ]] = matrix;
cmat]

eTotal[matrix_] :=
Sum[halo[matrix][[i, j]] (halo[matrix][[i, j + 1]] +
halo[matrix][[i + 1, j]]), {i, 1, Length[matrix]}, {j, 1,
Length[matrix]}];

Timing[
Clear[mdim, blist, bgrlist, binclasses, binreps, currbg, currmat,
currsym]

mdim = 4; (*matrix dimensions*)
blist = Tuples[{0, 1}, {mdim,
mdim}];  (* All square binary matrices of size mdim *)

(* blist grouped into sets of matrices with equal number of 1s, and \
bgrlist = GatherBy[blist, {sTotal, eTotal}];

(* Initialise list of equivalence classes *)
binclasses = {};
binreps = {};
]

Timing[
Do[
Do[
currbg = bgrlist[[i, j]];
While[Length[currbg] >= 1,
currmat = currbg[];
currsym = symlist[currmat];
binclasses = Append[binclasses, {currmat, Length[currsym]}];
binreps = Append[binreps, currmat];
If[i != 1 + (mdim^2 /2 ),
binreps =
Append[binreps, ConstantArray[1, {mdim, mdim}] - currmat];
];
currbg = DeleteCases[currbg, x_ /; MemberQ[currsym, x] == True];
]
, {j, 1, Length[bgrlist[[i]]]}];
, {i, 1, Ceiling[Length[bgrlist]/2]}];
Length[binreps]
]


EDIT Here is my code for the above algorithm. Note, that I have defined the symmetries of the square as functions, even though some of them (eg. reflection in diagonal) are simply transposing the matrix. I did this for clarity (for readers on here), rather than efficiency!

(*Symmetry functions.*)
vsym[matrix_] := Reverse[matrix, {2}];
hsym[matrix_] := Reverse[matrix];
d1sym[matrix_] := Transpose[matrix];
r1sym[matrix_] := vsym[d1sym[matrix]];
d2sym[matrix_] := r1sym[vsym[matrix]];

(*List of all distinct symmetries of given matrix *)
symlist[matrix_] :=
DeleteDuplicates[{matrix, vsym[matrix], hsym[matrix], d1sym[matrix],
d2sym[matrix], r1sym[matrix], r1sym[r1sym[matrix]],
r1sym[r1sym[r1sym[matrix]]]}];

(* All square binary matrices of size mdim *)
binlist[mdim_] := binlist = Tuples[{0, 1}, {mdim, mdim}];

(* Initialise list of equivalence classes *)
binclasses = {};

binlist;

While[ Length[binlist] >= 1,
(currmat = binlist[];
currsym = symlist[currmat];
binclasses = Append[binclasses, {currmat, Length[currsym]}];
binlist = DeleteCases[binlist, x_ /; MemberQ[currsym, x] == True]);
]

• Welcome to SE! You'll probably get more feedback if you also post the code you have so far. Aug 6 '13 at 8:35
• This may help, but I don't have time right now to write up a full answer here. Aug 6 '13 at 13:11
• @Owen Daniel Is there some constraint that my answer did not satisfy? Aug 13 '13 at 13:27
• @Hector. No you matched the constraints, but as I pointed out, my code runs significantly faster than yours did, I'm waiting to see whether somebody can actually provide a well optimized version. Aug 13 '13 at 16:31
• @OwenDaniel Faster than 4 seconds for the 4x4 case? I just cannot believe it. Aug 13 '13 at 18:31

SECOND ATTEMPT - 2nd Edition

The following functions are defined as pure functions for speed.

toMatrix = Partition[IntegerDigits[#, 2, 16], 4] &;
toInteger = FromDigits[Flatten[#], 2] &;
allSyms = Module[{s1, s3, s4, s6}, {#, s1 = Reverse[#, {2}], Reverse[#],
s3 = Transpose[#], s4 = Reverse[s3, {2}],
Reverse[Transpose[s1], {2}], s6 = Reverse[Transpose[s4], {2}],
Reverse[Transpose[s6], {2}]}] &;


You can check what they do using for example:

MatrixForm /@ allSyms[toMatrix]


Finally, the following code gives you the answer in 4 seconds 3.8 seconds!

casesToCheck = Range[0, 65535];
Length[#]} & /@ (Reap[
NestWhile[
Complement[#,
Sow[Union[toInteger /@ allSyms[toMatrix[#[]]]]]] &,
casesToCheck, (n = Length[#]) > 0 &]][[2, 1]]); Length[answer]]
(*{3.796875, 8548}*)


You can check the answers now. For example, the first 3 are answer[[1 ;; 3]]: The classes with only one element are Cases[answer, {_, 1}, {1}] : Old slow solution

Have you tried the following?

matrices = Partition[IntegerDigits[#, 2, 4], 2] & /@ Range[0, 15];;
vsym[matrix_] := Reverse[matrix, {2}];
hsym[matrix_] := Reverse[matrix];
d1sym[matrix_] := Transpose[matrix];
r1sym[matrix_] := vsym[d1sym[matrix]];
d2sym[matrix_] := r1sym[vsym[matrix]];
Tally[matrices,
vsym[#1] == #2 || hsym[#1] === #2 || d1sym[#1] == d1sym[#2] ||
r1sym[#1] == #2 || d2sym[#1] == #2&] • Thanks for the response. There seem to be a couple of errors in the code, and you're missing two of the symmetries (r2 and r3, which are r1 applied to themselves). This seems to be quite a bit slower than the code I posted above when doing the 4 x 4 case, I guess because Tally is comparing every pair. Thanks again. Aug 8 '13 at 8:08
• Tally does not make all possible comparisons. For example, Tally[Array[a[#, Mod[#, 2]] &, {4}], (Print[#1, #2];Last[#1] === Last[#2]) &] will print only 4 comparisons. However, Tally does not remove from the pool as you do. As for the missing symmetries, you can put them easily into the code. Aug 8 '13 at 18:40

Here is an implementation of the change I suggested to Hector's code:

toMatrix = Partition[IntegerDigits[#, 2, 16], 5] &;
toInteger = FromDigits[Flatten[#], 2] &;
allSyms =
Module[{s1, s3, s4, s6}, {#, s1 = Reverse[#, {2}], Reverse[#],
s3 = Transpose[#], s4 = Reverse[s3, {2}],
Reverse[Transpose[s1], {2}], s6 = Reverse[Transpose[s4], {2}],
Reverse[Transpose[s6], {2}]}] &;
casesToCheck = Range[0, 2^25];

Timing[answer = {MatrixForm@toMatrix@First[#], Length[#]} & /@ (Reap[
With[{syms = Union[toInteger /@ (allSyms@toMatrix[#])]},
If[# == syms[], Sow[syms]];
] & /@ casesToCheck

(*{6.593587, 8548}*)


For comparison: Four seconds on Hector's machine corresponds to about six seconds on mine. This solution is a little slower than Hector's, but it should scale better to the 5x5 case. I was surprised not to get a speed-up; this is a linear-time algorithm, but apparently the cost of calling Union[toInteger/@(allSyms@toMatrix[#])] for every single matrix is enough to nullify any advantage in the 4x4 case.

I thought it would be interesting to implement a few other (still linear-time) solutions using Hector's matrix functions and compare performance.

Here is the same solution, but with allSyms memoized (i.e. the integers corresponding to the matrices symmetric to the matrix corresponding to a given integer are remembered in the definition of the function sortedSymInts):

clearSortedSymInts[] := (
Clear[sortedSymInts];
sortedSymInts[matint_] :=
With[{foo = Union[toInteger /@ allSyms[toMatrix[matint]]]},
(sortedSymInts[#] = foo) & /@ foo;
foo
]
);

clearSortedSymInts[];
Length[#]} & /@ (Reap[
With[{syms = sortedSymInts[#]},
If[# == syms[], Sow[syms]];
] & /@ casesToCheck
]

(*{1.736454, 8548}*)


Now it's much faster, at the expense of requiring more memory.

Here are two reimplementations of my previous solution, using Hector's matrix functions. The first memoizes integers representing matrices symmetric to matrices which have already been seen. The second memoizes the matrices themselves, and only uses toMatrix to generate the initial list.

First:

clearNewSym[] := (
Clear[newSym];
newSym[
matint_] := ((newSym[#] = False) & /@ (sortedSymInts[matint]);
True)
);

clearSortedSymInts[];
clearNewSym[];
Timing[binclasses = {MatrixForm[toMatrix[#]],
Length[sortedSymInts[#]]} & /@ Select[casesToCheck, newSym];
Length[binclasses]]

(*{1.534411, 8548}*)


Second:

clearNewerSym[] := (
);
clearSortedSymInts[];
Timing[binclasses2 = {MatrixForm[#], Length[Union[allSyms[#]]]} & /@

(*{1.241104, 8548}*)


Surprisingly (to me), the last is the fastest. Apparently it is faster to memoize the matrices than to call toInteger and then memoize the integers. The hash function used by Mathematica's rule processing code must be very fast for lists (or perhaps slow for integers???).

• (1) The 3.98 seconds were in my little Surface Pro. I got 6 seconds too in my 6-core AMD 2.70Ghz desktop. (2) It is interesting that hashing a matrix and looking up the hashed matrix is faster than transforming the matrix, hashing the integer, and looking up the integer. (3) It would be interesting to see if compiling toInteger makes any difference. Aug 15 '13 at 22:20

an[u_] := Module[{n, perm, p, tp, ptp}, n = Length[u];
perm =
Flatten[Outer[List, Permutations@Range[n], Permutations@Range[n],
1], 1];
p = Map[u[[#[], #[]]] &, perm];
tp = Transpose /@ p;
ptp = Union[p, tp]
];

eqc = Union /@ Gather[an /@ test];


Displaying data:

{Map[MatrixForm, #, {2}], Length @@ #} & /@eqc // Grid • Hi, thanks for your answer. There is a small typo in display result, where you meant an2. Your code appears to be counting the number of distinct symmetries for each matrix, however is not sorting them into equivalence classes (I don't think!). I think that the sorting stage is the computationally heavy part! Thanks again. Aug 8 '13 at 8:24
• @OwenDaniel Thank you. Sorry for the typo. I misinterpreted the question but now it is a starting point for me. Sorry. Aug 8 '13 at 8:31
• Not right yet but getting closer Aug 8 '13 at 9:51
• @OwenDaniel let me know if the above revision works. Aug 8 '13 at 13:54
• Hi. Just tried it out on the 4 x 4 case and I got: "No more memory available. Mathematica kernel has shut down. Try quitting other applications and then retry."! Will try again on another computer later! Aug 8 '13 at 14:45

The following can replace your second Timing block. The memoizing function newsymrep returns True if no matrix equivalent to the input has previously been passed as input.

symlist[matrix_] := {matrix, vsym[matrix], hsym[matrix],
d1sym[matrix], d2sym[matrix], r1sym[matrix], r1sym[r1sym[matrix]],
r1sym[r1sym[r1sym[matrix]]]};

Clear[newsymrep];
newsymrep[matrix_] :=
With[{foo = symlist[matrix]}, (newsymrep[#] = False) & /@ foo;
True];

Timing[
binreps = Select[blist, newsymrep];
binclasses = {MatrixForm[#], Length@Union@symlist[#]} & /@ binreps;
Length[binreps]]

(*{1.513778, 8548}*)


I removed DeleteCases from symlist; as symlist is used here, it's faster to redundantly memoize.

• Your solution seems to lose information about the size of each class. Aug 14 '13 at 6:06
• Sorry. It should be correct now. Aug 14 '13 at 17:00