What matrices preserve the $L_1$ norm for positive, unit norm vectors?
Solution 1:
The matrices that preserve the set $P$ of probability vectors are those whose columns are members of $P$. This is obvious since if $x \in P$, $M x$ is a convex combination of the columns of $M$ with coefficients given by the entries of $x$. Each column of $M$ must be in $P$ (take $x$ to be a vector with a single $1$ and all else $0$), and $P$ is a convex set.
Solution 2:
Since you originally asked about $L^1$ spaces I dared to add this comment.
If one wants to preserve the integral in (finite-dimensional and with finite measure ) $L^1$ spaces rather than the norm of $\ell^p$, the matrices $M$ that do this are more general than the stochastic matrices.
One can define these matrices with two components, labeled $S$ (for the stochastic component) and $G$ (for the generalized permutation matrix component) such that that $M= S * G$, where * represents the Hadamard product.
The $S$ matrices are effectively stochastic matrices as shown by Robert Israel.
The $G$ matrix is given by the unique matrix resulting of the outer product $u_{\mu} \otimes \frac{1}{u_{\mu}} := | u_{\mu} \rangle \langle \frac{1}{u_{\mu}} |$ of the unique column vector $u_{\mu} :=\left(\begin{array}{c}\mu_1 \\ \mu_2 \\ \ldots \\ \mu_2\end{array}\right)$ and the also unique row vector $\frac{1}{u_{\mu}} :=\left(\frac{1}{\mu_1} \ \frac{1}{\mu_2} \ \ldots \ \frac{1}{\mu_n}\right)$:
$G:=u_{\mu} \otimes \frac{1}{u_{\mu}} = \left(\begin{array}{cccc} 1 & \frac{\mu_2}{\mu_1} & \ldots & \frac{\mu_n}{\mu_1} \\ \frac{\mu_1}{\mu_2} & 1 & \ldots & \frac{\mu_n}{\mu_2} \\ \ldots & \ldots & \ldots & \ldots \\ \frac{\mu_1}{\mu_n} & \frac{\mu_2}{\mu_n} & \ldots & 1 \end{array}\right)$
where $\mu_i$ are the measures of the generating family of subsets $\{ A_i \}$ of the underlying sigma algebra, i.e. $\mu_i := \mu(A_i)$ and $n = |\{ A_i \}|$.
To give you an example of where the stochastic component $S$ is absent, take the stochastic matrix $S$ to be simply a permutation matrix. In this case your $M$ that preserves the integral is a generalized permutation matrix whose non-zero elements are of the form $A_{i,j} =\frac{\mu_{j}}{\mu_i}$.
To see why the measure values $\mu_i$ are needed in the definition of $M$ recall that a $L^p$ space is defined given a measure space $(X,\Sigma,\mu)$. So if the $L^1$ space is finite dimensional then the vectors $v$ in $L^1$ are simple functions, whose integral is defined as the product $\langle u_{\mu}|v\rangle$. And if the measure is finite, then this integral is always well defined.
I hope I made myself clear.