Knowing that there is a probability matrix M (where all columns add to 1) which when applied to a given vector P produces the same vector P, what is the best solution to find M? I can get my head around it for a 2x2 example but cannot work out a general solution for larger matrices.
My simple example is where we know our steady state vector P:
\begin{bmatrix}1/3\\2/3\end{bmatrix}
So with eigen value = 1 Looking for matrix M where M.P = P
Therefore (M - I)P = 0
\begin{bmatrix}a-1&b\\c&d-1\end{bmatrix}
This gives equations:
(a-1) 1/3 + b 2/3 = 0
a = 0.5 b
c 1/3 + (d-1) 2/3 = 0
c = 2 - 2d
and because it is a probability matrix we know:
c = 1-a and d = 1-b
We can substitute to find the values of M:
\begin{bmatrix}0.2&0.4\\0.8&0.6\end{bmatrix}
But my question is how do I find the solution for bigger vectors such as the following where M would be 4x4 or even bigger where M is 10x10 etc. What is the fastest way to compute the solution?
\begin{bmatrix}0.1\\0.2\\0.4\\0.3\end{bmatrix}
The first equation, $\frac13 (a-1) + \frac23 b = 0,$ is actually equivalent to $a = 1 - 2b.$ It does not necessarily imply that $a = \frac12b.$ It just happens that by setting $a=\frac15$ and $b = \frac25,$ both the actual condition implied by $(M-I)P$ and the additional condition $a = \frac12b$ are satisfied. But that is not the only solution.
For example, try $$ M = M_1 = \begin{pmatrix} \frac13 & \frac13 \\ \frac23 & \frac23 \end{pmatrix}. $$
This satisfies both $\frac13 (a-1) + \frac23 b = 0$ and $\frac13 c + \frac23 (d-1) = 0,$ but most importantly, its columns each sum to $1$ and it satisfies $M_1P = P.$ In fact, it satisfies $M_1v = P$ for any vector $v$ whose entries have the sum $1.$
A trivial example, of course, is $M=I$, since $IP = P.$
Moreover, if $M_1P = P$ and $M_2P=P$ then $(qM_1)P = qP$ and $(rM_2P)=P,$ so $(qM_1 + rM_2)P = (q+r)P.$ If also $q+r=1$ then $(qM_1 + rM_2)P = P,$ that is, the linear combination of matrices $qM_1 + rM_2$ satisfies the conditions you set for $M.$ In other words, any weighted average of two matrices that satisfy the conditions for $M$ (that is, each matrix is multiplied by a scalar, the sum of the scalars is $1,$ and after scalar multiplication the resulting matrices are added together) will also satisfy the requirements of $M.$ Since you can choose $q$ to be anything you want and then set $r = 1 - q$ to satisfy $q+r=1,$ this gives you infinitely many choices as long as you can find two different solutions for $M$ (which we have done).
For example, the matrix you found in the question is $\frac65 M_1 - \frac15 I.$ To derive this linear combination of matrices you simply set $q M_1 + r I$ equal to the matrix in the question and see if you can solve for $q$ and $r$. As it turns out, you can, and the result is $q=\frac65,$ $r=-\frac15.$ You can check the result by writing out the results of the scalar multiplications and matrix addition in $\frac65 M_1 - \frac15 I$ and comparing the final result with your original matrix.
As we can see from these examples, the choice of a matrix $M$ is not uniquely determined by the requirements that $M$ be a probability matrix and that $MP = P.$
A general approach to finding an $n\times n$ matrix $M$ for a vector $P$ of $n$ entries is to identify a non-zero entry in $P$ (there must be at least one of these). Suppose $p_k \neq 0.$ Then construct a set of $n(n-1)$ matrices of dimension $n\times n$ as follows:
For each $i$ and $j$ such that $1\leq i \leq n,$ $1\leq j \leq n,$ and $i\neq k,$ construct a matrix $B(j,i)$ with entries $b_{ji} = 1$ and $b_{jk} = \frac{1}{p_k}(m_{ji} p_i),$ and set all other entries of this matrix to zero.
Each such matrix $B(j,i)$ satisfies the condition $B(j,i)\,P = 0.$ These matrices are all independent of each other and span a vector space $\mathcal L$ of $n(n-1)$ dimensions. If $L$ is any matrix in that vector space, then $LP = 0.$
We can construct an affine space of matrices $\mathcal A$ by adding the identity matrix $I$ to each matrix $L$ in $\mathcal L.$ If $A$ is a matrix in $\mathcal A$ then for some $L \in \mathcal L,$ $$ AP = (I + L)P = P + LP = P + 0 = P,$$ and therefore $M=A$ is a solution of $MP=P.$
Another way to put it is, given that $p_k\neq 0,$ you can construct an $n\times n$ matrix $A$ by putting any real numbers you like everywhere except in the $k$th column and then setting $a_{j,k}$ so that $$ a_{j,k} = \frac{1}{p_k} (p_j - a_{j,1}p_1 - \cdots - a_{j,k-1}p_{k-1} - a_{j,k+1}p_{k+1} - \cdots - a_{j,n}p_n). $$ The independent choices you make in each of the $n(n-1)$ entries other than the $k$th column give you a space of $n(n-1)$ dimensions. As long as $P\neq 0,$ however, $M=0$ is not a solution to $MP = P$ and therefore the space is an affine space rather than a vector space.
But this only shows that a matrix in this affine space is a solution to the equation $MP=P.$ Not all matrices in the space are probability matrices. The requirement that $M$ must be a probability matrix imposes more constraints.
To help discuss the implications of these constraints, suppose that a particular probability matrix $A$ can be written $$ A = I + r_{11} B(1,1) + \cdots + r_{1n} B(1,n) + \cdots + r_{n1} B(n,1) + \cdots + r_{nn} B(n,n), $$ that is, as the sum of $I$ and a linear combination of the matrices $B(j,i)$ previously described where $j$ runs from $1$ to $n$ and $i$ runs from $1$ to $n$ excluding $k.$
One set of constraints says that every entry of $A$ must be in $[0,1].$ For each $i \neq k,$ we have $a_{ii} = 1 + r_{ii},$ so the constraint implies that $-1 \leq r_{ii} \leq 0.$ But for $i\neq k$ and $i\neq j$ we have $a_{ji} = r_{ji},$ so the constraint implies that $0 \leq r_{ji} \leq 1.$
A further constraint is that the sum of each column must be $1.$ But the matrix $I$ has entry $1$ in each column, so the sum of entries in column $i$ of $A$ (where $i\neq k$) is $$ 1 = 1 + r_{1i} + \cdots + r_{ni}, $$ so the sum of all the coefficients $r_{ji}$ must be zero. This implies that the contributions of the matrices $B(j,i)$ to the sum of entries in column $k$ is also zero, so we find that the $k$th column also sums to $1$ without imposing any further constraint. However, the constraints on the sum of $r_{ji}$ for every $i$ such that $i\neq k$ reduces the dimension of the solution by $n - 1.$ The remaining solution is a convex subset (bounded by the conditions $-1 \leq r_{ii} \leq 0$ and $0 \leq r_{ji} \leq 1$ for $i\neq j$) within an $(n-1)^2$-dimensional affine space.
When $n=2$ the solution is a convex subset of a space of dimension $(2-1)^2 = 1.$ So it is no coincidence after all that the matrix found in the question is a linear combination of $I$ and the matrix $M_1$ found earlier in this answer. All solutions to that particular problem have that property.
For the vector $P = \begin{pmatrix} 1/3 \\ 2/3 \end{pmatrix},$ the solution space of $M$ consists of all matrices that can be written $$ \begin{pmatrix} 1&0 \\ 0&1 \end{pmatrix} + r \begin{pmatrix} -2 & 1 \\ 0&0 \end{pmatrix} - r \begin{pmatrix} 0&0 \\ -2 & 1 \end{pmatrix} $$ where $0 \leq r \leq 1.$
In particular, $$ \begin{pmatrix} 0.2 & 0.4 \\ 0.8 & 0.6 \end{pmatrix} = \begin{pmatrix} 1&0 \\ 0&1 \end{pmatrix} + \frac25 \begin{pmatrix} -2 & 1 \\ 0&0 \end{pmatrix} - \frac25 \begin{pmatrix} 0&0 \\ -2 & 1 \end{pmatrix} $$ and $$ \begin{pmatrix} \frac13 & \frac13 \\ \frac23 & \frac23 \end{pmatrix} = \begin{pmatrix} 1&0 \\ 0&1 \end{pmatrix} + \frac13 \begin{pmatrix} -2 & 1 \\ 0&0 \end{pmatrix} - \frac13 \begin{pmatrix} 0&0 \\ -2 & 1 \end{pmatrix}. $$