October 8, 2018

\(\DeclareMathOperator{\res}{Res}\)

This one’s another puzzle from work:

Consider a pigeon coop with \(n\) pigeonholes, arranged in a straight line. When a pigeon arrives at the coop, it will roost in a pigeonhole only if it is empty, and both neighboring pigeonholes are also empty. It selects such a pigeonhole uniformly at random, enters the pigeonhole, and does not leave. At some point, the coop will fill up, but not every pigeonhole will be occupied. What is the expected density of pigeons in the coop, as \(n\) grows large?

If you run a few simulations, you get that it’s about \(0.432332\ldots\). But this isn’t any easily recognizable number. What is it in closed form?

This problem illustrates one of the things I find really cool about math: the boundaries between different disciplines are essentially fictitious. This is a combinatorics problem, and so we might expect to be using arguments involving counting, bijections, and other finite tools. But instead we’ll sprint as fast as we can into the realm of analysis and solve the problem there.

Let \(a_n\) be the expected number of pigeons for a coop with \(n\) holes. Then we can come up with a recurrence relation for \(a_n\).

Consider what happens when the first pigeon arrives in an unoccupied coop. If it arrives in the first hole, then we can imagine deleting the first hole and its neighbor from the coop, leaving us with an unoccupied coop of size \(n - 2\). If it lands in the last hole, we have the same situation. Otherwise, it lands somewhere in the middle; when a pigeon comes to rest in the \(k\)th hole (I’m going to \(1\)-index, by the way), it splits the coop into two smaller coops, one with \(k - 2\) holes, and the other with \(n - k - 1\) holes. Since each hole is equally likely, we can average over all values of \(k\) to get a first draft of our recurrence relation:

$$ a_n = 1 + \frac{1}{n} \left( a_{n-2} + a_{n-2} + \sum_{k=2}^{n-1} (a_{k-2} + a_{n-k-1}) \right) $$

This can be prettied up with some mild re-indexing:

$$ a_n = 1 + \frac{2}{n} \sum_{k=0}^{n-2} a_k $$

We can do even better though! If we consider \(n a_n - (n-1) a_{n-1}\), we can collapse most of our terms:

$$
\begin{align*}
n a_n - (n-1) a_{n-1} &= \left( n + 2 \sum_{k=0}^{n-2} a_k \right) - \left( n-1 + 2 \sum_{k=0}^{n-1} a_k \right) \\
n a_n - (n-1) a_{n-1} &= 1 + 2 a_{n-2} \\
a_n &= \frac{1}{n} ( 1 + (n-1) a_{n-1} + 2 a_{n-2} )
\end{align*}
$$

This isn’t a linear recurrence relation, so we can’t apply linear algebra tricks to it. So we fall back on the Swiss Army knife of recurrence relations: the generating function.

Let \(G(z) = a_0 + a_1 z + a_2 z^2 + a_3 z^3 + \cdots\). We don’t know what this function is yet, but we can use the recurrence relation to pin down what it is.

\begin{align*}
G(z) &= \sum_{n=0}^\infty a_n z^n \\
G'(z) &= \sum_{n=1}^\infty n a_n z^{n-1} \\
&= a_1 + \sum_{n=2}^\infty n a_n z^{n-1} \\
&= a_1 + \sum_{n=2}^\infty \left( 1 + (n-1) a_{n-1} + 2 a_{n-2} \right) z^{n-1}
\end{align*}

Dealing with the three pieces separately makes this much easier to read (and also to write *wink*):

$$ \sum_{n=2}^\infty z^{n-1} = \frac{z}{1 - z} $$

$$ \sum_{n=2}^\infty (n-1) a_{n-1} z^{n-1} = \sum_{n=1}^\infty n a_n z^n = z G'(z) $$

$$ \sum_{n=2}^\infty 2 a_{n-2} z^{n-1} = 2 \sum_{n=0}^\infty a_n z^{n+1} = 2z G(z) $$

Putting it all together, we get a differential equation for \(G(z)\):

$$ G'(z) = 1 + \frac{z}{1 - z} + z G'(z) + 2z G(z) $$

Cleaning it up a little, we see that it’s first order and linear, so we can put those diff eq skills to use:

$$ G'(z) = \frac{2z}{1 - z} G(z) + \frac{1}{(1 - z)^2} $$

The details aren’t super important, but basically you use an integrating factor and get:

$$ G(z) = \frac{1 + C e^{-2z}}{2(z-1)^2} $$

What should \(C\) be? We’ll have to use our initial conditions, and one of them is particularly straightforward: \(G(0) = a_0\), which we know is \(0\), and so \(C = -1\).

At this point, let’s stop and recollect our thoughts. We’ve defined a function \(G(z)\) whose power series coefficients are \(a_n\), the average number of pigeons in a coop of size \(n\). Our solution is now encoded in quite a peculiar way: how fast do the coefficients of \(G(z)\) grow?

To figure this out, let’s put the “analytic” in “analytic combinatorics”, and consider some contour integrals. Fix some \(R > 1\), and define \(I_n\) to be the integral of \(G(z)/z^{n+1}\) around the circle of radius \(R\) at the origin (taken counter-clockwise).

What is \(I_n\)? We can evaluate it using the residue theorem. There are two poles, one at \(z = 0\), and the other at \(z = 1\). The former is easy to compute; the residue is the coefficient on the \(z^{-1}\) term, which is exactly \(a_n\). The second does not admit such a nice description, and so we compute it the usual way:

\begin{align*}
\res\left( \frac{G(z)}{z^{n+1}}, 1\right) &= \lim_{z \to 1} \frac{d}{dz} (z-1)^2 \frac{G(z)}{z^{n+1}} \\
&= \lim_{z \to 1} \frac{d}{dz} \frac{1 - e^{-2z}}{2 z^{n+1}} \\
&= \lim_{z \to 1} \frac{2 z e^{-2z} - (n+1)(1 - e^{-2z})}{2 z^{n+2}} \\
&= \frac{(n+3)e^{-2} - (n+1)}{2}
\end{align*}

So \(\frac{1}{2 \pi i} I_n = a_n + \frac{(n+3)e^{-2} - (n+1)}{2}\). What good does this do us?

If you’ve seen this trick before, you know that \(I_n\) drops exponentially to \(0\) as \(n\) increases, but if not, here’s the justification. Let \(M\) be the largest value (in terms of absolute value) that \(G\) attains on the circle \(|z| = R\). Then the triangle inequality tells us:

$$ | I_n | = \left| \int_{C_R} \frac{G(z)}{z^{n+1}}~dz \right| \le \int_{C_R} \left| \frac{G(z)}{z^{n+1}} \right|~dz \le \int_{C_R} \frac{M}{R^{n+1}}~dz = \frac{2 \pi M}{R^n} $$

So as \(n \to \infty\), \(I_n\) drops to \(0\), and so \(a_n\) approaches \(\frac{(n+1)-(n+3)e^{-2}}{2}\). Therefore, the expected density of pigeons, \(a_n/n\), approaches \((1 - e^{-2})/2\), or about \(0.432332\).

There were other solutions that people came up with for this problem, but what I really like about this one is that it demonstrates a way to approach these problems in general, and (at least IMO) it’s a pretty unexpected one. If someone asked me to figure out how fast the coefficients of a power series grow, the residue theorem would not be the first thing on my mind. And yet, not only does it get the job done, it works for many other similar problems, in essentially the same way. I’m not much of an analysis person, but my understanding is that this kind of trick is common in analytic combinatorics, and I think that’s pretty cool!