Randomized Strong Recursive Skeletonization: Simultaneous Compression and LU Factorization of Hierarchical Matrices using Matrix-Vector Products

TL;DR

RSRS addresses the challenge of invertible factorization for $H^2$ matrices under strong admissibility in a matrix-free setting. It builds a simultaneous compression and LU factorization using randomized sketches of $\mathbf{A}$ and $\mathbf{A}^*$, with structure imposed by post-processing (block nullification and extraction) rather than explicit sparsity patterns. The method achieves near-linear complexity in the problem size and produces an invertible factorization suitable as a direct solver or robust preconditioner, even for ill-conditioned or indefinite operators. Numerical results across 2D and 3D integral and differential problems demonstrate robustness, sample-efficiency, and scalable performance, highlighting RSRS as a practical tool for large-scale PDEs and integral equations.

Abstract

The hierarchical matrix framework partitions matrices into subblocks that are either small or of low numerical rank, enabling linear storage complexity and efficient matrix-vector multiplication. This work focuses on the $H^2$-matrix format constructed under the strong admissibility condition, which has two key properties: (1) a compressed representation that approximates far-field interactions with low-rank blocks while near-field interactions are stored densely, and (2) a nested basis structure that reuses basis matrices across hierarchy levels. Although these matrices support fast Cholesky and LU factorizations, implementing them - especially for 3D PDE discretizations - remains challenging due to the nested recursions and recompressions involved. This paper introduces an algorithm that simultaneously compresses and factorizes a general ${H}^{2}$ matrix, using only the action of the matrix and its adjoint on vectors. The number of required matrix-vector products is independent of the matrix size, and depends only on the problem geometry and a rank parameter. The resulting LU factorization is invertible and can serve as an approximate direct solver, with accuracy influenced by the spectral properties of the matrix. To achieve competitive sample complexity, the method employs dense Gaussian test matrices without explicitly encoding structured sparsity. Samples are drawn only once at the start of the algorithm; as the factorization proceeds, structure is dynamically introduced into the test matrices through efficient linear algebraic operations. Numerical experiments demonstrate robustness to indefiniteness and ill-conditioning, as well as the efficiency of the method for integral and differential equations in 2D and 3D.

Figures (8)

• Figure 1: Illustration of strong skeletonization for a single box $B$. Top row: the matrix structure at each stage (with modified entries in blue). Bottom row: corresponding pseudo-1D geometry and index ordering. "Sparsify" applies the ID to far-field interactions and identifies redundant indices $\mathsf R$. "Diagonalize" then eliminates the remaining interactions involving $\mathsf R$ via block elimination.
• Figure 2: For a pseudo-1D geometry, the active points as well as the corresponding matrix are shown at various stages of the computation. First, the redundant points of the boxes on the finest level are diagonalized. To continue the computation, the remaining active points are regrouped according to the next coarse level of the tree, which introduces low-rank subblocks that can be further diagonalized.
• Figure 3: The analog of Figure \ref{['fig:SRS1d']} when the computational domain is a unit box. Observe that many more blocks get updated for a true two-dimensional domain. For every diagonalization step, up to 35 pairwise interactions are introduced between non-neighboring boxes, which must be stored and later recompressed.
• Figure 4: Two types of test matrices are needed for RSRS. (\ref{['fig:setting_blocknull']}) To sketch far-field interactions between a box and its distant neighbors. (\ref{['fig:setting_blockextract']}) To extract near-field interactions between a subset of box points and neighboring points. Rather than designing structured test matrices directly, we apply block nullification and block extraction, which use linear transformations of dense Gaussian matrices to introduce the desired structure.
• Figure 5: Reconstruction time, memory usage, and solve time vs. degrees of freedom for matrix $\bm{\mathsf{A}}$ arising from the discretization of a boundary integral equation with entries given in (\ref{['eq:Aentries_exp1']}).

Theorems & Definitions (3)

• Remark 1: Well-conditioning of the sparsifying matrices
• Remark 2: Modified near-field interactions
• Remark 3: Nested bases