Staging Blocked Evaluation over Structured Sparse Matrices

The matrices used in many computational settings are naturally sparse, holding a small percentage of nonzero elements. Storing such matrices in specialized sparse formats enables algorithms that avoid wasting computation on zeros, significantly accelerating common matrix computations like sparse matrix-vector multiplication (SpMV) and sparse matrix-matrix multiplication (SpMM). In many real-world sparse matrices, however, nonzero elements are densely clustered in subregions of the matrix. For matrices that feature this sort of structured sparsity, hybrid formats can further improve performance by representing these subregions as dense blocks. Existing hybrid formats either fix the dimensions of dense blocks, padding irregular regions with zeros and wasting computation, or incur run-time overhead when iterating over variable-sized blocks. This paper presents SABLE, a framework for accelerating structured sparse matrix computations by using staging to achieve the best of both of these approaches. Ahead of execution, SABLE inspects the matrix to identify variable-sized dense subregions, which it stores using a new hybrid format. It then eliminates the overhead typically associated with variable-sized blocks by using staging to generate specialized code that is amenable to vectorization. We evaluate SABLE on SpMV and SpMM kernels using matrices from the popular SuiteSparse data set. SABLE outperforms the best available SpMV baseline by ${\sim}$10% on average, and SpMM baselines by ${\sim}$20%. When parallelized, SABLE achieves further speedups of up to ${\sim}7\times$ on SpMV and SpMM over the best fully-sparse baseline when using 8 threads.

💡 Research Summary

The paper introduces SABLE, a framework designed to accelerate computations on structured sparse matrices—matrices in which non‑zero entries tend to cluster into dense sub‑regions rather than being uniformly scattered. Traditional sparse storage formats such as CSR or COO store each non‑zero individually and rely on indirection arrays for traversal. While these formats reduce memory usage, they incur significant overhead when the ratio of zeros to non‑zeros is high, because each element still requires indirect addressing and branch‑heavy loops. Block‑based formats like BCSR or VBR mitigate some of this overhead by grouping rows and columns into fixed‑size or variable‑size dense blocks, improving cache locality and enabling vectorized kernels. However, fixed‑size blocks often require padding with zeros, wasting computation, and variable‑size blocks introduce runtime loop‑bound determination costs, which also degrade performance.

SABLE tackles these issues by combining an inspector phase (the “partitioner”) with a staging‑based code generation phase (the “compiler”). The partitioner first inspects a given matrix (provided in Matrix Market format) and automatically identifies high‑density sub‑matrices, called δ‑dense blocks, where at least δ = 50 % of the entries are non‑zero. Two heuristic thresholds guide the search: a minimum density of 0.5 and a minimum area A_min = 2500 (≈ 50 × 50). The algorithm works on both CSR and CSC representations, scanning possible outer‑dimension spans (rows for CSR, columns for CSC) and, for each start position, building a histogram of non‑zero counts across the inner dimension. A two‑pointer sliding window over this histogram finds the inner‑dimension interval that maximizes a scoring function: score = area × (density − 0.5)^1.5. The best candidate block is trimmed by iteratively removing boundary rows or columns that would improve overall density, then accepted if it meets the density and area thresholds and does not overlap previously accepted rows or columns. Accepted blocks are marked, and the algorithm recurses on the four surrounding sub‑regions (above, below, left, right). This recursive decomposition guarantees non‑overlapping blocks, a requirement for the storage format that follows. The partitioner discards blocks that are only one column wide (since the current dense kernels iterate over rows) and finally checks that at least 10 % of the matrix’s non‑zeros are covered by accepted blocks; otherwise the matrix is deemed unsuitable for SABLE.

The identified blocks are stored in a new hybrid format called VBR‑C (Variable Block Row – Compressed). VBR‑C extends the classic Variable Block Row (VBR) format: δ‑dense blocks are stored as dense sub‑matrices in a contiguous value array, together with indirection arrays that record block row and column boundaries and the start index of each block. The remaining non‑dense region is stored in a user‑chosen sparse format (e.g., CSR, CSC). This hybrid approach allows dense blocks to be processed with highly optimized dense kernels (BLAS routines or hand‑written shape‑specialized kernels) while preserving the ability to plug in any future sparse library for the rest of the matrix.

The staging compiler receives the VBR‑C metadata (block coordinates, dimensions, and density) and generates specialized C code for each block type. For blocks with density > 0.5, the compiler emits calls to BLAS (cblas_dgemv for SpMV, cblas_dgemm for SpMM) or to custom kernels that are tuned for the exact block shape. Because the loop bounds are known at compile time, the generated code is free of the indirect‑addressing overhead typical of generic sparse kernels and can be automatically vectorized by the compiler. For blocks with lower density, the compiler emits calls to the user‑provided sparse kernel. Empty blocks are omitted entirely. The generated code, together with the VBR‑C data and the user’s sparse kernel, forms the “executor” that runs at runtime. The executor dispatches each block to its appropriate kernel, aggregates the results, and produces the final output vector (SpMV) or matrix (SpMM).

Experimental evaluation uses a representative subset of the SuiteSparse matrix collection, covering a variety of real‑world sparsity patterns. On single‑threaded runs, SABLE’s SpMV implementation outperforms the best CSR‑based baseline by an average of 10 % and up to 2× on certain matrices. For SpMM, average speedup is 20 % with a peak of 2.6×. When parallelized with OpenMP across eight threads, the framework achieves up to 7× speedup over the fully‑sparse baseline, demonstrating that the block‑level parallelism scales well. The partitioner’s runtime is modest (seconds for most matrices) and respects a four‑hour timeout for the few extremely large or irregular cases, which still produce useful partial partitions.

The paper’s contributions are fourfold: (1) a heuristic partitioning algorithm that automatically discovers non‑overlapping high‑density blocks, (2) the VBR‑C hybrid storage format that cleanly separates dense and sparse regions, (3) the SABLE staging compiler that generates structure‑aware, vectorizable code, and (4) a thorough empirical study showing consistent performance gains on both SpMV and SpMM across diverse matrices. The authors argue that the approach is broadly applicable to domains such as pruned neural networks, graph analytics, and scientific simulations where structured sparsity is common. Future work could explore adaptive selection of δ, more sophisticated block‑shape kernels, or integration with auto‑tuning frameworks to further push the performance envelope.