Splitting methods based on the nonzero diagonal pattern for computing matrix functions

Abstract

We consider the task of approximating a matrix function $f(A)$, where $A$ is a matrix in which only a relatively small number of (not necessarily consecutive) sub- and superdiagonals contain nonzero entries. Approximating $f$ by a low-degree polynomial $p$ allows us to obtain sparse approximations to $f(A)$, which one can efficiently work with (while, in general, $f(A)$ is a dense matrix, even when $A$ is sparse). Our approach is based on carefully inspecting the locations where nonzeros can occur in $p(A)$, and identifying the entries in $A$ that influence them. In particular, we illustrate how this approach can be used for efficiently approximating the trace of $f(A)$ and identify how this approach is related to established (stochastic) probing methods for trace estimation. Another application area in which our approach works particularly well is the computation of functions of Toeplitz matrices. Here, studying the sparsity pattern of $p(A)$ allows us to reduce the computation of the whole matrix polynomial to that of a single small-scale submatrix, yielding an algorithm that scales exceptionally well to large problem sizes.