Counting non-isomorphic relations

On a set $X$ of $n$ elements, how many non-isomorphic relations are there? The number of relations on a set of $n$ elements is $|\mathcal{P}(X \times X)|=2^{n^2}$, but is there any way to give a general formula for how many of these are non-isomorphic? I am interested in this question for an application to modal logic. Thanks!

In particular, I consider two relations $R$ and $S$ to be isomorphic if there exists a bijection $\varphi: X \rightarrow X$ such that $\varphi(x) R \varphi(y) \Leftrightarrow x S y$

A relation on $X$ can be modeled as a graph on the vertex set $X$. Since it is a hard problem to determine the number of isomorphism types of graphs, I guess that your answer does not have a simple answer, and that there is no simple formula for the number of non-isomorphic relations.

ADDITION This is sequence 595 in the online encyclopedia of integer sequences. I'm a bit surprised that there is indeed a formula given: $$a(n) = \sum_{1s_1+2s_2+\ldots =n} \frac{\operatorname{fixA}[s_1, s_2, ...]}{1^{s_1}s_1!\cdot 2^{s_2}s_2!\cdot\ldots}$$ where $$\operatorname{fixA}[s_1, s_2, \ldots] = 2^{\sum_{i, j\geq 1} \gcd(i, j)s_i s_j}$$ I guess the formula arises as an application of the Cauchy-Frobenius lemma.

Here are some more details concerning this computation. We solve the problem using the Polya Enumeration Theorem. This computation is very similar to the enumeration of non-isomorphic graphs shown at this MSE link.

We need to compute the cycle index of the edge permutation group of the complete bipartite graph on $2n$ nodes. This is almost the same as the edge permutation group of a graph on $n$ nodes. The only difference is that the appropriate group is obtained from the action of the symmetric group $S_n$ on ordered pairs rather than unordered pairs as in the edge permutation group. This actually simplifies the analysis because a pair of vertices on an even cycle does not close the cycle after traversing $n/2$ elements but rather traverses all of the cycle.

The following Maple code computes this cycle index from the cycle index of the symmetric group, which, it must be pointed out, is faster than iterating over all partitions as in the other answer.

with(numtheory);
with(group):
with(combinat):

pet_cycleind_symm :=
proc(n)
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_varinto_cind :=
proc(poly, ind)
local subs1, subsl, polyvars, indvars, v, pot;

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

    subsl := [];

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

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

        subsl := [op(subsl), v=subs(subs1, poly)];
    od;

    subs(subsl, ind);
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_cycleind_rel :=
proc(n)
        option remember;
        local dsjc, flat, p, cyc1, cyc2, l1, l2, res;

        if n=0 then return 1; fi;
        if n=1 then return a[1] fi;

        res := 0;

        for dsjc in pet_cycleind_symm(n) do
            flat := pet_flatten_term(dsjc);

            p := 1;

            for cyc1 in flat[2] do
                l1 := op(1, cyc1);
                for cyc2 in flat[2] do
                    l2 := op(1, cyc2);

                    p := p *
                    a[lcm(l1,l2)]^(l1*l2/lcm(l1, l2));
                od;
            od;

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

        res;
end;

gf :=
proc(n)
option remember;
local gf;

    gf := pet_varinto_cind(1+z, pet_cycleind_rel(n));
    expand(gf);
end;

This is the cycle index for $n=3$: $$1/6\,{a_{{1}}}^{9}+1/2\,{a_{{2}}}^{4}a_{{1}}+1/3\,{a_{{3}}}^{3}$$ and here it is for $n=4$: $$1/24\,{a_{{1}}}^{16}+1/4\,{a_{{1}}}^{4}{a_{{2}}}^{6}+1/3\,{a_{{3}}}^{5}a_{{1}}+1/8\, {a_{{2}}}^{8}+1/4\,{a_{{4}}}^{4}. $$ The substituted cycle index for $n=4$ is $${z}^{16}+2\,{z}^{15}+9\,{z}^{14}+32\,{z}^{13}+95\,{z}^{12}+203\,{z}^{11}+373\,{z}^{ 10}+515\,{z}^{9}\\+584\,{z}^{8}+515\,{z}^{7}+373\,{z}^{6}+203\,{z}^{5}+95\,{z}^{4}+32 \,{z}^{3}+9\,{z}^{2}+2\,z+1 .$$ For $n=5$ we get $$ {\frac {1}{120}}\,{a_{{1}}}^{25}+1/12\,{a_{{1}}}^{9}{a_{{2}}}^{8}+1/6\,{a_{{3}}}^{7} {a_{{1}}}^{4}+1/8\,{a_{{2}}}^{12}a_{{1}} \\+1/4\,{a_{{4}}}^{6}a_{{1}}+1/6\,{a_{{2}}}^{2 }{a_{{6}}}^{2}{a_{{3}}}^{3}+1/5\,{a_{{5}}}^{5}.$$ The count of the number of non-isomorphic binary relations is $$2, 10, 104, 3044, 291968, 96928992, 112282908928, 458297100061728,\\ 6666621572153927936, 349390545493499839161856,\\66603421985078180758538636288,\\ 46557456482586989066031126651104256, \\ 120168591267113007604119117625289606148096,\\ 1152050155760474157553893461743236772303142428672,\\ 41233441401686067929518834693764864820973598636585435136,\ldots$$ This is indeed OEIS A000595 as noted by the other poster.

Addendum, September 4, 2018. Since we are only counting these relations as opposed to classifying them by the number of edges in the corresponding bipartite graph, it is not necessary to use PET. Burnside is perfectly sufficient. We can also streamline the computation of the cycle index. The modified program is shown below, data produced is the same as the first version.

with(numtheory);
with(group):
with(combinat):

pet_cycleind_symm :=
proc(n)
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_rel :=
proc(n)
option remember;
local dsjc, flat, p, cyc1, cyc2, l1, l2, res;

    if n=0 then return 1; fi;
    if n=1 then return a[1] fi;

    res := 0;

    for dsjc in pet_cycleind_symm(n) do

        p := 1;

        for cyc1 in indets(dsjc) do
            l1 := op(1, cyc1);
            for cyc2 in indets(dsjc) do
                l2 := op(1, cyc2);

                p := p *
                a[lcm(l1,l2)] ^
                (degree(dsjc, cyc1)*degree(dsjc, cyc2) *
                 l1*l2/lcm(l1, l2));
            od;
        od;

        res := res + lcoeff(dsjc)*p;
    od;

    res;
end;

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

    cind := pet_cycleind_rel(n);
    vars := indets(cind);

    subs({seq(v = 2, v in vars)}, cind);
end;

Note that the bipartite graph here is special, consisting of $2n$ nodes with $n$ on the left and $n$ on the right, both labeled $1$ to $n$. The cycle index represents the action on the edges by the simultaneous action of the symmetric group $S_n$ on the two sets of labels.