Counting directed bicliques using Burnside's lemma

Question

Let $b_{n}$ be the number of different directed $K_{n,n}$ graphs, assuming that $G$ and $H$ are considered identical when $G$ is isomorphic either with $H$ or with its transpose $H^T$ (i.e. same graph with all edges reversed).

I'm trying to use Burnside's lemma to calculate $b_3$. The symmetry group seems to be $S_3 \oplus S_3 \oplus \mathbb{Z}_2 \oplus \mathbb{Z}_2$ (permutations of vertices on each side, horizontal reflection, and reversing the edges).

It's easy to find $b_1 = 1$ and $b_2 = 4$. I've written a brute-force program that outputs $b_3 = 13$. However, when trying to calculate it with Burnside's lemma, I'm getting $\frac{2112}{144} = 14 \frac{2}{3}$. Actually, I've also written a program that uses Burnside's lemma to calculate it, and I've checked this result manually.

The calculation goes as follows:

$(2^9+2^6) + 3*(2^6+2^3)*2 + 2*(2^3+2^2)*2 + 3*3*(2^5+2^6) + 2*3*(2^2+2^3)*2 + 2*2*(2^3+2^2) = 2112$

The factors outside of the parentheses count different types of vertex permutations on both sides of the graph, and the numbers inside the parentheses correspond to fixed points for these permutations, accordingly, without or with both horizontal reflection and edge reversal (there are no fixed points when either of them is used alone, due to odd-length cycles).

Do you have any hint what might be wrong here?

Edit:

The sequence is similar to A091059, but it's not the same. My program can successfully calculate this sequence (up to the 14th element, then it's getting too slow) when I tell it not to take horizontal reflection into account. This suggests something might be wrong with the reflection, I just still can't see why.

Let's consider directed $K_{2,2}$ graphs. A091059 says there are $5$ such graphs, and we can draw them:

But the last two graphs are isomorphic – they are just mirror reflections of each other. So $b_2 = 4$.

Answer 1

What we have here is an instance of Power Group Enumeration as described by Harary and Palmer, Graphical Enumeration. The algorithm is documented at the following MSE-link I. We require the cycle index $Z(Q_n)$ of the action on the edges of the permutations of the two parts of the graph, possibly combined with a horizontal reflection. This is the slot permutation group. We distribute edges of one of $k$ colors into these slots, and the group acting on them is the symmetric group with cycle index $Z(S_k)$. The cycle index $Z(Q_n)$ was computed at the following MSE-link II. We have e.g.

$$Z(Q_3) = {\frac {{a_{{1}}}^{9}}{72}} +1/6\,{a_{{1}}}^{3}{a_{{2}}}^{3} +1/8\,a_{{1}}{a_{{2}}}^{4}+1/4\,a_{{1}}{a_{{4}}}^{2} +1/9\,{a_{{3}}}^{3}+1/3\,a_{{3}}a_{{6}}.$$

and

$$Z(Q_4) = {\frac {{a_{{1}}}^{16}}{1152}} +{\frac {{a_{{1}}}^{8}{a_{{2}}}^{4}}{96}} +{\frac {5\,{a_{{1}}}^{4}{a_{{2}}}^{6}}{96}} +{\frac {{a_{{1}}}^{4}{a_{{3}}}^{4}}{72}} +{\frac {17\,{a_{{2}}}^{8}}{384}} \\ +1/12\,{a_{{1}}}^{2}a_{{2}}{a_{{3}}}^{2}a_{{6}} +1/8\,{a_{{1}}}^{2}a_{{2}}{a_{{4}}}^{3} +1/18\,a_{{1}}{a_{{3}}}^{5} +1/6\,a_{{1}}a_{{3}}{a_{{6}}}^{2} \\ +1/24\,{a_{{2}}}^{2}{a_{{6}}}^{2} +{\frac {19\,{a_{{4}}}^{4}}{96}} +1/12\,a_{{4}}a_{{12}}+1/8\,{a_{{8}}}^{2}.$$

With these ingredients we are ready to run the PGE algorihm. We get for two swappable types of edges the sequence

$$1, 4, 13, 104, 1507, 64203, 8426875, 3671999389, 5366787092478, \\ 26433809041087192, 441089058039611200394, 25113998661290096278734134, \ldots$$

and for three types

$$1, 6, 84, 7946, 5413511, 25231086540, 800871112032930, \\ 177544715836044855636, 281653040526999655665449719, \\ 3266495639384107667257990172349726, \\ 282129919925994006382238965837655927175534, \\ 184379837924757642947198903200667422197524750679153, \ldots $$

The Maple code for this is quite compact and shown below.

with(combinat);

pet_cycleind_symm :=
proc(n)
local l;
option remember;

    if n=0 then return 1; fi;

    expand(1/n*add(a[l]*pet_cycleind_symm(n-l), l=1..n));
end;

pet_cycleind_knn :=
proc(n)
option remember;
local cindA, cindB, sind, t1, t2, term, res,
    cmb, len, l1, l2, cycs, uidx, vidx,
    u, v, inst1;

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

    cindA := 0;

    for t1 in sind do
        for t2 in sind do
            res := 1;

            for u in indets(t1) do
                l1 := op(1, u);

                for v in indets(t2) do
                    l2 := op(1, v);

                    len := lcm(l1, l2);
                    res := res *
                    a[len]^(degree(t1, u)*degree(t2, v)
                            *l1*l2/len);
                od;
            od;

            cindA := cindA + lcoeff(t1)*lcoeff(t2)*res;
        od;
    od;

    cindB := 0;

    for term in sind do
        res := 1;

        # edges on different cycles of different sizes
        for cmb in choose(indets(term), 2) do
            u := op(1, cmb); v := op(2, cmb);

            l1 := 2*op(1, u); l2 := 2*op(1, v);
            res := res *
            a[lcm(l1, l2)]^((l1*l2/2/lcm(l1, l2))*
                            degree(term, u)*degree(term, v));
        od;

        # edges on different cycles of the same size
        for u in indets(term) do
            l1 := 2*op(1, u); inst1 := degree(term, u);
            # a[l1]^(1/2*inst1*(inst1-1)*l1*l1/2/l1)
            res := res *
            a[l1]^(1/2*inst1*(inst1-1)*l1/2);
        od;

        # edges on identical cycles of some size
        for u in indets(term) do
            l1 := 2*op(1, u); inst1 := degree(term, u);
            if type(l1/2, even) then
                # a[l1]^((l1/2)^2/l1);
                res := res *
                (a[l1]^(l1/4))^inst1;
            else
                # a[l1/2]^(l1/2/(l1/2))*a[l1]^(((l1/2)^2-l1/2)/l1)
                res := res *
                (a[l1/2]*a[l1]^(l1/4-1/2))^inst1;
            fi;
        od;


        cindB := cindB + lcoeff(term)*res;
    od;

    (cindA+cindB)/2;
end;

knn_swap_edge_cols :=
proc(n,k)
option remember;
local idx_slots, idx_cols, res, term_a, term_b,
    v_a, v_b, inst_a, inst_b, len_a, len_b, p, q;

    if n = 1 then
        idx_slots := [a[1]];
    else
        idx_slots := pet_cycleind_knn(n);
    fi;

    if k = 1 then
        idx_cols := [a[1]];
    else
        idx_cols := pet_cycleind_symm(k);
    fi;

    res := 0;

    for term_a in idx_slots do
        for term_b in idx_cols do
            p := 1;

            for v_a in indets(term_a) do
                len_a := op(1, v_a);
                inst_a := degree(term_a, v_a);

                q := 0;

                for v_b in indets(term_b) do
                    len_b := op(1, v_b);
                    inst_b := degree(term_b, v_b);

                    if len_a mod len_b = 0 then
                        q := q + len_b*inst_b;
                    fi;
                od;

                p := p*q^inst_a;
            od;

            res := res +
            lcoeff(term_a)*lcoeff(term_b)*p;
        od;
    od;

    res;
end;

Answer 2

Your $b_n$ should be the same sequence as OEIS A091059, which counts the number of $n \times n$ matrices over symbol set $\{1,2\}$ equivalent under any permutation of rows, columns, or the symbol set. The correspondence is: \begin{align} \text{left node} &\leftrightarrow \text{row}\\ \text{right node} &\leftrightarrow \text{column}\\ \text{edge from left to right} &\leftrightarrow 1\\ \text{edge from right to left} &\leftrightarrow 2 \end{align}

But the OEIS entry has $b_2=5\not=4$ and $b_3=18\not=13$. Maybe the discrepancy is because of the horizontal reflection, which I think should not be included.