---
title: Sparse-Dense Matrix Multiplication
url: https://www.emergentmind.com/topics/sparse-dense-matrix-multiplication-spmm
type: topic
---

# Sparse-Dense Matrix Multiplication

Sparse-Dense Matrix Multiplication (SpMM) refers to the operation $C = A \cdot B$, where $A \in \mathbb{R}^{m \times k}$ is a sparse matrix (typically with $\text{nnz}(A) \ll mk$) and $B \in \mathbb{R}^{k \times n}$ is a dense matrix, producing a dense or semi-dense output $C \in \mathbb{R}^{m \times n}$. SpMM is a central primitive in scientific computing, graph analytics, deep learning, and numerous algorithmic domains, with extensive architecturally specialized and software-optimized implementations across CPUs, GPUs, FPGAs, and custom accelerators [2508.04077][2312.12766][2007.03179][2504.06443][2412.11007][2109.11081][2312.05639].

## 1. Mathematical Formulation and Computational Complexity

Let $A$ be sparse and $B$ dense. The canonical row-wise formula is:
\[
C_{i,j} = \sum_{\ell=1}^k A_{i,\ell} \cdot B_{\ell,j}
\]
Usually, only nonzero entries of $A$ are processed, so:
\[
C_{i,j} = \sum_{\ell \in \text{nz}(A_i)} A_{i,\ell} \cdot B_{\ell,j}
\]
where $\text{nz}(A_i)$ is the set of columns in row $i$ with nonzeros.

The total floating-point operation count is $2 \cdot \text{nnz}(A) \cdot n$ (each nonzero incurs $n$ multiplies and $n$ adds), with memory traffic governed by reading $\text{nnz}(A)$ values and indices, $k\times n$ dense elements from $B$ (possibly streamed/tiled), and writing $m\times n$ outputs [2508.04077][2412.11007][2312.12766][2007.03179].

Generalization to user-defined semirings is routine in many algebraic graph kernels, i.e.,
\[
C_{i,j} = \bigoplus_{\ell \in \text{nz}(A_i)} A_{i,\ell} \otimes B_{\ell,j}
\]
where $(S, \oplus, \otimes)$ satisfies monoid and distributive properties [2508.04077].

## 2. Storage Formats and Dataflow Implications

Most SpMM kernels assume $A$ is stored in Compressed Sparse Row (CSR) or Blocked extensions (e.g., BCSR, GCOO, ME-BCRS), while $B$ and $C$ are dense row-major. 

Common CSR representation:
- `rowPtr[0..m]`: index into nonzero lists per row,
- `colInd[0..nnz-1]`: column indices for nonzeros,
- `val[0..nnz-1]`: values.

Advanced formats optimize memory hierarchy and core alignment:
- Grouped COO (GCOO): partitions columns into coarse groups for shared-memory reuse [2005.14469].
- Blocked formats (ME-BCRS, BCSR): exploit 2D tile regularity to maximize contiguous copy and hardware acceleration [2412.11007][2412.08902][2503.01253].

Choice of sparse format strongly affects coalesced memory access, vectorization, and the mapping to accelerator MMA units (e.g., for GPUs and SME on Arm) [2412.11007][2504.06443][2005.14469][2511.08158].

## 3. Algorithmic Strategies on Modern Architectures

### 3.1 GPU/TPU-Tailored SpMM

- **CUDA Scalar Kernel:** Assigns one thread or warp per row, exploiting output (dense) parallelism [2007.03179][2412.08902]. To mitigate load imbalance and uncoalesced access, warp merging, coalesced row caching, or batched schemes are deployed [1903.11409].

- **Tensor/Core Unit (TCU) Exploitation:** 
  - Dense zero-filling: Tiles of sparse $A$ are zero-padded to TCU shape (e.g., $16\times16$), then processed via MMA [2504.06443].
  - Hybrid approaches (e.g., cuTeSpMM): Use a "TCU-synergy" metric (tile fill rate) to decide when to invoke TCU vs. CUDA core kernels [2504.06443][2603.08734].
  - Fine-grained methods (e.g., FlashSparse) minimize redundancy by swapping operands ($A B = (B^T A^T)^T$), aligning sparse granularity to the narrow TCU dimension, and removing superfluous zero padding [2412.11007].
  - Specialized schemes (RSH-SpMM): Employ adaptive row partitioning and row-structured tiling to balance dense-tile formation against residual irregular rows routed to CUDA cores [2603.08734].

### 3.2 CPU/SME/Other Accelerators

- **SIMD-Aware JIT and Autotuning:** Systems like JITSPMM generate architecture and matrix-shape-aware inner loops at runtime, aggressively unrolling vectorized accumulations, optimizing register allocation, and selecting between row/nnz/merge partitioning for multithreading [2312.05639].
- **Hybrid SME/NEON** (Armv9): LOOPS partitions A into CSR rows for NEON AXPY and narrow-vector BCSR tiles for SME's outer product units (fmopa), scheduling via a lightweight two-level performance model [2511.08158].
- **General-purpose Hardware:** FusedMM fuses SpMM with related kernels (e.g., SDDMM), exposing user-defined accumulation, load-balancing, and cache-blocked register allocation for best SIMD efficiency [2011.06391]. 

### 3.3 Custom/Flexible Accelerators

- **Streaming Dataflow (Sextans):** Implements II=1 (initiation interval) pipelines for sparse-tuple stream, N-dense memory windows, and on-chip accumulation with pointer-based task decomposition [2109.11081].
- **IOPS and Unified Inner/Outer Product:** Fuses inner-product (maximal output locality) and outer-product (maximal zero skipping) in a PE mesh, with on-chip bookkeeping of partial sums and adaptive tiling based on buffer budgets and input sparsity [2312.12766].
- **Distributed RDMA/SHMEM Contexts:** SpMM is tiled over process grids, often using asynchrony and dynamic replication to minimize communication, with communication-eliding fusion for SDDMM→SpMM chains [2203.07673][2311.18141].

### 3.4 Batched and Semi-external Techniques

- **Batched SpMM:** Aggregates many small SpMMs (e.g., multiple GCN minibatches), assigning warps/subwarps for high GPU occupancy and shared-memory staging, critical in small-graph scenarios [1903.11409].
- **Semi-external SpMM:** For out-of-core large-scale problems, the sparse matrix resides on SSDs, while B is tiling-streamed through memory, orchestrated via asynchronous threads and task queues, achieving near in-memory performance [1602.02864].

## 4. Performance Models and Bottleneck Analysis

A central analytic rubric is the "roofline model", bounding throughput by the minimum of compute peak and memory bandwidth, scaled by operational intensity (FLOPs / bytes moved). For naïve CSR SpMM:
\[
\text{OI} = \frac{2 ~ \text{nnz}(A) \cdot n}{[\text{nnz}(A) (n+1) + m n ]}
\]
Operational intensity increases by blocking, shared-memory reuse, or aggregating multiply-accumulates across repeated column or tile indices [2005.14469][2503.01253][2312.12766].

Metrics such as "TCU-synergy" (average tile density) [2504.06443], block physical fill (in block formats), and arithmetic-to-I/O ratios determine the practical kernel and hardware mapping.

At moderate-to-high sparsity, SpMM is typically memory bandwidth–bound; as sparsity increases further, compute-to-load ratio drops and the kernel can become memory-bound even with highly tuned implementation [2412.11007][2503.01253]. On accelerator arrays, communication costs can dominate, requiring topology-aware distribution and fusion [2203.07673].

## 5. Architecture-Specific Optimization Strategies

| Hardware    | Key Strategies                                                      | Representative Papers       |
|-------------|---------------------------------------------------------------------|----------------------------|
| CPUs        | SIMD blocking, JIT, cache tiling, balanced partition, semiring fusion| [2312.05639][2011.06391][2508.04077] |
| NVIDIA GPUs | Memory-coalesced loads, warp-merge, TCU tile-packing, hybrid cores  | [2412.08902][2412.11007][2504.06443][2007.03179][2603.08734] |
| Arm SME     | Partitioned CSR/BCSR, hybrid NEON/SME scheduler, model-guided split | [2511.08158]                |
| FPGAs       | Streaming II=1, pointer-based tile queuing, HBM pipelining          | [2109.11081]                |
| CS-3        | SELLPACK multi-channel streaming, chunked I/O/PE sync               | [2604.27985]                |
| Distributed | 1.5D/2.5D dense/sparse shifting, communication-eliding, FusedMM     | [2203.07673][2311.18141]    |

*Editor’s term: This table succinctly organizes the major architecture-specific strategies for SpMM using only information from referenced data.*

## 6. Applications and Case Studies

SpMM underpins:
- Graph neural network (GNN) operations (e.g., message-passing, pooling, feature propagation) [2007.03179][2412.08902][2011.06391][2109.11081][2412.11007]
- Scientific simulations (FEM/CFD), eigensolvers (LOBPCG, Arnoldi), low-rank matrix factorization, clustering, and recommendation algorithms [2508.04077][1602.02864]
- Block-sparse attention (transformers), CTQFT tensor networks, collaborative filtering, PCA/NMF, and sampling-based randomized linear algebra [2508.04077][2412.11007]
- Billion-node graph analytics, PageRank, block-Krylov eigensolving, NMF at semi-external scale [1602.02864]

Typical matrix shapes and sparsity regimes vary from extremely tall/skinny or power-law-structured graphs (bioinformatics, social networks) to blocky/supernode matrices in quantum simulations.

## 7. Directions, Open Challenges, and Best Practices

While the last five years have seen major advances, several avenues remain at the forefront:

- **Irregularity Mitigation:** Designs like RSH-SpMM and HC-SpMM highlight data-dependent adaptive partitioning, cross-kernel fusion, and row-structured reordering as essential to stable high throughput under real-world sparsity distributions [2603.08734][2412.08902].
- **Format Auto-Tuning:** ME-BCRS, GCOO, and hybrid CSR/BCSR partition require empirical profiling or online adaptation to maximize tile density and SIMD/TCU utilization per matrix [2412.11007][2511.08158][2005.14469].
- **Pipeline and Prefetch Tuning:** Fine-grained pipelining (double buffering, staged prefetches) is critical to saturate DRAM and on-chip memory bandwidth, especially on GPUs/SME [2503.01253][2412.11007][2603.08734].
- **Extending Beyond SpMM:** FusedMM and IOHP kernels extend ideas across SDDMM, sparse-sparse MM, and higher-order contractions, exploiting similar dataflow, accumulator, and reduction logic [2312.12766][2011.06391][2412.11007].
- **Energy Efficiency and Cross-Platform Comparison:** Recent works report up to $10$–$70\times$ GFLOPS/W edge for SME/CPU frameworks versus GPUs for selected workloads [2511.08158].
- **Out-of-Core and Distributed SpMM:** Semi-external strategies, distributed-memory RDMA kernels, and adaptive communication-eliding synthesis are essential for scaling to billion-edge graphs and multi-node ML systems [1602.02864][2311.18141][2203.07673].

**Best practice** (across hardware): orchestrate blocking, coalesced memory access, SIMD/TCU/Tile mapping, and adapted row/column partitioning guided by operational intensity, bandwidth, and (increasingly) data-dependent dynamic scheduling for practical, robust SpMM performance.

---

**Key references:**
- [2508.04077] "The Ubiquitous Sparse Matrix-Matrix Products"
- [2412.11007] "FlashSparse: Minimizing Computation Redundancy for Fast Sparse Matrix Multiplications on Tensor Cores"
- [2312.12766] "IOPS: An Unified SpMM Accelerator Based on Inner-Outer-Hybrid Product"
- [2603.08734] "RSH-SpMM: A Row-Structured Hybrid Kernel for Sparse Matrix-Matrix Multiplication on GPUs"
- [2511.08158] "LOw-cOst yet High-Performant Sparse Matrix-Matrix Multiplication on Arm SME Architectures"
- [2504.06443] "cuTeSpMM: Accelerating Sparse-Dense Matrix Multiplication using GPU Tensor Cores"
- [2007.03179] "GE-SpMM: General-purpose Sparse Matrix-Matrix Multiplication on GPUs for Graph Neural Networks"
- [2109.11081] "Sextans: A Streaming Accelerator for General-Purpose Sparse-Matrix Dense-Matrix Multiplication"
- [1602.02864] "Semi-External Memory Sparse Matrix Multiplication for Billion-Node Graphs"
- [2312.05639] "JITSPMM: Just-in-Time Instruction Generation for Accelerated Sparse Matrix-Matrix Multiplication"
- [2011.06391] "FusedMM: A Unified SDDMM-SpMM Kernel for Graph Embedding and Graph Neural Networks"
- [2412.08902] "HC-SpMM: Accelerating Sparse Matrix-Matrix Multiplication for Graphs with Hybrid GPU Cores"
- [1903.11409] "Batched Sparse Matrix Multiplication for Accelerating Graph Convolutional Networks"
- [2203.07673] "Distributed-Memory Sparse Kernels for Machine Learning"
- [2311.18141] "RDMA-Based Algorithms for Sparse Matrix Multiplication on GPUs"
- [2503.01253] "NM-SpMM: Accelerating Matrix Multiplication Using N:M Sparsity with GPGPU"
- [2005.14469] "Efficient Sparse-Dense Matrix-Matrix Multiplication on GPUs Using the Customized Sparse Storage Format"
- [2604.27985] "Exploring Sparse Matrix Multiplication Kernels on the Cerebras CS-3"

Source: https://www.emergentmind.com/topics/sparse-dense-matrix-multiplication-spmm