A Schur–Padé Algorithm for Fractional Powers of a Matrix
Nicholas John Higham, Lijing Lin · SIAM Journal on Matrix Analysis and Applications · 2011
A new algorithm is developed for computing arbitrary real powers [Formula: see text] of a matrix [Formula: see text]. The algorithm starts with a Schur decomposition, takes [Formula: see text] square roots of the triangular factor [Formula: see text], evaluates an [[Formula: see text]] Padé approximant of [Formula: see text] at [Formula: see text], and squares the result [Formula: see text] times. The parameters [Formula: see text] and [Formula: see text] are chosen to minimize the cost subject to achieving double precision accuracy in the evaluation of the Padé approximant, making use of a result that bounds the error in the matrix Padé approximant by the error in the scalar Padé approximant with argument the norm of the matrix. The Padé 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 [Formula: see text], yielding increased accuracy. Since the basic algorithm is designed for [Formula: see text], a criterion for reducing an arbitrary real [Formula: see text] to this range is developed, making use of bounds for the condition number of the [Formula: see text] problem. How best to compute [Formula: see text] for a negative integer [Formula: see text] 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 [Formula: see text].