H²-Matrix-Matrix Multiplication

• H²-matrix-matrix multiplication is a technique that leverages nested low-rank bases to compress matrices and reduce storage complexity.
• The method employs forward and backward transformations with BLAS-3 operations to handle dense blocks and improve arithmetic intensity and cache reuse.
• Adaptive basis construction and rigorous error control ensure the hierarchical structure is maintained, enabling efficient large-scale applications such as PDE solvers and electromagnetic simulations.

Searching arXiv for papers on $\mathcal H^2$-matrix-matrix multiplication and closely related algorithms. $\mathcal H^2$-matrix-matrix multiplication denotes a family of algorithms for multiplying hierarchical, data-sparse matrices that exploit nested low-rank bases. In the $\mathcal H^2$ format, admissible off-diagonal blocks are represented as

$$G_{\hat t \times \hat s} \approx V_t\, S_b\, W_s^*,$$

with local cluster bases $V_t$, $W_s$ and small coupling matrices $S_b$, while inadmissible blocks are stored as small dense blocks. This structure reduces storage and matrix-vector multiplication to $\mathcal{O}(nk)$ for an $n$-dimensional matrix with local rank $k$, and it also furnishes the basis for more advanced multiplication tasks: applying an $\mathcal H^2$0-matrix to a dense block of right-hand sides, and multiplying two $\mathcal H^2$1-matrices while preserving or reconstructing an $\mathcal H^2$2 representation of the product (Christophersen, 21 Sep 2025, Börm, 2024).

1. Representation and algebraic setting

$\mathcal H^2$3-matrices are hierarchical matrices designed for dense operators with local low-rank structure, including matrices arising from boundary element methods, solution operators of partial differential equations, and integral operators (Christophersen, 21 Sep 2025, Börm, 2024). Their defining feature is nestedness: bases at coarser clusters are built from child bases through transfer matrices, so basis information is not duplicated across the hierarchy. This nested structure distinguishes $\mathcal H^2$4-matrices from classical $\mathcal H^2$5-matrices and is the reason storage can be reduced to $\mathcal H^2$6 rather than $\mathcal H^2$7 (Börm, 2024).

For admissible leaves, the low-rank form

$\mathcal H^2$8

or, in equivalent notation,

$\mathcal H^2$9

encodes the block through nested row and column cluster bases and a small coupling matrix (Ma et al., 2019, Börm, 2023). Inadmissible or nearfield blocks remain dense. The resulting arithmetic is hierarchical: traversal over row and column cluster trees replaces direct dense multiplication, and small-rank operations replace large dense ones wherever admissibility holds.

This representation creates two distinct multiplication problems. The first is $\mathcal H^2$0-matrix times a dense matrix of $\mathcal H^2$1 right-hand sides, which is central in block iterative methods. The second is $\mathcal H^2$2, in which both factors are hierarchically compressed and the output must again be approximated in $\mathcal H^2$3 form. The latter is algorithmically harder because the product generally does not lie in the original cluster bases of the factors (Christophersen, 21 Sep 2025, Ma et al., 2019).

2. Multiplication with dense blocks of right-hand sides

When an $\mathcal H^2$4-matrix is multiplied by a dense matrix $\mathcal H^2$5, the vector algorithm of Börm and Hackbusch is promoted to a block algorithm with three stages: forward transformation, multiplication phase, and backward transformation (Christophersen, 21 Sep 2025).

In the forward transformation, the column cluster tree is traversed bottom-up. For a leaf cluster $\mathcal H^2$6, one computes

$\mathcal H^2$7

while non-leaf clusters combine children through transfer matrices. The effect is to project the $\mathcal H^2$8 right-hand sides into the nested column bases. In the multiplication phase, each admissible block $\mathcal H^2$9 contributes

$$G_{\hat t \times \hat s} \approx V_t\, S_b\, W_s^*,$$0

a $$G_{\hat t \times \hat s} \approx V_t\, S_b\, W_s^*,$$1 matrix-matrix multiplication on small coupling matrices, whereas each inadmissible block contributes a dense block multiplication

$$G_{\hat t \times \hat s} \approx V_t\, S_b\, W_s^*,$$2

Finally, the backward transformation traverses the row cluster tree top-down and reconstructs the output through the row bases and transfer matrices (Christophersen, 21 Sep 2025).

A practical implementation introduces a preparation phase that builds, for each row cluster, a block-list $$G_{\hat t \times \hat s} \approx V_t\, S_b\, W_s^*,$$3 of all relevant leaf blocks. This permits conflict-free parallel processing because all row clusters can be processed concurrently, and forward and backward transformations can be parallelized across clusters on the same tree level. The arithmetic is systematically cast into BLAS level-3 form: admissible blocks use GEMM on $$G_{\hat t \times \hat s} \approx V_t\, S_b\, W_s^*,$$4 and $$G_{\hat t \times \hat s} \approx V_t\, S_b\, W_s^*,$$5 operands, inadmissible dense blocks use GEMM on all $$G_{\hat t \times \hat s} \approx V_t\, S_b\, W_s^*,$$6 columns at once, and transfer-matrix applications are likewise promoted to batched or block BLAS-3 operations (Christophersen, 21 Sep 2025).

The significance of this formulation is hardware-oriented. Grouping many right-hand sides into a single dense matrix raises arithmetic intensity, improves cache reuse, and avoids the poor locality of repeated GEMV-based matrix-vector products. In the reported experiments, the block $$G_{\hat t \times \hat s} \approx V_t\, S_b\, W_s^*,$$7-matrix-matrix multiplication achieved speedups of 10–13x for $$G_{\hat t \times \hat s} \approx V_t\, S_b\, W_s^*,$$8, and up to 50x for large $$G_{\hat t \times \hat s} \approx V_t\, S_b\, W_s^*,$$9 and suitable chunk sizes; for sufficiently large $V_t$0 and sufficiently high-accuracy matrices, the implementation saturated available memory bandwidth (Christophersen, 21 Sep 2025).

3. Multiplication of two $V_t$1-matrices

The product of two $V_t$2-matrices is more delicate than multiplication by a dense block because the natural bases of the product are not, in general, the original bases of the factors. Recent work therefore centers on adaptive basis construction and recompression (Ma et al., 2019, Börm, 2024, Börm, 2023).

A 2023 algorithm presents a two-phase strategy. First, it constructs an intermediate representation of $V_t$3 on an induced block tree with induced cluster bases; second, it performs structural recompression so that the result matches a prescribed output block tree (Börm, 2023). The induced row basis for a cluster $V_t$4 is described as

$V_t$5

where the additional terms encode product contributions that cannot be represented in the original basis alone. These bases are then compressed recursively by SVD-based low-rank approximation while preserving nested transfer relations (Börm, 2023).

A 2024 algorithm emphasizes mathematically rigorous error control by constructing the exact product in an intermediate hierarchical representation before any truncation is applied (Börm, 2024). Its central data structures are cluster trees, block trees, and a 3-way product tree whose nodes are triples $V_t$6 tracking recursive products $V_t$7. For admissible blocks,

$V_t$8

with precomputed basis product matrices

$V_t$9

The paper also uses semi-uniform matrices of the form $W_s$0 as recursive accumulators (Börm, 2024).

A 2019 approach pursues explicit accuracy control while preserving the original $W_s$1 structure of the product (Ma et al., 2019). Instead of fixing the cluster bases beforehand, it instantaneously updates them during the matrix product computation according to a prescribed truncation accuracy $W_s$2. This is framed as an alternative to earlier formatted multiplication algorithms that inherited cluster bases directly from the input factors and therefore did not control the product accuracy (Ma et al., 2019).

4. Basis adaptation and error control

Accuracy control is the central differentiator among $W_s$3 algorithms. Earlier formatted multiplication schemes could retain linear complexity only by assuming that a suitable basis for the result was already known or by reusing the bases of the factors. Later work argues that this assumption is often invalid, because multiplication can generate new subspaces that are not captured by the original nested bases (Ma et al., 2019, Börm, 2023).

In the 2019 method, basis adaptation is performed during the product itself. For each product block, the algorithm collects contributions into Gram matrices, normalizes them, computes an SVD of the aggregated Gram matrix, and truncates singular vectors whose singular values fall below the prescribed threshold $W_s$4. This basis update proceeds recursively from leaves to the root, so that all additions and multiplications are either exact or computed based on prescribed accuracy, and the product remains in $W_s$5 form throughout (Ma et al., 2019).

The 2024 method takes a different route: it delays approximation. The exact product is represented in an intermediate hierarchical structure, and low-rank truncations are applied only during coarsening and recompression. For each block $W_s$6, the paper states the block-relative spectral-norm guarantee

$W_s$7

obtained by choosing a local tolerance

$W_s$8

where $W_s$9 and $S_b$0 are the maximal numbers of block rows and columns in the merge (Börm, 2024).

The 2023 linear-complexity algorithm also uses adaptive compression with localized error estimates, but its guarantee is expressed relative to the contributing input blocks:

$S_b$1

The paper notes that full block-relative error is not always guaranteed in this formulation, although the reported practical error is very close (Börm, 2023).

A recurrent misconception is that linear-time $S_b$2 multiplication and explicit accuracy control are automatically compatible. The literature summarized here distinguishes these objectives sharply: some algorithms are optimized for linear complexity with adaptive bases, while others accept an additional logarithmic factor in exchange for rigorous block-relative error control (Börm, 2024, Börm, 2023).

5. Complexity, hardware efficiency, and algorithmic trade-offs

The principal algorithmic variants can be organized by the type of right-hand factor and the form of accuracy control.

Setting	Core mechanism	Stated complexity or guarantee
$S_b$3 dense block	Forward transform, block multiplication, backward transform; GEMM throughout	$S_b$4; linear in $S_b$5, $S_b$6, and linearly scalable with $S_b$7 (Christophersen, 21 Sep 2025)
$S_b$8 with basis updates	Instantaneous change of cluster bases using SVD truncation	$S_b$9 for constant rank; $\mathcal{O}(nk)$0 time and $\mathcal{O}(nk)$1 memory for variable-rank cases described in the paper (Ma et al., 2019)
$\mathcal{O}(nk)$2 with exact intermediate product	Product trees, basis trees, coarsening, recompression	$\mathcal{O}(nk)$3 with block-relative error estimates for submatrices (Börm, 2024)
Adaptive linear-complexity $\mathcal{O}(nk)$4	Induced block structure followed by structural recompression	$\mathcal{O}(nk)$5 for fixed $\mathcal{O}(nk)$6 under the stated tree and sparsity conditions (Börm, 2023)

For multiplication by dense blocks, complexity is dominated by $\mathcal{O}(nk)$7 work in the forward and backward transformations and the block multiplication phase. Each admissible block costs $\mathcal{O}(nk)$8, each inadmissible block costs $\mathcal{O}(nk)$9 for small block sizes $n$0, $n$1, and the total cost remains linear in $n$2 and $n$3 while scaling linearly with the number of columns $n$4 (Christophersen, 21 Sep 2025). The key implementation point is not only asymptotic count but arithmetic intensity: for $n$5, the ratio of FLOP to data moved is much higher than in repeated matrix-vector multiplication, so cache and memory bandwidth can be used more effectively (Christophersen, 21 Sep 2025).

For $n$6, the complexity depends on how the product basis is managed. The 2019 method gives

$n$7

time and

$n$8

memory, implying $n$9 time and memory for constant rank and $k$0 time with $k$1 memory when rank grows linearly with electrical problem size (Ma et al., 2019). The 2024 method reports $k$2, with the logarithmic factor reflecting recursive product-tree and recompression operations (Börm, 2024). The 2023 method states $k$3 under three technical assumptions: leaf size proportional to $k$4, bounded sparsity, and a bounded refinement ratio between fine and coarse block trees (Börm, 2023).

These algorithms therefore expose a genuine trade-off. This suggests that the multiplication problem is not characterized by a single “best” complexity class; rather, the practical choice depends on whether the dominant requirement is BLAS-3 throughput for many right-hand sides, linear-time adaptive basis construction, or rigorous blockwise error guarantees.

6. Applications, performance evidence, and limitations

The most immediate application of $k$5-matrix-matrix multiplication by dense blocks is in block Krylov subspace methods. The 2025 implementation uses the block multiplication kernel to accelerate block-CG and block-GMRES with hierarchical preconditioners, and reports that matrix operations become a small part of total runtime, with practical speedups of 30–50x over a standard sequential approach (Christophersen, 21 Sep 2025). This is specifically relevant when multiple systems with the same hierarchical matrix must be solved.

For $k$6, the applications emphasized in the literature include efficient preconditioners, matrix functions, and fast solvers for large-scale electromagnetic analysis (Börm, 2024, Ma et al., 2019). The 2019 accuracy-controlled product is demonstrated on large-scale capacitance extraction and electromagnetic scattering problems involving millions of unknowns on a single core, with time and memory scaling linearly in constant-rank capacitance examples and as predicted by theory in variable-rank scattering examples (Ma et al., 2019). The 2023 adaptive linear-complexity method reports multiplication for boundary element matrices of the 3D Laplacian, including examples up to $k$7, with the main multiplication completed in 53.5 seconds for the final matrix formation (on a single core) (Börm, 2023). The 2024 method states that experiments confirm multiplication for matrices of size $k$8 while maintaining the prescribed block-relative spectral error bound (Börm, 2024).

Several limitations recur. In the block right-hand-side setting, the memory footprint grows with the number of columns $k$9 because auxiliary matrices scale as $\mathcal H^2$00; for very large $\mathcal H^2$01, chunked processing is recommended (Christophersen, 21 Sep 2025). In the $\mathcal H^2$02 setting, the intermediate representation can be substantially larger than the final output: the 2023 method notes that it is sometimes 3x larger, although only transiently (Börm, 2023). Another limitation is methodological: the only linear-time algorithms known before the recent adaptive approaches either required the very special structure of HSS-matrices or required a suitable basis for all submatrices in advance (Börm, 2023).

Taken together, these developments define $\mathcal H^2$03-matrix-matrix multiplication as a mature but internally differentiated area of hierarchical linear algebra. One branch focuses on exploiting modern hardware through GEMM-based block application to many right-hand sides; another focuses on preserving or rebuilding $\mathcal H^2$04 structure under matrix products with adaptive bases; and a third prioritizes mathematically rigorous block-relative error control through exact intermediate representations and subsequent recompression (Christophersen, 21 Sep 2025, Börm, 2024).