Optimal Solvers for Linear Systems with Fractional Powers of Sparse SPD Matrices

In this paper we consider efficient algorithms for solving the algebraic equation ${\mathcal A}^\alpha {\bf u}={\bf f}$, $0< \alpha <1$, where ${\mathcal A}$ is a symmetric and positive definite matrix obtained form finite difference or finite element approximations of second order elliptic problems in ${\mathbb R}^d$, $d=1,2,3$. The method is based on the best uniform rational approximation of the function $t^{\beta-\alpha}$ for $0 < t \le 1$ and natural $\beta$, and the assumption that one has at hand an efficient method (e.g. multigrid, multilevel, or other fast algorithm) for solving equations like $({\mathcal A} +c {\mathcal I}){\bf u}= {\bf f}$, $c \ge 0$. The provided numerical experiments on model problems with ${\mathcal A}$ obtained by finite element approximation of elliptic equations in one and three spacial dimensions confirm the efficiency of the proposed algorithms.