Hierarchical Clustering for SpGEMM

• Hierarchical clustering for SpGEMM is a technique that groups matrix rows with shared nonzero indices to improve B-row reuse and lower data movement.
• It leverages row reordering, Jaccard similarity–based clustering, and a row-clustered CSR format to convert structural overlap into temporal locality.
• Empirical evaluations show speedups up to 4.68× with reduced preprocessing overhead compared to graph partitioning, highlighting its cost-effective performance gains.

Searching arXiv for the cited papers and closely related work on hierarchical clustering and SpGEMM. Hierarchical clustering for sparse matrix–sparse matrix multiplication (SpGEMM) is a locality-oriented technique that groups rows of the first input matrix $A$ so that rows with overlapping nonzero column indices are processed together, thereby improving reuse of rows of the second input matrix $B$ during multiplication. In the formulation developed for Gustavson-style SpGEMM, the method combines row reordering, cluster formation, and a row-clustered sparse format to convert structural overlap in $A$ into temporal locality for $B$, with the explicit goal of reducing data movement in a kernel that is typically memory-bound rather than arithmetic-bound (Islam et al., 28 Jul 2025).

1. Computational setting and locality objective

SpGEMM is a key kernel in many scientific applications and graph workloads, yet it is bottlenecked by data movement due to irregular memory access patterns. In Gustavson’s row-wise algorithm, processing row $i$ of $A$ requires streaming the rows of $B$ indexed by all nonzero column positions $p$ in that row, and then scattering products into an irregular sparse accumulator for row $i$ of $C$. Two sources of irregularity dominate: the set of rows of $B$0 touched by adjacent rows of $B$1 changes rapidly and often overlaps, forcing repeated loads of the same $B$2 rows, and the accumulator accesses are sparse and unpredictable (Islam et al., 28 Jul 2025).

For $B$3, $B$4, and $B$5, let the set of nonzero column indices in row $B$6 of $B$7 be

$B$8

The scalar multiply-add count in Gustavson’s algorithm is

$B$9

where $A$0 is the set of nonzeros in row $A$1 of $A$2 (Islam et al., 28 Jul 2025).

The locality opportunity is cross-row reuse in $A$3. If rows $A$4 and $A$5 satisfy $A$6, then they both require the same rows $A$7 for $A$8. When the two rows are processed with sufficient temporal proximity, those $A$9 rows can remain hot in cache and avoid reloads from main memory. This is the key distinction from many SpMV-oriented reorderings, which mostly target bandwidth reduction or adjacency of nonzeros to a dense vector. For SpGEMM, simply reordering rows of $B$0 while retaining CSR row-major traversal does not ensure that $B$1’s hot rows remain in cache across adjacent $B$2 rows (Islam et al., 28 Jul 2025).

To formalize this dependence structure, the method defines a row-intersection graph $B$3 over rows of $B$4, with $B$5, and for $B$6, an undirected edge $B$7 weighted by

$B$8

A partition $B$9 of the rows induces a cut

$i$0

A locality-maximizing clustering then seeks to minimize

$i$1

subject to balance constraints $i$2 and optionally $i$3. Equivalently, it maximizes the total intra-cluster overlap (Islam et al., 28 Jul 2025).

The same objective can be expressed directly in terms of data movement on $i$4. For a cluster $i$5, define the set of unique $i$6 rows touched by the cluster as

$i$7

If $i$8 is the number of bytes needed to stream row $i$9 of $A$0, then the row-wise traffic model is

$A$1

whereas the cluster-wise model is

$A$2

The reduction factor is approximately

$A$3

If $A$4 fits in cache, the benefit approaches the ideal suggested by $A$5 (Islam et al., 28 Jul 2025).

2. Construction of hierarchical clusters

The hierarchical clustering method is a bottom-up merging scheme driven by overlap estimates obtained from one binary SpGEMM between $A$6 and $A$7. The procedure sets all values of $A$8 to $A$9, computes

$B$0

and interprets each entry $B$1 as the exact overlap count $B$2. For each row $B$3, it keeps the topK neighbors $B$4 by Jaccard similarity

$B$5

above a threshold $B$6, producing candidate pairs. The paper states that this replaces expensive LSH while yielding exact overlap counts (Islam et al., 28 Jul 2025).

Cluster growth proceeds by greedily merging disjoint-set representatives using a max-heap keyed by similarity, subject to a maximum cluster size $B$7 and a minimum similarity threshold $B$8. The method is hierarchical in the sense of progressively coarsening row groups guided by overlap. In the reported experiments, practical defaults were $B$9 and $p$0 (Islam et al., 28 Jul 2025).

A central design choice is the decoupling of reordering from the clustered matrix format. Any reordering can be applied to $p$1, or no reordering can be applied, and clusters can then be formed using one of three schemes: fixed-length clustering, variable-length clustering, or hierarchical clustering. Fixed-length clustering uses contiguous blocks of $p$2 rows. Variable-length clustering uses contiguous blocks determined by Jaccard similarity to a representative row, subject to $p$3 and $p$4. Hierarchical clustering merges noncontiguous rows via bottom-up union-find using $p$5-derived candidates, followed by relayout (Islam et al., 28 Jul 2025).

This decoupling is conceptually important. The row-clustered format is independent of how clusters were chosen, and its benefits accrue whenever intra-cluster overlap is high, regardless of the reordering heuristic used. A plausible implication is that hierarchical clustering should be viewed less as a single monolithic preprocessing step than as one element in a modular locality pipeline (Islam et al., 28 Jul 2025).

3. Row-clustered sparse format and cluster-wise execution

To realize the locality implied by clustering, the method introduces a row-clustered CSR layout, denoted CSR$p$6, that stores nonzeros grouped by column within a cluster, enabling column-wise traversal over merged rows. The format contains a cluster pointer array, and for each cluster $p$7, a column header array listing the union

$p$8

in ascending column order. For each column $p$9, it stores a compact list of $i$0 pairs for rows $i$1 that satisfy $i$2. Variable-length CSR$i$3 uses sparse per-column lists and an additional pointer array to values; a fixed-length variant can encode a dense row mask for the cluster (Islam et al., 28 Jul 2025).

Indices are remapped so that rows are local within a cluster, $i$4, allowing tight loops over cluster rows. Metadata per cluster includes the cluster size, a mapping from row-local indices to global row IDs or an implicit mapping via original IDs and cluster sizes, and offsets to each column’s sublists. The reported space overhead stems from padding when $i$5 is large but many rows in $i$6 do not have that column; empirically it is $i$7 in $i$8 of cases, and is often smaller than standard CSR because indices are shared across rows within a column of a cluster (Islam et al., 28 Jul 2025).

The associated access pattern is column-major within each cluster. For every cluster $i$9, the algorithm traverses each column $C$0 once, loads row $C$1 into cache, and then updates every row $C$2 for which $C$3 is present. In condensed form, the multiply phase is:

• parallelize over clusters $C$4;
• for each column $C$5, load row $C$6;
• for each $C$7 in the column-sublist of $C$8, recover the global row $C$9;
• for each $B$00, update the accumulator entry for $B$01.

This changes the reuse mechanism from accidental cache residency under row-wise traversal to explicit one-read-per-cluster reuse of $B$02 when overlap is present (Islam et al., 28 Jul 2025).

The arithmetic work remains

$B$03

but memory traffic for $B$04 approaches $B$05 rather than $B$06. Preprocessing complexity depends on the clustering method. Fixed-length clustering requires $B$07 conversion. Variable-length clustering requires $B$08 for similarity scans of consecutive rows, plus $B$09 conversion. Hierarchical clustering requires one $B$10 with unit values, topK filtering, union-find merges with complexity $B$11 where $B$12, followed by $B$13 conversion (Islam et al., 28 Jul 2025).

4. Relationship to matrix reordering

The study evaluates ten reordering algorithms: Original, Random (Shuffled), Reverse Cuthill–McKee (RCM), Approximate Minimum Degree (AMD), Nested Dissection (ND), Graph Partitioning (GP via METIS, edge cut), Hypergraph Partitioning (HP via PaToH, cut-net/quality), Gray code ordering, Rabbit, Degree-based, and SlashBurn. Their objectives differ. CM/RCM target bandwidth reduction to improve locality. AMD/ND target fill-in reduction. GP/HP minimize cut and thereby implicitly maximize intra-part overlap of $B$14, which aligns well with the cut objective on the row-intersection graph. Gray, Rabbit, Degree, and SlashBurn aim to group structurally similar rows or communities and to reduce random accesses (Islam et al., 28 Jul 2025).

The empirical finding is that reordering based on graph partitioning provides better SpGEMM performance than existing alternatives at the cost of high preprocessing time. The paper also states that GP/HP perform best for SpGEMM because their optimization more directly minimizes cross-part edges weighted by overlap, which is strongly correlated with reducing duplicate reads of $B$15 across parts. However, many GP/HP instances require $B$16 the cost of one SpGEMM to amortize, whereas hierarchical clustering achieves competitive speedups at $B$17 cost in $B$18 of inputs (Islam et al., 28 Jul 2025).

This creates a cost–benefit distinction between reordering and clustering. When preprocessing can be amortized over many multiplies, GP or HP may be justified. When preprocessing budget is tighter, hierarchical clustering offers a lower-cost route to exploiting overlap. The paper’s explicit recommendation is to treat reordering and clustering independently. If a domain already provides a good ordering, fixed-length clustering may suffice. If ordering is unknown or expensive, hierarchical clustering yields high-quality clusters, and an implied reordering, in one pass. Combining reordering and clustering can be synergistic but is matrix-dependent (Islam et al., 28 Jul 2025).

A common misconception is that row reordering alone is sufficient for SpGEMM locality. The study directly argues against this: simply reordering rows of $B$19 and keeping CSR row-major traversal does not ensure that frequently used rows of $B$20 remain in cache across adjacent rows of $B$21. Cluster-wise computation is therefore not merely an implementation detail of reordering, but the mechanism that converts structural overlap into realized reuse (Islam et al., 28 Jul 2025).

5. Empirical characterization

The evaluation uses 110 square matrices from SuiteSparse, including 26 from prior SpGEMM studies and 32 from recent graph-oriented suites. Selection ensures $B$22M, $B$23B, and reduced redundancy among grouped publishers, while full SNAP and DIMACS10 are kept. Results are averaged over 10 runs on Perlmutter CPU nodes with AMD EPYC 7763, 64 cores, DDR4, $B$24 MiB/core, OpenMP with 64 threads, Intel icpc 2024.1 $B$25, and a hash-table accumulator. The workloads are $B$26 and $B$27 tall-skinny frontier matrices from CombBLAS BC (Islam et al., 28 Jul 2025).

For hierarchical clustering, the reported average speedup is $B$28, with improvements on $B$29 of matrices. The distribution shows speedups up to $B$30, with most between $B$31 and $B$32. The preprocessing cost is $B$33 the cost of a single SpGEMM on $B$34 of inputs (Islam et al., 28 Jul 2025).

For reordering alone, HP achieves geomean $B$35 with positive speedups on $B$36 of inputs, and GP and RCM are also strong. Best-per-matrix reordering achieves up to $B$37 geomean across datasets, but with high overhead (Islam et al., 28 Jul 2025).

For clustering without reordering, fixed-length and variable-length clustering speed up $B$38 and $B$39 of cases, respectively. With HP or GP preprocessing, they improve $B$40 of inputs with $B$41 geomean. On tall-skinny SpGEMM, reorderings that help $B$42 also help $B$43, and hierarchical clustering is often beneficial across frontier iterations, which the paper interprets as confirming reuse benefits independent of the specific $B$44 operand (Islam et al., 28 Jul 2025).

The paper also supplies an illustrative example. If

$B$45

then $B$46, $B$47, $B$48, and so forth. For cluster $B$49, the touched set is $B$50. In row-wise SpGEMM, $B$51 is read three times across rows $B$52, and $B$53 twice. In cluster-wise traversal, $B$54 and $B$55 are each read once per cluster, reducing $B$56 to $B$57 by $B$58–$B$59 for those rows (Islam et al., 28 Jul 2025).

6. Limitations, applicability, and related uses of hierarchy

The method is explicitly most beneficial when rows of $B$60 share substantial overlap in $B$61, but similar rows are not adjacent. Low-overlap matrices, including diagonally dominant or random sparsity patterns, yield little reuse; clustering then adds overhead and padding. Fixed-length clustering can inflate memory because of padding columns with sparse participation. Excessively large clusters can expand $B$62 beyond cache, negating locality. The paper therefore recommends choosing $B$63 and $B$64 so that the expected $B$65 fits in $B$66, and stopping merging when adding a row increases $B$67 beyond cache-, NUMA-, or bandwidth-optimal bounds, or when $B$68 (Islam et al., 28 Jul 2025).

There are also interactions with parallelism. Union-find merging can produce uneven clusters, and load imbalance may hurt parallel efficiency. Variable or heterogeneous cluster sizes require dynamic scheduling. Per-row accumulators should be private or lock-free. Symbolic computation can be fused with numeric by first-touch marker arrays per row. The format primarily improves $B$69-row reuse; hash-based kernels benefit most, though heap- and SPA-based kernels can also benefit from reduced $B$70-row read duplication (Islam et al., 28 Jul 2025).

GPU-specific constraints are not addressed in the row-clustering work; the paper states that extending CSR$B$71 to GPUs requires careful tiling. This is distinct from the sense of hierarchy in "Communication-Avoiding SpGEMM via Trident Partitioning on Hierarchical GPU Interconnects," where “hierarchy-aware partitioning/clustering” refers to grouping computation and communication to match the hardware hierarchy of local and global interconnects, rather than clustering rows by overlap. That work introduces Trident, a hierarchy-aware 2D distributed SpGEMM algorithm for modern heterogeneous supercomputers, and explicitly clarifies that this usage of clustering is not graph or data clustering (Bellavita et al., 22 Mar 2026).

A broader antecedent appears in "Rapid Near-Neighbor Interaction of High-dimensional Data via Hierarchical Clustering," which introduced a method for obtaining a matrix permutation that renders a desirable sparsity profile, guided by the principle of a block-sparse matrix with dense blocks, using lower-dimensional embedding, hierarchical data clustering, multi-level matrix compression storage, and multi-level interaction computations (Pitsianis et al., 2017). This suggests a wider lineage of locality-oriented hierarchical organization for sparse computations. However, the row-overlap formulation for SpGEMM is more specific: it derives clusters directly from shared $B$72-row access requirements and couples them to a cluster-wise access pattern and a row-clustered CSR layout (Islam et al., 28 Jul 2025).