Computing the matrix fractional power with the double exponential formula

Abstract

Two quadrature-based algorithms for computing the matrix fractional power $A^\alpha$ are presented in this paper. These algorithms are based on the double exponential (DE) formula, which is well-known for its effectiveness in computing improper integrals as well as in treating nearly arbitrary endpoint singularities. The DE formula transforms a given integral into another integral that is suited for the trapezoidal rule; in this process, the integral interval is transformed into an infinite interval. Therefore, it is necessary to truncate the infinite interval to an appropriate finite interval. In this paper, a truncation method, which is based on a truncation error analysis specialized to the computation of $A^\alpha$, is proposed. Then, two algorithms are presented---one where $A^\alpha$ is computed with a fixed number of abscissa points and one with $A^\alpha$ computed adaptively. Subsequently, the convergence rate of the DE formula for Hermitian positive definite matrices is analyzed. The convergence rate analysis shows that the DE formula converges faster than Gaussian quadrature when $A$ is ill-conditioned and $\alpha$ is a non-unit fraction. Numerical results show that our algorithms achieve the required accuracy and are faster than other algorithms in several situations.