How many $n\times m$ binary matrices are there, up to row and column permutations?

I'm interested in the number of binary matrices of a given size that are distinct with regard to row and column permutations.

If $\sim$ is the equivalence relation on $n\times m$ binary matrices such that $A \sim B$ iff one can obtain B from applying a permutation matrix to A, I'm interested in the number of $\sim$-equivalence classes over all $n\times m$ binary matrices.

I know there are $2^{nm}$ binary matrices of size $n\times m$, and $n!m!$ possible permutations, but somehow I fail to get an intuition on what this implies for the equivalence classes.

-
Intuitively, the average matrix has trivial stabilizer, so there ought to be roughly 2^{nm}/n!m! equivalence classes. This is probably a very hard question in general. – Qiaochu Yuan Feb 15 '11 at 12:02
There is the set $S:=[m]\times[n]$ on which the group $G:=S_m\times S_n$ acts, and we have to color $S$ with colors $0$ and $1$. How many colorings are there when two colorings that differ by a $g\in G$ are considered the same? Now there is a famous theory that addresses exactly this kind of questions; it is called Polya counting theory. I could imagine that your problem is a standard example in the field. – Christian Blatter Feb 15 '11 at 13:24
If you view A as the incidence matrix of an unweighted undirected bipartite graph. Then I think the question you're asking is how many unique bipartite graphs up to isomorphism are there with the two vertex groups having n,m vertices respectively. – JSchlather Feb 15 '11 at 17:53

This is solved here using Pólya enumeration theory. For the square case ($n=m$), see this sequence.

Comment: I found these by searching for $1,2,7$ in the OEIS.

-
The first is broken now but I would love to know the answer. Is there an alternative source? – Raphael Feb 21 '14 at 20:38
@Raphael Updated the link. – Yuval Filmus Feb 21 '14 at 20:46

Here is a computational contribution that treats the case of a square matrix. As pointed out this problem can be solved using the Polya Enumeration Theorem. We just need to compute the cycle index of the group acting on the slots of the matrix, subsitute $1+z$ into the cycle index to get the generating function and finally set $z=1$ to obtain the count.

These cycle indices are easy to compute and we do not need to iterate over all $(n!)^2$ pairs of permutations (acting on rows and columns) but instead it is sufficient to iterate over pairs of terms from the cycle index $Z(S_n)$ of the symmetric group $S_n$ according to their multiplicities to obtain the cycle index $Z(Q_n)$ of the combined action on rows and columns. The number of terms here is the much better count of the number of partitions of $n$ squared (upper bound).

The algorithm now becomes very simple -- iterate over pairs of terms as described above, compute two representative permutations, let them permute rows and columns, factor the result into cycles and add it to the cycle index being computed.

This gives the following cycle indices (only the first five are shown): $$Z(Q_2) = 1/4\,{a_{{1}}}^{4}+3/4\,{a_{{2}}}^{2},$$ $$Z(Q_3) = 1/36\,{a_{{1}}}^{9}+1/6\,{a_{{1}}}^{3}{a_{{2}}}^{3} +2/9\,{a_{{3}}}^{3}+1/4\,a_{{1}}{a_{{2}}}^{4}+1/3\,a_{{3}}a_{{6}},$$ $$Z(Q_4) = {\frac {{a_{{1}}}^{16}}{576}}+1/48\,{a_{{1}}}^{8}{a_{{2}}}^ {4}+1/36\,{a_{{1}}}^{4}{a_{{3}}}^{4}+{\frac {17\,{a_{{2}}}^ {8}}{192}}\\+{\frac {13\,{a_{{4}}}^{4}}{48}}+1/16\,{a_{{1}}}^ {4}{a_{{2}}}^{6}+1/6\,{a_{{1}}}^{2}{a_{{3}}}^{2}a_{{2}}a_{{ 6}}\\+1/9\,a_{{1}}{a_{{3}}}^{5}+1/12\,{a_{{2}}}^{2}{a_{{6}}}^ {2}+1/6\,a_{{4}}a_{{12}}$$ and $$Z(Q_5) = 1/12\,{a_{{1}}}^{2}{a_{{4}}}^{2}a_{{3}}a_{{12}}+1/24\,{a_{{ 1}}}^{3}{a_{{4}}}^{5}a_{{2}}+1/36\,{a_{{2}}}^{5}{a_{{3}}}^{ 3}a_{{6}}\\+1/24\,{a_{{1}}}^{2}{a_{{2}}}^{4}a_{{3}}{a_{{6}}}^ {2}+{\frac {{a_{{1}}}^{5}{a_{{4}}}^{5}}{240}}+{\frac {{a_{{ 2}}}^{5}{a_{{3}}}^{5}}{360}}+1/36\,{a_{{1}}}^{4}{a_{{3}}}^{ 7}\\+1/36\,{a_{{1}}}^{6}{a_{{3}}}^{3}{a_{{2}}}^{2}a_{{6}}+{ \frac {a_{{1}}{a_{{2}}}^{12}}{64}}+1/48\,{a_{{1}}}^{3}{a_{{ 2}}}^{11}+{\frac {{a_{{1}}}^{10}{a_{{3}}}^{5}}{360}}\\+{ \frac {{a_{{1}}}^{5}{a_{{2}}}^{10}}{480}}+{\frac {{a_{{1}}} ^{9}{a_{{2}}}^{8}}{144}}+{\frac {{a_{{1}}}^{15}{a_{{2}}}^{5 }}{720}}+{\frac {{a_{{1}}}^{25}}{14400}}\\+1/10\,a_{{5}}a_{{ 20}}+1/15\,a_{{10}}a_{{15}}+1/18\,{a_{{2}}}^{2}{a_{{3}}}^{5 }a_{{6}}+1/16\,a_{{1}}{a_{{4}}}^{5}{a_{{2}}}^{2}\\+1/24\,{a_{ {2}}}^{5}a_{{3}}{a_{{6}}}^{2}+{\frac {13\,{a_{{5}}}^{5}}{ 300}}+1/12\,a_{{2}}a_{{3}}{a_{{4}}}^{2}a_{{12}}+1/30\,{a_{{ 5}}}^{3}a_{{10}}\\+1/15\,{a_{{5}}}^{2}a_{{15}}+1/20\,a_{{5}}{ a_{{10}}}^{2}+1/16\,a_{{1}}{a_{{4}}}^{6}+1/36\,{a_{{2}}}^{2 }{a_{{6}}}^{2}{a_{{3}}}^{3}.$$

Substituting $1+z$ into these cycle indices we quickly obtain the sequence (with a very reasonable time complexity) $$2, 7, 36, 317, 5624, 251610, 33642660, 14685630688,\\ 21467043671008, 105735224248507784,1764356230257807614296,\\ 100455994644460412263071692,19674097197480928600253198363072,\\ 13363679231028322645152300040033513414,\\ 31735555932041230032311939400670284689732948,\ldots$$ which is indeed OEIS A002724.

Note that the cycle indices make it possible to enumerate configurations with more than two possible entries or with entries having different weights. For example, with a 3x3 square and three colors $A,B$ and $C$ we get the generating function $$Z(Q_3)(A+B+C) = 1/36\, \left( A+B+C \right) ^{9}+1/6\, \left( A+B+C \right) ^{3} \left( {A}^{2}+{B}^{2}+{C}^{2} \right) ^{ 3}\\+2/9\, \left( {A}^{3}+{B}^{3}+{C}^{3} \right) ^{3}+1/ 4\, \left( A+B+C \right) \left( {A}^{2}+{B}^{2}+{C}^{2 } \right) ^{4}\\+1/3\, \left( {A}^{3}+{B}^{3}+{C}^{3} \right) \left( {A}^{6}+{B}^{6}+{C}^{6} \right)$$ which expands to $${A}^{9}+{A}^{8}B+{A}^{8}C+3\,{A}^{7}{B}^{2}+3\,{A}^{7}B C+3\,{A}^{7}{C}^{2}+6\,{A}^{6}{B}^{3}+10\,{A}^{6}{B}^{2 }C\\+10\,{A}^{6}B{C}^{2}+6\,{A}^{6}{C}^{3}+7\,{A}^{5}{B}^ {4}+17\,{A}^{5}{B}^{3}C+28\,{A}^{5}{B}^{2}{C}^{2}\\+17\,{ A}^{5}B{C}^{3}+7\,{A}^{5}{C}^{4}+7\,{A}^{4}{B}^{5}+22\, {A}^{4}{B}^{4}C+43\,{A}^{4}{B}^{3}{C}^{2}+43\,{A}^{4}{B }^{2}{C}^{3}\\+22\,{A}^{4}B{C}^{4}+7\,{A}^{4}{C}^{5}+6\,{ A}^{3}{B}^{6}+17\,{A}^{3}{B}^{5}C+43\,{A}^{3}{B}^{4}{C} ^{2}+54\,{A}^{3}{B}^{3}{C}^{3}\\+43\,{A}^{3}{B}^{2}{C}^{4 }+17\,{A}^{3}B{C}^{5}+6\,{A}^{3}{C}^{6}+3\,{A}^{2}{B}^{ 7}+10\,{A}^{2}{B}^{6}C+28\,{A}^{2}{B}^{5}{C}^{2}\\+43\,{A }^{2}{B}^{4}{C}^{3}+43\,{A}^{2}{B}^{3}{C}^{4}+28\,{A}^{ 2}{B}^{2}{C}^{5}+10\,{A}^{2}B{C}^{6}+3\,{A}^{2}{C}^{7}\\+ A{B}^{8}+3\,A{B}^{7}C+10\,A{B}^{6}{C}^{2}+17\,A{B}^{5}{ C}^{3}+22\,A{B}^{4}{C}^{4}+17\,A{B}^{3}{C}^{5}\\+10\,A{B} ^{2}{C}^{6}+3\,AB{C}^{7}+A{C}^{8}+{B}^{9}+{B}^{8}C+3\,{ B}^{7}{C}^{2}+6\,{B}^{6}{C}^{3}+7\,{B}^{5}{C}^{4}\\+7\,{B }^{4}{C}^{5}+6\,{B}^{3}{C}^{6}+3\,{B}^{2}{C}^{7}+B{C}^{ 8}+{C}^{9}.$$

This is the Maple code for this computation.


pet_cycleind_symm :=
proc(n)
local p, s;
option remember;

if n=0 then return 1; fi;

end;

pet_autom2cycles :=
proc(src, aut)
local numa, numsubs;
local marks, pos, cycs, cpos, clen;

numsubs := [seq(src[k]=k, k=1..nops(src))];
numa := subs(numsubs, aut);

marks := Array([seq(true, pos=1..nops(aut))]);

cycs := []; pos := 1;

while pos <= nops(aut) do
if marks[pos] then
clen := 0; cpos := pos;

while marks[cpos] do
marks[cpos] := false;
cpos := numa[cpos];
clen := clen+1;
od;

cycs := [op(cycs), clen];
fi;

pos := pos+1;
od;

return mul(a[cycs[k]], k=1..nops(cycs));
end;

pet_varinto_cind :=
proc(poly, ind)
local subs1, subs2, polyvars, indvars, v, pot, res;

res := ind;

polyvars := indets(poly);
indvars := indets(ind);

for v in indvars do
pot := op(1, v);

subs1 :=
[seq(polyvars[k]=polyvars[k]^pot,
k=1..nops(polyvars))];

subs2 := [v=subs(subs1, poly)];

res := subs(subs2, res);
od;

res;
end;

pet_flatten_term :=
proc(varp)
local terml, d, cf, v;

terml := [];

cf := varp;
for v in indets(varp) do
d := degree(varp, v);
terml := [op(terml), seq(v, k=1..d)];
cf := cf/v^d;
od;

[cf, terml];
end;

pet_multiset_term :=
proc(varp)
local terml, d, cf, v;

terml := [];

cf := varp;
for v in indets(varp) do
d := degree(varp, v);
terml := [op(terml), [v, d]];
cf := cf/v^d;
od;

[cf, terml];
end;

pet_flat2rep :=
proc(f)
local p, q, res, cyc, t, len;

q := 1; res := [];

for t in f do
len := op(1, t);
cyc := [seq(p, p=q+1..q+len-1), q];
res := [op(res), seq(cyc[p], p=1..nops(cyc))];
q := q+len;
od;

res;
end;

pet_cycleind_sqmat :=
proc(n)
option remember;
local cind, sind, t1, t2, src, autom, p, q,
flat1, flat2, rep1, rep2;

cind := 0;

if n=1 then
sind := [a[1]];
else
sind := pet_cycleind_symm(n);
fi;

src := [seq(seq([p, q], q=1..n), p=1..n)];

for t1 in sind do
flat1 := pet_flatten_term(t1);
rep1 := pet_flat2rep(flat1[2]);

for t2 in sind do
flat2 := pet_flatten_term(t2);
rep2 := pet_flat2rep(flat2[2]);

autom :=
[seq([rep1[src[q][1]],
rep2[src[q][2]]], q=1..n*n)];

cind := cind + flat1[1]*flat2[1]*
pet_autom2cycles(src, autom);
od;
od;

cind;
end;

pet_cycleind_sqmat2 :=
proc(n)
option remember;
local cind, sind, t1, t2, q, cyc1, cyc2,
flat1, flat2, len, len1, len2;

cind := 0;

if n=1 then
sind := [a[1]];
else
sind := pet_cycleind_symm(n);
fi;

for t1 in sind do
flat1 := pet_flatten_term(t1);

for t2 in sind do
flat2 := pet_flatten_term(t2);

q := 1;

for cyc1 in flat1[2] do
len1 := op(1, cyc1);

for cyc2 in flat2[2] do
len2 := op(1, cyc2);

len := lcm(len1, len2);
q := q * a[len]^(len1*len2/len);
od;
od;

cind := cind + flat1[1]*flat2[1]*q;
od;
od;

cind;
end;

v :=
proc(n)
option remember;
local cind, gf;

cind := pet_cycleind_sqmat(n);
gf := pet_varinto_cind(1+z, cind);

subs(z=1, gf);
end;

w :=
proc(n)
option remember;
local res, sind, t1, t2, p1, p2,
mset1, mset2, cf;

res := 0;

if n=1 then
sind := [a[1]];
else
sind := pet_cycleind_symm(n);
fi;

for t1 in sind do
mset1 := pet_multiset_term(t1);

for t2 in sind do
mset2 := pet_multiset_term(t2);

cf := 0;

for p1 in mset1[2] do
for p2 in mset2[2] do
cf := cf +
gcd(op(1, p1[1]), op(1, p2[1]))*
p1[2]*p2[2];
od;
od;

res := res + mset1[1]*mset2[1]*2^cf;
od;
od;

res;
end;


Remark. Observe that the function w provides an implementation of the closed form computation given in the OEIS entry. This is an important check that confirms the correctness of the generating function results. Note also the function pet_cycleind_sqmat2 which implements the classic approach to the calculation of these cycle indices based on the simple observation that for two cycles, one of length $l_1$ from a row permutation $\alpha$ and another of length $l_2$ from a column permutation $\beta$ their contribution to the disjoint cycle decomposition product for $(\alpha,\beta)$ in the cycle index $Z(Q_n)$ is by inspection $$a_{\mathrm{lcm}(l_1, l_2)}^{l_1 l_2 / \mathrm{lcm}(l_1, l_2)} = a_{\mathrm{lcm}(l_1, l_2)}^{\gcd(l_1, l_2)}.$$

This MSE Meta Link has many more PET computations by various users.

-