Logo
ScienceStack
Table of Contents
Fetching ...
Sign In

Computing Functions of Symmetric Hierarchically Semiseparable Matrices

Angelo A. Casulli, Daniel Kressner, Leonardo Robol

TL;DR

This work develops a fast, scalable framework for computing functions of symmetric HSS matrices by leveraging a generalized telescopic decomposition that exploits nested low-rank structure. The method reduces large-scale rational Krylov computations to operations on small matrices, yielding nearly linear or polylogarithmic complexity depending on the function, with rigorous error bounds tied to optimal rational approximations. A key advance is the integration of telescopic decompositions with Beckeermann-style low-rank updates to produce accurate approximations of $f(A)$ from a few random vectors and small subproblems. Numerical experiments across inverses, exponentials, inverse square roots, and the sign function demonstrate substantial speedups over existing HSS-based divide-and-conquer and other Krylov-based approaches, including for matrix inversion. The framework also provides practical tools for converting between HSS representations and telescopic forms, and points to future extensions to nonsymmetric matrices and deeper theoretical understanding of the sign function result.

Abstract

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.

Computing Functions of Symmetric Hierarchically Semiseparable Matrices

TL;DR

This work develops a fast, scalable framework for computing functions of symmetric HSS matrices by leveraging a generalized telescopic decomposition that exploits nested low-rank structure. The method reduces large-scale rational Krylov computations to operations on small matrices, yielding nearly linear or polylogarithmic complexity depending on the function, with rigorous error bounds tied to optimal rational approximations. A key advance is the integration of telescopic decompositions with Beckeermann-style low-rank updates to produce accurate approximations of from a few random vectors and small subproblems. Numerical experiments across inverses, exponentials, inverse square roots, and the sign function demonstrate substantial speedups over existing HSS-based divide-and-conquer and other Krylov-based approaches, including for matrix inversion. The framework also provides practical tools for converting between HSS representations and telescopic forms, and points to future extensions to nonsymmetric matrices and deeper theoretical understanding of the sign function result.

Abstract

The aim of this work is to develop a fast algorithm for approximating the matrix function of a square matrix 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 , 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 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.
Paper Structure (19 sections, 7 theorems, 64 equations, 7 tables, 4 algorithms)

This paper contains 19 sections, 7 theorems, 64 equations, 7 tables, 4 algorithms.

Table of Contents

  1. Introduction
  2. Rational Krylov subspaces
  3. Hierarchically semiseparable matrices
  4. Block diagonal matrices and depth reduction
  5. Telescopic decompositions
  6. A general telescopic decomposition
  7. From HSS matrices to telescopic decompositions
  8. Converting a general telescopic decomposition into a standard one
  9. Symmetric telescopic decompositions
  10. Computing telescopic decompositions for functions of symmetric HSS matrices
  11. Pole selection
  12. Exponential function
  13. Markov functions
  14. Numerical experiments
  15. Computation of the inverse
  16. ...and 4 more sections

Key Result

Proposition 2.2

Let $\boldsymbol{A}:=\mathrm{blkdiag}(A_i)$ and $\boldsymbol{B}:=\mathrm{blkdiag}(B_i)$ with $A_i\in \mathbb{R}^{n_i\times n_i}$ and $B_i\in \mathbb{R}^{n_i\times b_i}$ and let $\boldsymbol{\xi}_k$ be a list of poles. Then an orthonormal basis of $\mathcal{Q}(\boldsymbol{A},\boldsymbol{B},\boldsymbo

Theorems & Definitions (21)

  • Definition 2.1
  • Proposition 2.2
  • proof
  • Definition 3.1
  • Remark 3.2
  • Remark 3.3
  • Definition 3.4
  • Definition 4.1
  • Proposition 4.2
  • proof
  • ...and 11 more