Ocean: Fast Estimation-Based Sparse General Matrix-Matrix Multiplication on GPU Yifan Li
Giulia Guidi
Cornell University Ithaca, NY, USA [email protected]
Cornell University Ithaca, NY, USA [email protected]
arXiv:2604.19004v1 [cs.DC] 21 Apr 2026
Abstract In computational science and data analytics, many workloads involve irregular and sparse computations that are inherently difficult to optimize for modern hardware. A key kernel is Sparse General Matrix–Matrix Multiplication (SpGEMM), which underpins simulations, graph analytics, and machine learning applications. SpGEMM exhibits irregular memory access patterns and workload imbalance, making it challenging to achieve high performance on GPUs. Current GPU SpGEMM solutions typically rely on a two-pass workflow to address load imbalance and reduce memory access. The symbolic pass, which determines the number of output elements per row, accounts for roughly 28% of the total runtime on average. In this work, we question the necessity of exact symbolic computation and introduce an estimation-based SpGEMM workflow. Our approach replaces the costly symbolic step with lightweight HyperLogLog estimators, combined with an analysis strategy that dynamically selects the optimal workflow and guides accumulator configuration. In addition, we introduce a hybrid accumulator design, including a novel hash-based accumulator that leverages both shared and global memory. Our approach consistently outperforms leading GPU SpGEMM implementations across a wide range of both square and rectangular matrices, achieving speedups of 1.4×–2.8× on NVIDIA A100 and H100 architectures.
CCS Concepts • Theory of computation → Massively parallel algorithms; • Computing methodologies → Linear algebra algorithms. ACM Reference Format: Yifan Li and Giulia Guidi. 2026. Ocean: Fast Estimation-Based Sparse General Matrix-Matrix Multiplication on GPU. In 2026 International Conference on Supercomputing (ICS ’26), July 06–09, 2026, Belfast, United Kingdom. ACM, New York, NY, USA, 13 pages. https://doi.org/10.1145/3797905.3807868
1
Introduction
Data-driven, irregular, and sparse computations are central to modern data analytics and scientific computing, but they are notoriously hard to optimize on current hardware. A key kernel is Sparse General Matrix–Matrix Multiplication (SpGEMM), which underpins simulations, graph analytics, and machine learning applications. SpGEMM multiplies two sparse matrices, A and B, to produce
This work is licensed under a Creative Commons Attribution 4.0 International License. ICS ’26, Belfast, United Kingdom © 2026 Copyright held by the owner/author(s). ACM ISBN 979-8-4007-2522-7/2026/07 https://doi.org/10.1145/3797905.3807868
an output matrix C. This operation is fundamental to many applications, including algebraic multigrid, Markov clustering, and computational biology [3, 7, 16, 17, 31]. It is also used in machine learning, such as in graph neural networks [32, 33], and has the potential to be adopted in sparse attention mechanisms. GPUs have recently become the focus of SpGEMM research, as their high computational throughput and memory bandwidth make them more powerful than CPUs, especially for large matrices. The underlying mathematical operation is the same as general matrix multiplication, but the high sparsity of the input matrices makes conventional dense matrix multiplication inefficient. SpGEMM implementations [11, 20, 29, 34, 36, 38] operate only on the nonzero elements of the input matrices to avoid unnecessary computation on zero entries. This computation paradigm, however, poses challenges for parallelization, especially on massively parallel architectures such as GPUs. The often irregular nonzero patterns of input matrices result in unpredictable computation and memory access behavior. Moreover, the nonzero pattern of the output matrix and the number of intermediate products are unknown before execution. Therefore, efficient approaches depend on dynamic load balancing and fine-grained allocation of on-chip resources. Currently, most GPU-based SpGEMM solutions use variants of Gustavson’s algorithm to address these issues. Gustavson’s algorithm computes each row of the output matrix C independently, providing natural parallelization opportunities on GPUs, where a thread block is typically assigned to each output row. This reduces the scratchpad memory required per block and allows load balancing across rows to be partially managed by the hardware scheduler. In addition, state-of-the-art solutions typically adopt a two-pass approach to fully utilize hardware resources during numeric computation [11, 29, 38]. The symbolic pass is executed first to calculate the number of nonzero elements in each output row. The numeric pass then performs the numerical computations and writes the results. The output of the symbolic pass is critical for guiding the allocation of scratchpad memory during the numeric pass to ensure efficient resource utilization. Both the symbolic pass and numeric pass use accumulators to store temporary results. Different types of accumulators, such as hash-based and dense accumulators, have been studied extensively [10, 11, 15, 23, 29, 36, 38]. Each accumulator offers advantages for specific matrix types. The accumulator is stored in scratchpad memory, with global memory used as a fallback when the required capacity exceeds the scratchpad’s maximum limit [11, 27, 29]. Current strategies have addressed several major challenges in GPU-based SpGEMM, including dynamic accumulator selection [29], global and local load balancing [29, 36], and architecture-specific optimizations [11, 34]. However, computational efficiency remains
ICS ’26, July 06–09, 2026, Belfast, United Kingdom
well below hardware limits. For example, spECK and HSMU-SpGEMM achieve only about 3−4% of peak FP64 performance on an RTX 3090 Ti on average [38], largely due to memory access bottleneck and complex control flow. The two computation passes, often exceeding 85% of total runtime, constitute the main bottleneck. Challenges. Two notable limitations in the widely adopted computational paradigm for GPU-based SpGEMM can restrict its effi1 The cost of size prediction can outweigh the benefits of ciency. ○ load balancing and precise resource allocation. Our experimental results show that, for the state-of-the-art solution spECK [29], the symbolic step accounts for ≈ 28% of the total execution time, even though it produces only a single value per output row. A lightweight solution for size prediction has the potential to significantly reduce this overhead. Current output size estimation approaches often lack accuracy and granularity, and can be more costly than 2 Current accumulator symbolic computation for some matrices. ○ designs may not perform efficiently for both short and long output rows. Short rows, defined as having fewer nonzeros than a warp (typically 32), underutilize hardware due to their limited size. In contrast, long rows, with more nonzeros than the capacity of a hash accumulator on a single streaming multiprocessor, revert to nonadaptive global-memory kernels and incur substantial performance penalties on large matrices. Contributions. In this paper, we address these challenges and propose an estimation-based SpGEMM solution. In particular, we replace the symbolic pass with fast HyperLogLog-based estimation. HyperLogLog [13] provides efficiency and sufficient accuracy for per-row output size prediction. Because the estimation-based workflow may incur higher computation and memory access costs, we introduce an analysis step that uses input statistics and sampling to dynamically select between estimation-based and symbolic-based workflows. The information gathered during this step is also used to further optimize the accumulation kernels. In addition, we optimize the accumulation of short and long output rows, often overlooked by existing solutions, through a hybrid accumulator design and cooperation between shared and global memory. Our solution, Ocean, is implemented in CUDA C++1 . Our experiments, conducted on an NVIDIA A100 platform using 337 square matrices, show that Ocean achieves speedups of 1.4×, 2.6×, 3.5×, and 2.0× over spECK, opSparse, TileSpGEMM, and HSMUSpGEMM, respectively. Consistent speedup is also observed on the NVIDIA H100 platform and on a dataset of 64 rectangular matrices. In summary, our main contributions include: (1) A new paradigm for GPU SpGEMM that uses HyperLogLog estimation instead of exact symbolic computation. To the best of our knowledge, this is the first work to use HyperLogLog to accelerate sparse linear algebra, highlighting the potential of estimators for sparse primitives; (2) An analysis step that accurately predicts the cost of estimation and selects the best workflow with minimal overhead; (3) A hybrid accumulator design that leverages three types of accumulators and a specialized configuration that exploits cooperation between shared and global memory; (4) Ocean, an open-source GPU implementation, that achieved speedups of 1.4× to 2.8× over the state-of-the-art. 1 Code available at https://github.com/CornellHPC/Ocean-SpGEMM
Li and Guidi
2 2
1
1
1
4
1 3
5 1
2
0
2
1
9
2
1
5
1
9
1
2
0
Figure 1: Computation of a single output row using Gustavson’s algorithm with a hash-based and dense accumulator. Nonzero elements of the input row are multiplied with their corresponding rows in matrix B to generate intermediate products. These intermediate products are stored in accumulators until all computations are completed.
2
Background
In this section, we provide background information on the basics of SpGEMM, state-of-the-art approaches for SpGEMM on GPUs, and HyperLogLog estimator.
2.1
CSR format and Gustavson’s Algorithm
Given sparse matrices A ∈ R𝑚×𝑘 and B ∈ R𝑘 ×𝑛 , sparse matrix– matrix multiplication computes C = AB, where C ∈ R𝑚×𝑛 . For SpGEMM computation, the input matrices are typically stored in Compressed Sparse Row (CSR) format, where two arrays store the column indices and values of the nonzero entries. The nonzeros in each row are stored contiguously, and a sorted array of row offsets indicates the starting position of each row. This storage format reduces space usage by eliminating redundant row index storage. The CSR format is particularly well suited for SpGEMM computation using Gustavson’s algorithm [18]. Gustavson’s algorithm computes each row of the output matrix C independently, with each row depending only on a single row of the input matrix A. Each row is computed using two nested loops. The outer loop traverses the nonzeros in the corresponding row of A, and the inner loop traverses the row of B associated with each of these nonzeros. During this accumulation process, intermediate results are stored in memory using accumulators. Because Gustavson’s algorithm accesses nonzeros within the same row of both input matrices contiguously, the CSR representation enables efficient traversal. The combination of the CSR format and Gustavson’s algorithm is used in most SpGEMM approaches [11, 23, 25, 27, 29, 38].
2.2
Accumulators
Gustavson’s algorithm uses accumulators to store intermediate results until an entire output row is computed. Figure 1 illustrates a simple example of computing a single output row using Gustavson’s algorithm with two different accumulator designs. GPUs consist of multiple independent compute units called streaming multiprocessors (SMs). Each SM has a fixed amount of hardware resources, including compute cores, registers, and fast on-chip scratchpad memory. Scratchpad memory not explicitly allocated to shared memory is automatically configured as L1 cache, while SMs share
Ocean: Fast Estimation-Based Sparse General Matrix-Matrix Multiplication on GPU
Analysis Binning
AB
Symbolic Computation Binning
C
Numeric Computation Post-Processing
Figure 2: Common GPU SpGEMM Workflow. a unified L2 cache. Global memory is shared by SMs and is much larger but significantly slower. Because accumulators are accessed frequently, they are typically placed in scratchpad memory. The simplest design is the dense accumulator, in which the entire output row is represented as a dense array [18]. This approach provides fast access but is impractical for many matrices because the output row may exceed the limited capacity of scratchpad memory. Hash-based accumulators are more space-efficient and are commonly used in GPU SpGEMM implementations [11, 14, 27, 29]. They use a hash table to store the column indices and values of nonzero output elements. During numeric computation, the hash table size is typically set to ∼ 1.5× the number of nonzero elements in the output row to reduce hash collisions [10, 11, 29]. Compared to dense accumulators, hash-based designs are generally more memoryefficient because of the sparsity of output rows, even though they require storing element indices. Their access latency is higher, as each update involves at least one additional lookup for the index and sometimes requires atomic exchanges on the hash table. ESC accumulators (Expand–Sort–Compact) use a different strategy [4]. Rather than merging intermediate products with identical column indices during row traversal, they preserve the intermediate products. Once all products are generated, a sorting step groups products with the same column indices. Finally, a compaction step aggregates these groups to produce the output row. ESC accumulators require more temporary memory and are generally faster for matrices with small numbers of intermediate products.
2.3
ICS ’26, July 06–09, 2026, Belfast, United Kingdom
not known in advance, the computed indices and values cannot be written directly to their final positions in the output CSR matrix. As a result, an additional memory copy step is required to reorganize the results into CSR format. Gustavson’s algorithm is also commonly used for symbolic computation. During symbolic accumulation, numerical values are discarded, and only the column indices of the output elements are retained. This reduces the scratchpad memory requirement by 67% for hash-based accumulators and by 98% for dense accumulators, assuming column indices are stored in uint32 and values in fp64. Thus, symbolic computation is generally more efficient than numeric computation. Fine-grained accumulator selection is another key factor in achieving good performance [11, 29, 38]. Because allocating the exact scratchpad memory needed for each row is impractical in GPU programming, a common strategy is to predefine multiple kernel configurations with fixed scratchpad sizes and thread counts. Rows are assigned to different kernels based on their output size before computation kernels are launched. This process is typically referred to as binning. Certain implementations further extend this approach by supporting multiple types of accumulators [23, 29]; for example, spECK uses both hash-based and dense accumulators. Hash-based accumulators are generally more space-efficient for sparse rows, while dense accumulators provide lower latency and better space efficiency for rows whose nonzero entries are narrowly distributed—that is, the distance between the first and last nonzero is small. To enable this selection, the binning or computation steps must choose between accumulator types, which requires additional information, such as the nonzero span of each output row. For this purpose, some solutions introduce an analysis step as the first stage of the workflow to collect relevant data [11, 25, 29, 38]. Figure 2 illustrates the common workflow. The symbolic step computes the number of output elements for each row, while the numeric step performs the actual computation. Binning is applied before both the symbolic and numeric steps to maximize performance. The numeric binning is guided by the number of output nonzeros per row, and the symbolic binning is typically guided by the number of intermediate products, serving as a safe upper bound. A post-processing step finalizes the output to comply with the CSR format (e.g., if the numeric computation produces unsorted rows, the post-processing step sorts the entries in each row).
GPU SpGEMM Workflow
Because they efficiently use scratchpad memory and offer great flexibility, hash-based accumulators are the most commonly used accumulators for GPU SpGEMM [2, 11, 14, 25, 29]. To achieve optimal efficiency, the size of the hash table must be determined for each row. However, on GPUs, shared memory is allocated to a thread block before execution. As a result, achieving high utilization requires knowing the number of output elements before numeric accumulation begins. This size estimation phase is commonly called the symbolic pass, while the subsequent computation is known as the numeric pass. Most recent GPU SpGEMM approaches use a symbolic pass for accurate size prediction. The prediction is also motivated by the need to determine the output matrix structure and minimize memory movement. Due to the constraints of the CSR format, output indices and values must be stored in contiguous memory. If the exact per-row output size is
2.4
HyperLogLog
HyperLogLog (HLL) is a probabilistic cardinality estimation algorithm that provides near-optimal space efficiency for counting distinct elements in data streams [13]. It is based on the observation that the distribution of leading zeros in uniformly hashed values provides information about the number of distinct elements. Each HLL sketch maintains a compact array of 𝑚 registers. Each element is hashed and assigned to a register using a fixed number of lower hash bits. Each register records the maximum number of leading zeros observed among the hash values assigned to it. The collection of registers captures the distribution of these maxima, reflecting the probability distribution of the dataset’s cardinality. By computing the harmonic mean across the 𝑚 registers and applying a statistical correction, HLL provides an unbiased estimate with a
ICS ’26, July 06–09, 2026, Belfast, United Kingdom
Li and Guidi
x
x x
x
x
x
x
x
x
Sketch
x
x
Sketch
x
Sketch
Sketch
Sketch
Figure 3: The estimation of per-row nonzeros using HyperLogLog with construct-and-merge. The sketches are built for each row of matrix B, merged according to each row of matrix A, and the number of output nonzeros are estimated from the merged sketch.
√ relative error of 1.04/ 𝑚. As the number of registers in the sketch increases, the precision of the HLL estimator improves. Due to its small memory footprint and high efficiency, HLL has been widely adopted in large-scale data processing for query optimization, network monitoring, and analytics over massive datasets [19]. To the best of our knowledge, it has not been applied to SpGEMM or any sparse linear algebra primitives.
3
Methods
In this section, we describe our SpGEMM solution, Ocean, which stands for Optimizing SpGEMM with Cardinality Estimation and Hybrid Accumulator. Ocean optimizes the symbolic and numeric steps, which together take over 85% of runtime. Ocean uses Gustavson’s method for numeric computation but differs from the traditional two-pass approach by using HyperLogLog to predict per-row output sizes and guide kernel selection. Our redesigned workflow introduces a lightweight analysis step for cost prediction and workflow selection, along with an overflow handling mechanism, enabling robust performance across diverse matrices. Finally, Ocean employs specialized numeric accumulators for short and long rows through a hybrid design.
3.1
Prediction with HyperLogLog
HyperLogLog (HLL) is adopted for per-row output size prediction in SpGEMM and is well suited to this task for three key reasons: (1) Controlled error: The estimation error remains within an acceptable range with a reasonable number of registers per sketch (see Section 5.3). (2) Constant memory footprint: HLL uses a constant amount of memory regardless of the number of nonzeros, making it effective for long rows where traditional accumulators require large kernels or global buffers. (3) Parallel updates: The sketch updates avoid compare-and-swap operations and can be implemented using atomicMax. These properties allow HLL to serve as a high-performance substitute for the accumulators used in the symbolic step. A key advantage of HyperLogLog is its ease of merging, which further increases the efficiency of size prediction. Combining multiple HLL
Figure 4: Overview of the Ocean SpGEMM workflow, with estimation-based components highlighted in darker backgrounds. The analysis step gathers metrics for workflow selection and load balancing. Then, size prediction uses either symbolic accumulation or estimation to predict per-row output sizes. Finally, numeric computation performs the multiplication using the collected statistics, followed by postprocessing to produce the CSR output.
sketches only requires taking the element-wise maximum across corresponding registers, provided the same hash function is used to construct them. This property is especially useful for per-row size prediction. Because the rows of matrix B are typically accessed multiple times, a sketch can be constructed for each row of B as a preprocessing step. Consequently, for each row of A, the corresponding sketches can be merged to produce the sketch of the corresponding row in C. The nonzeros in each row of C can then be estimated from the merged sketch, as shown in Figure 3. This construct-and-merge approach further reduces computation and memory access as long as the number of registers is smaller than the average number of elements per row of B. Moreover, the irregular symbolic accumulation is transformed into two regular, cachefriendly stages: sketch construction and sketch merging. Overall, HLL-based estimation avoids the high cost of symbolic computation while providing sufficiently accurate per-row size prediction.
3.2
Estimation-Based Workflow
HLL-based estimation cannot be integrated into existing SpGEMM solutions by simply replacing the symbolic component. Because the HLL estimator does not guarantee a strict error bound, a fallback mechanism must be provided to handle overflow. In addition, when matrix B is highly sparse, the cost of estimation may exceed that of the traditional symbolic step. In such cases, it is preferable to use the standard symbolic approach. In this work, we redesign the
Ocean: Fast Estimation-Based Sparse General Matrix-Matrix Multiplication on GPU
x x
x
x
x x
x x
x
x x
x
x
x
x
x
x
x
x
x
x
Figure 5: Computation of Input Expansion Ratio (𝐸𝑅) and Output Compression Ratio (𝐶𝑅). 𝐸𝑅 denotes the ratio between the number of intermediate products and the number of nonzeros in input matrix A. 𝐶𝑅 denotes the ratio between the number of intermediate products and the number of nonzeros in output matrix C.
SpGEMM workflow and its key components to support estimationdriven execution. Figure 4 illustrates the resulting Ocean pipeline. The discussion below focuses on the two central challenges of this approach: (i) overflow handling and (ii) cost prediction. Overflow Handling. Overflow during numeric accumulation occurs when HLL underestimates the number of output elements in a row, which can severely degrade performance. Fortunately, common accumulation kernel designs can tolerate estimation errors within a certain range. The hash table expansion factor, typically set to 1.5×, tolerates some underestimation without noticeable performance impact. Moreover, the binning process rounds up the expanded hash table size to a predefined bin size, further absorbing estimation errors. These two mechanisms greatly reduce the risk of overflow. However, numeric accumulation may still overflow in rare cases. Both hash-based and dense accumulators are susceptible to overflow. For hash-based accumulators, accumulation cannot proceed when the allocated shared memory is insufficient. For dense accumulators, accumulation proceeds safely because shared memory is allocated according to the output index range. Under the current design with pre-allocated global output memory, however, the final results may exceed the available space. Ocean handles overflow with a single fallback kernel, launched after all normal accumulation kernels complete. This kernel uses the largest dense configuration, iterates over the entire row, and can handle arbitrarily long rows. Output memory is allocated based on the number of intermediate products, which serves as an upper bound. This approach uses conservative resource allocation and therefore reduces efficiency; however, the overall impact is limited because overflow occurs in only a small fraction of rows. Cost Prediction and Workflow Choice. HLL-based estimation is generally efficient, but, in some cases, it can be more expensive than symbolic accumulation. An effective yet lightweight analysis step is required to predict the relative cost of estimation and revert to the symbolic approach when necessary. Here, we propose the use of the Input Expansion Ratio (ER) and Output Compression Ratio (CR) for this purpose. The input expansion ratio is used to estimate the relative cost of the symbolic step compared to the estimation step. HLL estimators require a fixed number of registers to guarantee a target precision,
ICS ’26, July 06–09, 2026, Belfast, United Kingdom
regardless of the number of processed elements. For some input matrices, merging HLL sketches requires more computation and memory traffic than a full accumulation pass, making the constructand-merge approach inefficient. Here, we define the Expansion Ratio (𝐸𝑅) as the ratio of intermediate products to the number of nonzeros in matrix A. 𝐸𝑅 estimates the relative memory traffic of symbolic accumulation compared to HLL-based estimation. Because size prediction is typically memory-bound, if this ratio falls below a predefined threshold (e.g., one quarter of the number of registers), the matrix is considered very sparse and the symbolicbased workflow is enabled. These statistics are collected through an O (𝑛𝑛𝑧 A ) analysis of matrix A, similar to the analysis step in spECK, where 𝑛𝑛𝑧 A denotes the number of nonzeros in A. The output Compression Ratio 𝐶𝑅 has been introduced in prior work [11, 24, 29] to characterize different matrices. It is defined as the ratio of the number of intermediate products to the number of nonzeros in the output matrix C. In this context, we use 𝐶𝑅 to estimate the relative cost of the additional post-processing required by the estimation-based workflow. Recall that the CSR format requires indices and values to be stored in contiguous memory. Unlike symbolic-based approaches that determine the output structure in advance, estimation-based approaches obtain this information only after numeric computation is complete. A compaction step at the end is required to relocate the output nonzeros and produce a valid CSR representation, which introduces additional memory movement overhead. 𝐶𝑅 measures the relative cost between symbolic accumulation and compaction overhead in terms of memory traffic; lower values indicate higher compaction memory traffic and favor a non-estimation-based workflow. A precise or approximate value of 𝐶𝑅 serves as a useful indicator, but it cannot be determined until a complete symbolic or estimation step has been performed. Relying on the estimated 𝐶𝑅 computed from a full estimation pass would introduce unnecessary overhead for matrices with high compaction costs. To address this limitation, we introduce the Sampled Compression Ratio. During the analysis step, an HLL sketch is constructed for each row of matrix B. A subset of rows from matrix A is randomly selected as samples, and the corresponding HLL sketches are merged and evaluated. The resulting sampled 𝐶𝑅 serves as an indicator for selecting either the estimation-based workflow or the precise symbolic step. The cost of computing the sampled 𝐶𝑅 is kept low to ensure the analysis step remains lightweight. HLL sketch construction has a complexity of O (𝑛𝑛𝑧 B ). Only a small fraction of rows in matrix A (3%) are sampled and evaluated. The sketch construction and merging process accounts for only 3% of the total runtime on average. To further optimize performance for a small number of extremely sparse matrices, both symbolic accumulation and HLL-based estimation are disabled when the average number of intermediate products per row falls below a lower threshold. In this case, the estimation-based workflow is selected, using the number of intermediate products as an upper-bound estimate. Overall, our analysis uses multiple metrics to choose between estimation-based workflows (including upper-bound estimation) and the symbolic-based workflow. From matrices A and B, we obtain 𝐸𝑅 and the number of intermediate products per row, while row sampling estimates 𝐶𝑅. Figure 5 illustrates these statistics.
ICS ’26, July 06–09, 2026, Belfast, United Kingdom
// Shared Memory FP 64 Atomic Add
1 2
. L _x _8: LDS .64 R4 , [ RZ ] ; DADD R6, R 4 , 1 ; ATOMS . CAST . SPIN .64 R 6 , [ RZ ] , R 4 , R 6 ; ISETP . EQ . U 32. AND P 0 , PT , R 6 , 0 x 1 , PT ; ISETP . EQ . U 32. AND . EX P 0 , PT , R 7 , RZ , PT , P 0 ; @!P0 BRA `(. L_ x _8) ;
3 4 5 6 7 8 9
// Global Memory FP 64 Atomic Add RED .E. ADD . F 64. RN . STRONG . GPU [ R 2.64] , R 8 ;
10 11
Listing 1: FP64 Atomic Add on Shared and Global Memory.
Together, these statistics characterize the relative cost of estimationbased workflows and can be computed with minimal overhead. The workflow selection criteria based on these statistics are presented in Section 4.3 and Table 1.
3.3
Hybrid Accumulator
Then, we redesign the hybrid accumulator scheme to better handle very short and very long rows. Prior hybrid design has combined hash-based and dense accumulators with multiple kernel configurations [29]. Building on this approach, we introduce: (1) An enhanced hash-based accumulator that leverages both shared and global memory to efficiently handle long rows, and (2) ESC accumulators for short rows with a limited number of intermediate products. Here, long rows are defined as rows with many nonzeros, typically exceeding the maximum capacity of a hash accumulator on a single streaming multiprocessor, while short rows are rows with fewer nonzeros than the size of a warp, which is 32 on NVIDIA GPUs. Enhanced Hash-Based Accumulator. Hash-based accumulators must store both the column indices and the values of the output elements. Column indices are typically 32 bits, while values are usually 64 bits. Because scratchpad memory has limited capacity, rows that exceed a certain size must switch to global memory, which makes hash-based accumulators inefficient for long rows. Our experiments show that we can store values in global memory without impacting performance under certain conditions described below, contrary to the common belief that the entire hash table must reside in scratchpad memory for fast accumulation. This allows hash-based kernels to accumulate rows that are 3× longer without significant performance penalty. Our design is motivated by two key observations: (1) The index access patterns and operation frequency are considerably more complex than those for value access. The indices serve as element identifiers and must be read, compared, and swapped using atomic operations during insertion, while values are only involved in atomicAdd operations. (2) In practice, storing values in global memory does not significantly degrade performance because of GPU architectural features. Hash-based accumulators are currently limited more by memory access latency than by bandwidth. Prior studies show that the latency of atomic operations in global memory is only 2–4× higher than that of shared memory atomics in general [21]. Moreover, because the return value of the atomic
Li and Guidi
operation is not needed in this context, the compiler uses a fire-andforget access pattern, further reducing the impact of latency [35]. Finally, atomic operations on shared memory for FP64 data types are not natively supported on NVIDIA GPUs and are compiled into compare-and-swap loops, reducing the advantage of shared memory placement. Listing 1 shows the SASS code for an atomicAdd operation compiled for sm80. For shared memory, a loop with compareand-swap and branching instructions is required. For global memory, the operation compiles into a single hardware-supported SASS instruction that uses a fire-and-forget access pattern. ESC Accumulator. Certain matrices exhibit extremely high sparsity, and their average number of intermediate products falls below the smallest block size (typically 64, chosen to ensure good GPU occupancy). In such cases, the ESC accumulator, which first expands the intermediate products and then sorts and compacts them into the final results, can be selected to take advantage of its specialization for short rows. The ESC accumulator is particularly well suited to the upper-bound estimation workflow, as its configuration does not depend on the number of output elements. Instead, the number of intermediate products per row serves as the selection criterion. The kernel is further optimized to process multiple rows concurrently, enabling efficient accumulation for short rows. Accumulator Selection. In total, three types of accumulators– hash-based, dense, and ESC–are used, each with multiple configurations for different block sizes. The accumulators are selected on a per-row basis to ensure proper load balance and best performance. As in previous work, the configuration that requires the fewest resources for a given row is selected, as a smaller kernel generally helps reduce load imbalance and synchronization costs. Dense kernels are preferred over hash-based kernels when both require the same amount of resources, except for the largest hash-based kernel, which uses our hybrid-memory design. The ESC accumulator is selected only in the upper-bound estimation workflow.
4
Kernel Configuration and Optimization
In practice, achieving high performance in GPU SpGEMM requires careful kernel design. Ocean leverages information from the analysis step to accelerate accumulation and applies architecture-aware optimizations with carefully tuned parameters.
4.1
Assisted Kernels
The estimation-based workflow is not chosen for matrices that would lead to high estimation and compaction costs. In these cases, symbolic computation is required, and information gathered during the analysis step can be used to accelerate symbolic accumulation. Prior work uses the number of intermediate products as the accumulator selection metric for the symbolic step [11, 29]. This approach often overestimates the required scratchpad memory, resulting in underutilization of hardware resources. In contrast, Ocean computes the average and standard deviation of the compression ratio from sampled rows and derives a conservative estimate of the compression ratio for the entire matrix. The number of intermediate products for each row is then divided by this estimated ratio and used to select the accumulator.
Ocean: Fast Estimation-Based Sparse General Matrix-Matrix Multiplication on GPU
Ocean uses the analysis step to speed up dense accumulators in both symbolic and numeric phases, akin to assisted symbolic computation. Dense accumulators maintain a bitmap to track nonempty elements, which is frequently updated because each insertion requires at least one bitmap access. Ocean uses 𝐶𝑅 to determine whether to query the bitmap before writing. Querying the bitmap increases read latency slightly but can significantly reduce write operations, which are often a bottleneck in shared memory. Therefore, the bitmap query is enabled only when the estimated compression ratio is above a predefined threshold.
4.2
Method
Selection Criteria
Upper-bound Estimation HLL Estimation Symbolic Accumulation
𝑛𝑝𝑟𝑜𝑑𝑢𝑐𝑡𝑠 avg < 64 𝑛𝑝𝑟𝑜𝑑𝑢𝑐𝑡𝑠 avg ≥ 64 ∧ 𝐸𝑅 ≥ 8 ∧ 𝐶𝑅 ≥ 8 Remaining cases
Table 1: Ocean’s selection criteria for different SpGEMM workflows: 𝐸𝑅 is the Input Expansion Ratio, and 𝐶𝑅 is the (Sampled) Output Compression Ratio.
Optimization
For hash-based accumulators, output rows must be sorted after accumulation to comply with the CSR format. The column index serves as the sorting key, and the corresponding accumulated value is permuted along with the key during the sorting process. This sorting step can incur significant overhead, especially for long rows. To accelerate the sorting of key–value pairs, we use indirect sorting [22]. Rather than sorting key–value pairs directly, key– ptr pairs are generated and sorted. The ptr element points to the corresponding value element and occupies fewer bits because of the small address space of on-chip memory. Once sorting is complete, values are written to the destination according to the sorted pointers. This approach reduces memory movement and lowers register pressure during sorting. In addition, we observe that most matrices have a column index range that does not require the full 32-bit representation. Because ptr is bounded and occupies no more than 14 bits, we concatenate the key and ptr into a single 32-bit integer when possible. Radix sort is used, and the bits corresponding to the ptr field are ignored during sorting [22]. This optimization further reduces register pressure, decreases memory traffic, and accelerates sorting. Independent Thread Scheduling (ITS), introduced in the Volta architecture, allows threads within a warp to diverge [28]. Our profiling shows that ITS increases memory access requests and warp stall samples inside the core accumulation loop, leading to performance degradation. This effect can be mitigated by compiling for an earlier architecture or by manually calling __syncwarp(). Ocean uses the former approach.
4.3
ICS ’26, July 06–09, 2026, Belfast, United Kingdom
Parameter Configuration
Parameter configuration is key for high performance and, for Ocean, primarily involves HyperLogLog settings, workflow selection criteria, and accumulation kernel configurations. HLL parameters determine estimation precision. In Ocean, we use 32 registers per HLL sketch when the input expansion ratio is less than 48, and 64 registers otherwise. This configuration balances estimation accuracy and computational efficiency. An empirical analysis of estimation precision is provided in Section 5.3. The sampling ratio in the analysis step is set to 0.03, with a minimum of 600 and a maximum of 10, 000 sampled rows. The threshold for enabling bitmap queries is set to 2. The relative variance of 1/𝐶𝑅 is: 1 𝑛𝑠𝑎𝑚𝑝𝑙𝑒𝑑
2
2
2
[𝜖 + 𝐶𝑉 (1 + 𝜖 )]
where 𝐶𝑉 = 𝜎𝐶 /𝜇𝐶 is the coefficient of variation of the density of output rows, and 𝜖 = 1.04/√𝑚 is the relative error of the HLL estimator with 𝑚 registers. Then, applying Chebyshev’s inequality [30] provides a probabilistic error bound on the sampled 𝐶𝑅. For a matrix with 200,000 rows, a sampling rate of 0.03, and 64 HLL registers, the relative error is below 3%, 6%, and 17% at 95% confidence for matrices with 𝐶𝑉 = 0.5, 1, and 3, corresponding to balanced, exponential, and heavily skewed output row distributions. Given the sampled 𝐶𝑅 is used for workflow selection, these error levels are acceptable even for highly skewed matrices. Table 1 summarizes the criteria for workflow selection. The upper-bound method, which skips both symbolic accumulation and HLL-based estimation, is used when the average number of intermediate products per row is below 64. This threshold matches the minimum number of threads per block required for full occupancy. For the remaining matrices, the HLL-based estimation workflow is selected when both 𝐸𝑅 and sampled 𝐶𝑅 are greater than or equal to 8. The threshold for 𝐸𝑅 is set to 8 because each HLL register occupies one byte, while each column index requires four bytes; at this point, HLL sketch merging incurs a memory access cost comparable to symbolic accumulation. The threshold for 𝐶𝑅 is set to 8, chosen empirically to amortize the cost of output compaction. For accumulation kernel configuration, Ocean uses the same profile for all architectures. The maximum shared memory allocated per SM for accumulation is capped at 128 KB, with the remaining capacity reserved for L1 cache. Our kernel configurations follow prior work [11, 29]. In total, five normal kernels and two specialized kernels are used for both hash-based and dense accumulators. The largest normal kernel uses half of the available shared memory and 1024 threads to achieve full occupancy. Each subsequent normal kernel uses half the resources of the previous one, resulting in progressively smaller accumulator sizes for shorter rows while achieving full occupancy. One specialized kernel uses the same resource configuration as the largest normal kernel and processes long rows using global memory. The other specialized kernel matches the smallest normal kernel and processes multiple short rows concurrently. For ESC accumulators, two configurations are provided, processing two or four rows concurrently with 64 threads.
5
Evaluation
This section evaluates Ocean through comparisons with state-ofthe-art solutions, analysis of HLL estimation and sampling accuracy, and a detailed breakdown of its core methods.
ICS ’26, July 06–09, 2026, Belfast, United Kingdom
Experiment Setup
Ocean is evaluated on two platforms: Perlmutter at NERSC and DeltaAI at NCSA. Each Perlmutter GPU node has one AMD EPYC 7763 CPU and four NVIDIA A100 GPUs with 40 GB of memory each. The system runs SUSE Linux Enterprise Server 15, and all experiments use CUDA 12.9. DeltaAI is an ARM-based supercomputer. Each compute node contains four NVIDIA GH200 superchips, each integrating a 72-core NVIDIA Grace CPU and an NVIDIA H100 GPU with 96 GB of memory. DeltaAI also runs SUSE Linux Enterprise Server 15, and CUDA 12.4 is used. Ocean is evaluated against four state-of-the-art GPU SpGEMM implementations: spECK [29], opSparse [11], TileSpGEMM [26], and HSMU-SpGEMM [38]. The code base of MOSparse [34] is not open source, so we do not include a comparison with it. For completeness, we include cuSPARSE in one experimental setting, although prior studies show it is generally less performant than the state of the art [14]. For each implementation, we follow the authors’ recommended guidelines to configure, compile, and run on our platforms. In some cases, minor code modifications are required to ensure a fair comparison. In particular, we modify spECK to include output memory allocation in the reported execution time, and we adjust configurations in HSMU-SpGEMM to support rectangular matrices. In our experiments, a few matrices fail the internal correctness checks of HSMU-SpGEMM. Our attempts to contact the authors were unsuccessful. Ocean and other SpGEMM approaches are evaluated on a collection of over 400 matrices, divided into two datasets for clarity: one set of square matrices and one set of rectangular matrices. The first dataset consists of 337 square matrices from the SuiteSparse Matrix Collection. This dataset was originally used in HSMU-SpGEMM [38] and includes matrices from the SuiteSparse matrix collection that require at least 100 million floating-point operations (FLOPs) to compute 𝐴𝐴 and 𝐴𝐴𝑇 . FLOPs are counted as twice the number of intermediate products. The second dataset consists of rectangular matrices from the SuiteSparse Matrix Collection, filtered to have no more than 10 billion nonzeros and at least 100 million FLOPs for computing 𝐴𝐴𝑇 , resulting in 64 matrices. The matrices from the two datasets cover a wide range of applications and display diverse structural properties, such as matrices that produce dense output and symmetric matrices. For the square dataset, the evaluation focuses on computing 𝐴𝐴. For the rectangular dataset, we evaluate the computation of 𝐴𝐴𝑇 . For each test, 10 warm-up runs are followed by 10 measured runs, and the average runtime of the measured runs is reported. A timeout of 30 seconds is imposed for all runs. If a tool fails on a matrix due to a runtime error or timeout, its runtime is replaced with the slowest valid runtime observed for that matrix across all tools. This penalizes failures, prevents artificially inflated performance, and ensures comparable results.
5.2
Overall Performance
The first section of Table 2 summarizes the overall performance statistics for the square-matrix dataset on NVIDIA A100 GPUs. Figure 6 illustrates the corresponding performance trends, with the number of intermediate products on the x-axis and the average achieved GFLOPS on the y-axis. Out of the 337 matrices, Ocean achieves
GPU (dataset)
Tool
#best
#2nd best
#fail
GFLOPS
A100 (square)
cuSparse spECK opSparse Tile HSMU Ocean
0 22 16 5 0 294
0 185 86 9 19 38
82 0 20 82 0 0
3.39 46.2 24.2 18.1 32.1 63.7
H100 (square)
spECK opSparse Ocean
21 4 312
204 110 23
0 18 0
66.2 40.1 108.0
A100 (rect.)
spECK opSparse HSMU Ocean
7 15 0 42
39 3 6 16
0 4 15 1
27.3 10.8 15.3 36.3
H100 (rect.)
spECK opSparse Ocean
8 8 48
38 12 14
0 4 0
42.6 17.1 60.7
Table 2: Comparison of SpGEMM approaches on different datasets and architectures. The square dataset includes 337 matrices. The rectangular dataset includes 64 matrices. The average GFLOPS is the geometric mean across matrices.
102
Ocean spECK opSparse
Tile HSMU
GFLOPS
5.1
Li and Guidi
101
108
Products
109
Figure 6: The smoothed GFLOPS achieved over the square matrices dataset on the A100. Ordered by the number of intermediate products. The line shade indicates the deviation.
the best performance on 86% of the matrices and ranks second on an additional 11%. cuSparse, spECK, opSparse, TileSpGEMM, and HSMU-SpGEMM achieve the best performance on 0%, 7.8%, 3.9%, 1.7%, and 0% matrices, respectively. Ocean achieves a geometric mean performance of 63.23 GFLOPS. On average, Ocean outperforms cuSparse, spECK, opSparse, TileSpGEMM, and HSMUSpGEMM by a factor of 18.8×, 1.4×, 2.6×, 3.5×, and 2×, respectively. The trend plot illustrates the performance of different SpGEMM approaches across matrices with varying numbers of intermediate products. As shown in the figure, Ocean consistently achieves the highest performance across the entire dataset, and its performance advantage increases as the number of intermediate products grows. This trend can be attributed to two main factors:
Ocean: Fast Estimation-Based Sparse General Matrix-Matrix Multiplication on GPU
8
spECK (baseline) Ocean opSparse HSMU
4 Speedup vs spECK
Speedup vs spECK
8
spECK (baseline) Ocean opSparse
4 2 1 1/2 1/4 1/8
ICS ’26, July 06–09, 2026, Belfast, United Kingdom
2 1 1/2 1/4
0
50
100
150 200 Matrix ID
250
300
350
1/8
0
10
20
30 Matrix ID
40
50
60
Figure 7: Relative speedup of SpGEMM solutions compared to spECK. The speedups below 1/8 are clipped to 1/8 in the plot. Results for square matrices (left) are shown on H100, and results for rectangular matrices (right) are shown on A100.
(1) The relative overhead of the analysis components decreases as overall computation increases; (2) The estimation-based workflow and enhanced hash-based accumulators are effective at handling longer rows, making Ocean more efficient for larger matrices. Of the six evaluated implementations, spECK achieves the secondbest overall performance. Its lightweight analysis, load-balancing strategy, and hybrid use of hash and dense accumulators enable consistently strong performance across a wide range of matrices. In addition, spECK disables global load balancing when the analysis results are uniform across all rows. This reduces overhead for certain matrices and allows spECK to outperform Ocean in these cases. HSMU-SpGEMM ranks third in overall performance. Our results suggest that the binary search-based design of its numeric accumulator has higher complexity than hash-based accumulators, which outweighs its greater utilization of shared memory. opSparse generally ranks fourth in performance. It outperforms Ocean and spECK on a small number of matrices, most of which are extremely sparse, due to its overhead optimizations and specialized design for the smallest accumulator kernel. However, its design prevents it from achieving consistent performance across matrices, a limitation also reflected in the high variance in performance. TileSpGEMM exhibits the second-lowest overall performance, although it achieves the best results on six matrices. Its tile-based computation strategy is effective primarily for matrices with block-sparse structure, which limits its general applicability. The square dataset is also evaluated on DeltaAI, equipped with NVIDIA H100 GPUs. Because DeltaAI is an ARM-based system, TileSpGEMM and HSMU-SpGEMM fail to compile as they rely on x86 intrinsics and are therefore excluded from this evaluation. The results are presented in the second section of Table 2. Figure 7 (left) further illustrates the relative speedup of Ocean and opSparse, using spECK as the baseline. spECK is chosen as the baseline because it successfully runs on all matrices and exhibits stable performance. Ocean achieves consistent speedup across the dataset, being the best-performing solution for 92% of the matrices, with average speedups of 1.6× over spECK and 2.7× over opSparse. In addition, we evaluate the rectangular dataset on NVIDIA A100 GPUs to demonstrate the generality of Ocean. TileSpGEMM does not support rectangular matrices and is therefore excluded from
this evaluation. Figure 7 (right) shows the relative speedup of the four remaining SpGEMM implementations. Ocean demonstrates consistently strong performance on this dataset, achieving the best performance on 42 matrices and the second-best performance on 18 matrices. On average, Ocean achieves speedups of 1.3×, 3.4×, and 2.4× over spECK, opSparse, and HSMU-SpGEMM, respectively. The evaluation on H100 is provided in Table 2. On the A100, Ocean fails on one matrix, JP, due to memory constraints. The estimation-based approaches inherently incur a peak memory usage of about 2.2× because of GPU memory management limitations. Specifically, results produced by the numeric step are not stored in consecutive memory locations and must be compacted to satisfy the CSR format. A small auxiliary buffer can compact non-consecutive outputs back into the original array, but the array remains over-allocated because in-place shrinking is not supported on GPUs. In situations where GPU memory is limited, this issue can potentially be mitigated by falling back to the symbolic step if insufficient memory is detected. The input expansion ratio and the sampled output compression ratio can be used to estimate the memory footprint of the estimation-based workflow.
5.3
Estimation Precision
To evaluate the precision of HyperLogLog under different register counts, we conduct experiments on the square dataset, excluding matrices that do not use the HLL-based estimation approach. The number of registers per HLL sketch is set to 32, 64, and 128, and the accuracy of per-row nonzero estimation is measured. Figure 8 (left) shows the cumulative distribution function (CDF) of the average relative estimation error. For each matrix, the relative estimation error is calculated as the arithmetic mean of the per-row relative errors. The average errors across all matrices with 32, 64, and 128 registers are 0.13, 0.10, and 0.07, respectively, with corresponding RMS averages of the standard deviation of 0.17, 0.15, and 0.14. These results indicate that increasing HLL precision consistently reduces estimation errors. They also show that HLL can provide sufficient accuracy for per-row nonzero estimation with a reasonable number of registers. Our empirical results show that fixing the register count at 32 or 64 causes an average overall slowdown of 1.4% and 1%, respectively, compared to the dynamic configuration described in Section 4.3.
5.4
Contribution Breakdown
Finally, we assess Ocean’s core innovation through an ablation study on the square dataset using an A100, evaluating four incremental versions of Ocean. The baseline (V1) disables the estimationbased workflow, the assisted kernels, and the hybrid accumulators (including the enhanced hash-based accumulator and ESC accumulators). V2 introduces the estimation-based workflow (labeled E), V3 adds the assisted kernels (labeled AS), and V4 incorporates hybrid accumulator support (labeled HA). Table 3 summarizes the results. Each version is compared with its immediate predecessor to report incremental speedup, while the final row labeled Overall compares V4 with V1. Results are reported separately for matrices that use symbolic-based prediction and those that use estimation.
100%
80%
80%
Cumulative Probability
100%
60% 40% 20% 0% 0.0
num_reg = 32 num_reg = 64 num_reg = 128
60% 40% 20% 0%
0.1 0.2 0.3 Relative Estimation Error
num_reg = 32 num_reg = 64 num_reg = 128
0
2 4 Overflow Rows (%)
Figure 8: CDF of the average relative estimation error of HLL Estimators (left) and ratio of overflow rows (right) under different number of HLL registers on the Square dataset. Analysis
Symbolic
Estimation
Numeric
Post-Processing
1.0 Normalized Runtime
0.8 0.6 0.4 0.2
pa V3 ra-7 pa V4 ra-7
ex d V2ata_ 1 ex da t V3 a_ 1
roa V1d_us a roa V2d_us a
C_ O V PF_ TS 1 1047 C_ O V2PF_1 04 7
0.0
TS
To compare the precision of Ocean’s HLL estimator with Cohen’s estimator [8], both estimators are configured to use the same amount of memory per output row, as estimation kernels are memory-bound. Cohen’s estimator has a memory access pattern similar to HLL on matrix A but requires substantially more memory accesses on matrix B. Therefore, this configuration favors Cohen’s estimator. Using 64 bytes per output row, HLL achieves a smaller average relative error for per-row estimation on all 148 matrices in the HLL estimation-based workflow, with an average error 2.1× smaller than Cohen’s. Even when Cohen’s estimator is given 4× more memory, HLL still outperforms it on 116 of the 148 matrices, consistent with theoretical expectations. These results show that HLL provides more accurate per-row output size predictions for SpGEMM than existing estimation approaches. In addition, we evaluate the fraction of overflow rows using HLL-based estimation. The HLL estimates are treated the same as symbolic results for binning. Each estimate is first multiplied by a coefficient (1.5 by default) and then rounded up to the nearest accumulator size. A row is considered overflowing when the actual number of output elements exceeds 80% and 100% of the allocated capacity for hash-based accumulators and dense accumulators, respectively. In this setting, the average overflow ratios are 1.2%, 0.3%, and well below 0.1% for 32, 64, and 128 registers, respectively. The corresponding maximum overflow ratios are 7%, 2.7%, and 0.5%. Figure 8 (right) shows the cumulative distribution of overflow ratios. These results indicate that overflows caused by estimation are rare and have limited impact on numeric accumulation. To further mitigate overflows at lower precision, we increase the expansion coefficient to 2.0 when using 32 registers. The sampling strategy proposed in Section 3.2 enables accurate estimation of the output compression ratio 𝐶𝑅. On the square dataset, the average relative sampling errors are 0.05, 0.04, and 0.03 for 32, 64, and 128 registers, respectively. For the three register configurations, only 2, 1, and 1 matrices are assigned to different workflow-selection categories when comparing the sampled 𝐶𝑅 with the ground-truth 𝐶𝑅. These error levels are sufficient for reliable workflow selection. The sampling rate of 0.03 is chosen empirically to balance accuracy and overhead: reducing it to 0.01 causes a 2.1% average slowdown, as the less accurate 𝐶𝑅 estimate leads to suboptimal workflow selection despite benefiting 21 individual matrices; increasing it to 0.05 causes a 3.5% slowdown due to the higher sampling cost.
Li and Guidi
Cumulative Probability
ICS ’26, July 06–09, 2026, Belfast, United Kingdom
Figure 9: Runtime breakdown of individual components across different versions. Runtimes are normalized to the first version for each matrix for clarity.
Of the 337 matrices, 148 are assigned to the HLL-based estimation workflow according to the criteria in Table 1, 26 use upper-bound estimation, and the remaining 163 use precise prediction. Ocean’s performance is generally consistent across input matrices of varying density throughout the dataset. The estimation-based workflow provides an average speedup of 1.3× for the selected matrices by significantly reducing the cost of the symbolic step. The speedup is greater for matrices with density in [10−2, 10−1 ) and for those below 10−5 . The maximum speedups are observed on matrices TSC_OPF_1047 for HLL-based estimation and road_usa for the upper-bound method, respectively. Figure 9 further illustrates the breakdown of runtime for individual stages. For TSC_OPF_1047, HLL-based estimation substantially reduces the cost of size prediction compared to the symbolic step, but slightly increases numeric kernel time due to suboptimal binning. For road_usa, the upper-bound method eliminates the cost of output prediction and improves numeric performance through more effective binning guided by the upper bound. To validate the cost-prediction model and the accuracy of workflow selection, we conduct an additional evaluation using 𝑉 4, forcing matrices to use either the estimation-based workflow or the symbolic-based workflow. Of the 163 matrices that initially selected the symbolic-based workflow, only 10.8% benefit from switching to the estimation-based workflow with a speedup of at least 1.05×.
Ocean: Fast Estimation-Based Sparse General Matrix-Matrix Multiplication on GPU
ICS ’26, July 06–09, 2026, Belfast, United Kingdom
Symbolic-based Speedup
Estimation-based Speedup
Overall Avg.
Version
Method
Avg.
Min.
Max.
Avg.
Min.
Max.
Speedup
GFLOPS
V1 V2 V3 V4
Baseline (V1)+E (V2)+AS (V3)+HA
1.04× 1.07×
0.96× 0.93×
1.65× 1.63×
1.30× 1.03× 1.04×
0.95× 0.99× 0.55×
2.77× 2.19× 2.53×
1.14× 1.03× 1.06×
50.8 58.3 60.2 63.7
1.12×
0.92×
1.66×
1.40×
0.99×
3.77×
1.25×
63.7
Overall
Table 3: Comparison of Ocean with different optimization strategies.
In contrast, among the 174 matrices that originally selected the estimation-based workflow, 92% achieve a speedup greater than 1.05× compared to reverting to the symbolic-based workflow. These results demonstrate that our cost-prediction model successfully selects the optimal workflow for most matrices. The assisted kernels provide an additional average speedup of 1.03×, with a maximum speedup of 2.19× on the matrix exdata_1. This matrix produces many dense output rows with a very high compression ratio, where frequent updates to the bitmap cause substantial memory traffic. As shown in Figure 9, the assisted kernels reduce the numeric step runtime by more than 50%. Overall, the performance gains from assisted kernels, particularly the assisted symbolic step, are modest. Our results show that, although an estimated compression ratio allows for allocating less shared memory and selecting smaller kernels, this does not always improve performance. Kernel efficiency depends on multiple factors, including intermediate products and load balancing. These results indicate that smaller kernel sizes do not inherently yield better performance, and additional metrics should guide kernel assignment. The hybrid accumulator provides an average speedup of 1.06×. The enhanced hash-based accumulator benefits three times as many matrices as the ESC accumulator, though each method yields an average speedup of 1.2× on its respective set of affected matrices. The largest performance improvement is observed on matrix SiO2, which has a small number of scattered rows requiring multiple dense-accumulator iterations, leading to a tail effect. Figure 9 illustrates the runtime breakdown for another matrix, para-7, which also benefits from the hybrid accumulator. Rows that would otherwise require multiple dense accumulator iterations are instead handled by the enhanced hash-based accumulator. The time spent on post-processing, including output sorting, increases, but overall performance improves, achieving a 1.6× speedup. A notable exception is torso1, which experiences a slowdown of 0.55× and is the only matrix with a slowdown greater than 0.9×. For torso1, 1,215 out of 116,158 rows are assigned to the enhanced hash-based accumulator. These rows generate a large number of intermediate products, resulting in high global memory traffic and performance degradation. The enhanced hash-based accumulator effectively hides memory access latency; however, adding a mechanism to detect and avoid throughput-bound cases would further improve its robustness. Finally, we analyze the runtimes of different components. For Ocean, analysis takes 7% of runtime, with HLL sampling just 2%, confirming the analysis is lightweight. For matrices processed with HLL-based estimation, the estimation step takes an average of 4% of
the runtime, and post-processing accounts for about 8%. In contrast, with the baseline version, the symbolic step for such matrices consumes 30% of the runtime, while post-processing accounts for about 4%. As a reference, the symbolic step takes an average of 28% of the total runtime in spECK. These results show that, although memory compaction cost increases, estimation incurs much less overhead than symbolic computation, resulting in a net performance gain.
6
Related Work
cuSparse [27] is one of the earliest SpGEMM solutions to use a two-pass (symbolic and numeric) approach. It uses hash-based accumulators for accumulation. spECK [29] follows the classic two-pass design. It uses both hash-based and dense accumulators. A lightweight analysis is introduced to dynamically select among different accumulators and algorithms. Its design aims to deliver robust performance across a wide range of input matrices. opSparse [11] also follows the two-pass approach and proposes multiple low-level, architecture-specific optimizations. TileSpGEMM [26] uses tiles as the minimum unit of computation. This novel computation strategy eliminates the need for symbolic computation, but a similar preprocessing step is still required to determine the output tile pattern. HSMU-SpGEMM [38] replaces hash-based accumulation with binary search on a pre-generated column-index array, eliminating the need for a load factor and improving shared memory utilization by about 1.5×. HSMU-SpGEMM also follows the two-step design, but the symbolic phase is modified to generate the column-index array. MOSparse [34] uses a LightGBM-based model to select between the two-pass “precise” method and the symbolic-free “upper-bound” method, while also incorporating architecture-specific optimizations and parameter tuning. The code is not publicly available. AC-SpGEMM [37] emphasizes its use of ESC accumulators with adaptive chunking. It performs well on very sparse matrices but is generally slower than newer approaches such as spECK [14, 29]. Therefore, we focus on more recent GPU SpGEMM approaches. Orthogonal to single-GPU approaches, recent work has addressed distributed multi-GPU SpGEMM using RDMA-based communication and hierarchy-aware decomposition [5, 6]. In particular, Bellavita et al. introduced trident partitioning, a 2D-1D hybrid scheme that reduces internode communication by leveraging fast intranode GPU interconnect [5]. These distributed approaches rely on a local single-GPU SpGEMM kernel; Ocean could serve as such a kernel to improve end-to-end performance. Overall, existing GPU SpGEMM solutions prioritize accumulation in shared memory. If shared memory is insufficient, they fall
ICS ’26, July 06–09, 2026, Belfast, United Kingdom
back to global memory. To our knowledge, no previous work has proposed hash-based accumulators that utilize both shared and global memory, or combined hash-based, dense, and ESC accumulators within a single design. Prior GPU SpGEMM solutions do not use estimation-based techniques to replace per-row size prediction. MOSparse sometimes avoids symbolic computation using the “upper-bound” method, but it lacks probabilistic estimation and faces two limitations: (i) potentially high GPU memory usage depending on output compression, and (ii) reliance on training data, which may not generalize to unseen matrices. Prior work has explored output size estimation for SpGEMM, but they do not target GPUs. Cohen et al. [8, 9] formulated the problem as a graph problem and applied Monte Carlo-based algorithms for fast estimation. Their approach can be applied to a chain of sparse matrix multiplications. More recently, Amossen et al. [1] improved the algorithm using pairwise independent hash functions. Du et al. [12] proposed a sampling-based estimation method for CPUs. The latter two approaches are not directly applicable to fast GPU-based SpGEMM for two main reasons: (i) their estimation procedures cannot be efficiently parallelized on GPUs, and (ii) they do not provide direct per-row output estimation, which is critical on GPUs due to scratchpad memory constraints. Dividing the per-row intermediate product count by a global compression ratio gives only a rough estimate and can lead to significant errors in computations with non-uniform compression ratios. Cohen’s estimator supports per-row output estimation and has access patterns for matrix A similar to Ocean’s HLL approach. It, however, requires multiple rounds of random access to the nodes corresponding to columns of matrix B, making it less efficient than Ocean. In addition, it does not support fast sampling or 𝐶𝑅 estimation.
7
Conclusion
SpGEMM is a critical kernel for simulations, graph analytics, and machine learning, yet most GPU implementations rely on a twopass variant of Gustavson’s algorithm to optimize accumulation. This work questions the necessity of the initial GPU pass for precise size prediction and instead proposes a fast, estimation-driven approach. In particular, we introduce HyperLogLog as a lightweight estimator for SpGEMM, demonstrating that it offers sufficient accuracy and is highly parallelizable and efficient on GPUs. In addition, we design a complete workflow, Ocean, that incorporates sampling and method selection to enable effective integration of HyperLogLog into GPU SpGEMM. Building on this estimationcentered design, we optimize accumulation kernels using sampled information and further propose a hybrid accumulator scheme. Our extensive performance comparisons and ablation studies across multiple versions of Ocean validate the effectiveness of our design choices. Overall, Ocean consistently outperforms existing implementations on more than 400 matrices. Key directions for future work include predicting accumulation kernel performance across configurations, since shared memory requirements alone provide only a lower bound on efficiency. Understanding the structure of input matrices and adopting a more fine-grained workflow selection can also lead to higher performance. Beyond SpGEMM, HyperLogLog also shows promise for tasks such as row reordering and workload characterization. Output structure
Li and Guidi
prediction is unique to SpGEMM, but HLL could also be applied to other sparse linear algebra primitives and tasks, such as guiding row reordering or characterizing workloads for SpMM.
Acknowledgments This research used resources from the National Energy Research Scientific Computing Center, a DOE Office of Science User Facility supported by the Office of Science of the U.S. Department of Energy under Contract No. DE-AC02-05CH11231, using NERSC award ASCR-ERCAP0030076. This work used DeltaAI at the National Center for Supercomputing Applications (NCSA) through allocation CIS251351 from the Advanced Cyberinfrastructure Coordination Ecosystem: Services & Support (ACCESS) program, which is supported by U.S. National Science Foundation grants #2138259, #2138286, #2138307, #2137603, and #2138296. This material is based upon work supported by the National Science Foundation under Grant IIS-2435801. The authors disclose that generative AI and editing assistants were used for grammar checking and to improve the clarity of the writing. The authors thank Julian Bellavita for reading the manuscript and providing feedback, and Qingyao Sun for contributing to the error analysis of 𝐶𝑅 sampling.
References [1] Rasmus Resen Amossen, Andrea Campagna, and Rasmus Pagh. 2014. Better Size Estimation for Sparse Matrix Products. Algorithmica 69, 3 (July 2014), 741–757. doi:10.1007/s00453-012-9692-9 [2] Pham Nguyen Quang Anh, Rui Fan, and Yonggang Wen. 2016. Balanced hashing and efficient gpu sparse general matrix-matrix multiplication. In Proceedings of the 2016 International Conference on Supercomputing. 1–12. [3] Ariful Azad, Georgios A Pavlopoulos, Christos A Ouzounis, Nikos C Kyrpides, and Aydin Buluç. 2018. HipMCL: a high-performance parallel implementation of the Markov clustering algorithm for large-scale networks. Nucleic acids research 46, 6 (2018), e33–e33. [4] Nathan Bell, Steven Dalton, and Luke N. Olson. 2012. Exposing Fine-Grained Parallelism in Algebraic Multigrid Methods. SIAM Journal on Scientific Computing 34, 4 (2012), C123–C152. doi:10.1137/110838844 [5] Julian Bellavita, Lorenzo Pichetti, Thomas Pasquali, Flavio Vella, and Giulia Guidi. 2026. Communication-Avoiding SpGEMM via Trident Partitioning on Hierarchical GPU Interconnects. In Proceedings of the 2026 International Conference on Supercomputing (ICS ’26) (Belfast, United Kingdom). ACM, New York, NY, USA, 1–13. doi:10.1145/3797905.3800543 [6] Benjamin Brock, Aydın Buluç, and Katherine Yelick. 2024. RDMA-based Algorithms for Sparse Matrix Multiplication on GPUs. In Proceedings of the 38th ACM International Conference on Supercomputing. 225–235. [7] Aydın Buluç. 2025. The Ubiquitous Sparse Matrix-Matrix Products. (2025). arXiv:2508.04077 [math.NA] https://arxiv.org/abs/2508.04077 [8] Edith Cohen. 1996. On optimizing multiplications of sparse matrices. In International Conference on Integer Programming and Combinatorial Optimization. Springer, 219–233. [9] Edith Cohen. 1997. Size-Estimation Framework with Applications to Transitive Closure and Reachability. J. Comput. System Sci. 55, 3 (Dec. 1997), 441–453. doi:10.1006/jcss.1997.1534 [10] Steven Dalton, Sean Baxter, Duane Merrill, Luke Olson, and Michael Garland. 2015. Optimizing sparse matrix operations on gpus using merge path. In 2015 IEEE International Parallel and Distributed Processing Symposium. IEEE, 407–416. [11] Zhaoyang Du, Yijin Guan, Tianchan Guan, Dimin Niu, Linyong Huang, Hongzhong Zheng, and Yuan Xie. 2022. OpSparse: A Highly Optimized Framework for Sparse General Matrix Multiplication on GPUs. IEEE Access 10 (2022), 85960–85974. doi:10.1109/ACCESS.2022.3196940 [12] Zhaoyang Du, Yijin Guan, Tianchan Guan, Dimin Niu, Nianxiong Tan, Xiaopeng Yu, Hongzhong Zheng, Jianyi Meng, Xiaolang Yan, and Yuan Xie. 2022. Predicting the Output Structure of Sparse Matrix Multiplication with Sampled Compression Ratio. doi:10.48550/arXiv.2207.13848 arXiv:2207.13848 [cs]. [13] Philippe Flajolet, Éric Fusy, Olivier Gandouet, and Frédéric Meunier. 2007. HyperLogLog: the analysis of a near-optimal cardinality estimation algorithm. Discrete Mathematics & Theoretical Computer Science DMTCS Proceedings vol. AH,..., Proceedings (Jan. 2007). doi:10.46298/dmtcs.3545 Publisher: Centre pour la Communication Scientifique Directe (CCSD).
Ocean: Fast Estimation-Based Sparse General Matrix-Matrix Multiplication on GPU
ICS ’26, July 06–09, 2026, Belfast, United Kingdom
[14] Jianhua Gao, Weixing Ji, Fangli Chang, Shiyu Han, Bingxin Wei, Zeming Liu, and Yizhuo Wang. 2023. A Systematic Survey of General Sparse Matrix-matrix Multiplication. ACM Comput. Surv. 55, 12 (Dec. 2023), 1–36. doi:10.1145/3571157 [15] John R Gilbert, Cleve Moler, and Robert Schreiber. 1992. Sparse matrices in MATLAB: Design and implementation. SIAM journal on matrix analysis and applications 13, 1 (1992), 333–356. [16] Giulia Guidi, Marquita Ellis, Daniel Rokhsar, Katherine Yelick, and Aydın Buluç. 2021. BELLA: Berkeley efficient long-read to long-read aligner and overlapper. In SIAM Conference on Applied and Computational Discrete Algorithms (ACDA21). SIAM, 123–134. [17] Giulia Guidi, Oguz Selvitopi, Marquita Ellis, Leonid Oliker, Katherine Yelick, and Aydın Buluç. 2021. Parallel string graph construction and transitive reduction for de novo genome assembly. In 2021 IEEE International Parallel and Distributed Processing Symposium (IPDPS). IEEE, 517–526. [18] Fred G Gustavson. 1978. Two fast algorithms for sparse matrices: Multiplication and permuted transposition. ACM Transactions on Mathematical Software (TOMS) 4, 3 (1978), 250–269. [19] Stefan Heule, Marc Nunkesser, and Alexander Hall. 2013. HyperLogLog in practice: algorithmic engineering of a state of the art cardinality estimation algorithm. In Proceedings of the 16th International Conference on Extending Database Technology. ACM, Genoa Italy, 683–692. doi:10.1145/2452376.2452456 [20] Abdullah Al Raqibul Islam, Helen Xu, Dong Dai, and Aydın Buluç. 2025. Improving SpGEMM Performance Through Matrix Reordering and Cluster-wise Computation. doi:10.48550/arXiv.2507.21253 arXiv:2507.21253 [cs]. [21] Zhe Jia, Marco Maggioni, Benjamin Staiger, and Daniele P. Scarpazza. 2018. Dissecting the NVIDIA Volta GPU Architecture via Microbenchmarking. doi:10. 48550/arXiv.1804.06826 arXiv:1804.06826 [cs]. [22] Donald E Knuth. 1998. The Art of Computer Programming: Sorting and Searching, Volume 3. Addison-Wesley Professional. [23] Weifeng Liu and Brian Vinter. 2014. An efficient GPU general sparse matrixmatrix multiplication for irregular data. In 2014 IEEE 28th international parallel and distributed processing symposium. IEEE, 370–381. [24] Yusuke Nagasaka, Satoshi Matsuoka, Ariful Azad, and Aydın Buluç. 2019. Performance optimization, modeling and analysis of sparse matrix-matrix products on multi-core and many-core processors. Parallel Comput. 90 (2019), 102545. [25] Yusuke Nagasaka, Akira Nukada, and Satoshi Matsuoka. 2017. High-Performance and Memory-Saving Sparse General Matrix-Matrix Multiplication for NVIDIA Pascal GPU. In 2017 46th International Conference on Parallel Processing (ICPP). IEEE, Bristol, United Kingdom, 101–110. doi:10.1109/ICPP.2017.19 [26] Yuyao Niu, Zhengyang Lu, Haonan Ji, Shuhui Song, Zhou Jin, and Weifeng Liu. 2022. TileSpGEMM: a tiled algorithm for parallel sparse general matrix-matrix multiplication on GPUs. In Proceedings of the 27th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming. ACM, Seoul Republic of Korea, 90–106. doi:10.1145/3503221.3508431 [27] NVIDIA Corporation. 2026. cuSPARSE Library. https://docs.nvidia.com/cuda/ cusparse/. Version 13.1. [28] NVIDIA Corporation. 2026. Volta Tuning Guide Release 13.1. NVIDIA Corporation. https://docs.nvidia.com/cuda/volta-tuning-guide/index.html. [29] Mathias Parger, Martin Winter, Daniel Mlakar, and Markus Steinberger. 2020. spECK: accelerating GPU sparse matrix-matrix multiplication through lightweight analysis. In Proceedings of the 25th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming. ACM, San Diego California, 362–375. doi:10.1145/3332466.3374521 [30] John G Saw, Mark CK Yang, and Tse Chin Mo. 1984. Chebyshev inequality with estimated mean and variance. The American Statistician 38, 2 (1984), 130–132. [31] Oguz Selvitopi, Saliya Ekanayake, Giulia Guidi, Georgios A Pavlopoulos, Ariful Azad, and Aydın Buluç. 2020. Distributed many-to-many protein sequence alignment using sparse matrices. In SC20: International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE, 1–14. [32] Alok Tripathy, Katherine Yelick, and Aydın Buluç. 2020. Reducing communication in graph neural network training. In SC20: International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE, 1–14. [33] Minjie Wang, Da Zheng, Zihao Ye, Quan Gan, Mufei Li, Xiang Song, Jinjing Zhou, Chao Ma, Lingfan Yu, Yu Gai, Tianjun Xiao, Tong He, George Karypis, Jinyang Li, and Zheng Zhang. 2019. Deep Graph Library: A Graph-Centric, HighlyPerformant Package for Graph Neural Networks. arXiv preprint arXiv:1909.01315 (2019). [34] Yizhuo Wang, Hongpeng Lin, Bingxin Wei, Jianhua Gao, and Weixing Ji. 2025. Optimizing General Sparse Matrix-Matrix Multiplication on the GPU. ACM Trans. Archit. Code Optim. (Nov. 2025), 3774654. doi:10.1145/3774654 [35] Nicholas Wilt. 2013. The cuda handbook: A comprehensive guide to gpu programming. Pearson Education. [36] Martin Winter, Daniel Mlakar, Rhaleb Zayer, Hans-Peter Seidel, and Markus Steinberger. 2019. Adaptive sparse matrix-matrix multiplication on the GPU. In Proceedings of the 24th Symposium on Principles and Practice of Parallel Programming. ACM, Washington District of Columbia, 68–81. doi:10.1145/3293883.3295701 [37] Martin Winter, Daniel Mlakar, Rhaleb Zayer, Hans-Peter Seidel, and Markus Steinberger. 2019. Adaptive sparse matrix-matrix multiplication on the GPU. In
Proceedings of the 24th symposium on principles and practice of parallel programming. 68–81. [38] Min Wu, Huizhang Luo, Fenfang Li, Yiran Zhang, Zhuo Tang, Kenli Li, Jeff Zhang, and Chubo Liu. 2025. HSMU-SpGEMM: Achieving High Shared Memory Utilization for Parallel Sparse General Matrix-Matrix Multiplication on Modern GPUs. In 2025 IEEE International Symposium on High Performance Computer Architecture (HPCA). IEEE, Las Vegas, NV, USA, 1452–1466. doi:10.1109/HPCA61900.2025. 00109