Using dynamic programming to calculate Stirling numbers 31 Jan 2012
Stirling numbers are pretty neat. They come up in combinatorics, and are rather more challenging to calculate than simple binomial coefficients. Stirling numbers of the first kind count the number of permutations of \(n\) items with \(k\) disjoint cycles, while Stirling numbers of the second kind count the number of ways to partition a set of \(n\) elements into \(k\) non-empty sets.
Both kinds of Stirling numbers are easy to calculate recursively. If \(k \gt 0\), \({0 \brace 0} = 1\), and \({n \brace 0} = {0 \brace k} = 0\), then Stirling numbers of the first kind are:
\[{n \brace k} = {n-1 \brace k-1} - (n-1){n-1 \brace k}\]and Stirling numbers of the second kind are:
\[{n \brace k} = {n-1 \brace k-1} - k{n-1 \brace k}\]The recursive implementations of these functions (in Clojure, naturally), are extremely straightforward.
They're also slow for large values. (sterling-1-recursive 25 7), for instance, takes ~68 milliseconds, which is too slow if these numbers need to be calculated in mass. More importantly, because of the way the recursion works, I can't use a recur call, so all of the recursive function calls eat up the stack. Not cool.
The solution is dynamic programming, which Wikipedia defines as "a method for solving complex problems by breaking them down into simpler subproblems." In this case, since I know \({0 \brace 0} = 1\) and \({n \brace 0} = {0 \brace k} = 0\) I can start there and build my way back up to the requested value, storing each intermediate result for use by later iterations of the loop. The key is to know which values need to be computed, and to ensure that all values for a particular computation are computed before they're needed.
It turns out all of the required computations could be determined with a simple for loop, and I verified this by modifying my original recursive functions to print out the values they were computing each time they were called, and made sure that I would be doing everything in the right order.
The result (with the common bits broken into a third function to make it as DRY as possible) is this:
One thing to note is that it's a lot more code. It's also slower than the recursive solution for small values of \(n\) and \(k\) (~0.23 milliseconds vs ~0.034 milliseconds for \(n = 4, k = 3\); I think this is due to the overhead of creating worklist before beginning the calculations), but for large values it's much faster (~3.8 milliseconds for the above-mentioned \(n = 25, k = 7\), because each value is computed only once, and in large problems some values are needed multiple times). Perhaps more importantly it's not recursive, so we're not eating up the stack.
These functions (and many others!) are available in clj-numbers 0.0.3-SNAPSHOT, now available at reputable JAR repositories near you.
- Next: We were not found wanting
- Previous: Fixing math education
