# Expected number of runs

Let $S[16]$ be a binary array i.e, elements of $S$ are 0/1 with elements $S[i]$ are taken uniformly and independently form $\{0,1\}$. Let $k$ be a random element taken uniformly from $\{0,1\}$. I have run the following process.

1. Take $a,b,c$ randomly and independently from $[1,16]$.
2. Change $k=k+S[a]$
3. $S[b]=k+S[c]$.

Run the whole process (1,2,3) until all elements of $S$ are filled with 0. What is expected number of runs for this event?

All additions are binary i.e, in $Z_2$.

-
Are you implying that this is independent of the initial content of $S$? Or did you forget to specify that? Or are the $S[i]$ also initially uniformly and independently distributed? Also, I presume you mean that $k$ is uniformly distributed over $\{0,1\}$? And that addition is carried out in $\mathbb Z_2$? – joriki Jul 24 '11 at 11:47
$S[i]$ are initially uniformly and independently distributed over $\{0,1\}$. $K$ is also uniformly distributed over $\{0,1\}$. And additions are in $Z_2$ – diamond kukur Jul 24 '11 at 12:07
Best to include that information in the question so people don't have to read through the comments to understand the question. – joriki Jul 24 '11 at 12:10
Note that it doesn't make any difference to your process if the array is shuffled after each iteration of the steps. Thus, it should be sufficient to track the number $n$ of elements of $S$ that are 1. In particular, it should be sufficient to analyze the Markov chain with state $(n,k) \in \{0,\ldots,16\} \times \{0,1\}$. – Ilmari Karonen Jul 24 '11 at 12:10

Note that it doesn't make any difference to your process if the array is shuffled after each iteration of the steps. Thus, it should be sufficient to track the number $n$ of elements of $S$ that are 1.

In particular, it should be sufficient to analyze the Markov chain with state $(n,k) \in \{0,\ldots,16\} \times \{0,1\}$. The transitions of this chain are given by

$$k_{t+1} = \begin{cases} 1-k_t & \text{with probability }\tfrac{n_t}{16} \\ k_t & \text{otherwise} \end{cases}$$

and

$$n_{t+1} = \begin{cases} n_t+1 & \text{with probability }\tfrac{n_t}{16}\left(1-\tfrac{n_t}{16}\right)\text{ if }k_{t+1}=0 \\ n_t-1 & \text{with probability }\tfrac{n_t}{16}\left(1-\tfrac{n_t}{16}\right)\text{ if }k_{t+1}=0 \\ n_t+1 & \text{with probability }\left(1-\tfrac{n_t}{16}\right)^2\text{ if }k_{t+1}=1 \\ n_t-1 & \text{with probability }\left(\tfrac{n_t}{16}\right)^2\text{ if }k_{t+1}=1 \\ n_t & \text{otherwise.} \end{cases}$$

If $n$ hits $0$, the process ends; all we need is the expected time to hit $(0,0)$ or $(0,1)$ from the given initial state ($n \sim \operatorname{Bin}(16, \frac 1 2)$, $k$ uniform). Let me write out the transition matrix and get back on that...

OK, since the transition matrices are $2N \times 2N$, writing down the matrix for $N=16$ here is rather impractical. Let me show instead a cut-down example with $N = 4$.

In general, it's convenient to write the transition matrix as $M = M_n M_k$, where $M_k$ gives the transition probabilities for step 2 and $M_n$ for step 3 in your process. Encoding the state $(n,k)$ as index $i = 2n+k-1$ and working with column probability vectors, the transition matrices for $N=4$ are

$$M_k = \frac 1 4 \begin{bmatrix} 3 & 1 & 0 & 0 & 0 & 0 & 0 & 0 \\ 1 & 3 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 2 & 2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 2 & 2 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 1 & 3 & 0 & 0 \\ 0 & 0 & 0 & 0 & 3 & 1 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 4 \\ 0 & 0 & 0 & 0 & 0 & 0 & 4 & 0 \end{bmatrix} \text{ and } M_n = \frac 1 {16} \begin{bmatrix} 10 & 0 & 4 & 0 & 0 & 0 & 0 & 0 \\ 0 & 6 & 0 & 4 & 0 & 0 & 0 & 0 \\ 3 & 0 & 8 & 0 & 3 & 0 & 0 & 0 \\ 0 & 9 & 0 & 8 & 0 & 9 & 0 & 0 \\ 0 & 0 & 4 & 0 & 10 & 0 & 0 & 0 \\ 0 & 0 & 0 & 4 & 0 & 6 & 0 & 16 \\ 0 & 0 & 0 & 0 & 3 & 0 & 16 & 0 \\ 0 & 0 & 0 & 0 & 0 & 1 & 0 & 0 \end{bmatrix},$$

giving

$$M = \frac 1 {64} \begin{bmatrix} 30 & 10 & 8 & 8 & 0 & 0 & 0 & 0 \\ 6 & 18 & 8 & 8 & 0 & 0 & 0 & 0 \\ 9 & 3 & 16 & 16 & 3 & 9 & 0 & 0 \\ 9 & 27 & 16 & 16 & 27 & 9 & 0 & 0 \\ 0 & 0 & 8 & 8 & 10 & 30 & 0 & 0 \\ 0 & 0 & 8 & 8 & 18 & 6 & 64 & 0 \\ 0 & 0 & 0 & 0 & 3 & 9 & 0 & 64 \\ 0 & 0 & 0 & 0 & 3 & 1 & 0 & 0 \end{bmatrix}.$$

To calculate the expected survival time of the process, we need the fundamental matrix $F = I + M + M^2 + \cdots = (I - M)^{-1}$. The elements $F_{ij}$ give the expected number times the process visits state $i$ if started from state $j$. The calculation of $F$ is straightforward, but yields somewhat ungainly fractions even for $N = 4$; rounded to 2 decimal digits,

$$F \approx \begin{bmatrix} 4.96 & 4.00 & 4.48 & 4.48 & 4.48 & 4.48 & 4.48 & 4.48 \\ 2.40 & 4.00 & 3.20 & 3.20 & 3.20 & 3.20 & 3.20 & 3.20 \\ 3.96 & 4.50 & 6.48 & 5.48 & 5.67 & 5.79 & 5.79 & 5.79 \\ 6.12 & 7.50 & 8.56 & 9.56 & 9.37 & 9.25 & 9.25 & 9.25 \\ 3.36 & 4.00 & 5.01 & 5.01 & 6.65 & 6.04 & 6.04 & 6.04 \\ 3.36 & 4.00 & 5.01 & 5.01 & 5.83 & 6.86 & 6.86 & 6.86 \\ 0.84 & 1.00 & 1.25 & 1.25 & 1.54 & 1.64 & 2.64 & 2.64 \\ 0.21 & 0.25 & 0.31 & 0.31 & 0.40 & 0.39 & 0.39 & 1.39 \end{bmatrix}.$$

Summing over the columns of $F$, we get the expected total survival time for each initial state (except, of course, those with $n=0$, for which the survival time is $0$)

$$E = \mathbf 1^T \; F \approx \begin{bmatrix} 25.21 & 29.25 & 34.31 & 34.31 & 37.15 & 37.65 & 38.65 & 39.65 \end{bmatrix},$$

which, multiplied by the initial state distribution vector $P_{2n+k-1} = {N \choose n} \frac 1 {2^{N+1}}$, gives a mean expected survival time of $\frac{7553}{240} \approx 31.47$ steps.

For $N = 16$, the corresponding calculation gives an expected survival time of about $7083$ steps. I can also give the expected survival times for other values of $N$ if you want.

Update: Here are the expected survival times for various $N$ and initial choices of $k$:

$$\begin{array}{r||c|c|c|} N & k_0 = 0 & k_0\text{ random} & k_0 = 1 \\ \hline 1 & 0.500000000 & 0.750000000 & 1.000000000 \\ 2 & 4.375000000 & 4.500000000 & 4.625000000 \\ 3 & 13.74927043 & 14.02225195 & 14.29523346 \\ 4 & 30.87208333 & 31.47083333 & 32.06958333 \\ 5 & 59.12594334 & 60.12456357 & 61.12318380 \\ 6 & 102.9786630 & 104.3672328 & 105.7558025 \\ 7 & 169.3232530 & 171.0531510 & 172.7830490 \\ 8 & 268.6316739 & 270.6468948 & 272.6621157 \\ 9 & 416.5658801 & 418.8161295 & 421.0663789 \\ 10 & 636.3579461 & 638.8020774 & 641.2462088 \\ 11 & 962.3500552 & 964.9556906 & 967.5613261 \\ 12 & 1445.231835 & 1447.973893 & 1450.715950 \\ 13 & 2159.757981 & 2162.617119 & 2165.476257 \\ 14 & 3216.099739 & 3219.061030 & 3222.022322 \\ 15 & 4776.525496 & 4779.577369 & 4782.629242 \\ 16 & 7079.897276 & 7083.030703 & 7086.164130 \\ 17 & 10477.62774 & 10480.83562 & 10484.04350 \\ 18 & 15486.43610 & 15489.71279 & 15492.98948 \\ 19 & 22865.71932 & 22869.06029 & 22872.40127 \\ 20 & 33730.97949 & 33734.38105 & 33737.78262 \\ \end{array}$$

Note that the middle column ($k_0$ random) is simply the average of the $k_0 = 0$ and $k_0 = 1$ columns; I've included it for completeness.

Ps. I ran the calculations up to $N = 40$, and found that for large $N$, $\exp(2.875+0.3804 N)$ gives a pretty good fit for the middle column. (Specifically, it's approximately the logarithmic least squares fit for $N = 35, \ldots, 40$.)

Addendum: For large $N$, there are approximate numerical methods that can be use to avoid inverting $I-M$. Generally, they're based on the fact that you don't really need $F = (I-M)^{-1}$, but just the vector $FP = (I-M)^{-1}P = P + MP + M^2P + \cdots$ (and that, further, you really only care about the sum of the elements of $FP$, not of $FP$ itself).

For example, you can approximate the series $FP = \sum_{m=0}^\infty M^mP$ by starting with the initial state distribution vector $P$ and repeatedly multiplying it with $M$. Of course, since the mean time to absorption is exponential in $N$, and since you have to sum several times that many terms before they become negligible, that naive approach isn't really practical.

However, if you start to compute that sum, you'll eventually notice that, for large enough $m$, $M^{m+1}P \approx \lambda M^mP$, where $\lambda$ is a constant slightly less than $1$. In fact, $\lambda$ will be the dominant eigenvalue of $M$, and the terms $M^mP$ will approach the corresponding eigenvector $\tilde P$. Thus, there will be some $m_0$ such that, for $m > m_0$, the remaining terms can be well approximated as $M^mP \approx \lambda^{m-m_0} M^{m_0}P$, and thus $\sum_{m=m_0}^\infty M^mP \approx (1-\lambda)^{-1} M^{m_0}P$.

Of course, that trick only really helps if the dominant eigenvalue is significantly closer to 1 than the others, but that's quite typical for chains like these, where the absorbing state ($n = 0$) is only reachable through a small number of relatively rare states ($n = 1$). For examples, I calculated the eigenvalues explicitly for $N=32$: the dominant eigenvalue is about 0.9999997076, while the second-largest eigenvalue is only 0.9692222154. Thus, $M^mP$ will converge to the quasi-stationary distribution $\tilde P$ about 10000 times faster than it will decay towards zero.

-
Thanks. Waiting for the final answer. – diamond kukur Jul 24 '11 at 13:00
Or $\frac{191915616959698822140954370121074028938435014457429042003}{27095127072672‌​200907704180806809457040439579268546560}$ steps, to be precise ;-) – joriki Jul 24 '11 at 16:31
Ilmari: Will you kindly point out me more preciously about the order and entries of $M_n$ and $M_k$? Thanks for your effort. – diamond kukur Jul 24 '11 at 17:45
The entries of $M_n$ and $M_k$ are just the transition probabilities given above. In particular, the diagonal entries of $M_k$ equal $1-\frac{n}{N}$ (where $n = 1+\lfloor\frac{j-1}{2}\rfloor$ and $j$ is the column index), while the off-diagonal entries (one per column, below the diagonal if $j$ is odd, above if $j$ is even) equal $\frac{n}{N}$. The off-diagonal elements of $M_n$, two rows above and below the diagonal, equal $\frac nN(1-\frac nN)$ for odd columns ($k=0$) and $(\frac nN)^2$ and $(1-\frac nN)^2$ respectively for even columns ($k=1$). – Ilmari Karonen Jul 24 '11 at 19:01
...and the diagonal elements of $M_n$ equal $1-2 \frac nN (1-\frac nN)$ for odd columns and $1-(\frac nN)^2-(1-\frac nN)^2$ for even columns. In particular, note that $M_k$ is a (right) stochastic matrix, i.e. its columns sum to one, while $M_n$ (and thus also $M$) is stochastic except for its first two columns, whose sum is one minus the probability of transitioning to an $n=0$ state. – Ilmari Karonen Jul 24 '11 at 19:11