  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 matrix-vector
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. Even
for matrix inversion, our algorithm exhibits superior performance, even if not
specifically designed for this task.
