3
$\begingroup$

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

  • 2
    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\}$.2011-07-24

1 Answers 1

7

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.

  • 0
    @Ilmari: A very polite and humble request. Will you kindly look at http://math.stackexchange.com/questions/54722/expectation-of-an-event ?2011-07-31