Computing Functions of Symmetric Hierarchically Semiseparable Matrices

The aim of this work is to develop a fast algorithm for approximating the matrix function f(A) of a square matrix A that is symmetric and has hierarchically semiseparable (HSS) structure. Appearing in a wide variety of applications, often in the context of discretized (fractional) differential and integral operators, HSS matrices have a number of attractive properties facilitating the development of fast algorithms. In this work, we use an unconventional telescopic decomposition of A, inspired by recent work of Levitt and Martinsson on approximating an HSS matrix from matrixvector products with a few random vectors. This telescopic decomposition allows us to approximate f(A) by recursively performing low-rank updates with rational Krylov subspaces while keeping the size of the matrices involved in the rational Krylov subspaces small. In particular, no large-scale linear system needs to be solved, which yields favorable complexity estimates and reduced execution times compared to existing methods, including an existing divide-and-conquer strategy. The advantages of our newly proposed algorithms are demonstrated for a number of examples from the literature, featuring the exponential, the inverse square root, and the sign function of a matrix. For the special case of matrix inversion, our algorithm reduces to a procedure previously proposed by Gillman, Young, and Martinsson [Front. Math. China, 7 (2012), pp. 217-247].