# A composable abstraction of hierarchical methods for matrix-vector product acceleration

Antoine Gicquel, Olivier Coulaud, Bérenger Bramas, Emmanuel Agullo

► To cite this version:

Antoine Gicquel, Olivier Coulaud, Bérenger Bramas, Emmanuel Agullo. A composable abstraction of hierarchical methods for matrix-vector product acceleration. Workshop on Approximate Computing in Numerical Linear Algebra 2025, Oct 2025, Paris, France. hal-05397428

HAL Id: hal-05397428

https://inria.hal.science/hal-05397428v1

Submitted on 4 Dec 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A composable abstraction of hierarchical methods for matrix-vector product acceleration

Antoine Gicquel

Joint-work with: Olivier Coulaud, Bérenger Bramas, Emmanuel Agullo

Inria Bordeaux (Concace team)

Workshop on Approximate Computing in Numerical Linear Algebra

7-10 October 2025

<figure><img src="image_2.png" /></figure>

SORBONNE
UNIVERSITÉ
# Introduction

## Matrix-vector product

### Context

Focus on the **unification of matrix-vector product** algorithms

$$ y = Ax \tag{1} $$

where $A \in \mathbb{R}^{N \times N}$ and $x, y \in \mathbb{R}^N$.
# Introduction

## Matrix-vector product

### Context

Focus on the **unification of matrix-vector product** algorithms

$$y = Ax \quad (1)$$

where $A \in \mathbb{R}^{N \times N}$ and $x, y \in \mathbb{R}^N$.

* If $A$ is a **rank-structured** matrix, time and memory complexity of (1) can be reduced from $O(N^2)$ to $O(N \log N)$ or $O(N)$ with a prescribed accuracy $\epsilon$ using **hierarchical methods**.

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>
# Introduction

## Model problem

### 1D N-body problem

<figure><img src="image_7.png" /></figure>

$$ \varphi_j = \sum_{i=1, i \neq j}^{N} \frac{1}{|s_i - s_j|} q_i \quad j \le N $$

$$ \mathbf{y} = \mathbf{A}\mathbf{x} \qquad \mathbf{A}_{ij} = \delta_{ij} \frac{1}{|s_i - s_j|} $$

<figure><img src="image_8.png" /><figcaption>Matrix A</figcaption></figure>
# Introduction

## Model problem

### 1D N-body problem

<figure><img src="image_9.png" /></figure>

$$ \varphi_j = \sum_{i=1, i \neq j}^{N} \frac{1}{|s_i - s_j|} q_i \quad j \le N $$

$$ \mathbf{y} = \mathbf{A}\mathbf{x} \qquad \mathbf{A}_{ij} = \delta_{ij} \frac{1}{|s_i - s_j|} $$

<figure><img src="image_10.png" /><figcaption>Matrix A</figcaption></figure>
# Introduction

## Model problem

### 1D N-body problem

<figure><img src="image_11.png" /></figure>

$$ \varphi_j = \sum_{i=1, i \neq j}^{N} \frac{1}{|s_i - s_j|} q_i \quad j \le N $$

$$ \mathbf{y} = \mathbf{A}\mathbf{x} \qquad \mathbf{A}_{ij} = \delta_{ij} \frac{1}{|s_i - s_j|} $$

<figure><img src="image_12.png" /><figcaption>Matrix A</figcaption></figure>

### Low-rank approximation

* Tiles with low numerical rank (green tiles)

* Factorization of the form

$$ A_{IJ} \approx \hat{F}\hat{G}^* $$

<figure><img src="image_13.png" /></figure>
# Introduction

## Model problem

<figure><img src="image_14.png" /></figure>

Flat partitioning
⇒ single-level tiling

<figure><img src="image_15.png" /></figure>

Hierarchical partitioning
⇒ multi-level tiling

## Applications (rank-structured matrices)

N-body-like problems, integral equations (BEM, VIM, ...)
# Introduction

## Taxonomy (Hierarchical methods) [Ashcraft et al. 2021]

<figure><img src="image_16.png" /></figure>

* Some pre-existing work on the unification of the FMM, $H$, $H^2$ and their variants [Yokota et al. 2017]

* Our goal: **Composable abstraction + generic matrix-vector product algorithm**
# Abstraction of a hierarchical method

## Component 1 - reordering and partitioning strategy

<figure><img src="image_17.png" /></figure>

Without reordering

<figure><img src="image_18.png" /></figure>

With reordering

<figure><img src="image_19.png" /></figure>

* **Geometric:** space-filling curves (e.g., Morton, Hilbert indexing, ...)

* **Algebraic:** graph partitioning
# Abstraction of a hierarchical method

## Component 1 - reordering and partitioning strategy

### Partitioning tree $\mathcal{T}$

* Partition of $\{1, ..., N\}$

* Based on the geometry: 1 node $l_i \leftrightarrow 1$ cluster $\Omega_i \leftrightarrow 1$ subset of $\{1, ..., N\}$

* **n-level partitioning**

* $n = 1$: **single-level partitioning**

* $n \ge 2$: **multi-level partitioning**

* $\mathcal{T}_r$: **row partitioning tree**

* $\mathcal{T}_c$: **column partitioning tree**

<figure><img src="image_20.png" /></figure>
# Abstraction of a hierarchical method

## Component 2 - admissibility criterion

* A pair of nodes $(I, J) \in T_r \times T_c$ defines a tile $A_{IJ}$

<figure><img src="image_21.png" /></figure>

## Admissibility criterion

Binary decision $\lambda : (I, J) \in T_r \times T_c \mapsto \{\text{false, true}\}$ that determines whether a tile $A_{IJ}$ can be **approximated** (i.e., is **admissible**) or not
# Abstraction of a hierarchical method

## Component 2 - admissibility criterion

### Examples

* Weak admissibility: $A_{IJ}$ is approximated if $I \neq J$

* Strong admissibility: $A_{IJ}$ is approximated if $\Omega_I$ and $\Omega_J$ are "far enough"

<figure><img src="image_22.png" /><figcaption>Single-level part.
+ Weak adm.</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Multi-level part.
+ Weak adm.</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>Single-level part.
+ Strong adm.</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Multi-level part.
+ Strong adm.</figcaption></figure>
# Abstraction of a hierarchical method

## Component 3 - approximation strategy

* For an admissible tile $A_{IJ}$:

<figure><img src="image_26.png" /></figure>

<figure><img src="image_27.png" /></figure>
# Abstraction of a hierarchical method

## Component 3 - approximation strategy

* For an **admissible** tile $A_{IJ}$:

<figure><img src="image_28.png" /></figure>

<figure><img src="image_29.png" /></figure>

* **Algebraic (low-rank) approximation** (randomized SVD, ACA, ID, ...)

$$A_{IJ} \approx B_{IJ} C_{IJ}^*$$

$$A_{IJ} \approx V_I S_{IJ} W_J^*$$
# Abstraction of a hierarchical method

## Component 3 - approximation strategy

* For an **admissible** tile $A_{IJ}$:

<figure><img src="image_30.png" /></figure>

<figure><img src="image_31.png" /></figure>

* **Algebraic (low-rank) approximation (randomized SVD, ACA, ID, ...)**

$$A_{IJ} \approx B_{IJ} C_{IJ}^*$$

$$A_{IJ} \approx V_I S_{IJ} W_J^*$$

* **Analytical (low-rank) approximation (interpolation, spherical harmonics, Taylor series, ...). Assuming that $A_{ij} = k(s_i, s_j)$:**

$$k(x, y) \approx \sum_{\lambda} \varphi_{\lambda}(x) \gamma_{\lambda}(y)$$

$$k(x, y) \approx \sum_{\alpha} \Phi_{\alpha}(x) \sum_{\beta} S_{\alpha, \beta} \Psi_{\beta}(y)$$
# Abstraction of a hierarchical method

## Component 3 - approximation strategy

* For an **admissible** tile $A_{IJ}$:

<figure><img src="image_32.png" /></figure>

<figure><img src="image_33.png" /></figure>

* **Algebraic (low-rank) approximation** (randomized SVD, ACA, $ID, \dots$)

$$A_{IJ} \approx B_{IJ} C_{IJ}^*$$

$$A_{IJ} \approx V_I S_{IJ} W_J^*$$

* **Analytical (low-rank) approximation** (interpolation, spherical harmonics, Taylor series...). Assuming that $A_{ij} = k(s_i, s_j)$:

$$k(x, y) \approx \sum_{\lambda} \varphi_{\lambda}(x) \gamma_{\lambda}(y)$$

$$k(x, y) \approx \sum_{\alpha} \Phi_{\alpha}(x) \sum_{\beta} S_{\alpha, \beta} \Psi_{\beta}(y)$$

* **Arithmetic approximation** (mixed-precision, adaptive accuracy, ...) [cf. C. Peaucelle's talk]
# Abstraction of a hierarchical method

## Component 4 - basis mutualization

* **Independent bases**

$$A_{IJ} \approx B_{IJ} C_{IJ}^*$$

<figure><img src="image_34.png" /></figure>
# Abstraction of a hierarchical method

## Component 4 - basis mutualization

* **Independent bases**

$$A_{IJ} \approx B_{IJ} C_{IJ}^*$$

<figure><img src="image_35.png" /></figure>

* **Shared bases**

$$A_{IJ} \approx V_I S_{IJ} W_J^*$$

<figure><img src="image_36.png" /></figure>
# Abstraction of a hierarchical method

## Component 4 - basis mutualization

<figure><img src="image_37.png" /></figure>

<figure><img src="image_38.png" /></figure>

$$A_{l_0 J_0} \approx V_{l_0} S_{l_0 J_0} W_{J_0}^*$$
# Abstraction of a hierarchical method

## Component 4 - basis mutualization

<figure><img src="image_39.png" /></figure>

## Nested bases [Börm 2010; Greengard and Rokhlin 1987]

* We can build the bases of any level from the bases of the leaf level

* We only need the translation operators $\{B_{J_k J_0}\}_{k=1,2}$ and $\{E_{l_k l_0}\}_{k=1,2}$
# Abstraction of a hierarchical method

## Component 4 - basis mutualization

### Shared bases

<figure><img src="image_40.png" /></figure>

### Nested bases

<figure><img src="image_41.png" /></figure>

### Matrix-vector product

Asymptotic cost: $O(N \log N)$

### Matrix-vector product

Asymptotic cost: $O(N)^1$

<footnote n="1">¹Under certain assumptions</footnote>
# Abstraction of a hierarchical method

Abstraction of a
hierarchical method
Abstraction of a hierarchical method

<figure><img src="image_42.png" /></figure>

• Reordering and partitioning:

► **Subdivide matrix** A into **tiles** {$A_{ij}^{(\ell)}$} at different levels $\ell$
Abstraction of a hierarchical method

<figure><img src="image_43.png" /></figure>

* **Reordering and partitioning:**

* **Subdivide matrix A into tiles {$A_{ij}^{(\ell)}$} at different levels $\ell$**

* **Admissibility criterion:**

* **Determines which tiles $A_{ij}^{(\ell)}$ are approximated (or admissible)**
Abstraction of a hierarchical method

<figure><img src="image_44.png" /></figure>

* Reordering and partitioning:

* Subdivide matrix $A$ into tiles $\{A_{IJ}^{(\ell)}\}$ at different levels $\ell$

* Admissibility criterion:

* Determines which tiles $A_{IJ}^{(\ell)}$ are approximated (or admissible)

* Approximation strategy:

* Approximates the admissible tiles $A_{IJ}^{(\ell)} \approx \hat{V}_I^{(\ell)} \hat{A}_{IJ}^{(\ell)} \hat{W}_J^{(\ell)}$*
Abstraction of a hierarchical method

<figure><img src="image_45.png" /></figure>

• Reordering and partitioning:

► **Subdivide matrix** $\mathbf{A}$ into **tiles** $\{A_{ij}^{(\ell)}\}$ at different levels $\ell$

• Admissibility criterion:

► **Determines which tiles** $A_{ij}^{(\ell)}$ **are approximated (or admissible)**

• Approximation strategy:

► **Approximates the admissible tiles** $A_{ij}^{(\ell)} \approx \hat{V}_I^{(\ell)} \hat{A}_{IJ}^{(\ell)} \hat{W}_J^{(\ell)}$*

• Basis mutualization:

► **Factorizes the approximation bases** $\hat{V}_I^{(\ell)}$ and $\hat{W}_J^{(\ell)}$
# Generic algorithm for the hierarchical matrix-vector product

## Prior to this work:

* Each method had its *own matrix-vector product algorithm*

* Each algorithm leverages the *method's specific features*
# Generic algorithm for the hierarchical matrix-vector product

## Prior to this work:

* Each method had its **own matrix-vector product algorithm**

* Each algorithm leverages the **method's specific features**

<caption>Algorithm 1: Generic H matrix-vector *</caption>

1 Procedure H-mv(A, x, y):
2     /* Compression stage */
3     ($\hat{x}^{(1)},...,\hat{x}^{(L)}$) = compress(x);
4     /* Multiplication stage */
5     for $\ell = L,...,1$ do
6         $\hat{y}_{\ell}^{(L)}$ = $\hat{A}_{\ell J}^{(L)} \hat{x}_{J}^{(L)}$;
7         end
8     end
9     for each tile $(I,J)$ at level 0 do
10        $y_I$ += $A_{IJ}^{(0)} x_J$;
11 end
12
13 /* Decompression stage */
14 y += decompress($\hat{y}^{(1)},...,\hat{y}^{(L)}$);

<figure><img src="image_46.png" /></figure>
### Particle interaction (3D)

* $N$ particles located at $s_i \in \Omega$ (prolate ellipsoid) with charge $q_i$

* **Goal:** Compute the potential $\varphi_i$ induced by pairwise interactions

$$ \varphi_j = \sum_{i=1, i \neq j}^{N} \frac{1}{\|s_i - s_j\|} q_i \quad \forall j = 1, \dots, N $$

$$ y = Ax \quad \text{with} \quad \begin{cases} A_{ij} = \delta_{ij} \frac{1}{\|s_i - s_j\|} \\ x_j = q_j \\ y_i = \varphi_i \end{cases} $$

<figure><img src="image_47.png" /><figcaption>Prolate ellipsoid</figcaption></figure>

### Implementation

* Sequential prototype in composyx (C++ / Comosability)

* hardware: 2x18-core Cascade Lake Intel Xeon Skylake Gold 6240 @ 2.6GHz

### Reference software

* ScalFMM (C++ / FMM), hmat-oss (C++ / $\mathcal{H}$-matrix)
# Numerical results

## Numerical validation

### Tested configurations

<table><thead><tr><th></th><th>Partitioning</th><th>Admissibility</th><th>Approximation</th><th>Bases</th></tr></thead><tbody><tr><th>FMM</th><td>multi-level<br>(adaptive-octree)</td><td>strong</td><td>analytical<br>(Chebyshev interpolation)</td><td>nested</td></tr><tr><th>H</th><td>multi-level<br>(binary tree)</td><td>strong</td><td>algebraic (ACA)</td><td>independent</td></tr></tbody></table>
# Numerical results

## Numerical validation

### Tested configurations

<table><thead><tr><th></th><th>Partitioning</th><th>Admissibility</th><th>Approximation</th><th>Bases</th></tr></thead><tbody><tr><th>FMM</th><td>multi-level<br>(adaptive-octree)</td><td>strong</td><td>analytical<br>(Chebyshev interpolation)</td><td>nested</td></tr><tr><th>H</th><td>multi-level<br>(binary tree)</td><td>strong</td><td>algebraic (ACA)</td><td>independent</td></tr></tbody></table>

### Relative norm-2 error

* Direct computation: $y_{exact} = Ax$

* Generic algorithm: $y_{approx} = \hat{A}x$

$$ Error = \frac{\|y_{exact} - y_{approx}\|_2}{\|y_{exact}\|_2} $$

* Parameter controlling the accuracy: $\varepsilon$

<figure><img src="image_48.png" /></figure>
# Numerical results

## Preliminary benchmarks

### Tested configurations

<table><thead><tr><th></th><th>Partitioning</th><th>Admissibility</th><th>Approximation</th><th>Bases</th></tr></thead><tbody><tr><th>FMM</th><td>multi-level<br>(adaptive-octree)</td><td>strong</td><td>analytical<br>(Chebyshev interpolation)</td><td>nested</td></tr><tr><th>H</th><td>multi-level<br>(binary tree)</td><td>strong</td><td>algebraic (ACA)</td><td>independent</td></tr></tbody></table>
# Numerical results

## Preliminary benchmarks

### Tested configurations

<table><thead><tr><th></th><th>Partitioning</th><th>Admissibility</th><th>Approximation</th><th>Bases</th></tr></thead><tbody><tr><th>FMM</th><td>multi-level<br>(adaptive-octree)</td><td>strong</td><td>analytical<br>(Chebyshev interpolation)</td><td>nested</td></tr><tr><th>H</th><td>multi-level<br>(binary tree)</td><td>strong</td><td>algebraic (ACA)</td><td>independent</td></tr></tbody></table>

<figure><img src="image_49.png" /></figure>

<figure><img src="image_50.png" /></figure>
# Numerical results

## Global comparison

### Tested configurations

<table><thead><tr><th></th><th>Partitioning</th><th>Admissibility</th><th>Approximation</th><th>Bases</th></tr></thead><tbody><tr><td>FMM</td><td>multi-level<br>(adaptive-octree)</td><td>strong</td><td>analytical<br>(Chebyshev interpolation)</td><td>nested</td></tr><tr><td>H</td><td>multi-level<br>(binary tree)</td><td>strong</td><td>algebraic (ACA)</td><td>independent</td></tr><tr><td>HODLR</td><td>multi-level<br>(binary tree)</td><td>weak</td><td>algebraic (ACA)</td><td>independent</td></tr><tr><td>BLR</td><td>single-level</td><td>strong</td><td>algebraic (ACA)</td><td>independent</td></tr></tbody></table>
# Numerical results

## Global comparison

### Tested configurations

<table><thead><tr><th></th><th>Partitioning</th><th>Admissibility</th><th>Approximation</th><th>Bases</th></tr></thead><tbody><tr><td>FMM</td><td>multi-level<br>(adaptive-octree)</td><td>strong</td><td>analytical<br>(Chebyshev interpolation)</td><td>nested</td></tr><tr><td>H</td><td>multi-level<br>(binary tree)</td><td>strong</td><td>algebraic (ACA)</td><td>independent</td></tr><tr><td>HODLR</td><td>multi-level<br>(binary tree)</td><td>weak</td><td>algebraic (ACA)</td><td>independent</td></tr><tr><td>BLR</td><td>single-level</td><td>strong</td><td>algebraic (ACA)</td><td>independent</td></tr></tbody></table>

<figure><img src="image_51.png" /></figure>
# Conclusion

## Summary

* A **composable abstraction** unifying hierarchical methods

* 4 core components: partitioning, admissibility, approximation, bases

* Composability at the algorithmic level

* A **generic algorithm** for the matrix-vector product

* Bridges the gap between FMM, $H$, $H^2$, HODLR, BLR, ...

* Implementation of a **sequential prototype** in the composyx library

* Consistent results with existing software
# Conclusion

## Summary

* A **composable abstraction** unifying hierarchical methods

* 4 core components: partitioning, admissibility, approximation, bases

* Composability at the algorithmic level

* A **generic algorithm** for the matrix-vector product

* Bridges the gap between FMM, $H$, $H^2$, HODLR, BLR, ...

* Implementation of a **sequential prototype** in the composyx library

* Consistent results with existing software

## Perspectives

* Implementation of some **optimizations** (e.g., M2L operator for the FMM + symmetries) to increase performance

* **Parallelization** of the generic algorithm

* Shift focus to **factorization** algorithms

* More applications to "real-world" problems
References

<figure><img src="image_52.png" /></figure>

Amestoy, Patrick, Cleve Ashcraft, Olivier Boiteau, Alfredo Buttari, Jean-Yves l'Excellent, and Clément Weisbecker (2015). "Improving multifrontal methods by means of block low-rank representations". In: *SIAM Journal on Scientific Computing* 37.3, A1451–A1474.

<figure><img src="image_53.png" /></figure>

Ashcraft, Cleve, Alfredo Buttari, and Theo Mary (2021). "Block Low-Rank Matrices with Shared Bases: Potential and Limitations of the BLR ^2 Format". In: *SIAM Journal on Matrix Analysis and Applications* 42.2, pp. 990–1010.

<figure><img src="image_54.png" /></figure>

Börm, Steffen (2010). *Efficient numerical methods for non-local operators: H2-matrix compression, algorithms and analysis*. Vol. 14. European Mathematical Society.

<figure><img src="image_55.png" /></figure>

Börm, Steffen, Lars Grasedyck, and Wolfgang Hackbusch (2003). "Introduction to hierarchical matrices with applications". In: *Engineering analysis with boundary elements* 27.5, pp. 405–422.

<figure><img src="image_56.png" /></figure>

Greengard, Leslie and Vladimir Rokhlin (1987). "A fast algorithm for particle simulations". In: *Journal of computational physics* 73.2, pp. 325–348.

<figure><img src="image_57.png" /></figure>

Yokota, Rio, Huda Ibeid, and David Keyes (2017). "Fast multipole method as a matrix-free hierarchical low-rank approximation". In: *Eigenvalue Problems: Algorithms, Software and Applications in Petascale Computing: EPASA 2015, Tsukuba, Japan, September 2015*. Springer, pp. 267–286.