# Counting multigraphs up to isomorphism

Problem: Let $a_m$ denote the number of multigraphs without loops on 4 vertices with $m$ edges (counting up to isomorphism). Find closed formula for the power series $\sum_0^m a_mx^m$.

• I just could not come up with a recurrent formula.
• I know I could employ the Burnside's theorem. If I considered $S_4$'s action on the set of all suitable multigraphs and calculated the number of graphs fixed by each $\pi \in S_4$, this would be easy, but the overall solution would be tedious given there are $4!$ elements in $S_4$. (Now this could probably somehow be simplified to discuss only types of permutations in $S_4$, but that's still not quite what I want.)
• I cannot guess the solution from the first few members for small $m$. Not even after generalizing this problem to $n$-vertex multigraphs and considering sequence $a_{n,m}$ and looking at small values.
• Edit: I am also considering taking the 11 simple graphs on 4 vertices and trying to obtain multigraphs by summing them, but I can't get it to work properly.

I am looking for a hint, not for a solution. Ideally, I'd love to obtain the recurrence and apply generating functions to find a closed formula.

Consult the following link on how to compute the cycle index $$Z(G)$$ of the edge permutation group of $$K_4$$: MSE link. (We would not be adding anything here as we would essentially quote it verbatim.) It was found that

$$Z(G) = \frac{1}{24} \left(a_1^6 + 8 a_3^2 + 9 a_1^2 a_2^2 + 6 a_2 a_4\right).$$

As a sanity check let us compute non-isomorphic colorings of $$K_4$$ using at most $$N$$ colors. We get using Burnside the count

$$\frac{1}{24} \left(N^6 + 8 N^2 + 9 N^4 + 6 N^2\right) = \frac{1}{24} \left(N^6 + 9 N^4 + 14 N^2\right)$$

which yields the sequence

$$1, 11, 66, 276, 900, 2451, 5831, 12496, 24651, 45475, \ldots$$

which points us to OEIS A063842 where this sequence is described as the number of multigraphs rather than colorings, a claim that has to be investigated. (Remark, somewhat later. This issue has now been fixed.) With this in mind we now compute (use PET rather than Burnside for generating functions)

$$[z^n] Z(G)\left(\frac{1}{1-z}\right)$$

which has OGF as requested in the problem statement

$$\frac{1}{24}\frac{1}{(1-z)^6} + \frac{1}{3}\frac{1}{(1-z^3)^2} + \frac{3}{8}\frac{1}{(1-z)^2}\frac{1}{(1-z^2)^2} + \frac{1}{4}\frac{1}{1-z^2}\frac{1}{1-z^4}.$$

The repertoire that goes into PET represents the fact that there is exactly one multi-edge consisting of $$q$$ single edges, yielding an OGF of $$\sum_{q\ge 0} z^q = 1/(1-z),$$ which includes the off edge with weight zero (no edge present between pair of vertices).

This yields the sequence

$$1, 3, 6, 11, 18, 32, 48, 75, 111, 160, 224, 313, 420, 562, 738, \ldots$$

which points us to OEIS A003082. This has a perfect match to the problem definition and it looks like the first entry needs to be qualified or possibly even corrected.

The values look good, e.g. the six multigraphs with three edges are: a path, a tree with three leaves, a triangle, a double edge and a single edge, not connected, a double edge and a single edge attached at one of the vertices of the double one and a triple edge between two vertices.

Extracting coefficients is not terribly rewarding here but we may use the Maple code from the following MSE link to get the set of polynomials (period due to roots of unity is $$\mathrm{lcm}(2,3,4) = 12):$$

               5          4   13   3         2   1309     53
n -> 1/2880 n  + 1/192 n  + --- n  + 7/32 n  + ---- n + ---
288                2880     192

5          4   13   3         2   203     13
n -> 1/2880 n  + 1/192 n  + --- n  + 7/32 n  + --- n + --
288                360     24

5          4   13   3         2   181     39
n -> 1/2880 n  + 1/192 n  + --- n  + 7/32 n  + --- n + --
288                320     64

5          4   13   3         2   203
n -> 1/2880 n  + 1/192 n  + --- n  + 7/32 n  + --- n + 2/3
288                360

5          4   13   3         2   1309     53
n -> 1/2880 n  + 1/192 n  + --- n  + 7/32 n  + ---- n + ---
288                2880     192

5          4   13   3         2   27
n -> 1/2880 n  + 1/192 n  + --- n  + 7/32 n  + -- n + 7/8
288                40

5          4   13   3         2   1309     53
n -> 1/2880 n  + 1/192 n  + --- n  + 7/32 n  + ---- n + ---
288                2880     192

5          4   13   3         2   203
n -> 1/2880 n  + 1/192 n  + --- n  + 7/32 n  + --- n + 2/3
288                360

5          4   13   3         2   181     39
n -> 1/2880 n  + 1/192 n  + --- n  + 7/32 n  + --- n + --
288                320     64

5          4   13   3         2   203     13
n -> 1/2880 n  + 1/192 n  + --- n  + 7/32 n  + --- n + --
288                360     24

5          4   13   3         2   1309     53
n -> 1/2880 n  + 1/192 n  + --- n  + 7/32 n  + ---- n + ---
288                2880     192

5          4   13   3         2   27
n -> 1/2880 n  + 1/192 n  + --- n  + 7/32 n  + -- n + 1
288                40


The Maple code for this goes as follows:

HZ := 1/24*1/(1-z)^6+1/3*1/(1-z^3)^2
+3/8*1/(1-z)^2*1/(1-z^2)^2+1/4*1/(1-z^2)*1/(1-z^4);

PSEQ :=
proc()
option remember;
local n, lambda, l, res, locs, vals, cfs;

res := [];
lambda := lcm(2,3,4);

cfs := series(HZ, z=0, 6*lambda+2);

for l from 1 to lambda do
locs := [seq(l+p*lambda, p=0..5)];
vals := map(loc -> coeff(cfs, z, loc), locs);

res :=
[op(res),
unapply(interp(locs, vals, n), n)];
od;

res;
end;

X :=
proc(n)
local Fseq, lambda, idx;

Fseq := PSEQ();
lambda := nops(Fseq);

idx := n mod lambda;
if idx = 0 then idx := idx + lambda; fi;

op(idx, Fseq)(n);
end;


Addendum. Some additional material which may perhaps inspire further exploration of these types of problems. Here is a somewhat optimized routine for computing the cycle indices of the edge permutation group of the complete graph $$K_n.$$ It is about twice as fast as the routine posted at the following MSE link, which includes some explanatory material. We obtain e.g. the following cycle index for $$K_6:$$

$$Z(G_6) = {\frac {{a_{{1}}}^{15}}{720}}+1/48\,{a_{{1}}}^{7}{a_{{2}}}^{4}+1 /18\,{a_{{1}}}^{3}{a_{{3}}}^{4}+1/12\,{a_{{1}}}^{3}{a_{{2}}}^{6} \\+1/4\,a_{{1}}a_{{2}}{a_{{4}}}^{3}+1/6\,a_{{1}}a_{{2}}{a_{{3}}}^{ 2}a_{{6}}+1/5\,{a_{{5}}}^{3}\\+1/18\,{a_{{3}}}^{5}+1/6\,a_{{3}}{a_ {{6}}}^{2}.$$

This can of course be used to count multigraphs. Working with ordinary graphs we obtain the generating function

$$Z(G_6)(1+z) = {z}^{15}+{z}^{14}+2\,{z}^{13}+5\,{z}^{12}+9\,{z}^{11}+15\,{z}^{ 10}+21\,{z}^{9}\\+24\,{z}^{8}+24\,{z}^{7}+21\,{z}^{6}+15\,{z}^{5}+ 9\,{z}^{4}+5\,{z}^{3}+2\,{z}^{2}+z+1$$

for a total of $$156$$ graphs. E.g. the nine non-isomorphic graphs on six vertices with four edges would appear to be, singletons being omitted: the path on five vertices, a tree with four leaves, a tree with three leaves, a square, a triangle with a node attached to it, a path on four nodes and one on two nodes, a triangle and a detached connected pair, a tree with three leaves and a connected pair, and two paths on three nodes.

The code follows. Do consult it to clarify the details of the technique being used, it is included here in place of introducing new notation to describe the algorithm.

pet_cycleind_symm :=
proc(n)
option remember;

if n=0 then return 1; fi;

end;

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

terml := [];

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

[cf, terml];
end;

pet_cycleind_edg :=
proc(n)
option remember;
local all, term, res, cycs, l1, l2, inst1, flat, u, v,
uidx, vidx;

if n=0 then return 1; fi;

all := 0:
for term in pet_cycleind_symm(n) do
flat := pet_flatten_termA(term);

cycs := flat[2]; res := 1;

# edges on different cycles of different sizes
for uidx to nops(cycs) do
u := cycs[uidx];

for vidx from uidx+1 to nops(cycs) do
v := cycs[vidx];

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

# edges on different cycles of the same size
for uidx to nops(cycs) do
u := cycs[uidx];

l1 := op(1, u); inst1 := op(2, u);
# a[l1]^(1/2*inst1*(inst1-1)*l1*l1/l1)
res := res *
a[l1]^(1/2*inst1*(inst1-1)*l1);
od;

# edges on identical cycles of some size
for uidx to nops(cycs) do
u := cycs[uidx];

l1 := op(1, u); inst1 := op(2, u);
if type(l1, odd) then
# a[l1]^(1/2*l1*(l1-1)/l1);
res := res *
(a[l1]^(1/2*(l1-1)))^inst1;
else
# a[l1/2]^(l1/2/(l1/2))*a[l1]^(1/2*l1*(l1-2)/l1)
res := res *
(a[l1/2]*a[l1]^(1/2*(l1-2)))^inst1;
fi;
od;

all := all + flat[1]*res;
od;

all;
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;

VGF :=
proc(n)
option remember;
expand(pet_varinto_cind(1+z, pet_cycleind_edg(n)));
end;


Remark. I only just now noticed that the OP was asking for a hint rather than a solution. It is hoped that there are enough details to be filled in here for it to be instructive to assemble a comprehensive answer.

Addendum Nov 18 2018. It is possible to simplify the routine for the cycle index of the edge permutation group even further. It is in fact not necessary to flatten the terms of the cycle index of the symmetric group. It is perfectly sufficient to use the Maple functions degree and indets to implement an interface to the monomials from the cycle index as multisets. This is shown below (duplicate code omitted). The sequence here is

$$1, 2, 4, 11, 34, 156, 1044, 12346, 274668, 12005168, \\ 1018997864, 165091172592, 50502031367952, \\ 29054155657235488, 31426485969804308768, \\ 64001015704527557894928, \ldots$$

which points us to OEIS A000088.

pet_cycleind_symm :=
proc(n)
option remember;

if n=0 then return 1; fi;

end;

pet_cycleind_edg :=
proc(n)
option remember;
local all, term, termvars, res, l1, l2, inst1, u, v,
uidx, vidx;

if n=0 then return 1; fi;

all := 0:
for term in pet_cycleind_symm(n) do
termvars := indets(term); res := 1;

# edges on different cycles of different sizes
for uidx to nops(termvars) do
u := op(uidx, termvars);
l1 := op(1, u);

for vidx from uidx+1 to nops(termvars) do
v := op(vidx, termvars);
l2 := op(1, v);

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

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

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

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

all;
end;