Multivariable recurrence relations

I’m going to ramble on about what I found out about multivariable recurrence relations today, and digress here and there. One thing that I found to be particularly interesting was a way to show that Pascal’s triangle has binomial coefficients as entries.

1. An example of methods

The following is from a 6.006 problem set that Colleen showed me (which I paraphrase, keeping some of the silly flavor of the original).

Right when an upcoming concert is announced, some extraordinary fans of the band dutifully line up in front of the ticket office. When it’s time, the ticket master opens his door and says, “Hey! So you think you’re all fans? I’ll show you what a real fan is!” and then points a large fan at the them, turning it on. Each person is either sitting in a flimsy chair or standing; the fan blows away those that are sitting in a chair and those that are standing in front of a sitting person, but the line is so long that the last person in line is standing in front of a brick wall and so is not blown away. How many ways can n fans be sitting or standing so that exactly k fans remain? (and write a program to compute this.)

The most remarkable term is T_u(n − 1, k − 1): the k − 1 is there because we assume the last person is standing, so they know they will remain; hence, one fewer required fan in the remainder of the line.

These are sufficient for writing a O(nk)-time algorithm to compute T_u(n, k).^[1] But, I’m more interested right now in understanding these mutually-recursive multivariable functions.

From the second, we see that T_d(n − 1, k) = T_d(n, k) − T_u(n − 1, k), and, combining this with the first, we get T_d in terms of T_u:

which we can put back into the first to get the following recurrence relation:

From now on, we’ll write T = T_u.

Just to get a sense of the function, here is a table of T(n, k) with the the rows and columns being n and k.^[2]

My motivating question is whether there’s a closed-form solution to compute T. While I didn’t actually find one, I learned about multivariable generating functions in the process, and this is mainly what I’m wanting to write about. (However, see Update 6/28/17.) The multivariable case of a generating function is similar to the single variable case, except that there is a c_ijxⁱy^j for every term c_ij (in what might be called a bi-sequence). Hence, we let

be our generating function.

I managed to compute g at first as a generating function of generating functions, fixing each variable one at a time, using Mathematica to do some of the hairy symbolic manipulation. However, it turns out there is a much simpler way. To see how we’d do this, consider that

by defining T of negative numbers to be zero. That is, multiplication of g by x^ay^b gives a generating function g ′ of

Hence, the right hand side of the recurrence relation has the generating function

and when we consider the boundary conditions (namely that T(0, 0) = 1), we have

and so

Aside: if y = 0 (which is equivalent to saying k = 0 since y = 0 kills all terms involving non-zero k), this reduces to the recurrence relation of the Fibonacci sequence. In fact, a combinatorial proof is that k = 0 corresponds to having each person stand only when there is a sitting person behind them; that is, tile an n × 1 board with 2 × 1 and 1 × 1 pieces, representing a stander-sitter and a sitter, respectively. I’ve seen that this problem is counted by the Fibonacci sequence, so it’s enough of a proof for me at least!

Unfortunately, this polynomial is irreducible, so there’s not really anything I know we can do to find a closed form. Perhaps interestingly, using the geometric series 1/(1 − t) = 1 + t² + t³ + …, we can rewrite g in the form

so to compute T(n, k), we need only keep track of which of these terms contribute to xⁿy^k, which can only happen for exponents (roughly) from (n + k)/3 to n + k.

2. The binomial coefficient

Another function which is conducive to study using multivariable recurrences is the binomial coefficient. Let’s say we start with Pascal’s triangle:

and the closed form for the generating function is

Similarly rewriting using the geometric series, we have

Importantly, each monomial in f comes from exactly one of the (x + xy)ⁿ terms, which are each composed of monomials xⁿy^{n − k} (and xⁿ can only come from this particular term). What are the coefficients of these monomials? They’re the binomial coefficients of course!^[3] Therefore, the entries in Pascal’s triangle are all binomial coefficients!

The identity involving binomial coefficients (which we’ll denote by B), namely

follows from the recurrence relation.

What about going the other way? These x + xy terms aren’t as nice as x + y, and (x + y)ⁿ gives binomial coeffients just as well. So, consider the generating function

which has h(x, y) − xh(x, y) − yh(x, y) = 1 and a corresponding recurrence relation

with H(0, 0) = 1. This is actually Pascal’s triangle again, except transformed into a right triangle.

This suggests a notation I sometimes would rather have for the Binomial coefficient:^[4] We define [a, b] to be the number of ways of ordering a of one kind of indistinguishable object among b of another (that is, it’s (a + b) choose a). With this notation, the Binomial theorem would look like

The recurrence relation provides the following identity:

What is the intuition for this? The first object in any permutation is either an object of the first kind or of the second kind; add up each case.

Using similar combinatorial reasoning, the notation can be extended to a multinomial coefficient with a similar identity:

which has the generating function

3. Stirling’s numbers

One reason for keeping the binomial coefficients with the same notation is that there is a striking parallel between it and what are called Stirling’s numbers. Define S(n, k) and C(n, k), where S(n, k) is the number of ways to partition n objects into k non-empty subsets, and where C(n, k) is the number of elements in the symmetric group Sⁿ which have k permutation cycles.^[5] These (including B for parallelism) have the following recurrence relations:

Each considers taking a single object from n objects, and then handling two cases for the object:

• For B: whether it is not chosen or chosen to be among the k chosen objects;
• For S: whether it is in one of the k subsets with others or in its own singleton; or
• For C: whether it is in one of the cycles with others (in which case it can be between any of the n − 1 objects in the permutation) or in its own 1-cycle.

Is there similar notation for S and C like [a, b] for B(a + b, b) which preserves the parallelism? If so, can they be extended into a multi-S and a multi-C? If not, then it seems the parallelism between B and the other two might be something of a red herring.

Like how [a, b] is B(a + b, a), suppose we define [a, b]_S = S(a + b, a). This gives [a, b]_s = a[a, b − 1]_S + [a − 1, b]_S using the above recurrence. I’m far from convinced that this has given us anything.

4. Conclusion

You can use deal with generating functions even if the recurrence relations is over multiple variables. Going the other way, this suggests that polynomial long division can be done by first determining the corresponding recurrence relation.

Also, the seemingly innocuous polynomial

is hiding all of Pascal’s triangle.

5. Update 6/28/17

Chris Orr, an e-mail correspondent, was able to come up with a closed form from the generating function using elementary methods (and surely a good deal of ingenuity). He also thought to put the table into the Sloane On-Line Encyclopedia of Integer Sequences and discovered it is sequence A053538. For 0 ≤ k ≤ n, T(n, k) is the sum over all p in 0 ≤ p ≤ n of the number of ways to place p balls into n slots with k in the rightmost p slots. The Encyclopedia gives the closed form

which one may verify gives the same table as above using the Mathematica command

TableForm[
 Table[Sum[Binomial[p, k] Binomial[n - p, p - k], {p, 0, n}],
   {n, 0, 8}, {k, 0, n}]]

Mathematica is able to figure out SeriesCoefficient[1/(1 - x - x y), {x, 0, n}, {y, 0, k}], but interestingly it gives a very complicated result for SeriesCoefficient[1/(1 - x - x^2 - x y + x^2 y), {x, 0, n}, {y, 0, k}].

def T(n, k) :
    if n < 0 or k < 0 :
        return 0
    valsu = [0, 1] + [0]*k
    valsd = [0, 1] + [0]*k
    for i in xrange(1, n+1) :
        valsu2 = [0]
        valsd2 = [0]
        for j in xrange(0, k + 1) :
            valsu2.append(valsu[j] + valsd[j+1])
            valsd2.append(valsu[j+1] + valsd[j+1])
        valsu = valsu2
        valsd = valsd2
    return valsu[k+1]