Limit of power of Markov matrix

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.

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

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

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

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 $