## Roots of Stochastic Matrices and Fractional Matrix Powers

In Markov chain models in finance and healthcare a transition matrix over a certain time interval is needed but only a transition matrix over a longer time interval may be available. The problem arises of determining a stochastic $p$th root of astochastic matrix (the given transition matrix). By exploiting the theory of functions of matrices, we develop results on the existence and characterization of stochastic $p$th roots. Our contributions include characterization of when a real matrix hasa real $p$th root, a classification of $p$th roots of a possibly singular matrix,a sufficient condition for a $p$th root of a stochastic matrix to have unit row sums,and the identification of two classes of stochastic matrices that have stochastic $p$th roots for all $p$. We also delineate a wide variety of possible configurationsas regards existence, nature (primary or nonprimary), and number of stochastic roots,and develop a necessary condition for existence of a stochastic root in terms of the spectrum of the given matrix. On the computational side, we emphasize finding an approximate stochastic root: perturb the principal root $A^{1/p}$ or the principal logarithm $\log(A)$ to the nearest stochastic matrix or the nearest intensity matrix, respectively, if they are not valid ones;minimize the residual $\normF{X^p-A}$ over all stochastic matrices $X$ and also over stochastic matrices that are primary functions of $A$. For the first two nearness problems, the global minimizers are found in the Frobenius norm. For the last two nonlinear programming problems, we derive explicit formulae for the gradient and Hessian of the objective function $\normF{X^p-A}^2$ and investigate Newton's method, a spectral projected gradient method (SPGM) and the sequential quadratic programming method to solve the problem as well as various matrices to start the iteration. Numerical experiments show that SPGM starting with the perturbed $A^{1/p}$to minimize $\normF{X^p-A}$ over all stochastic matrices is method of choice.Finally, a new algorithm is developed for computing arbitrary real powers $A^\a$ of a matrix $A\in\mathbb{C}^{n\times n}$. The algorithm starts with a Schur decomposition,takes $k$ square roots of the triangular factor $T$, evaluates an $[m/m]$ Pad\'e approximant of $(1-x)^\a$ at $I - T^{1/2^k}$, and squares the result $k$ times. The parameters $k$ and $m$ are chosen to minimize the cost subject to achieving double precision accuracy in the evaluation of the Pad\'e approximant, making use of a result that bounds the error in the matrix Pad\'e approximant by the error in the scalar Pad\'e approximant with argument the norm of the matrix. The Pad\'e approximant is evaluated from the continued fraction representation in bottom-up fashion, which is shown to be numerically stable. In the squaring phase the diagonal and first superdiagonal are computed from explicit formulae for $T^{\a/2^j}$, yielding increased accuracy. Since the basic algorithm is designed for $\a\in(-1,1)$, a criterion for reducing an arbitrary real $\a$ to this range is developed, making use of bounds for the condition number of the $A^\a$ problem. How best to compute $A^k$ for a negative integer $k$ is also investigated. In numerical experiments the new algorithm is found to be superior in accuracy and stability to several alternatives,including the use of an eigendecomposition, a method based on the Schur--Parlett\alg\ with our new algorithm applied to the diagonal blocks and approaches based on the formula $A^\a = \exp(\a\log(A))$.