<s>[INST] How to evaluate an infinite product of idempotent stochastic matrices? How to find the infinite products of the following matrices?
I have a set of $n\times n$ doubly stochastic matrices, all of them are idempotent and when they act on a vector, what they do is to take the average of two consecutive elements in a vector. That is, each matrix is in the form of $I_k\oplus\pmatrix{\frac12&\frac12\\ \frac12&\frac12}\oplus I_{n-k-2}$. For example, when $n=4$, the set comprises of the following members:
$$
\begin{pmatrix}
\frac12 & \frac12 & 0 & 0 \\
\frac12 & \frac12 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1 \\
\end{pmatrix},
\ \ \begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & \frac12 & \frac12 & 0 \\
0 & \frac12 & \frac12 & 0 \\
0 & 0 & 0 & 1 \\
\end{pmatrix}
\ \text{ and }\
\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & \frac12 & \frac12 \\
0 & 0 & \frac12 & \frac12 \\
\end{pmatrix}.
$$
When they are drawn randomly to form an infinite sequence, the infinite product converges to:
$$
\frac14
\begin{pmatrix}
1&1&1&1\\ 1&1&1&1\\ 1&1&1&1\\ 1&1&1&1
\end{pmatrix} = \frac{1}{4}J_4.
$$
This is for $4$-by-$4$ matrices. By numerical experiments, I saw that for $n$-by-$n$ matrices, we get the same result, i.e. the product is: 
$$\frac{1}{n}J_n.$$
Is there any way to prove that? (Although all the matrices are diagonalizable, they are multiplied in random order infinite times, so it is not straightforward how to prove that the result will converge.)
 [/INST] it's a actually a very simple problem.  As a gut check note that your infinite product is bounded -- the product of $k$ stochastic matrices (n x n) has a Frobenius norm $\leq n$ for all natural numbers $k$.  
It's enough to prove that
$\mathbf e_j^T\cdot \text{"infinite product"} = \frac{1}{n}\mathbf 1^T$
WP1, for each standard basis vector. $\mathbf e_j$ 
There's no need to reinvent the wheel here: your problem is just a standard finite state time homogenous markov chain, in disguise.  
I.e. for your specific example (with the obvious generalization the the n x n case)  what you have is
$P =  
\frac{1}{3}\begin{pmatrix}
\frac12 & \frac12 & 0 & 0 \\
\frac12 & \frac12 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1 \\
\end{pmatrix} + \frac{1}{3}\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & \frac12 & \frac12 & 0 \\
0 & \frac12 & \frac12 & 0 \\
0 & 0 & 0 & 1 \\
\end{pmatrix}
+ \frac{1}{3}\begin{pmatrix}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & \frac12 & \frac12 \\
0 & 0 & \frac12 & \frac12 \\
\end{pmatrix} = \begin{pmatrix}
\frac{5}{6} & \frac{1}{6} & 0 & 0 \\
\frac{1}{6} & \frac{2}{3} & \frac{1}{6} & 0 \\
0 & \frac{1}{6} & \frac{2}{3} & \frac{1}{6} \\
0 & 0 & \frac{1}{6} & \frac{5}{6} \\
\end{pmatrix}
$
(I took 'they are drawn randomly' to mean uniform at random but this can be tweaked without much complication so long as each probability is $\in (0,1)$)  
In effect you just applied total probability by conditioning on some unrelated event (say a roll of dice).  To clarify: each transition matrix is applied with probability $\frac{1}{3}$ based on some event (roll of dice) that is independent of past selection and more to the point, it is independent of the current state we are in -- i.e. your selection mechanism preserves the markov property.  
$P$ is doubly stochastic, has a single communicating class and is aperiodict (e.g. because there is at least one positive number on the diagonal).  
so
$\mathbf e_j^T\cdot \prod_{k=1}^\infty P = \frac{1}{n}\mathbf 1^T$
by standard results from markov chains, renewal theory or Perron Frobenius theory  
which completes the proof.  
</s>