Limit of power of Markov matrix Recently, a friend gave me the following problem:

Let $$M = \begin{pmatrix} p & q & r \\ q & r & p \\ r & p & q \end{pmatrix}$$ where $p, q, r > 0$ and $p + q + r = 1$. Prove that $$\lim_{n \rightarrow \infty} M^n = \begin{pmatrix} 1/3 & 1/3 & 1/3 \\ 1/3 & 1/3 & 1/3 \\ 1/3 & 1/3 & 1/3 \end{pmatrix}$$

Some of my observations include that $M^{2^n}$ remains cyclic, so I tried to bound (with no luck) the maximum and minimum to see if they both converged to $1/3$. 
It would be nice if someone could post an elementary solution to this problem.  
 A: Let
$$ \sigma = p^2 + q^2 + r^2 - qr - rp - pq = \frac{1}{2} \left( (q-r)^2 + (r-p)^2 + (p-q)^2 \right) \geq 0.  $$
Note $ \sigma = (p+q+r) \sigma = p^3 + q^3 + r^3 - 3 pqr.$
Which is why the characteristic polynomial of your matrix is
$$ \lambda^3 - \lambda^2 - \sigma \lambda + \sigma = (\lambda - 1)(\lambda^2 - \sigma) $$
On the other hand
$$ 1 = p^2 + q^2 + r^2 +2 qr +2 rp +2 pq, $$ so
$$ \sigma = 1 - 3 (qr + rp + pq) < 1. $$
The eigenvalues are $1, \sqrt \sigma, -\sqrt \sigma.$ The diagonalized matrix $ D = P^T MP$ (where $P$ is orthogonal) raised to an arbitrarily large power, converges to the diagonal matrix with entries $1,0,0.$ Reverse the diagonalization, we get a symmetric rank one matrix.. give me a minute
Where was I?
$$
P =
\left(
\begin{array}{rr}
\frac{1}{\sqrt 3} & a & b \\
\frac{1}{\sqrt 3} & c & d \\
\frac{1}{\sqrt 3} & e & f 
\end{array}
\right)
$$
$$
D_\infty =
\left(
\begin{array}{rr}
1 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 
\end{array}
\right)
$$
We then just multiply out
$$
M_\infty = P D_\infty P^T =
\left(
\begin{array}{rr}
\frac{1}{ 3} & \frac{1}{ 3} & \frac{1}{ 3} \\
\frac{1}{ 3} & \frac{1}{ 3} & \frac{1}{ 3} \\
\frac{1}{ 3} & \frac{1}{ 3} & \frac{1}{ 3} 
\end{array}
\right)
$$
A: firstly we note that all even powers of $M$ lie within a two-dimensional subalgebra whose members we denote $A(x,y)$ with $x,y \in \mathbb{R}$, since
$$
M^2 = \begin{pmatrix} x_2 & y_2 & y_2 \\ y_2 & x_2 & y_2 \\ y_2 & y_2 & x_2 \end{pmatrix} = A(x_2,y_2)
$$
where
$$
x_2=p^2+q^2+r^2
$$
and
$$
y_2 = pq+qr+rp
$$
with
$$
x_2+2y_2 = (p+q+r)^2 =1
$$
for subsequent squarings (of $M$) we will have:
$$
x_{2^{n+1}} = x_{2^n}^2+2y_{2^n}^2 \\
y_{2^{n+1}} = 2x_{2^n}y_n+y_{2^n}^2
$$
from which we have:
$$
x_{2^{n+1}}+ 2y_{2^{n+1}} = (x_{2^n}+ 2y_{2^n})^2 = 1\tag{1}
$$
and
$$
x_{2^{n+1}}-y_{2^{n+1}} = (x_{2^n}-y_{2^n})^2  = \lambda^{2^n} \tag{2}
$$
where $0 \le \lambda = x_2-y_2 \lt 1$
it is easily verified that if
$$
A(X,Y) = A(x,y)A(x',y')
$$
then 
$$
X-Y = (x-y)(x'-y')
$$
any even number $2N$ may be written
$$
2N = \sum_{k=1} a_k2^k
$$
with $a_k \in \{0,1\}$ so
$$
M^{2N} = \prod_{a_k=1} M^{2k}
$$
so, using (2)
$$
x_{2N}-y_{2N} = \lambda^N
$$
this shows that the subsequence $M^2,M^4,M^6,\dots$ converges to $A(\frac13,\frac13)$#
to deal with the odd powers we note that
$$
M^{2N+1} = M(M^{2N}) = M\bigg(A(\frac13,\frac13) + \frac13\lambda^N A(2,-1)\bigg) = A(\frac13,\frac13) +O(\lambda^N)
$$
A: Here is a slightly streamlined approach. Let $E$ the matrix with all entries equal to $1$ and let $A=M^2$. To prove that $M^n$ converges to $\frac13E$, it suffices to show that $A^n$ converges to $\frac13E$, because such convergence implies also that $M^{2n+1}=A^nM$ converges to $\frac13EM=\frac13E$.
As another answer has pointed out, $A=M^2$ is in the form of
$$
A=M^2=\pmatrix{x&y&y\\ y&x&y\\ y&y&x},
$$
where $y=pq+qr+rp$. Since all row sums of $M$ are equal to $1$, row sums of $A$ must be equal to $1$ too, i.e. $x+2y=1$. Furthermore, as $p+q+r=1$, if we substitute the equality $p^2+q^2+r^2=1-2(pq+qr+rp)$ into the inequality $(p-q)^2+(q-r)^2+(r-p)^2\ge0$, we get $y=pq+qr+pr\le\frac13$. Thus
$$
A=(1-3y)I+yE,
$$
where $0\le 1-3y<1$.
Let $U$ be any real orthogonal matrix whose first column is the unit vector $\frac1{\sqrt{3}}(1,1,1)^T$. Then $E=3U\operatorname{diag}(1,0,0)U^T$. Consequently,
$$
\begin{align}
\lim_{n\to\infty}A^n
&=\lim_{n\to\infty}U\,[\,(1-3y)I+3y\operatorname{diag}(1,0,0)\,]^n\,U^T\\
&=\lim_{n\to\infty}U\operatorname{diag}(1,\,1-3y,\,1-3y)^nU^T\\
&=U\operatorname{diag}(1,0,0)U^T=\frac13E.
\end{align}
$$
A: Except for the last line, this is a purely algebraic answer, as there is a lot of special structure in this problem.  (Note Gerschgorin discs are an elementary and powerful but underutilized tool to bound eigenvalues.  They also have a nice picture associated with them.)  
you have a real valued symmetric transition matrix.  Hence it is diagonalizable with real eigenvalues, and mutually orthonormal eigenvectors.  
$M = U \Lambda U^T$
It's rows and columns both sum to one i.e. it is double stochastic.  Hence we know $\lambda_1 = 1$ and $\mathbf u_1 \propto \mathbf 1$. (In fact $\mathbf u_1 = \frac{1}{\sqrt{3}}\mathbf 1$)
$trace(M) = p + q + r = 1 = \lambda_1 + \lambda_2 + \lambda_3 = 1 +\lambda_2 + \lambda_3$ 
Hence we know $\lambda_2 + \lambda_3 = 0$. Put differently: $\lambda_2 = - \lambda_3$.
Use Gerschgorin discs to find that the minimum eigenvalue $\gt -1$.  
I.e. each disc has a center point $\gt 0$ and a radius $\lt 1$.  (It's really more of a line segment than a disc since we know the eigenvalues are real, but we can still draw the disc as an upper bound.)
Hence $\vert\lambda_2\vert \lt 1$ and $\vert\lambda_3\vert \lt 1$
$M^n = U \Lambda^n U^T = \big(\lambda_1^n \mathbf u_1 \mathbf u_1^T + \lambda_2^n \mathbf u_2 \mathbf u_2^T + \lambda_3^n \mathbf u_3 \mathbf u_3^T \big) = \big(1^n \mathbf u_1 \mathbf u_1^T + \lambda_2^n \mathbf u_2 \mathbf u_2^T + \lambda_3^n \mathbf u_3 \mathbf u_3^T \big)$
now at the end, take a limit:
$\lim_{n \to \infty} M^n = \lim_{n \to \infty}\big(1^n \mathbf u_1 \mathbf u_1^T + \lambda_2^n \mathbf u_2 \mathbf u_2^T + \lambda_3^n \mathbf u_3 \mathbf u_3^T\big) = 1 \mathbf u_1 \mathbf u_1^T + 0 \mathbf u_2 \mathbf u_2^T + 0 \mathbf u_3 \mathbf u_3^T = \mathbf u_1 \mathbf u_1^T $
A: This a problem which is probably intended to be solved using Markov chains. I will assume you have a basic knowledge of the subject.

Suppose $(X_n)_{n\in\mathbb{N}}$ is a homogenous markov chain with transition matrix given by 
$$P=\begin{pmatrix}
p & q & r\\ 
q & r & p\\ 
r & p & q
\end{pmatrix}$$
where $p+q+r=1$ and $p,q,r>0$. Then, clearly, the chain is irreducible and ergodic. 

Ergodic theorem: If $(X_n)_{n\in\mathbb{N}}$ is an irreducible and ergodic Markov chain, then, it admits a unique stationnary distribution which satisfies
  $$\pi_j=\lim_{n\to \infty} p_{i,j}^{(n)}\qquad \forall i,j.$$

Solving the equation $\pi P=\pi$ for your chain gives 
$$\pi=(1/3,1/3,1/3)$$
and thus, by the cited theorem,
$$\lim_{n\to \infty} P^{(n)}=\begin{pmatrix}
\pi_1 & \pi_2 & \pi_3\\ 
\pi_1 & \pi_2 & \pi_3\\ 
\pi_1 & \pi_2 & \pi_3
\end{pmatrix}=\begin{pmatrix}
1/3 & 1/3 & 1/3\\ 
1/3 & 1/3 & 1/3\\ 
1/3 & 1/3 & 1/3
\end{pmatrix}$$
which is the result we were looking for.
