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.)

Let’s just go through with finding the recurrence relation — it doesn’t take too long. We define Tu(n, k) to be the number of configurations of n fans so that k remain, given that, behind the last fan (i.e., after these n fans), there is either a standing fan or a brick wall; the idea is that we know from the given that if the last fan is standing, that fan will remain. Let Td(n, k) be the same, except it supposes that the last person is in front of a fan who’s sitting, so, sitting or standing, there is no way the last fan can remain. From these, we get the following two equations, where the left-hand side of the addition corresponds to the last person in line standing, and the right-hand side, sitting:

Tu(n, k) = Tu(n − 1, k − 1) + Td(n − 1, k)
Td(n, k) = Tu(n − 1, k) + Td(n − 1, k)

The most remarkable term is Tu(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 Tu(n, k).[1] But, I’m more interested right now in understanding these mutually-recursive multivariable functions.

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

Td(n, k) = Tu(n, k) + Tu(n − 1, k) − Tu(n − 1, k − 1), 

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

Tu(n, k) = Tu(n − 1, k − 1) + Tu(n − 1, k) + Tu(n − 2, k) − Tu(n − 2, k − 1)

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

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]

1 1
2 1 1
3 3 1 1
5 5 4 1 1
8 10 7 5 1 1
13 18 16 9 6 1 1
21 33 31 23 11 7 1 1
34 59 62 47 31 13 8 1 1

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 cijxiyj for every term cij (in what might be called a bi-sequence). Hence, we let

g(x, y) = Σi, jT(i, j)xiyj

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

xaybg(x, y) = Σi, jT(i, j)xi + ayj + a = Σi, jT(ia, jb)xiyj

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

T ′(n, k) = T(na, kb).

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

xyg(x, y) + xg(x, y) + x2g(x, y) − x2yg(x, y), 

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

g(x, y) = 1 + xyg(x, y) + xg(x, y) + x2g(x, y) − x2yg(x, y), 

and so

g(x, y) = 
1 − xx2xy + x2y

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 + t2 + t3 + …, we can rewrite g in the form

g(x, y) = 1 + (x + x2 + xyx2y) + (x + x2 + xyx2y)2 + (x + x2 + xyx2y)3 + …, 

so to compute T(n, k), we need only keep track of which of these terms contribute to xnyk, 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:

1 1
1 2 1
1 3 3 1
1 4 6 4 1
which is generated by the recurrence P(i, j) = P(i − 1, j) + P(i − 1, j − 1), where P is 0 everywhere outside the triangle, and P(0, 0) = 1, which is the tip of the triangle. Using similar techniques, and letting f(x, y) be the generating function of P, we have
f(x, y) = 1 + xf(x, y) + xyf(x, y), 

and the closed form for the generating function is

f(x, y) = 
1 − xxy

Similarly rewriting using the geometric series, we have

f(x, y) = 1 + (x + xy) + (x + xy)2 + (x + xy)3 + ….

Importantly, each monomial in f comes from exactly one of the (x + xy)n terms, which are each composed of monomials xnynk (and xn 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

B(n, k) = B(n − 1, k) + B(n − 1, k − 1), 

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)n gives binomial coeffients just as well. So, consider the generating function

h(x, y) = 1 + (x + y) + (x + y)2 + … = 
1 − xy

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

H(a, b) = H(a − 1, b) + H(a, b − 1)

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

(x + y)n = Σa + b = n [a, b]xayb.

The recurrence relation provides the following identity:

[a, b] = [a − 1, b] + [a, b − 1].

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:

[a1, a2, …, an] = [a1 − 1, a2, …, an] + [a1, a2 − 1, …, an] + … + [a1, a2, …, an − 1], 

which has the generating function

M(x1, x2, …, xn) = 
1 − x1x2 − … − xn

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 Sn which have k permutation cycles.[5] These (including B for parallelism) have the following recurrence relations:

B(n, k) = B(n − 1, k) + B(n − 1, k − 1)
S(n, k) = kS(n − 1, k) + S(n − 1, k − 1)
C(n, k) = (n − 1)C(n − 1, k) + C(n − 1, k − 1)

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

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

1 − xy

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

T(n, k) = Σ0 ≤ p ≤ n B(p, k)B(np, pk), 

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

 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}].

[1] For completeness, here’s a short Python program to compute T(n, k) in O(nk) time. This performs some extraneous computations out of simplicity, but we’re not trying to break a speed record here.
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]

[2] Like Pascal’s triangle, the data in each diagonal parallel to the major diagonal is described by a polynomial. I couldn’t find any reasonable pattern in these polynomials, though.

[3] Where we’re defining binomial coefficients to the the coefficients of polynomials (x + y)n, for varying n. The “of course” refers to how they are binomial coefficients by definition.

[4] I haven’t seen this notation anywhere else; I’d definitely like to see where others have used similar notation if it’s been done.

[5] Knuth has better notation for these, but, with some amount of not-unnoticed irony, I’m not able to typeset them with my LaTeX-like typesetting software which is used this website. Anyway, Knuth describes these, which are known as Stirling’s numbers of the second and first kind, respectively, much better than I do in Two notes on notation.