This week’s Riddler Classic is about random walks on a lattice:

Two delirious ducks are having a difficult time finding each other in their pond. The pond happens to contain a 3×3 grid of rocks.

Every minute, each duck randomly swims, independently of the other duck, from one rock to a neighboring rock in the 3×3 grid — up, down, left or right, but not diagonally. So if a duck is at the middle rock, it will next swim to one of the four side rocks with probability 1/4. From a side rock, it will swim to one of the two adjacent corner rocks or back to the middle rock, each with probability 1/3. And from a corner rock, it will swim to one of the two adjacent side rocks with probability 1/2.

If the ducks both start at the middle rock, then on average, how long will it take until they’re at the same rock again? (Of course, there’s a 1/4 chance that they’ll swim in the same direction after the first minute, in which case it would only take one minute for them to be at the same rock again. But it could take much longer, if they happen to keep missing each other.)

Extra credit: What if there are three or more ducks? If they all start in the middle rock, on average, how long will it take until they are all at the same rock again?

Here is my solution:

[Show Solution]

### One duck

We will solve this problem by starting with a simpler one and building up to the full problem. Let’s consider a single duck on the grid of rocks. At each step, the duck moves randomly to one of the other nearby rocks. We can describe this process using a Markov chain. Each state is a different rock, and the directed edges are labeled with the transition probabilities. Here is a diagram of what the Markov chain looks like:

If we number the states in the same way as above, then the transition matrix for this Markov chain is:

\[

P = \begin{bmatrix}

0 & \tfrac12 & 0 & \tfrac12 & 0 & 0 & 0 & 0 & 0 \\

\tfrac13 & 0 & \tfrac13 & 0 & \tfrac13 & 0 & 0 & 0 & 0 \\

0 & \tfrac12 & 0 & 0 & 0 & \tfrac12 & 0 & 0 & 0 \\

\tfrac13 & 0 & 0 & 0 & \tfrac13 & 0 & \tfrac13 & 0 & 0 \\

0 & \tfrac14 & 0 & \tfrac14 & 0 & \tfrac14 & 0 & \tfrac14 & 0 \\

0 & 0 & \tfrac13 & 0 & \tfrac13 & 0 & 0 & 0 & \tfrac13 \\

0 & 0 & 0 & \tfrac12 & 0 & 0 & 0 & \tfrac12 & 0 \\

0 & 0 & 0 & 0 & \tfrac13 & 0 & \tfrac13 & 0 & \tfrac13 \\

0 & 0 & 0 & 0 & 0 & \tfrac12 & 0 & \tfrac12 & 0 \\

\end{bmatrix}

\]The way to interpret this matrix is that $P_{ij}$ is the probability that we will transition from $i$ to $j$. This explains why the rows sum to 1 (the matrix is right-stochastic). We can express this fact neatly using matrix multiplication as: $P \mathbf{1} = \mathbf{0}$, where $\mathbf{0}$ and $\mathbf{1}$ are column vectors of all-zeros and all-ones, repectively. Mathematically, we can propagate probability distributions through this Markov chain via matrix multiplication. If we have some initial distribution $\mathbf{a}$ across the states, for example

\[

\mathbf{a}^\mathsf{T} = \begin{bmatrix} \tfrac12 & \tfrac12 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}

\]i.e. we are equally likely to be in states 1 or 2. Then at the next step, the probability distribution will transition $\mathbf{a}\to\mathbf{b}$ with

\[

\mathbf{b}^\mathsf{T} = \mathbf{a}^\mathsf{T}P = \begin{bmatrix} \tfrac16 & \tfrac14 & \tfrac16 & \tfrac14 & \tfrac16 & 0 & 0 & 0 & 0 \end{bmatrix}

\]So in the next step, we might be in states 1, 2, 3, 4, 5 with the associated probabilities above.

It’s important to distinguish *states* e.g. $\{1,2,\dots,9\}$ from *distributions over states*, which are the vectors of probabilities such as the $\mathbf{a}$ and $\mathbf{b}$ used above. In the case where a distribution across states is degenerate, i.e. concentrated entirely on a particular state $s$, then I’ll use the following notation to denote the corresponding distribution over states:

\[

\left( \mathbf{e}_s\right)_i = \begin{cases} 1 & \text{if } i=s \\ 0 & \text{otherwise} \end{cases}

\]So, for example, $\mathbf{e}_2^\mathsf{T} = \begin{bmatrix} 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \end{bmatrix}$.

### Stopping time

We are interested in the expected hitting time, which is the number of steps it will take on average to travel from some initial state $s$ to some target state $t \in \mathcal{T}$. To keep this general, I’m assuming there can be multiple target states, indicated by the set $\mathcal{T}$. Let’s define $\mathbf{q}_s$ to be the expected hitting time to get from state $s$ to any state in our terminal set $\mathcal{T}$. It turns out that $\mathbf{q}$ satisfies the recursive relationship:

\[

\mathbf{q}_s = \begin{cases}

0 & \text{if }s \in \mathcal{T} \\

1 + \sum_{i} P_{si} \mathbf{q}_i & \text{otherwise}

\end{cases}

\]The first case is clear: if we start in the terminal set, then we have already arrived, so the hitting time is zero. If we are outside the terminal set, then the expected hitting time will be $1$ plus the weighted sum of the hitting times of wherever we end up after the next transition. The recursion above is essentially the Bellman equation from dynamic programming.

Define the vector $\mathbf{t}_i = \begin{cases} 0 & \text{if }i \in \mathcal{T}\\ 1 & \text{otherwise} \end{cases}$.

We can rewrite the equation above in concise vector form as:

\[

\mathbf{q} = \textrm{diag}(\mathbf{t}) \left( \mathbf{1} + P \mathbf{q} \right)

\]Using the fact that $\textrm{diag}(\mathbf{t})\mathbf{1} = \mathbf{t}$, we can further simplify and obtain:

\[

\left( I-\textrm{diag}(\mathbf{t}) P \right) \mathbf{q} = \mathbf{t}

\]For the case of one duck, if we set the terminal set to be $\mathcal{T} = \{5\}$, then we substitute $\mathbf{t} = \begin{bmatrix} 1 & 1 & 1 & 1 & 0 & 1 & 1 & 1 & 1 \end{bmatrix}^\mathsf{T}$ above and obtain:

\[

\begin{bmatrix}

1 & -\tfrac12 & 0 & -\tfrac12 & 0 & 0 & 0 & 0 & 0 \\

-\tfrac13 & 1 & -\tfrac13 & 0 & -\tfrac13 & 0 & 0 & 0 & 0 \\

0 & -\tfrac12 & 1 & 0 & 0 & -\tfrac12 & 0 & 0 & 0 \\

-\tfrac13 & 0 & 0 & 1 & -\tfrac13 & 0 & -\tfrac13 & 0 & 0 \\

0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\

0 & 0 & -\tfrac13 & 0 & -\tfrac13 & 1 & 0 & 0 & -\tfrac13 \\

0 & 0 & 0 & -\tfrac12 & 0 & 0 & 1 & -\tfrac12 & 0 \\

0 & 0 & 0 & 0 & -\tfrac13 & 0 & -\tfrac13 & 1 & -\tfrac13 \\

0 & 0 & 0 & 0 & 0 & -\tfrac12 & 0 & -\tfrac12 & 1 \\

\end{bmatrix}

\begin{bmatrix}\mathbf{q}_1\vphantom{\tfrac12} \\ \mathbf{q}_2\vphantom{\tfrac12} \\ \mathbf{q}_3\vphantom{\tfrac12} \\

\mathbf{q}_4\vphantom{\tfrac12} \\ \mathbf{q}_5\vphantom{\tfrac12} \\ \mathbf{q}_6\vphantom{\tfrac12} \\

\mathbf{q}_7\vphantom{\tfrac12} \\ \mathbf{q}_8\vphantom{\tfrac12} \\ \mathbf{q}_9\vphantom{\tfrac12} \end{bmatrix}

= \begin{bmatrix}

1\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\

0\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\1\vphantom{\tfrac12} \\1\vphantom{\tfrac12}

\end{bmatrix}

\]Inverting this matrix and solving for $\mathbf{q}$, here is the result in array form:

\[

\begin{array}{|c|c|c|} \hline

\mathbf{q}_1 = 6 & \mathbf{q}_2 = 5 & \mathbf{q}_3 = 6 \\ \hline

\mathbf{q}_4 = 5 & \mathbf{q}_5 = 0 & \mathbf{q}_6 = 5 \\ \hline

\mathbf{q}_7 = 6 & \mathbf{q}_8 = 5 & \mathbf{q}_9 = 6 \\ \hline

\end{array}

\]So if we start at node 1, it will take us on average 6 moves to reach node 5. It is by no means obvious that the expected hitting times should be integers, since they represent an average path length over an infinite number of possible paths. They just turn out to be integers in this case. The symmetry observed here also makes sense, and we could leverage it to simplify the problem to a Markov chain with only 3 states, but I’ll leave that as an exercise to the reader! (the reason I didn’t do this from the start is that I want to generalize easily to the multi-duck case, and the symmetries there are more complicated).

What about the variant where we start at the origin but we want to know how much time it will take us to return to the origin? In this case, we simply start counting from our second move, and we add 1 to the answer to account for the move we skipped. If we start at node 5, i.e. our initial distribution is $\mathbf{e}_5^\mathsf{T}$, then our distribution on the second turn is $\mathbf{e}_5^\mathsf{T} P$, and so the expected number of moves to start at 5 and return to 5 is: $\mathbf{e}_5^\mathsf{T} P \mathbf{q} + 1 = 6$.

### Two ducks

At first glance, the two-duck version might seem much more challenging than the on-duck version. As we shall see, the problem is essentially the same once we look at it the right way. The key is to imagine a Markov chain where instead of the states being $\{1,2,\dots,9\}$, the states are $\{(1,1), (1,2), \dots, (9,9)\}$. In other words, there are 81 states, consisting of all possible pairs $(s_1,s_2)$, where $s_i \in \{1,2,\dots,9\}$ is the position of duck $i$.

What would be the transition matrix for such a Markov chain? The transition probability $(a_1, a_2) \to (b_1, b_2)$ is simply $P_{a_1,b_1} P_{a_2,b_2}$, i.e. the product of the respective transition probabilities for each duck. This means that if we order our states according to the first duck first (lexicographical ordering), then the transition matrix for this augmented Markov chain will be $P\otimes P$, where $\otimes$ is the Kronecker product.

So we can solve the problem in a very similar fashion. This time, our terminal set consists of all the nodes $(s_1, s_2)$ where $s_1=s_2$. There are 9 such nodes, and we can form the associated $\mathbf{t}$ by evaluating $\mathbf{t} = \mathbf{1}-\textrm{vec}(I)$, where $\mathrm{vec}$ is the *vectorized* version of the matrix, which you obtain by enumerating the columns. For example:

\[

\mathrm{vec} \begin{bmatrix}a_{11} & a_{12} & a_{13} \\

a_{21} & a_{22} & a_{23} \\

a_{31} & a_{32} & a_{33}

\end{bmatrix} =

\begin{bmatrix}

a_{11} \\ a_{21} \\ a_{31} \\ a_{12} \\ a_{22} \\ a_{32} \\ a_{13} \\ a_{23} \\ a_{33} \end{bmatrix}

\]the new equation for $\mathcal{q}$ is $\bigl( I-\textrm{diag}(\mathbf{t}) (P\otimes P) \bigr) \mathbf{q} = \mathbf{t}$ and the solution to the problem where you start at $(5,5)$ and want to know the expect time it takes for us to reach the terminal set (where the two ducks are on the same rock) can be computed following the same logic as in the one-duck case. That is, $\mathbf{e}_{(5,5)}^\mathsf{T} (P\otimes P) \mathbf{q} + 1$.

### A cautionary note

While it may seem like this is straightforward algebra, there is a major issue I neglected to discuss: invertibility. When we solve the equation for $\mathbf{q}$ in the one-duck case, the equation looks like $A\mathbf{q} = \mathbf{t}$, and the solution is $\mathbf{q} = A^{-1}\mathbf{t}$. But what if $A$ isn’t invertible? This corresponds to cases where the Markov chain is not connected. So if you start on one island and the terminal set is on another island, you’ll never get there; the stopping time is infinite. That doesn’t happen with one duck, but it *never* meet.

From a linear algebra standpoint, it just means that even though the matrix is not invertible, if we restrict our attention to the “island” (read: subspace) where the two ducks have the same parity, then the associated transition matrix for that island will be invertible and everything will be fine. From a practical standpoint, the equation above can still be solved; we just need to replace the inverse by a pseudoinverse, and ignore all the components of $\mathbf{q}$ for which the two components have different parity. More on this later…

### Computation

I computed the result using the following Matlab code:

% transition matrix for one duck P = [ 0 1/2 0 1/2 0 0 0 0 0 1/3 0 1/3 0 1/3 0 0 0 0 0 1/2 0 0 0 1/2 0 0 0 1/3 0 0 0 1/3 0 1/3 0 0 0 1/4 0 1/4 0 1/4 0 1/4 0 0 0 1/3 0 1/3 0 0 0 1/3 0 0 0 1/2 0 0 0 1/2 0 0 0 0 0 1/3 0 1/3 0 1/3 0 0 0 0 0 1/2 0 1/2 0 ]; T = kron(P,P); % transition matrix t = 1-vec(eye(9)); % terminal states % starting distribution s = kron([0 0 0 0 1 0 0 0 0],[0 0 0 0 1 0 0 0 0])'; % compute stopping time as a rational number stop_time = s'*T*( (eye(81)-diag(t)*T)\t ) + 1 rats(stop_time)

The resulting expected stopping time was $\frac{363}{74}$, or approximately $4.905$ steps. We can also find an approximate answer via simulation. Here is what you get by simulating one million trials:

The staggered decrease of the probabilities is a real effect; it’s not a result of approximation error!

### Many ducks

If we have $n$ ducks, we can clearly just generalize the approach used above. This time, the transition matrix will be $\underbrace{P\otimes \cdots \otimes P}_{n\text{ times}}$. The issue here is that our transition matrix will be quite large: $9^n \times 9^n$, to be precise. One way to shrink the number of states is to return to the idea of parity. The reason our transition matrix is large is because it describes *all* transitions between all possible duck configurations. Since we know all ducks start on the same node and that parity will be preserved, we only need to worry about the subset of states for which all ducks have the same parity. Namely, there will be $5^n$ odd states and $4^n$ even states. So we can restrict ourselves to a smaller $(5^n+4^n)\times(5^n+4^n)$ transition matrix. This is still large, but it’s an improvement. We will also leverage the fact that the transition matrix will be sparse, which will help with computation.

After all this work, we end up with the following result (for up to 6 ducks):

The values are as follows:

Number of ducks | Expected rendez-vous time |
---|---|

2 | 4.9054 |

3 | 18.4360 |

4 | 66.7420 |

5 | 237.3955 |

6 | 825.3364 |

It’s clear the expected rendez-vous time grows exponentially, but I haven’t found a practical way to approximate it or bound it. Unfortunately, it’s pretty difficult to extend my approach beyond $n=6$ ducks. To give you an idea of scale, our naive transition matrix would have had $9^6=531,441$ rows and columns. After our reduction procedure, we’re down to $5^6+4^6=19,721$. The transition matrix is sparse (about 98.5% of its entries are zeros), but that still leaves about 6 million nonzero entries, which makes computing $\mathbf{q}$ a challenge.

Further reductions are possible if we leverage symmetries. Let’s say we have 2 ducks. If we don’t do any reductions, there are $9^2=81$ states. If we reduce using even/odd parity, we drop to $5^2+4^2 = 41$.

- From a stopping time standpoint, the only thing that matters is how many ducks are on each stone; it doesn’t matter which ducks are on which stones. This immediately allows us to reduce $5^n+4^n$ to ${n+4 \choose 4} + {n+3\choose 3}$, which is a dramatic improvement. In the two-duck case, this brings us to $25$.
- Deeper symmetries occur as well. For example, a pair of ducks on 1 and 3 is equivalent to a pair of ducks on 7 and 9 by rotational symmetry. In the two-duck case, this brings us to $8$ states. Counting and bookkeeping these states might require some group theory…

Ultimately, these reductions would have a dramatic effect, reducing the scaling from exponential to polynomial. I didn’t implement any of these additional reductions because at this point (up to $n=6$) the exponential trend was clear. I didn’t think that extending the trend further would be particularly enlightening.

Hi Laurent,

Very nice. Your approach is much more elegant than what I did.

I approached it a little differently, building a Fokker-Planck equation to evolve from a starting configuration through the probabilities of each successive configuration, removing cases where the ducks join up again, and I get the same answers.

Regarding an approximate answer:

For large numbers of ducks, you might imagine that, before a move, the ducks are “randomly” distributed between the points of a given parity, and then compute the chance that they all converge to the same spot.

When the ducks are on the 4 edge-states (2,4,6,8 in your notation), each duck has a 1/3 chance to go to the middle (5) and a 1/6 chance each to go to a given corner (1,3,7,9). Therefore the chance that they all go to the middle is 1/3^N and the chance that they all go to the same corner is 4/6^N (4 for the 4 corners).

When the ducks are on an odd state, each duck has a 1/4 chance to go to a given edge, so there is a 1/4^N chance that they all go to (2) and the same for each of (4,6,8) for a total of 4/4^N chance that they will converge on the next move.

Averaging over these two cases, the chance of “finishing” on a given move is

( 1/3^N + 4/6^N + 4/4^N)/2

and the average number of moves is

2/( 1/3^N + 4/6^N + 4/4^N)

which is:

4.235 for N=2

16.94 for N=3

64.40 for N=4

234.31 for N=5

821.68 for N=6

2794.56 for N=7

etc (asymptotically, adding a duck lengthens the process by a factor of 3)

Of course this is not exact: the location of each duck is correlated and not completely random. But for large N the correlation is very small, and the approximation gets better. I was surprised how well it already worked for N=2…

I did exactly what you did, Guy. At side-to-non-side transitions, there’s always a 1/3^N chance of all ducks meeting in the center, so the expected first meeting of all ducks in the center is 2*3^N. For large N, all ducks meeting elsewhere becomes vanishingly improbable compared to meeting in the center, so that seems to be the limiting expectation, and it’s a good approximation already for N > 10. (See plot at link.)

Glad to see some new postings — I almost thought I saw some Julia code in here 😉

Also, good catch on the issues with periodic behavior. It’s a nuisance, and I’d be inclined to tweak the problem to be a lazy chain ($\delta \in (0,1)$ so $P’:= \delta I + (1-\delta)P$) or possibly look at an ensemble average. When I first read the problem I was hoping for a Kronecker product structure, but I thought it was too good to be true.

A couple other things — this is a *reversible* markov chain and in this particular case the steady probability vector entry for a given state is proportional to the number of connections it has — which is $\{\frac{1}{12}, \frac{1}{8}, \frac{1}{6}\}$ for the single duck case.

One final thing, since $\pi_i = \frac{1}{\bar{X_i}}$ i.e. the perron / steady state vector has probability given by the inverse of the expected time to return to state i given a start in state i, *and* the Kronecker product has nice properties with respect to eigenvalues and eigenvectors we can say that in the $n$ duck case, if they all start in the same location, the expected time until they all meet again at that same location is $\{12^n, 8^n, 6^n\}$. I suppose this gives a very crude upper bound on expected time until rendez-vous for the more general case being contemplated in this Riddler (i.e. meeting in some arbitrary box)

This problem is conceptually nice when we use $P’:= \delta I + (1-\delta)P$ for any $\delta \in (0,1)$. But the original problem with $\delta := 0$ has quite a bit of odd behavior. Algebraically speaking the oddities come from having an eigenvalue of +1 and -1, so when we do a kron product we have maximal eigenvalues given by $1\cdot 1=1$ and $-1\cdot-1=1$ and when we kron prod 3 matrices we have $1\cdot 1 \cdot 1 =1$ but also $1 \cdot -1 \cdot -1 = 1$, as well as $-1 \cdot 1 \cdot -1 = 1$, and $-1 \cdot -1 \cdot 1 = 1$, and so on for higher degrees of kronecker products. Having $k$ copies of the eigenvalue 1 tells us there are k disjoint, recurrent, graphs in here (perron frobenius). Each disjoint graph has its own real non-negative eigenvector that is disjoint (orthogonal) to the other eigenvectors associated with other graphs/ recurrent classes. So when we make some state (say the ‘middle’ one) an absorbing state, this doesn’t ‘hit’ k-1 of the eigenvectors with eigenvalue 1. This problem also makes a very good illustration of how the pseudo inverse is not in general continuous, but the actual inverse (when $\delta \in (0,1)$) is continuous. The non-lazy chain (i.e. stated problem) really is pecuiliar due to what happens when we mix periodicity (more than one eigenvalue with modulus one for a communicating class) and kronecker products.