  We propose two approaches, based on Riemannian optimization, for computing a
stochastic approximation of the $p$th root of a stochastic matrix $A$. In the
first approach, the approximation is found in the Riemannian manifold of
positive stochastic matrices. In the second approach, we introduce the
Riemannian manifold of positive stochastic matrices sharing with $A$ the Perron
eigenvector and we compute the approximation of the $p$th root of $A$ in such a
manifold. This way, differently from the available methods based on constrained
optimization, $A$ and its $p$th root approximation share the Perron
eigenvector. Such a property is relevant, from a modelling point of view, in
the embedding problem for Markov chains. The extended numerical experimentation
shows that, in the first approach, the Riemannian optimization methods are
generally faster and more accurate than the available methods based on
constrained optimization. In the second approach, even though the stochastic
approximation of the $p$th root is found in a smaller set, the approximation is
generally more accurate than the one obtained by standard constrained
optimization.
