💻 computer science

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

This paper presents a randomized algorithm that simultaneously compresses and factorizes $H^2$ -matrices using only matrix-vector products, achieving a sample complexity independent of matrix size while providing a robust, invertible approximate direct solver for integral and differential equations in 2D and 3D.

Original authors: Anna Yesypenko, Per-Gunnar Martinsson

Published 2026-02-02

📖 6 min read🧠 Deep dive

CC BY 4.0

Original authors: Anna Yesypenko, Per-Gunnar Martinsson

Original paper licensed under CC BY 4.0 (http://creativecommons.org/licenses/by/4.0/). ✨ This is an AI-generated explanation of the paper below. It is not written or endorsed by the authors. For technical accuracy, refer to the original paper. Read full disclaimer

Imagine you have a massive, incredibly complex puzzle. In the world of mathematics and physics, this puzzle is a giant "matrix"—a grid of numbers representing a problem like how heat spreads through a metal block or how sound waves bounce off a sphere.

Usually, solving this puzzle requires looking at every single number in the grid. If the puzzle has a million pieces, looking at every single one takes forever and requires a computer with a massive memory.

This paper introduces a new, clever way to solve these puzzles called Randomized Strong Recursive Skeletonization (RSRS). Here is how it works, explained through simple analogies:

1. The Problem: The "Too Big to Hold" Puzzle

In many scientific problems, the matrix is "dense," meaning almost every number is connected to every other number.

The Old Way: To solve the puzzle, you usually need to write down every single number on a giant sheet of paper. This is slow and uses up all your memory.
The H2-Matrix Idea: Scientists realized that while the puzzle looks messy, it actually has hidden patterns. If you look at two parts of the puzzle that are far apart, they interact in a very simple, predictable way (like a low-rank pattern). You don't need to write down every number for those distant parts; you just need a few "summary notes." This is called compression.

2. The Challenge: The "Black Box"

The tricky part is that in many real-world scenarios, we don't have the "sheet of paper" with all the numbers. We only have a Black Box.

You can put a list of numbers (a vector) into the box, and it spits out a new list of numbers (the result of the matrix acting on that vector).
But you can't peek inside to see the individual numbers.
Previous methods for solving these puzzles required peeking inside or using very specific, complicated test inputs. If you couldn't see the numbers, you were stuck.

3. The Solution: The "Magic Sketch"

The authors created a method to solve the puzzle using only the Black Box, without ever seeing the individual numbers. They call this RSRS.

Here is the step-by-step magic trick:

Step A: The Random "Splat"

Instead of trying to guess the structure of the puzzle, the researchers throw a bunch of random "darts" (random numbers) at the Black Box.

Think of this like spraying a wall with a hose. You don't know the shape of the wall, but the water hits it and splashes back.
By analyzing how the water splashes back (the output), they can start to figure out the shape of the wall.
Crucially, they only need to do this a fixed number of times, regardless of how huge the puzzle is. Whether the puzzle has 1,000 pieces or 1,000,000, the number of "splats" needed stays the same.

Step B: The "Skeleton" (The Bones of the Puzzle)

Once they have the splashes, they use a technique called Skeletonization.

Imagine the puzzle is a human body. You don't need to know the exact shape of every muscle and skin cell to understand how the body moves. You just need the skeleton (the bones).
The algorithm finds the "bones" of the matrix—the most important numbers that hold everything together. It ignores the "flesh" (the less important details) because the distant parts of the puzzle are simple enough to be summarized by these bones.

Step C: The Recursive "Russian Doll"

The puzzle is organized like a set of Russian nesting dolls (a hierarchy).

Start Small: They solve the puzzle for the tiniest dolls (the smallest groups of numbers).
Build Up: They take the "bones" they found in the small dolls and use them to build the solution for the slightly larger dolls.
Repeat: They keep doing this, moving from the smallest groups to the largest, until they have solved the whole thing.

Because they are building on the work they just did, they don't have to start over every time. This makes the process incredibly fast.

Step D: The "Magic Filter" (Block Nullification)

One of the paper's biggest innovations is how they handle the "Black Box" limitation.

Normally, to isolate a specific part of the puzzle, you'd need to tell the Black Box, "Ignore these numbers, only look at these." But you can't do that if you can't see the numbers.
The authors invented a "Magic Filter." They take their random "splats" and mathematically twist them so that they act as if they are ignoring the wrong parts and focusing only on the right parts.
It's like taking a photo of a crowd and using software to blur out everyone except the person you are interested in, without ever having to ask the crowd to stand still.

4. The Result: A Fast, Accurate Solver

By combining these steps, the algorithm produces a factorization.

Think of the original puzzle as a locked safe.
The algorithm doesn't just guess the combination; it builds a master key (an approximate inverse) that can open the safe almost instantly.
This key works even if the safe is rusty or broken (ill-conditioned), which usually makes other methods fail.

Why This Matters (According to the Paper)

No Peeking Required: You can solve these massive problems even if you can't see the individual numbers, only how they react to inputs.
Efficiency: The time it takes to solve the problem grows linearly with the size of the problem. If you double the size of the puzzle, it takes roughly double the time, not a million times longer.
Robustness: It works well for difficult 3D problems, like simulating sound waves (Helmholtz equation) or heat flow, where other methods often get stuck or take too long.

In short, the paper presents a way to take a giant, invisible, complex mathematical puzzle, throw some random darts at it, and use the splashes to build a skeleton key that solves the puzzle quickly and accurately, without ever needing to see the puzzle pieces themselves.

Technical Summary: Randomized Strong Recursive Skeletonization

Problem Statement
Hierarchical matrices ( $H^2$ -matrices) under the strong admissibility condition offer a framework for compressing dense matrices arising from partial differential equations (PDEs) and integral equations. These matrices partition interactions into near-field (stored densely) and far-field (approximated by low-rank blocks) components, utilizing a nested basis structure to achieve linear storage complexity and efficient matrix-vector multiplication. While fast matrix-vector products are well-established, the direct inversion (LU factorization) of $H^2$ -matrices remains challenging, particularly for three-dimensional problems. Existing inversion methods typically require explicit access to matrix entries and rely on repeated recompression of off-diagonal blocks during the recursive factorization process. Furthermore, obtaining an $H^2$ representation for matrices available only through their action on vectors (matrix-free settings) has been an open problem, with previous randomized approaches suffering from large prefactors or requiring structured test matrices that are difficult to design dynamically.

Methodology
The paper introduces Randomized Strong Recursive Skeletonization (RSRS), an algorithm that simultaneously compresses and factorizes general $H^2$ -matrices using only the action of the matrix $A$ and its adjoint $A^*$ on vectors. The method operates under the strong admissibility condition and avoids direct access to matrix entries.

The core of the methodology relies on two main components:

Randomized Sketching with Dense Gaussian Matrices: Instead of designing structured test matrices with specific zero patterns (which increases sample complexity and requires prior knowledge of skeleton indices), RSRS employs dense Gaussian test matrices $\Omega$ and $\Psi$ . The algorithm generates initial sketches $Y = A\Omega$ and $Z = A^*\Psi$ .
Dynamic Structural Imposition via Linear Transformations: As the factorization proceeds recursively through a hierarchical tree (from finest to coarsest levels), the algorithm introduces necessary structure into the test matrices via post-processing. Two key operations are applied to the dense sketches:
- Block Nullification: Computes a nullspace basis to annihilate contributions from near-field indices, effectively isolating far-field interactions for low-rank compression (Interpolative Decomposition).
- Block Extraction: Uses pseudoinverse transformations to extract specific dense near-field subblocks required for the elimination steps.

The algorithm recursively applies a "sparsify-then-diagonalize" procedure. At each box in the hierarchy, it computes an interpolative decomposition (ID) of the far-field interactions to identify redundant indices (skeletons). It then applies block elimination matrices to decouple these indices from the rest of the system. Crucially, the sketch matrices ( $Y, Z$ ) and test matrices ( $\Omega, \Psi$ ) are updated at each step using the derived elimination and sparsification matrices ( $L, U, E, F$ ). This allows the algorithm to maintain sketches of the evolving diagonalized system without ever accessing individual matrix entries or redrawing random samples.

Key Contributions

Matrix-Free Inversion: RSRS extends $H^2$ inversion to settings where matrix entries are inaccessible, relying solely on matrix-vector products. This broadens applicability to integral operators solvable via Fast Multipole Methods (FMM) or FFTs, and to dense blocks arising in sparse PDE solvers.
Simultaneous Compression and Factorization: Unlike traditional methods that separate compression and inversion, RSRS integrates these steps. This avoids the storage-intensive data structures and repeated recompressions required by previous approaches, maintaining updated off-diagonal blocks in compressed form.
Sample Efficiency: The number of required matrix-vector products ( $s$ ) is independent of the matrix size $N$ and depends linearly only on the maximal rank of off-diagonal blocks and the problem geometry. The method achieves this by using generic dense Gaussian test matrices, avoiding the high sample costs associated with structured test matrices.
Robustness: The algorithm is designed to handle indefinite and ill-conditioned matrices, such as those arising from 3D Helmholtz equations, without requiring the matrix to be positive definite.

Numerical Results
The authors validate RSRS through three numerical experiments:

2D Volume Integral Equation: A first-kind Laplace integral equation on a uniform grid. Despite the matrix becoming increasingly ill-conditioned as $N$ grows (condition number $\sim 10^6$ ), RSRS maintained low relative compression errors ( $10^{-7}$ to $10^{-6}$ ) and reduced GMRES iterations from thousands to fewer than 5 when used as a preconditioner.
3D Boundary Integral Equation: A second-kind Fredholm equation on a sphere. The method demonstrated stable accuracy and linear scaling in reconstruction time, memory, and solve time with respect to $N$ .
3D Helmholtz Schur Complement: A dense Schur complement arising from a sparse direct solver for the 3D Helmholtz equation. The matrix is indefinite. RSRS successfully compressed the Schur complement, reducing GMRES iterations from non-convergence (over 10,000) to single digits (3–9) for various problem sizes.

In all cases, the reconstruction time ( $T_{rec}$ ), solve time ( $T_{sol}$ ), and memory usage scaled linearly with $N$ , consistent with the theoretical complexity analysis of $O(k^2 N)$ and $O(k N)$ .

Significance and Claims
The paper claims that RSRS represents a significant departure from traditional $H^2$ inversion techniques by enabling direct factorization in matrix-free environments. By combining compression and inversion into a single recursive process driven by randomized sketching, the method offers competitive computational efficiency and reduced storage requirements. The authors emphasize that the use of dense Gaussian test matrices, coupled with dynamic linear transformations, eliminates the need for complex, problem-specific test matrix designs. The method is presented as a robust approximate direct solver suitable for a wide range of applications, including boundary integral equations, PDE discretizations, and uncertainty quantification, particularly where matrix entries are expensive or unavailable. The work suggests that the feedback effect of using sketches to define the transformations that update those sketches is benign in practice, allowing for accurate factorization without redrawing samples.