If you’re using a recent version of SciPY you might be using this idea already and not even know.

Bringing order to partition counts

Subsets have many varieties and they have some aspects that are distinguishable while other aspects may not be distinguishable.

Let’s start simple, we have 2 distinct items and we want to partition them into non-empty subsets. Allowing empty subsets can mess with the recursion options in case you’re wondering why they must be non-empty. Now for the two items, we can either put them into a separate set each or we can put them into the same set. The first is a partition into 2 non-empty subsets and the second is a partition into 1 non-empty subset. There is 1 way to do each (for now we don’t care about the ordering for the 2 subset case).

Now to go to 3 distinct items we recognize that we can either do everything from scratch or we can build up using what we have already in the 2 item case. So let’s use the 2 items case as a building block.

In so doing we consider how many ways to take either of the two items and add the third?

There a decision to be made when adding the new item; first do we make a new set which contains only this new item or do we add the new item to one of the existing non-empty subsets? For the first we have only 1 way to do this and for the second option we have \(k\) existing subsets (say) so we can pick any one of these \(k\).

Now if we arrange these numbers as a list for a given number of distinct items starting with 2 and moving to 3 we see that we have

n=1 -> [0,1,0,0]
n=2 -> [0,1,1,0]
n=3 -> [0,1,3,1]

many ways where the ith entry in the nth list tells us how many ways to make i non-empty subsets of n distinct items.

Now if you’re wondering where the 3 came from it is from choosing to add the 3rd item to either of the two singleton sets (2 ways) plus creating a new set with only the new item and adding that to the singleton set that contains the two previous items.

In general for \(n\) distinct items into \(i\) non-empty subsets we count these using this recursion, these subset numbers are listed in GKP’s concrete math book as Stirling numbers of the second kind where this recursion is given.

Two interesting things, the first is that there is no commonly agreed upon notation for these numbers. I follow GKP and use the bracket/brace notation,

\({n \brace k} = k{n-1 \brace k} + {n-1 \brace k-1}\),

for \(0<k \leq n\). The above equation also codifies the recursion relation for these numbers. Now you can use this to compute any value you like in a computer via what is called a dynamic programm or DP.

Now this is all well known, but one thing that is less well known is that you can do the DP in-place instead of using 2 adjacent rows of integers in separate memory locations, one storing the values for \(n-1\) and the other for \(n\) that is being populated in the current iteration.

Code implementations

To do this in place you need to traverse the row of integers right to left, also called traversing fro the back in some circles.

The steps are:

1. append a new entry to the list of integers at the end
2. multiply that value by the index in place
3. add the previous index value to the current indexed value

stop after index 1, assuming your list is 0-indexed.

This might not seem like much, it is only 1 list vs 2 after all but each item in the list becomes a large integer fairly quickly so this adds up for large n.

Here is a brief python code snippet to illustrate

def stirling2(n: int, k: int) -> int:
    # stirling numbers second kind ...
    if n==1 or k == n:
        return 1
    s2 = [0,1]
    for i in range(1, n-1):
        s2.append(1)
        for j in range(len(s2)-1, 1, -1):
            s2[j] *= j
            s2[j] += s2[j-1]
    s2.append(1)
    return s2[k]

Or javascript if you prefer…

/**
 * @param {number} n
 * @param {number} k
 * @return {number}
 */
var stirling2= function(n, k) {
    let row = [BigInt(0),BigInt(1)];
    // Note: you cal also do the tiling via a parallelogram
    // to reduce computation...
    for(let r=1; r<n-1; r++){
        if(row.length <= k){
            row.push(BigInt(1));
        }
        for(let c=Math.min(row.length-1, k); c>1; c--){
            row[c] *= BigInt(c);
            row[c] += row[c-1];
        }
    }
    row.push(BigInt(1));
    return Number(row[k]);
};

The above logic is roughly how I implemented the stirling2 function is SciPy for the exact=True branch.

Traverse in Reverse and Ordering partitions

The \(k\) value in the recursion formula above tells us that we need to multiply in place on index \(k\) before adding the value from \(k-1\). This means we must traverse from the back.

There is a related set of numbers which are a variant of these numbers, those don’t really have a name so I’ll call them ordered Bell numbers of order \(k\). You’ll see why in just a moment.

The ordered bell numbers are given by the formula:

\(b(n) = \sum_{k=1}^n b(n,k) = \sum_{i=1}^n k!{n \brace k}\).

Now I’m referring to the \(k\) component of this sum which is \(k!{n \brace k}\).

To calculate these the recursion is easily seen as

\(b(n,k) = kb(n-1, k) + kb(n-1,k-1)\).

In other words we’re multiplying both items by \(k\).

Now because the value \(k\) appears on both the recursion can go from either side e.g. you can traverse from the front, e.g. the usual way i=0, i=1, etc. and then divide by \(i!\) if you really want but you would in either case need to be concerned about numeric precision given how large the numbers become so quickly (in n).

Note that for the ordered Bell numbers of order \(k\), you can still do the computation in place similar to how you did the Stirling numbers of the second kind. The primary difference beside what the numbers actually count is that for ordered Bell numbers of order \(k\) you can iterate over a row in either direction and still do the computation in place.

Applications

For ordered bell numbers of order \(k\), we can interpret this as each non-empty subset being distinguished in some way, e.g. by different colors or distinct numbers on the subsets. This means 1 and then 2 as two non-empty subsets is distinct from 2 and then 1 as two non-empty subsets.

Another way to interpret is as the number of ties in a ranking of \(n\) elements where we have \(k\) distinct sets ranked. If there are no ties in a set then the set is a singleton.

The calculation of the number of possible ways to win a race with ties can come up when calculating summary metrics on search ranking results. This is especially true when you are working with human annotation teams and they are allowed to rank the values with ties in the ground truth search results.

Implementation Note

In the SciPy implementation the Stirling number computation is vectorized and given the nature of the partially overlapping subproblems I used a heap to schedule the computation for the individual cells of the matrices. There isn’t a better way to do this that I’m aware of and still be exact.

For the approximate branch-the default-the code is using the Lambert W function and uses the second order approximation from Temme see section 4 of the linked paper for details.

Cite this post

@misc{LucasRoberts2024subsetsandordering,
  author = {Lucas Roberts},
  title = {On subsets and ordering},
  year = 2024,
  howpublished = {\url{https://rlucas7.github.io/posts/2024/10/blog-post-1/}},
  note = {Accessed: 2024-10-15}
}