Sparsity-Aware Roofline Models for Sparse Matrix-Matrix Multiplication Matthew Qian∗ , Yahia Ramadan† , Suhita Anubha‡ , Ariful Azad§ ∗ Texas A&M University, TX, USA ([email protected]) † Texas A&M University, TX, USA ([email protected]) ‡ Indiana University, Bloomington, IN, USA ([email protected])
arXiv:2604.06637v1 [cs.DC] 8 Apr 2026
§ Texas A&M University, TX, USA ([email protected])
Abstract—Sparse matrix–dense matrix multiplication (SpMM) is a critical kernel in scientific computing, graph analytics, and machine learning, whose performance is often constrained by memory bandwidth. In this work, we investigate the applicability and limitations of roofline modeling for SpMM by explicitly accounting for the impact of matrix sparsity structure on arithmetic intensity and attainable performance. We evaluate three SpMM implementations: Compressed Sparse Row (CSR), Compressed Sparse Blocks (CSB), and Intel’s Math Kernel Library (MKL). Each implementation was tested using largescale matrices from the SuiteSparse collection and grouped by sparsity pattern, including block-structured, banded (diagonal), scale-free, and uniformly random matrices. We derive sparsityaware roofline models that incorporate memory traffic, cache locality, and blocking behavior, and demonstrate that a single model is insufficient to accurately predict performance across diverse structures. Experiments were conducted on an AMDbased Perlmutter compute node with a varying number of columns in the dense matrix. In particular, blocking and structured sparsity significantly alter effective arithmetic intensity. The results show that accurate roofline-based performance analysis of SpMM requires sparsity-aware modeling, and that data layout and blocking strategies must be evaluated in the context of matrix structure rather than through a single unified model.
I. I NTRODUCTION Sparse matrix–dense matrix multiplication (SpMM) is a fundamental operation in sparse linear algebra [1], [2]. It plays a central role in Graph Neural Networks (GNNs), supporting both forward and backward propagation during training and inference [3]. Beyond GNNs, SpMM underpins a wide range of graph analytics, including graph embedding [4], graph visualization [5], and batched PageRank computations [6]. In scientific computing, SpMM is widely used in modal analysis for the Finite Element Method (FEM) [7] and electronic structure simulations based on Density Functional Theory (DFT) [8]. More recently, SpMM has emerged as a key primitive in traditionally dense domains, such as sparsified Transformers for large language models (LLMs) [9] and sparsified convolutional filters in deep neural networks [10]. Given the broad and growing reliance on SpMM, accurate performance models are needed to predict its attainable performance on modern multicore architectures. The Roofline model [11] is one of the most widely used frameworks for performance analysis that relates achievable performance to arithmetic intensity, defined as the ratio of floating-point operations to memory traffic. Operations with low arithmetic
intensity are classified as memory-bound, with performance primarily limited by memory bandwidth. While the Roofline model has proven effective for dense computations and simpler sparse kernels such as sparse matrix–vector multiplication (SpMV), its application to SpMM has received comparatively little attention. SpMM is widely regarded as a memory-bound operation whose performance is strongly influenced by the sparsity structure of the matrix. For example, matrices with block or clustered sparsity often achieve significantly higher performance than matrices with uniformly random nonzeros. Consequently, a meaningful Roofline analysis of SpMM must explicitly account for sparsity structure. However, most existing Roofline-based models for SpMM either ignore sparsity or are tailored to specific application contexts, such as nuclear reactor simulations [12] or sparse attention mechanisms in LLMs [13], and thus fail to capture the wide performance variability observed across different sparsity patterns. In this paper, we argue that no single Roofline model can adequately characterize SpMM performance across diverse sparse matrices. Instead, we propose a set of sparsity-aware Roofline models that reflect distinct structural regimes. We develop four models by deriving arithmetic-intensity bounds for random sparsity, block-diagonal sparsity, general blockstructured sparsity, and scale-free sparsity. We show that random sparsity represents a worst-case scenario, providing a lower bound on SpMM performance, while structured sparsity enables significantly higher attainable performance. We validate our models using four classes of matrices corresponding to each sparsity structure and evaluate three SpMM implementations based on compressed sparse row (CSR), compressed sparse block (CSB) [14], and the Intel MKL library. Experimental results on an AMD multicore processor demonstrate that sparsity-aware Roofline models closely track observed performance. These results confirm that incorporating sparsity structure is essential for accurately predicting SpMM performance. Overall, the paper makes the following contributions. 1) We introduce sparsity-aware Roofline models for SpMM that explicitly capture the impact of sparsity structure. We derive the arithmetic intensity and modeling for random, diagonal, scale-free, and block-structured matrices.
TABLE I: Notations used in the paper. Symbol
Description
A B C n d nnz t D z
Sparse matrix of size n × n Dense matrix of size n × d (tall-and-skinny) Output matrix of size n × d (tall-and-skinny) Number of rows/columns in A Number of columns in B (and C) Number of nonzeros in A Block dimension for A Average number of nonzeros in a block of A Average number of nonempty columns in a block of A
FLOPs AI β
Total floating point operations in SpMM (≈ 2 · nnz · d) Arithmetic intensity (AI = FLOPs/Bytes) Peak memory bandwidth (122.6 GB/s)
2) We validate the proposed models by showing that sparsity-aware Roofline predictions closely match measured performance across CSR, CSB, and MKL SpMM implementations. II. BACKGROUND AND R ELATED W ORK SpMM multiplies a sparse matrix A ∈ Rn×n with a dense matrix B ∈ Rn×d to produce a dense matrix C ∈ Rn×d . Table I outlines the notation used throughout this paper. The performance of SpMM depends on multiple factors, including matrix dimensions, sparsity structure, data layout, implementation strategy, and the underlying hardware architecture. In this section, we discuss these factors to provide the necessary context and motivate the design of our Roofline performance models. A. Matrix Characteristics The sizes of matrices affect cache usage and thereby impact the model and performance. Table II shows five applications of SpMM in scientific computing and machine learning. All of these applications use a square sparse matrix and a talland-skinny dense matrix where d≪n. In addition, the distribution of nonzeros strongly influences SpMM performance. Random sparsity leads to irregular memory accesses and poor cache reuse, while banded structures improve locality by reusing portions of B across multiple multiplications. Higher average density increases arithmetic intensity and increases theoretical performance, whereas very sparse rows or columns increase bandwidth pressure. Modeling Assumptions. Based on the above discussion, we restrict our analysis to tall-and-skinny dense matrices. To account for the impact of sparsity structure on performance, we develop models for random, banded, block, and scale-free sparsity patterns. B. Algorithm and Data Layout Matrix format significantly affects SpMM performance. Formats like CSR, CSC, and CSB differ in how they store and access nonzeros, impacting cache locality and memory traffic. Blocking or tiling strategies, as in CSB, enhance cache reuse by working on submatrices, improving performance
on modern memory hierarchies. Additionally, the choice of multiplication algorithms such as column-by-column, inner product, or outer product affects data access patterns and parallelism, influencing efficiency and scalability. Data Layout Options. We consider SpMM implementations based on the CSR, CSC, and CSB data structures to capture commonly used sparse formats and enable fair and representative performance modeling and comparisons. C. The Roofline model The Roofline model [11] is a visually intuitive performance framework that relates a kernel’s attainable performance to the architectural limits of a computing platform. By plotting performance (GFLOP/s) against arithmetic intensity (AI), the model identifies whether a kernel is “memory-bound,” limited by the hardware’s peak bandwidth, or “compute-bound,” limited by its peak computational throughput. More precisely, arithmetic intensity is defined as the ratio of the total number of floating-point operations performed by a kernel to the total number of bytes of data transferred to and from the main memory during its execution (AI = FLOPs Bytes ). This metric quantifies the “computational density” of an algorithm; a higher AI suggests greater data reuse within local caches, allowing a kernel to potentially reach the hardware’s maximum processing speed rather than being throttled by the slower memory subsystem. The two bounding sections are defined by P = min(β · AI, π), where π is the peak compute throughput, and β is the peak memory bandwidth. For SpMM, the kernel is predominantly memory-bound. This is because each nonzero in the sparse matrix contributes relatively few floating point operations (on the order of 2d FLOPs) while requiring multiple memory accesses. As a result, the arithmetic intensity remains low, placing SpMM in the memory-bound region of the roofline model. In this regime, performance scales linearly with arithmetic intensity and is therefore bounded by the sloped line P = β · AI. Only when the arithmetic intensity exceeds the ridge point (AI > π/β) does the kernel become compute-bound, at which point performance is capped by π. D. Related Work To better reflect modern architectures, several extensions incorporate hierarchical bandwidth limits across cache levels, yielding cache-aware or hierarchical Roofline models that improve diagnostic power for locality-sensitive kernels [15]. Developing a Roofline model that captures the observed performance of sparse operations is more challenging than for dense kernels. This is because the standard Roofline model does not account for locality, sparsity structure, or indirect memory accesses. Prior work has used Rooflinestyle reasoning to analyze sparse linear algebra kernels such as SpMV [16], SpGEMM [17], and FusedMM [18], often concluding that irregular access patterns and limited locality reduce effective bandwidth well below hardware peaks [12].
TABLE II: Common matrix shapes and sparsity patterns arising in SpMM applications across machine learning and scientific computing. Applications
Type
A Shape
Common Sparsity
B Type
Shape
GNN and graph embedding Finite Element Method (FEM) Density Functional Theory (DFT) Transformer Sparsification Recommendation Systems
adjacency matrix stiffness/mass matrices Hamiltonian sparse attention map user-item matrix
square square square square rectangular
scale-free banded, mesh-local banded top-k/row, banded irregular
node features eigenvectors Kohn–Sham orbitals dense value matrix item embeddings
tall-and-skinny tall-and-skinny tall-and-skinny tall-and-skinny tall-and-skinny
SpMM occupies an intermediate position between SpMV and dense GEMM: while it inherits indirect accesses from the sparse operand, it can exploit reuse through the dense matrix, making its performance highly sensitive to matrix structure and dense-block shape. Several application-driven studies analyze SpMM and related kernels using Roofline-style models, both in scientific computing [12] and emerging machine-learning workloads such as sparse attention [13]. Optimized SpMM is a widely studied problem, with algorithms developed for multicore processors [19], GPUs [3], [20], [21], and distributed systems [22]–[24]. Works such as WISE [25] and DDB [26] attempt to optimize sparse computations by using machine learning techniques to predict the performance of various configurations. However, the performance of these algorithms is rarely modeled analytically, making it difficult to assess whether the observed performance approaches the peak achievable performance of the target platform. III. S PARSITY-AWARE ROOFLINE MODELS FOR S P MM To develop sparsity-aware roofline models, we consider n× n square sparse matrix and n×d tall-and-skinny dense matrix. Since we only perform a multiply-accumulate for each nonzero element in the sparse matrix A for every column (d) in the dense matrix B: FLOP = 2d · nnz.
(1)
The factor of 2 accounts for one multiplication and one addition per non-zero. To estimate memory traffic, we assume that all matrix values are stored in double-precision floating-point format, while indices in the sparse matrix are stored as 32-bit integers. The total memory traffic is computed as the sum of the costs of reading A, reading B, and writing C. For A in CSR format, we must read the values, column indices, and row pointers. Thus, total data movement for A is: TrafficA = nnz · sizeof(val) + nnz · sizeof(col idx) + (n + 1) · sizeof(row ptr) ≈ 12nnz bytes Since C is dense, it is typically written once, giving TrafficC = nd · sizeof(val) = 8nd bytes. The remaining question, therefore, is how B is accessed under different sparsity patterns of the sparse matrix. In the following, we
consider four representative sparsity patterns and develop a corresponding access model for B in each case. A. Random Sparsity Under a random sparsity pattern, the column indices of A are assumed to be uniformly distributed. If the dense matrix B does not fit in cache, each nonzero of A incurs a memory access to load the corresponding row of B from main memory, resulting in no effective data reuse for B. Hence, in the worst case, we access nnz rows of B for the entire multiplication. Thus, total data movement for B in the random sparsity model is: TrafficB = nnz · d · sizeof(val) = 8d · nnz bytes Using the previously defined FLOP count, and TrafficA and TrafficC , we can derive the arithmetic intensity for the random sparsity: AI (Random) =
2d · nnz (12 + 8d) · nnz + 8nd
(2)
B. Diagonal Sparsity For Diagonal sparsity, we assume that the nonzeros of A are concentrated along the main diagonal. More concretely, we assume that the number of nonzeros per row is small and bounded, and that the nonzeros have strong spatial locality between consecutive rows. This assumption allows access to B to have high temporal locality, such that rows that are accessed once in B are accessed repeatedly before being evicted from cache. Thus, with an optimal data structure, the memory traffic for B is amortized across all nonzeros in the diagonal region and does not scale with the number of nonzeros in A. Consequently, the dominant memory traffic for this diagonal case is dependent on reading from A, writing to C, and a single load of B into cache, so we have: AI (Diagonal) =
2d · nnz 12 · nnz + 16nd
(3)
C. Blocked Sparsity Assume A ∈ Rn×n is divided into t × t blocks, with a total of nnz nonzero elements and N total nonzero blocks. Let D = nnz/N be the average number of nonzero entries per nonzero block. We assume that it takes 8 bytes to read the values, and 4 bytes to read an index for values within a block. Since we
access each value in A once, the total memory traffic for A is simply 12 nnz. We also assume that we write 8 bytes to the dense matrix C, so the memory traffic for C is 8nd. To model the memory access for B, we assume that each t × t block of A must access B once for each of its t columns that contain a nonzero. To estimate the average number of such columns (denoted z) per block, we can assume that nonzeros within a single block are distributed randomly among its t columns. We define an indicator random variable Ii for each column i within a block. ( Ii =
1, 0,
if col i contains at least one nonzero if col i is empty
TABLE III: Sparse matrices used for SpMM performance evaluation. er 22 10 It is generated using the Erdős–Rényi model with 222 rows and columns and an average of 10 nonzeros per row. Pattern
Matrix Name
Rows
Cols
Nonzeros
Blocking
road usa hugebubbles-00010 asia osm 333SP
23,947,347 19,458,087 11,950,757 3,712,815
23,947,347 19,458,087 11,950,757 3,712,815
57,708,624 58,359,528 25,423,206 22,217,266
Scale-free
com-Orkut com-LiveJournal uk-2002
3,072,441 3,997,962 18,520,486
3,072,441 3,997,962 18,520,486
234,370,166 69,362,378 298,113,762
Diagonal
rajat31 ideal diagonal 22
4,690,002 4,194,304
4,690,002 4,194,304
20,316,253 4,194,304
Random
er 22 1 er 22 10 er 22 20
4,194,304 4,194,304 4,194,304
4,194,304 4,194,304 4,194,304
4,194,304 41,942,990 83,885,880
Thus, the expected number of nonzeros in a block is given by E[z] =
t X
TABLE IV: Test system CPU specifications (Perlmutter) E[Ii ] = t · Pr(col i is nonempty)
i=0
A column i is empty if all D nonzeros are outside of i. For a single nonzero, Pr(miss col i) = 1− 1t . So for D independent nonzeros, we have 1 Pr(col i is nonempty) = 1−Pr(col i is empty) = 1−(1− )D , t When t (block size) and D (number of nonzeros in a block) are large, we can use the following approximation: 1 (1 − )D ≈ e−D/t , t because for large t and D, the binomial distribution converges to a Poisson distribution. Finally, we have z ≈ E[z] = t(1−e−D/t ) occupied columns per block [27]. A row of B is accessed for each block z times, and thus the memory traffic for B is given by the expression 8dN z. Due to the cache tiling that the blocked structure provides, in practice, repeatedly accessing nearby rows or columns of B allows for effective cache reuse rather than main memory access. To roughly account for this reuse, we scale the memory traffic from B by a factor of 14 . This heuristic reflects that due to tiling, only part of a single theoretical memory access actually results in a main memory access. The exact fraction depends on hardware, cache behavior, and implementation, so we choose 14 as an estimate based on observed experimental results. Including this factor in our full equation gives us:
Property
Value
Architecture CPU model Sockets Cores / Threads L1 cache L2 cache L3 cache Vector extensions Memory NUMA domains
x86 64 AMD EPYC 7763 (Milan) 2 (testing restricted to 1 socket) 64 cores / 64 threads (used) 64 KiB per core (32 KiB data, 32 KiB instruction) 512 KiB per core 256 MiB per socket AVX2, FMA 512 GB DDR4 total 4 per socket (NPS=4)
a small number of high-degree hub nodes can account for a disproportionately large fraction of nonzeros. Let nhub be the number of hub nodes and nnzhub be nonzeros (edges) associated with hub nodes. In the appendix, we derived the following estimation for nnzhub : nnzhub = nnz · f (α−2)/(α−1) ,
where f is the highest fraction of nodes by degree. In our experiments, we set the hub fraction to f =0.1% of the nodes. In our model, we assume that the rows of B associated with hub nodes can be retained in fast memory (e.g., cache) and therefore are not repeatedly loaded, unlike in the random sparsity model. For the remaining non-hub entries of the sparse matrix, accesses to B are assumed to be random. Under these assumptions, the arithmetic intensity for scale-free sparsity is defined as: AI =
AI (Blocked) =
2d · nnz 8 nnz + 2dN z + 8nd
(4)
D. Scale-free Sparsity For scale-free graphs, we assume that the degree distribution p(k) of a graph follows a power law distributation p(k) ∝ k −α , where α is the power-law exponent with 2 < α < 3 for real-world networks [28]. For scale-free graphs,
(5)
2d · nnz 12nnz + 8d · (nnz − nnzhub ) + 8d · nhub + 8nd
(6)
IV. E VALUATING S P MM I MPLEMENTATIONS USING S PARSITY-AWARE ROOFLINE M ODELS A. Dataset For data collection, matrices were collected and categorized according to four structures: blocking, scale-free, diagonal, and uniform random. All matrices were obtained from the
TABLE V: SpMM performance (GFLOP/s) for multiple k values across different formats. Pattern
Blocking
d=1
Matrix
d=4
d = 16
d = 64
CSR
MKL
CSB
CSR
MKL
CSB
CSR
MKL
CSB
CSR
MKL
CSB
road usa 9.468 hugebubbles-00010 5.875 asia osm 7.301 333SP 5.284
11.0924 7.146 10.078 8.692
14.240 9.696 10.668 13.057
17.528 14.358 20.455 12.258
17.289 13.490 21.481 23.625
25.423 32.768 15.853 21.743 14.027 33.975 24.875 28.784
32.652 22.975 34.568 28.893
36.234 28.322 35.093 35.227
41.316 21.743 38.345 29.729
38.567 22.975 38.450 30.106
43.006 28.322 33.479 39.596
Scale-free
com-Orkut com-LiveJournal uk-2002
8.402 18.340 26.894 14.505 30.560 38.501 21.037 29.053 34.403 12.256 22.460 32.017 11.536 15.010 26.984 35.687 44.851 72.008 66.266 76.981 92.091 41.683 53.544 61.322 16.701 24.139 16.204 55.851 78.538 67.526 146.583 167.960 148.299 226.757 205.945 164.359
Diagonal
ideal diagonal 22 rajat31
1.988 7.266
1.167 9.565
5.886 23.546 10.558 6.840 8.5888 9.390 26.944 29.348 22.601 56.978
9.039 59.644
14.202 39.275
10.902 74.064
11.023 69.266
17.294 53.911
1.586 6.194 8.091
1.634 7.833 10.906
3.998 4.957 5.446 6.226 12.832 13.921 15.225 12.373 16.283 14.979 16.249 15.453
8.194 12.374 14.169
10.216 13.456 13.483
8.547 10.0322 11.564
5.320 11.185 10.429
11.509 17.036 17.001
Uniform Random er 22 1 er 22 10 er 22 20
SuiteSparse Matrix Collection [29], except for the uniform random matrices, which were generated using the Erdős–Rényi model, and ideal_diagonal_22, which represents an ideal diagonal sparsity pattern. Table III lists the details of each matrix used. All matrices were selected to exceed the capacity of on-chip caches. This choice allows us to isolate the main weakness of the memory-bound operation by analyzing cache utilization across different implementations.
7.841 12.284 13.575
CSR
MKL
CSB
80
12 10
60
8
40
6 4
20
2
0
0 1
4
16
1
64
(a) Uniform Random
16
64
(b) Diagonal
B. Experimental Settings To evaluate performance, we used written benchmarks that directly measure the execution time of the SpMM operation. All code was compiled using the GNU C++ compiler, with Intel oneAPI to link MKL and OpenMP for parallelization. Experiments were conducted on a single node of the Perlmutter system, and benchmarks were executed using 64 threads. Memory allocation was controlled using numactl --interleave=all, distributing memory pages across available NUMA domains on the node. Memory bandwidth was measured using the STREAM [30] benchmark, where we measured a peak bandwidth β of 122.6 GB/s. Table IV summarizes the CPU specifications of the system used for evaluation. For measuring, only the actual SpMM operation was recorded, disregarding file loading times and variable initialization. This ensures arithmetic intensity, and the performance for the roofline model includes only floating point operations from the matrix multiplication calculation. We compared each matrix with CSR, MKL, and CSB implementations, as well as by d, the number of columns in B. Specifically, the values used are d = 1, 4, 16, 64. The CSB implementation evaluated in this work is based on the original Cilk Plus implementation [14] and was modified and updated to support OpenMP in order to enable direct comparison for SpMM.
4
50
100
40
75
30
50
20
25
10
0
0 1
4
16
(c) Blocking
64
1
4
16
64
(d) Scale-Free
Fig. 1: SpMM Performance (GFLOP/s) vs Number of Columns (d) for various sparsity patterns. Insets are visualizations of the specific sparsity pattern . Matrices used are (a) er 22 1, (b) rajat31, (c) road usa, and (d) com-LiveJournal.
C. Performance of SpMM implementations We first evaluate the performance of three SpMM implementations based on CSR and CSC data structures, as well as Intel MKL, for all matrices listed in Table III. Experiments are conducted with dense matrix column dimensions d = 1, 4, 16, and 64, and the complete results are reported in Table V. We then summarize the performance using four representative matrices, one from each sparsity pattern, and present these results in Figure 1. We observe that uniformly random matrices consistently achieve the lowest performance, while scale-free matrices attain the highest performance across all three implementations