A Schur–Padé Algorithm for Fractional Powers of a Matrix

• Published in: SIAM journal on matrix analysis and applications Vol. 32; no. 3; pp. 1056 - 1078
• Main Authors: Higham, Nicholas J., Lin, Lijing
• Format: Journal Article
• Language: English
• Published: Philadelphia, PA Society for Industrial and Applied Mathematics 01.07.2011
• Subjects: Accuracy; Algebra; Algorithms; Approximation; Combinatorics; Combinatorics. Ordered structures; Decomposition; Designs and configurations; Eigenvalues; Exact sciences and technology; Experimental design; Linear and multilinear algebra, matrix theory; Mathematics; Numerical analysis; Numerical analysis. Scientific computation; Numerical linear algebra; Probability and statistics; Sciences and techniques of general use; Statistics; primary matrix function; Costs; matrix exponential; Error estimation; matrix root; Computing; MATLAB; Continued fractions; 65F30; Padé approximation; Numerical stability; Numerical algorithm; Real matrix; Evaluation; Square root; Error bound; Representation; matrix power; Matrix function; Integer; matrix logarithm; Numerical analysis; Experimental design; Linear algebra; Fractional power; Condition number; Schur decomposition; Padé approximant
• ISSN: 0895-4798; 1095-7162
• DOI: 10.1137/10081232X

A new algorithm is developed for computing arbitrary real powers Ap of a matrix A∈Cn×n. The algorithm starts with a Schur decomposition, takes k square roots of the triangular factor T, evaluates an [m/m] Padé approximant of (1-x)p at I-T1/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 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 Tp/2j, yielding increased accuracy. Since the basic algorithm is designed for p∈(-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 Ap problem. How best to compute Ak 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 Ap=exp(plog(A)). [PUBLICATION ABSTRACT]