Set Partitions and Bell Numbers

An equivalence relation on a set determines a partition of it, and vice versa — this is intuitively not hard to understand. But if I ask how many different partitions there are of a finite set with $n$ elements, that used to feel like a very simple question I never bothered to think carefully about. Then the day before yesterday, my abstract algebra teacher mentioned that this is actually a fairly hard problem, so I looked into it, and found there's quite a lot going on. Here I'll share a brief account of my investigation process, so readers can see what a "from nothing to something" journey looks like.

Below, let a finite set with $n$ elements be $\{1,2,\dots,n\}$, and denote the number of its partitions by $B(n)$.

Early stage: brute-force computation

The case $n=3$ is easy to enumerate:

$$\begin{aligned}&\{\{1,2,3\}\},\{\{1,2\},\{3\}\},\{\{1,3\},\{2\}\},\\ &\{\{2,3\},\{1\}\},\{\{1\},\{2\},\{3\}\}\end{aligned}$$more

So $B(3)=5$. $B(4)$ is already fairly large and not so easy to list by hand. Working from integer partitions, I derived a fairly involved formula for computing set partition numbers. Take $S(4)$ as an example — it has the following decompositions:

$$\begin{aligned}&4=1+1+1+1\\ &4=1+1+2\\ &4=1+3\\ &4=4\\ &4=2+2\end{aligned}$$

Each decomposition corresponds to a way of partitioning $\{1,2,3,4\}$: it tells us how many blocks there are and how many elements each block has. For instance, 4=1+1+1+1 means splitting $\{1,2,3,4\}$ into 4 blocks, one element per block, and the number of possible partitions of this type is

$$\frac{1}{4!}\binom{4}{1}\binom{3}{1}\binom{2}{1}=1$$

We divide by $4!$ because all the "1"s here play equivalent roles — different orderings of them count as the same partition. Similarly we have

$$\begin{aligned}&4=1+1+2:\quad \frac{1}{2!}\binom{4}{1}\binom{3}{1}=6\\ &4=1+3:\quad \binom{4}{1}=4\\ &4=4:\quad \binom{4}{4}=1\\ &4=2+2:\quad \frac{1}{2!}\binom{4}{2}=3\end{aligned}$$

So $B(4)=1+6+4+1+3=15$. Using the same approach we can get (I computed this by hand) $B(5)=52,B(6)=203$; fully enumerating all decompositions for 5 and 6 is still manageable by hand.

Middle stage: casting a wide net

At this point we already have part of the sequence: 5, 15, 52, 203. It doesn't seem to follow any obvious pattern, and intuitively this sequence seems to grow even faster than exponentially. Feeding this sequence into OEIS?, we find it (the "number of ways to partition a set of n labeled elements"):

1, 1, 2, 5, 15, 52, 203, 877, 4140, 21147, 115975, 678570, 4213597, 27644437, 190899322, ...

The FORMULA section there lists many curious formulas. But since only the formulas are given, with no derivation, I couldn't fully make sense of them. So I had to look elsewhere. Searching Google for keywords like "number of partitions of a set" turned up a large batch of similar results. These gave me some help — for example, introducing $S(n,k)$ to compute $B(n)$, where $S(n,k)$ denotes the number of ways to partition a set of $n$ elements into $k$ blocks. By this definition, clearly

$$B(n)=\sum_{k=1}^n S(n,k)$$

Computing $S(n,k)$ is then done recursively. Consider adding one element to a set of $n-1$ elements — there are two possibilities:

1. The new element forms its own block by itself, and the remaining $n-1$ elements are partitioned into $k-1$ blocks; the number of cases here is $S(n-1,k-1)$;
2. The $n-1$ (original) elements are partitioned into $k$ blocks, and the new element is placed into any one of them; the number of cases here is $k S(n-1,k)$.

This gives us the recurrence

$$S(n,k)=S(n-1,k-1)+k S(n-1,k)$$

This is probably the most convenient formula for programming, written in Python as

def sp(n,m):
    if m>n:
        return 0
    elif m==1 or m==n:
        return 1
    else:
        return sp(n-1,m-1)+m*sp(n-1,m)

def setpart(n):
    return sum([sp(n,m) for m in range(1,n+1)])

print(setpart(7))

But code like this is actually extremely slow to compute. Even on PyPy, computing $B(30)$ is already pushing it. I spent a long time yesterday trying to derive the generating function for $B(n)$ from the recurrence, but I didn't have enough conditions to get a complete result.

Middle-late stage: analytical computation

In my middle-stage search, most of the articles I found were similar, only introducing the concept and the recurrence without going into more detail. But one website mentioned at the end of its article that these numbers are called "Bell numbers," and following that lead, I found the Wikipedia page:

http://zh.wikipedia.org/zh-cn/贝尔数

Wikipedia had a lot of valuable material, and I also learned that the standard notation is $B_n$, while the $S(n,k)$ introduced above is called the "Stirling number of the second kind." Wikipedia gives the recurrence

$$B_{n+1}=\sum_{k=0}^{n}{{n \choose k}B_k}$$

The combinatorial interpretation of this formula can be found directly on Wikipedia, so I won't repeat it here. I spent a long time last night unable to extract any further information from this formula (turns out I'd misread the formula — oops), such as computing the generating function. So I searched the English Wikipedia for "Bell Number" and got more detailed results than the Chinese version, though still without a derivation. Today, following hints from the English Wikipedia, I worked through it again and obtained the generating function. Let

$$f(x)=\sum_{n=0}^{\infty}\frac{B_n}{n!}x^n$$

Differentiating it, substituting into the recurrence, and swapping the order of summation:

$$\begin{aligned}f'(x)&=\sum_{n=1}^{\infty}\frac{B_n}{(n-1)!}x^{n-1}=\sum_{n=0}^{\infty}\frac{B_{n+1}}{n!}x^n\\ &=\sum_{n=0}^{\infty}\sum_{k=0}^{n}\frac{1}{n!}{n \choose k}B_k x^n\\ &=\sum_{n=0}^{\infty}\sum_{k=0}^{n}\frac{1}{(n-k)!k!}B_k x^k x^{n-k}\\ &=\sum_{k=0}^{\infty}\sum_{n=k}^{\infty}\frac{1}{(n-k)!k!}B_k x^k x^{n-k}\\ &=\sum_{k=0}^{\infty}\frac{B_k x^k }{k!}\sum_{n=k}^{\infty}\frac{1}{(n-k)!}x^{n-k}\\ &=\sum_{k=0}^{\infty}\frac{B_k x^k }{k!}e^x\\ &=f(x)e^x\\ \end{aligned}$$

This gives us $f(x)=e^{e^x-c}$, and from $f(0)=1$ we get $c=1$, hence

$$\sum_{n=0}^{\infty}\frac{B_n}{n!}x^n=f(x)=e^{e^x-1}$$

This is a neat, elegant generating function — an "exponential of an exponential" — which shows just how fast $B_n$ grows. From the generating function, we can derive an infinite series formula for $B_n$. Since

$$\begin{aligned}e^{e^x}&=\sum_{m=0}^{\infty}\frac{1}{m!}e^{mx}\\ &=\sum_{m=0}^{\infty}\frac{1}{m!}\sum_{n=0}^{\infty}\frac{1}{n!}m^n x^n\\ &=\sum_{n=0}^{\infty}\frac{1}{n!}\sum_{m=0}^{\infty}\frac{1}{m!}m^n x^n \end{aligned}$$

comparing both sides of $f(x)$, we obtain

$$B_n=\frac{1}{e}\sum_{m=0}^{\infty}\frac{m^n}{m!}$$

This is the Dobinski formula.

Later stage?

Learning never ends, so there is no "later stage." Let's keep encouraging each other on this journey~

English translation of a post from 科学空间 | Scientific Spaces by 苏剑林. Original: https://kexue.fm/archives/2985
Translated automatically with claude-sonnet-5; all equations are reproduced verbatim from the source. Copyright remains with the original author.