A SCHUR-PADE ALGORITHM FOR FRACTIONAL POWERS OF A MATRIX

A new algorithm is developed for computing arbitrary real powers A(p) of a matrix A is an element of C-nxn. The algorithm starts with a Schur decomposition, takes k square roots of the triangular factor T, evaluates an [m/m] Pade approximant of (1 - x)(p) at I - T-1/2k, 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 Pade approximant, making use of a result that bounds the error in the matrix Pade approximant by the error in the scalar Pade approximant with argument the norm of the matrix. The Pade 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-p/2j, yielding increased accuracy. Since the basic algorithm is designed for p is an element of (-1, 1), a criterion for reducing an arbitrary real p to this range is developed, making use of bounds for the condition number of the A(p) 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 and approaches based on the formula A(p) = exp(p log(A)).