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

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 $\mathcal{H}^2$-matrix format, whose defining feature is the nested basis property which allows basis matrices to be reused across different levels of the hierarchy. While $\mathcal{H}^2$-matrices support fast Cholesky and LU factorizations, implementing these methods is challenging -- especially for 3D PDE discretizations -- due to the complexity of nested recursions and recompressions. Moreover, compressing $\mathcal{H}^2$-matrices becomes particularly difficult when only matrix-vector multiplication operations are available. This paper introduces an algorithm that simultaneously compresses and factorizes a general $\mathcal{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, depending only on the problem geometry and a rank parameter that captures low-rank interactions between well-separated boxes. The resulting LU factorization is invertible and can serve as an approximate direct solver, with its accuracy influenced by the spectral properties of the matrix. To achieve competitive sample complexity, the method uses dense Gaussian test matrices without explicitly encoding structured sparsity patterns. 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 the algorithm's robustness to indefiniteness and ill-conditioning, as well as its efficiency in terms of sample cost for challenging problems arising from both integral and differential equations in 2D and 3D.