CSE 102 – Tutorial 9 – Combinatorial generation

After you are finished with the tutorial, please upload your submission to Moodle by the indicated deadline.

Generating integer partitions

A partition of a positive number is a way of writing it as a sum of positive numbers, where by convention we order the parts from largest to smallest. Formally, a partition of an integer \(n\) is a sequence of positive integers \(p_1,\dots,p_k\) such that \(p_1 \ge \dots \ge p_k\) and \(n = p_1 + \dots + p_k\). For example, the number 4 has five distinct partitions:

\[\begin{align*} 4 &= 4 \\ 4 &= 3 + 1 \\ 4 &= 2 + 2 \\ 4 &= 2 + 1 + 1 \\ 4 &= 1 + 1 + 1 + 1 \end{align*}\]

We will represent partitions in Python as reverse-sorted lists of positive integers. Thus the five partitions of 4 are represented as the following lists:

p1 = [4]
p2 = [3, 1]
p3 = [2, 2]
p4 = [2, 1, 1]
p5 = [1, 1, 1, 1]

Each of these lists is reverse-sorted (all (p[i] >= p[i+1] for i in range(len(p)-1))) and sums to 4 (sum(p) == 4). Incidentally, observe that the above listing of partitions is also reverse-sorted itself, in the sense that p1 > p2 > p3 > p4 > p5 relative to the lexicographic order on lists.

There is a simple iterative algorithm for enumerating partitions exhaustively in reverse lexicographic order, starting with the singleton partition and ending with the all-1s partition.

• Start by setting p = [n], and enter into a loop.
• At each iteration, yield the current partition p and then try to find the largest index i such that p[i] > 1.
• If there is no such index, that means that p is the all-1s partition [1, ..., 1], and we are done.
• Otherwise, p consists of some prefix p[:i], followed by the value p[i], followed by a suffix p[i+1:] of trailing 1s. Thus we replace p[i] by p[i]-1, and replace the trailing 1s by as many possible repetitions of p[i]-1 as fit, followed by a remainder term if necessary.

For example, if we apply this procedure to enumerate the partitions of 7, then after the initial partition [7], the next four partitions are [6, 1], [5, 2], [5, 1, 1], and [4, 3].

The goal of the next series of exercises is to write a function that generates partitions of an integer uniformly at random. We’ll begin by considering how to enumerate partitions recursively.

The number of partitions of \(n\) is denoted \(p(n)\), also known as the partition function. For example, \(p(4) = 5\) as we saw above. To compute \(p(n)\), it is convenient to introduce more refined partition numbers \(p(n,k)\) counting the number of partitions of \(n\) into \(k\) parts. For example, looking at the above listing, we see that \(p(4,2) = 2\) since there are two partitions of four into two parts ([3, 1] and [2, 2]), and \(p(4,3) = 1\) since there is only one partition of four into three parts ([2, 1, 1]). The numbers \(p(n,k)\) can be computed recursively via the equation

• \(p(n,k) = p(n-1,k-1) + p(n-k,k)\)

for all \(n>0\) and \(1 \le k \le n\), with the equations

• \(p(n,k) = 0\) if \(k > n\)
• \(p(n,0) = 0\) if \(n > 0\)
• \(p(0,0) = 1\)

covering the remaining cases. The total number of partitions of \(n\) can then be computed by the formula

\[p(n) = \sum_{k=0}^n p(n,k)\]

Your next task will be to write a generator function that goes beyond counting to explicitly enumerate all partitions of an integer recursively. For this, it is important to understand where the recurrence relation and boundary cases for the partition numbers \(p(n,k)\) come from. When \(n > 0\) and \(1 \le k \le n\), the recurrence relation

\[p(n,k) = p(n-1,k-1) + p(n-k,k)\]

can be established by the following argument: a partition of \(n\) into \(k\) parts either has smallest part 1, in which case the rest is a partition of \(n-1\) into \(k-1\) parts, or else its smallest part is at least 2, and we can subtract 1 from each part to obtain a partition of \(n-k\) into \(k\) parts. In other words, assuming \(n > 0\) and \(1\le k\le n\), we can obtain any partition \(n = a_1 + \dots + a_k\) of \(n\) into \(k\) parts by either:

• starting with a partition \(n-1 = a_1 + \dots + a_{k-1}\), and adding a part \(a_k = 1\) to get \(n = a_1 + \dots + a_{k-1} + 1\); or

• starting with a partition \(n-k = b_1 + \dots + b_k\), and incrementing each of the parts to get \(n = (1+b_1) + \dots + (1+b_k)\).

In the boundary cases, we have

\[p(n,k) = 0\text{ if }k > n \qquad \text{and} \qquad p(n,0) = 0\text{ if }n > 0\]

because there is no way to express \(n\) as a sum of more than \(n\) positive integers, or to express it as an empty sum if \(n>0\). Finally, we do have

\[p(0,0) = 1\]

because the empty partition, i.e., the empty sum, is a valid partition of 0 into 0 parts.

Now, as discussed in Lecture 9, it is possible to turn a recursive enumerator into a random generator by replacing exhaustion with random choice. For example, the generator function perms_rec(n) below implements a recursive strategy for enumerating all permutations of the numbers \(1\dots n\), and randperm(n) uses a corresponding strategy to compute a single random permutation.

def perms_rec(n):
    if n == 0:
        yield []
        return
    for pi in perms_rec(n-1):            # iterate over all permutation of 1..n-1
        for i in range(n):               # iterate over all indices 0 <= i < n
            yield pi[:i] + [n] + pi[i:]  # build a permutation of 1..n

def randperm(n):
    if n == 0:
        return []
    pi = randperm(n-1)                   # pick a random permutation of 1..n-1
    i = random.randrange(n)              # pick a random index 0 <= i < n
    return pi[:i] + [n] + pi[i:]         # build a permutation of 1..n

This example is a bit special in that we immediately obtained a uniform random generator without any extra work: one can verify that every permutation of \(1\dots n\) has equal probability of being returned by randperm(n).

In the case of permutations we could use uniform choice when choosing the values of the random variables pi and i, but more generally, we need to use biased sampling to turn a recursive enumerator into a uniform random generator, weighing the different possible values for each random variable by the probability that a random object has that value. (Question: so why was it okay to select pi and i using uniform sampling?)

For the concrete case of partitions it works as follows. In the recursive construction of a partition of \(n\) with \(k\) parts, to generate a random partition we need to decide whether to construct one from a partition of \(n-1\) into \(k-1\) parts, or from a partition of \(n-k\) into \(k\) parts. But we should not decide by flipping a fair coin! Instead we should flip a biased coin, choosing the first construction with probability \(\frac{p(n-1,k-1)}{p(n,k)}\). Similarly, to generate a random partition of \(n\) with an arbitrary number of parts, we can pick \(k\) at random between \(0\) and \(n\), and then generate a random partition of \(n\) into \(k\) parts. But to get a uniform generator, we must not choose \(k\) uniformly within this range. Rather, we should choose the value \(k\) with probability \(\frac{p(n,k)}{p(n)}\). (For instance, when \(n>0\), we should never pick \(k=0\) since \(p(n,0) = 0\), and we should only pick \(k=1\) with probability \(1/p(n)\). The probability that a random partition of \(n=100\) has exactly one part is about 1 in 190 million.)

In order to apply this strategy to generate large random integer partitions (say, with \(n=1000\)), we need to be able to compute partition numbers efficiently, and for that it’s important to use dynamic programming. Moreover, since we may want to sample random partitions repeatedly, it makes sense to precompute a table of partition numbers in advance. To achieve this, we’ll bundle the random generator into a class RandPartsGen equipped with a constructor RandPartsGen(maxn) that takes an upper bound on the size of the partitions and precomputes a table of partition numbers. The constructed RandPartsGen object is equipped with a method randpart(n, k=None) that returns a uniformly random partition of the given integer n (optionally, with a given number of parts k), assuming that n <= maxn.

For the next series of optional problems, we consider generation of partitions with some additional restrictions. A partition is said to be “odd” if it contains only odd integer parts. For example, 7 has five odd partitions: \[\begin{align*} 7 &= 7 \\ 7 &= 5 + 1 + 1 \\ 7 &= 3 + 3 + 1 \\ 7 &= 3 + 1 + 1 + 1 + 1 \\ 7 &= 1 + 1 + 1 + 1 + 1 + 1 + 1 \end{align*}\]

Rejection sampling of odd partitions is inefficient, because as \(n\) grows larger a random partition of \(n\) is very unlikely to contain only odd parts. A much better approach is to generate them recursively, following a similar strategy to the one we used for unrestricted partitions, but now defining odd partitions in mutual recursion with even partitions (i.e., partitions containing only even integer parts). Indeed, if we let \(p_o(n,k)\) and \(p_e(n,k)\) stand for the numbers of odd and even partitions of \(n\) into \(k\) parts, respectively, then these partition counts satisfy the recurrence relations

\[\begin{align*} p_o(n,k) &= p_o(n-1,k-1) + p_e(n-k,k) \\ p_e(n,k) &= p_o(n-k,k) \end{align*}\]

for all \(n>0\) and \(1 \le k \le n\). We can argue for the recurrence relations as follows. A partition of \(n\) into \(k\) odd parts either has smallest part 1, in which case the rest is a partition of \(n-1\) into \(k-1\) odd parts, or else its smallest part is at least 3, and we can subtract 1 from each part to obtain a partition of \(n-k\) into \(k\) even parts. Hence \(p_o(n,k) = p_o(n-1,k-1) + p_e(n-k,k)\). On the other hand, a partition of \(n\) into \(k\) even parts, for \(1 \le k \le n\), necessarily has smallest part at least 2, and so we can subtract 1 from each part to obtain a partition of \(n-k\) into \(k\) odd parts. Hence \(p_e(n,k) = p_o(n-k,k)\). Finally, the odd and even partition numbers satisfy the same boundary cases as for the ordinary partition numbers \(p(n,k)\), that is \(p_o(n,k) = p_e(n,k) = 0\) if \(k > n\) or if \(n > 0\) and \(k = 0\) (since there are no partitions in these cases, and hence no odd/even partitions), and \(p_o(0,0) = p_e(0,0) = 1\) (since, vacuously, all of the parts in the empty partition are both even and odd!).

Another natural restriction that we may place on integer partitions is that all of the parts be distinct. For example, there are five ways of writing 7 as a sum of distinct positive numbers: \[\begin{align*} 7 &= 7 \\ 7 &= 6 + 1 \\ 7 &= 5 + 2 \\ 7 &= 4 + 3 \\ 7 &= 4 + 2 + 1 \end{align*}\] A partition \(n = p_1 + \dots + p_k\) into distinct parts \(p_1 > \dots > p_k\) is also called a strict partition of \(n\). Once again, we could try to use rejection sampling to generate strict partitions uniformly at random, but this is inefficient because a large random partition is very likely to have repeated parts.

Fortunately, it is possible to turn a uniform generator of odd partitions into a uniform generator of strict partitions! Back in the mid-1700s, Euler proved that the number of partitions of \(n\) with odd parts is always equal to the number of partitions of \(n\) with distinct parts. Although Euler’s proof used purely algebraic manipulation, later J. J. Sylvester gave a proof of a (strengthening of) Euler’s theorem by describing an explicit bijection between partitions with odd parts and partitions with distinct parts. In general, if two sets of objects \(X\) and \(Y\) are in bijection \(X \rightleftarrows Y\), then we can turn an efficient uniform generator of \(X\)s into an efficient uniform generator of \(Y\)s and vice versa — at least assuming that the bijection is efficiently computable.

Sylvester’s bijection is best described using a graphical representation of partitions known as Young tableau or Ferrers diagrams. The idea is to represent each part by a row of boxes (or dots), stacked one on top of the other in descending order. For example, here is a representation of the (odd) partition \(21 = 9 + 5 + 3 + 3 + 1\) as a tableau:

To describe the bijection, let’s begin by placing a black square in every other box of the first row, as well as all along the first column:

The first part of the new partition will be 9, corresponding to the number of black squares. Next we remove the first column, so that what remains is an even partition (\(8 + 4 + 2 + 2 = 16\)):

This time we place a red square in every other box of the first row and all along the first column:

The second part of the new partition will be 7, corresponding to the number of red squares. We remove the first column and the first row and repeat the process until there are no more boxes, as illustrated below:

In the end we’ve obtained a partition of 21 into distinct parts \(21 = 9 + 7 + 4 + 1\). Observe though that the number of parts is not preserved!

Now that you’ve implemented the bijection, it is easy to get a uniform generator of strict partitions, simply by composing with your uniform generator of odd partitions:

def randspart(rpg, n):
    return odd_to_strict(rpg.randopart(n))

Some typical outputs:

In [29]: randspart(rpg, 100)
Out[29]: [37, 22, 12, 9, 8, 6, 5, 1]
In [30]: randspart(rpg, 100)
Out[30]: [23, 22, 20, 19, 13, 3]
In [31]: randspart(rpg, 1000)
Out[31]: [106, 80, 74, 70, 64, 63, 62, 58, 50, 47, 42, 40, 39, 37, 35, 30, 25, 22, 21, 20, 15]

Again, you should write some tests to confirm that the output of randspart is truly uniform over strict partitions! It will be, assuming that you’ve implemented odd_to_strict and randopart correctly.

Finally, as we remarked above, Sylvester’s bijection does not preserve the number of parts, and so composing the uniform generator randopart(n, k) with odd_to_strict will not produce a uniform generator of strict partitions of n with k parts. However, it is still possible to generate such partitions uniformly at random, by another approach.