{"id":2739,"date":"2020-01-18T12:43:56","date_gmt":"2020-01-18T18:43:56","guid":{"rendered":"https:\/\/laurentlessard.com\/bookproofs\/?p=2739"},"modified":"2020-01-20T00:01:54","modified_gmt":"2020-01-20T06:01:54","slug":"delirious-ducks","status":"publish","type":"post","link":"https:\/\/laurentlessard.com\/bookproofs\/delirious-ducks\/","title":{"rendered":"Delirious ducks"},"content":{"rendered":"

This week’s Riddler Classic<\/a> is about random walks on a lattice:<\/p>\n

\nTwo delirious ducks are having a difficult time finding each other in their pond. The pond happens to contain a 3\u00d73 grid of rocks.<\/p>\n

Every minute, each duck randomly swims, independently of the other duck, from one rock to a neighboring rock in the 3\u00d73 grid \u2014 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.<\/p>\n

If the ducks both start at the middle rock, then on average, how long will it take until they\u2019re at the same rock again? (Of course, there\u2019s a 1\/4 chance that they\u2019ll 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.)<\/p>\n

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?\n<\/p><\/blockquote>\n

Here is my solution:
\n[Show Solution]<\/a><\/p>\n

\n

One duck<\/h3>\n

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<\/a>. 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:<\/p>\n

\"\"<\/p>\n

If we number the states in the same way as above, then the transition matrix for this Markov chain is:
\n\\[
\nP = \\begin{bmatrix}
\n0 & \\tfrac12 & 0 & \\tfrac12 & 0 & 0 & 0 & 0 & 0 \\\\
\n\\tfrac13 & 0 & \\tfrac13 & 0 & \\tfrac13 & 0 & 0 & 0 & 0 \\\\
\n0 & \\tfrac12 & 0 & 0 & 0 & \\tfrac12 & 0 & 0 & 0 \\\\
\n\\tfrac13 & 0 & 0 & 0 & \\tfrac13 & 0 & \\tfrac13 & 0 & 0 \\\\
\n0 & \\tfrac14 & 0 & \\tfrac14 & 0 & \\tfrac14 & 0 & \\tfrac14 & 0 \\\\
\n0 & 0 & \\tfrac13 & 0 & \\tfrac13 & 0 & 0 & 0 & \\tfrac13 \\\\
\n0 & 0 & 0 & \\tfrac12 & 0 & 0 & 0 & \\tfrac12 & 0 \\\\
\n0 & 0 & 0 & 0 & \\tfrac13 & 0 & \\tfrac13 & 0 & \\tfrac13 \\\\
\n0 & 0 & 0 & 0 & 0 & \\tfrac12 & 0 & \\tfrac12 & 0 \\\\
\n\\end{bmatrix}
\n\\]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<\/a>). 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
\n\\[
\n\\mathbf{a}^\\mathsf{T} = \\begin{bmatrix} \\tfrac12 & \\tfrac12 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\end{bmatrix}
\n\\]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
\n\\[
\n\\mathbf{b}^\\mathsf{T} = \\mathbf{a}^\\mathsf{T}P = \\begin{bmatrix} \\tfrac16 & \\tfrac14 & \\tfrac16 & \\tfrac14 & \\tfrac16 & 0 & 0 & 0 & 0 \\end{bmatrix}
\n\\]So in the next step, we might be in states 1, 2, 3, 4, 5 with the associated probabilities above.<\/p>\n

It’s important to distinguish states<\/em> e.g. $\\{1,2,\\dots,9\\}$ from distributions over states<\/em>, 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<\/a>, i.e. concentrated entirely on a particular state $s$, then I’ll use the following notation to denote the corresponding distribution over states:
\n\\[
\n\\left( \\mathbf{e}_s\\right)_i = \\begin{cases} 1 & \\text{if } i=s \\\\ 0 & \\text{otherwise} \\end{cases}
\n\\]So, for example, $\\mathbf{e}_2^\\mathsf{T} = \\begin{bmatrix} 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\end{bmatrix}$.<\/p>\n

Stopping time<\/h3>\n

We are interested in the expected hitting time<\/a>, 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:
\n\\[
\n\\mathbf{q}_s = \\begin{cases}
\n0 & \\text{if }s \\in \\mathcal{T} \\\\
\n1 + \\sum_{i} P_{si} \\mathbf{q}_i & \\text{otherwise}
\n\\end{cases}
\n\\]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<\/a> from dynamic programming<\/a>.<\/p>\n

Define the vector $\\mathbf{t}_i = \\begin{cases} 0 & \\text{if }i \\in \\mathcal{T}\\\\ 1 & \\text{otherwise} \\end{cases}$.
\nWe can rewrite the equation above in concise vector form as:
\n\\[
\n\\mathbf{q} = \\textrm{diag}(\\mathbf{t}) \\left( \\mathbf{1} + P \\mathbf{q} \\right)
\n\\]Using the fact that $\\textrm{diag}(\\mathbf{t})\\mathbf{1} = \\mathbf{t}$, we can further simplify and obtain:
\n\\[
\n\\left( I-\\textrm{diag}(\\mathbf{t}) P \\right) \\mathbf{q} = \\mathbf{t}
\n\\]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:
\n\\[
\n\\begin{bmatrix}
\n1 & -\\tfrac12 & 0 & -\\tfrac12 & 0 & 0 & 0 & 0 & 0 \\\\
\n-\\tfrac13 & 1 & -\\tfrac13 & 0 & -\\tfrac13 & 0 & 0 & 0 & 0 \\\\
\n0 & -\\tfrac12 & 1 & 0 & 0 & -\\tfrac12 & 0 & 0 & 0 \\\\
\n-\\tfrac13 & 0 & 0 & 1 & -\\tfrac13 & 0 & -\\tfrac13 & 0 & 0 \\\\
\n0 & 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\\\
\n0 & 0 & -\\tfrac13 & 0 & -\\tfrac13 & 1 & 0 & 0 & -\\tfrac13 \\\\
\n0 & 0 & 0 & -\\tfrac12 & 0 & 0 & 1 & -\\tfrac12 & 0 \\\\
\n0 & 0 & 0 & 0 & -\\tfrac13 & 0 & -\\tfrac13 & 1 & -\\tfrac13 \\\\
\n0 & 0 & 0 & 0 & 0 & -\\tfrac12 & 0 & -\\tfrac12 & 1 \\\\
\n\\end{bmatrix}
\n\\begin{bmatrix}\\mathbf{q}_1\\vphantom{\\tfrac12} \\\\ \\mathbf{q}_2\\vphantom{\\tfrac12} \\\\ \\mathbf{q}_3\\vphantom{\\tfrac12} \\\\
\n\\mathbf{q}_4\\vphantom{\\tfrac12} \\\\ \\mathbf{q}_5\\vphantom{\\tfrac12} \\\\ \\mathbf{q}_6\\vphantom{\\tfrac12} \\\\
\n\\mathbf{q}_7\\vphantom{\\tfrac12} \\\\ \\mathbf{q}_8\\vphantom{\\tfrac12} \\\\ \\mathbf{q}_9\\vphantom{\\tfrac12} \\end{bmatrix}
\n= \\begin{bmatrix}
\n1\\vphantom{\\tfrac12} \\\\1\\vphantom{\\tfrac12} \\\\1\\vphantom{\\tfrac12} \\\\1\\vphantom{\\tfrac12} \\\\
\n0\\vphantom{\\tfrac12} \\\\1\\vphantom{\\tfrac12} \\\\1\\vphantom{\\tfrac12} \\\\1\\vphantom{\\tfrac12} \\\\1\\vphantom{\\tfrac12}
\n\\end{bmatrix}
\n\\]Inverting this matrix and solving for $\\mathbf{q}$, here is the result in array form:
\n\\[
\n\\begin{array}{|c|c|c|} \\hline
\n\\mathbf{q}_1 = 6 & \\mathbf{q}_2 = 5 & \\mathbf{q}_3 = 6 \\\\ \\hline
\n\\mathbf{q}_4 = 5 & \\mathbf{q}_5 = 0 & \\mathbf{q}_6 = 5 \\\\ \\hline
\n\\mathbf{q}_7 = 6 & \\mathbf{q}_8 = 5 & \\mathbf{q}_9 = 6 \\\\ \\hline
\n\\end{array}
\n\\]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).<\/p>\n

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$.<\/p>\n

Two ducks<\/h3>\n

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$.<\/p>\n

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<\/a>), then the transition matrix for this augmented Markov chain will be $P\\otimes P$, where $\\otimes$ is the Kronecker product<\/a>.<\/p>\n

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<\/em> version of the matrix, which you obtain by enumerating the columns. For example:
\n\\[
\n\\mathrm{vec} \\begin{bmatrix}a_{11} & a_{12} & a_{13} \\\\
\na_{21} & a_{22} & a_{23} \\\\
\na_{31} & a_{32} & a_{33}
\n\\end{bmatrix} =
\n\\begin{bmatrix}
\na_{11} \\\\ a_{21} \\\\ a_{31} \\\\ a_{12} \\\\ a_{22} \\\\ a_{32} \\\\ a_{13} \\\\ a_{23} \\\\ a_{33} \\end{bmatrix}
\n\\]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$.<\/p>\n

A cautionary note<\/h3>\n

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 does<\/e> happen with two ducks! This is because every time a duck moves to a new stone, the parity of the stone changes from odd to even or vice versa. So if one duck started on rock 1 and the other started on rock 2, they would never<\/em> meet.<\/p>\n

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<\/a>, and ignore all the components of $\\mathbf{q}$ for which the two components have different parity. More on this later…<\/p>\n

Computation<\/h3>\n

I computed the result using the following Matlab code:<\/p>\n

\r\n% transition matrix for one duck\r\nP = [ 0 1\/2 0 1\/2 0  0  0  0  0\r\n     1\/3 0 1\/3 0 1\/3 0  0  0  0\r\n      0 1\/2 0  0  0 1\/2 0  0  0\r\n     1\/3 0  0  0 1\/3 0 1\/3 0  0\r\n      0 1\/4 0 1\/4 0 1\/4 0 1\/4 0\r\n      0  0 1\/3 0 1\/3 0  0  0 1\/3\r\n      0  0  0 1\/2 0  0  0 1\/2 0\r\n      0  0  0  0 1\/3 0 1\/3 0 1\/3\r\n      0  0  0  0  0 1\/2 0 1\/2 0 ];\r\n\r\nT = kron(P,P);     % transition matrix\r\nt = 1-vec(eye(9)); % terminal states\r\n\r\n% starting distribution\r\ns = kron([0 0 0 0 1 0 0 0 0],[0 0 0 0 1 0 0 0 0])';\r\n\r\n% compute stopping time as a rational number\r\nstop_time = s'*T*( (eye(81)-diag(t)*T)\\t ) + 1\r\nrats(stop_time)\r\n<\/pre>\n
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:<\/p>\n
\"\"<\/a><\/p>\n
The staggered decrease of the probabilities is a real effect; it’s not a result of approximation error!<\/p>\n
Many ducks<\/h3>\n
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<\/em> 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.<\/p>\n
After all this work, we end up with the following result (for up to 6 ducks):<\/p>\n
\"\"<\/a><\/p>\n
The values are as follows:<\/p>\n\n\n\n\n\n\n\n
Number of ducks<\/th>\nExpected rendez-vous time<\/th>\n<\/tr>\n
2<\/td>\n4.9054<\/td>\n<\/tr>\n
3<\/td>\n18.4360<\/td>\n<\/tr>\n
4<\/td>\n66.7420<\/td>\n<\/tr>\n
5<\/td>\n237.3955<\/td>\n<\/tr>\n
6<\/td>\n825.3364<\/td>\n<\/tr>\n<\/table>\n
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.<\/p>\n
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$.<\/p>\n