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.)
 A: The infinite product converges to $\frac1nJ_n$ with probability $1$.
Let $\mathcal S=\{S_1,S_2,\ldots,S_{n-1}\}$ be your set of doubly stochastic matrices (I choose the symbol $S$ to mean "stochastic"), where $S_k=I_{k-1}\oplus(\frac12J_2)\oplus I_{n-k-1}$. The sample space is $\mathcal S^{\mathbb N}$, the set of all infinite sequences $\mathbf X=(X_1,X_2,\ldots)$ such that $X_i\in\mathcal S$ for each $i$. That is, $X_i$ is the outcome of the $i$-th draw from $\mathcal S$. The measurable events are defined by the smallest $\sigma$-algebra that contains $\{\mathbf X\in\mathcal S^{\mathbb N}: X_i=S\}$ as a member for each $i\in\mathbb N$ and each $S\in\mathcal S$. We assume that the draws are i.i.d. uniform.
Fix any $S\in\mathcal S$. Let $E_{rm}$ denotes the event that $S$ occurs at most $r$ times in the first $m$ draws. Then
$$
P\left(\bigcap_{m\in\mathbb N}E_{rm}\right)=\lim_{m\to\infty}P(E_{rm})=0
\ \text{ and hence }
\ P\left(\bigcup_{r=0}^\infty\bigcap_{m\in\mathbb N}E_{rm}\right)\le\sum_{r=0}^\infty P\left(\bigcap_{m\in\mathbb N}E_{rm}\right)=0,
$$
i.e. the probability that $S$ occurs only finitely many times is zero. It follows from the finiteness of $\mathcal S$ that with probability $1$, all matrices in $\mathcal S$ occur infinitely many times in the sequence $\mathbf X$.
Now, given any fixed sequence $\mathbf X\in\mathcal S^{\mathbb N}$ in which every matrix in $\mathcal S$ shows up infinitely many times, we want to show that the infinite product $\prod_{i=1}^\infty X_i$ converges to $\frac1nJ_n$. Let
$$
A=\pmatrix{1&\frac12&\frac13&\cdots&\frac1n\\ &\frac12&\frac13&\cdots&\frac1n\\ &&\frac13&\cdots&\frac1n\\ &&&\ddots&\vdots\\ &&&&\frac1n}.
$$
Note that $J_nA=J_n$. Therefore $\frac1nJ_n=\frac1nJ_nA^{-1}$. So, if we can prove that
$$
\lim_{k\to\infty}\prod_{i=1}^kX_iA=\frac1nJ_n,\tag{1}
$$
then
$$
\lim_{k\to\infty}\prod_{i=1}^kX_i=\left(\lim_{k\to\infty}\prod_{i=1}^kX_iA\right)A^{-1}=\frac1nJ_nA^{-1}=\frac1nJ_n.
$$
For convenience, let us call a vector $v=(v_1,v_2,\ldots,v_n)^T\in\mathbb R^n$ a down vector if it is a probability vector whose entries are arranged in decreasing order, i.e. if $v_1\ge v_2\ge\cdots\ge0$ and $\sum_{i=1}^nv_i=1$. Since every column of $A$ is a down vector, in order to prove $(1)$, it suffices to prove that $\prod_{i=1}^kX_iv$ converges to $\frac1ne=\frac1n(1,1,\ldots,1)^T$ for every down vector $v$.
Let $v$ be a down vector. For any $S\in\mathcal S$, $u=Sv$ is also a down vector and $v_1\ge u_1$. So, if we iteratively define $v^{(0)}=v$ and $v^{(k)}=X_kv^{(k-1)}=\prod_{i=1}^kX_iv$, then $\{v_1^{(k)}\}_{k\in\mathbb N}$ is monotonic decreasing and converges to some limit $c$.
We now prove by mathematical induction that $\lim_{k\to\infty}v_i^{(k)}=c$ for every $i$. The base case has been settled in the above. In the induction step, suppose that $\lim_{k\to\infty}v_i^{(k)}=c$ for some $i$. Then for any $\epsilon>0$, there exists an integer $K$ such that $c+\epsilon\ge v_i^{(k)}\ge c-\epsilon$ for all $k\ge K$. Since $S_i$ shows up infinitely many times, we may assume that $X_K=S_i$.
Suppose that $S_i$ occurs at some time $k_0\ge K$ and suppose its next occurrence time is $k_1$. Since $S_i$ doesn't occur between these two time points, the iterates of $v_{i+1}$ must be decreasing from time $k_0$ to time $k_1-1$:
$$
v_{i+1}^{(k_0)}\ge v_{i+1}^{(k_0+1)}\ge v_{i+1}^{(k_0+2)}\ge\cdots\ge v_{i+1}^{(k_1-1)}.
$$
However, as $S_i$ occurs at time $k_i$, we also have $v_i^{(k_1)}=\frac12(v_i^{(k_1-1)}+v_{i+1}^{(k_1-1)})\ge c-\epsilon$. Therefore
$$
v_{i+1}^{(k_1-1)}\ge2(c-\epsilon)-v_i^{(k_1-1)}\ge2(c-\epsilon)-(c+\epsilon)=c-3\epsilon.
$$
It follows that $v_{i+1}^{(k_0)}\ge v_{i+1}^{(k_0+1)}\ge v_{i+1}^{(k_0+2)}\ge\cdots\ge v_{i+1}^{(k_1-1)}\ge c-3\epsilon$. Since $k_0$ can be any occurrence time of $S_i$ since time $K$, we get $v_{i+1}^{(k)}\ge c-3\epsilon$ for all $k\ge K$. In turn, $c+\epsilon\ge v_i^{(k)}\ge v_{i+1}^{(k)}\ge c-3\epsilon$ for all $k\ge K$. As $\epsilon$ is arbitrary, $\lim_{k\to\infty}v_{i+1}^{(k)}=c$. The induction step is now completed.
Thus $\lim_{k\to\infty}v_i^{(k)}=c$ for every $i$, i.e. $\lim_{k\to\infty}v^{(k)}=ce$. However, as every iterate $v^{(k)}$ is a probability vector, so is $ce$. Therefore $c=\frac1n$ and $\lim_{k\to\infty}v^{(k)}=\frac1ne$. Now we are done.
A: 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.  
