Sparse Matrix-Matrix Product (SpMM)

• SpMM is defined as multiplying a sparse matrix with a dense or sparse matrix, serving as a core operation in graph algorithms, neural networks, and scientific computing.
• It supports generalized algebraic semirings for custom arithmetic operations, making it adaptable for shortest path algorithms and aggregations in GNNs.
• Optimized SpMM algorithms use data structures like CSR with merge-based load balancing and GPU enhancements such as Tensor Cores to improve performance.

Sparse Matrix-Matrix Product (SpMM) denotes the multiplication of a sparse matrix with another matrix, either dense or sparse, and functions as a core computational pattern in scientific computing, graph algorithms, sparsely connected neural networks, graph neural networks, clustering, and many-to-many comparisons of biological sequencing data (Buluç, 6 Aug 2025). In the sparse-dense setting, a standard formulation is $Y_{M \times N} = A_{M \times K} X_{K \times N}$, where $A$ is sparse and $X,Y$ are dense, with elementwise computation

$Y_{i,n} = \sum_k A_{i,k}\cdot X_{k,n}.$

In the broader literature, the term is used alongside closely related variants such as SpGEMM, where both operands are sparse, and masked sparse products, where only a subset of output entries is required (Huang et al., 2021, Buluç, 6 Aug 2025, Milaković et al., 2021).

1. Scope, notation, and algebraic generalization

SpMM is commonly written as $C = A \times B$ or $Y \leftarrow AX$, with a sparse left operand and a dense or sparse right operand (Xiang et al., 8 Apr 2025, Buluç, 6 Aug 2025). The survey treatment in "The Ubiquitous Sparse Matrix-Matrix Products" emphasizes that the operation is not limited to classical arithmetic over the reals: it also arises on arbitrary algebraic semirings where scalar addition and multiplication are overloaded with user-defined functions, and on more general heterogeneous algebras where the domains of the input matrices can differ (Buluç, 6 Aug 2025).

This broader viewpoint places several established operators under a single conceptual umbrella. In that treatment, SpMM refers to the case where $A$ is sparse and $X,Y$ are dense; SpGEMM refers to sparse-sparse multiplication; masked products selectively compute only entries permitted by a mask; and degenerate cases include SpMV and SpMSV (Buluç, 6 Aug 2025, Milaković et al., 2021). For masked sparse-sparse multiplication, the mask can be integrated into the multiplication rather than applied as a post-processing step, so that masked-out entries are skipped during accumulation itself (Milaković et al., 2021).

The algebraic generalization is not merely notational. It is directly tied to applications such as shortest paths under min-plus or max-plus algebra, randomized algorithms using custom scalar operations, and GNN variants where aggregation and combination are not simple arithmetic reductions (Buluç, 6 Aug 2025). This suggests that implementation generality and performance portability are often in tension: supporting arbitrary semirings or user-defined functions can complicate low-level vectorization, code generation, and the use of specialized hardware.

2. Data structures and algorithmic design principles

A persistent theme in SpMM research is that sparsity structure and dense-matrix width jointly determine the effective kernel design. CSR remains the dominant baseline because it is the most common format in software stacks and avoids expensive preprocessing or format conversion in many workloads (Yang et al., 2018, Huang et al., 2020). GE-SpMM explicitly targeted CSR to remain compatible with GNN frameworks and to support SpMM-like operations without preprocessing overheads (Huang et al., 2020).

Early GPU design principles emphasized two ingredients: merge-based load-balancing and row-major coalesced memory access (Yang et al., 2018). The merge-based formulation distributes work by balanced nonzero ranges rather than by whole rows, which mitigates the severe imbalance caused by irregular row lengths; row-major coalesced access makes accesses to both the dense input and dense output more efficient on GPUs (Yang et al., 2018). That work also argued that SpMM differs materially from SpMV because latency hiding via instruction-level parallelism and the dense output dimension become central (Yang et al., 2018).

Subsequent GPU kernels refined this picture. "Efficient Sparse Matrix Kernels based on Adaptive Workload-Balancing and Parallel-Reduction" introduced Vectorized Segment Reduction, Vector-type Dense-row Loading, and Coalesced Sparse-row Caching, and showed that the best performing SpMM kernel depends on both the sparsity pattern of $A$ and the width $N$ of dense $A$0 (Huang et al., 2021). In that study, parallel-reduction is preferred for $A$1, while sequential reduction is preferred for $A$2, and row-length statistics determine when workload-balancing should be enabled (Huang et al., 2021).

GE-SpMM addressed the GNN case with Coalesced Row Caching and Coarse-grained Warp Merging, both aimed at reusing sparse-row data among threads while maintaining coalesced global-memory access (Huang et al., 2020). ParamSpMM generalized this adaptivity with the Parameterized Compressed Sparse Row (PCSR) format, whose parameters encode blocking, workload balancing, thread coarsening, and thread block size, and paired it with an ML-based SpMM-decider that predicts configurations from graph size, degree distribution, and data-locality features (Zhang et al., 15 May 2026). Its reported average speedup over cuSPARSE is 1.92x, and the decider achieves more than 99% of hand-tuned performance on test graphs (Zhang et al., 15 May 2026).

For sparse-sparse and masked settings, the data-structure question extends to accumulation. Masked SpGEMM work classified push-based and pull-based algorithms and developed Masked Sparse Accumulator, Hash Accumulator, Masked Compressed Accumulator, and Heap Accumulator, with performance depending strongly on matrix density, mask density, mask structure, and cache behavior (Milaković et al., 2021).

3. GPU kernels: CUDA cores, Tensor Cores, and hybrid execution

Modern GPU SpMM research has increasingly focused on matrix engines such as Tensor Cores, but the central difficulty is that these units are optimized for dense, regular matrix products rather than unstructured sparse operands. The abstract of cuTeSpMM notes that, for a substantial majority of the sparse matrices in the SuiteSparse collection, the performance of TC-GNN falls significantly short of state-of-the-art SpMM kernels that only utilize scalar cores (Xiang et al., 8 Apr 2025). Its response was to define TCU-Synergy, based on the non-zero structure and a modeled Operational Intensity, and to show that matrices with high TCU-synergy can outperform scalar-core SpMM implementations, while low-TCU-Synergy cases incur only slightly lower performance (Xiang et al., 8 Apr 2025).

FlashSparse attacked the same problem by reducing sparse granularity through a swap-and-transpose strategy using

$A$3

thereby turning the effective granularity from $A$4 to $A$5 for Tensor Core MMA and reducing computation redundancy (Shi et al., 2024). It combined this with a memory-efficient thread mapping strategy and the ME-BCRS storage format, reporting a geometric mean 5.5x speedup over DTC-SpMM and 3.22x over RoDe on H100 and RTX 4090 GPUs (Shi et al., 2024). Acc-SpMM pursued a broader Tensor-Core library design through data-affinity-based reordering, the BitTCF compressed format, a high-throughput pipeline, and adaptive sparsity-aware load balancing, with average speedups of 2.52x on RTX 4090, 1.91x on A800, and 1.58x on H100 over cuSPARSE (Zhao et al., 16 Jan 2025).

A more recent line of work uses hybrid or asynchronous GPU execution. AsyncSparse is explicitly co-designed for NVIDIA asynchronous features such as Tensor Memory Accelerator and warp specialization; it provides a BCSR kernel for structured sparsity and a WCSR kernel for irregular sparsity (Liu et al., 20 Apr 2026). On SuiteSparse matrices, the WCSR kernel outperforms prior SpMM kernels by 1.47x over AccSpMM and 6.24x over cuSPARSE, while the BCSR kernel yields a combined 2.66x end-to-end speedup on Qwen2.5-7B