Ocean: Fast Estimation-Based Sparse General Matrix-Matrix Multiplication on GPU

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.4x-2.8x on NVIDIA A100 and H100 architectures.

• The paper introduces an estimation-driven SpGEMM framework using HyperLogLog to replace the traditional symbolic pass, reducing runtime overhead by up to 1.3×.
• It employs a dynamically adaptive workflow and hybrid accumulator design that selects optimal kernel configurations for various row lengths.
• Experimental results demonstrate significant performance gains, with 1.4×–2.8× speedups on NVIDIA A100 and H100 GPUs across diverse sparse matrices.

Ocean: Estimation-Driven SpGEMM on GPUs via HyperLogLog and Hybrid Accumulation

Introduction

Sparse General Matrix-Matrix Multiplication (SpGEMM) is a fundamental kernel in scientific computing, graph analytics, and machine learning. The challenge in high-performance SpGEMM on GPU architectures arises from highly irregular data access patterns and unpredictable output structure due to sparsity in input matrices. State-of-the-art GPU SpGEMM algorithms rely on a two-pass workflow—symbolic and numeric—where the symbolic step computes per-row output sizes and directly impacts load balancing and scratchpad allocation. However, the symbolic pass constitutes, on average, 28% of the total runtime in leading implementations, despite only producing a single value per row. This work introduces "Ocean," an estimation-driven SpGEMM framework that replaces the exact symbolic computation phase with scalable, probabilistically accurate per-row output size estimation using HyperLogLog (HLL) estimators. Ocean also proposes an analysis-driven, dynamically adaptive workflow, and a hardware-conscious hybrid accumulation design. The approach consistently outperforms state-of-the-art GPU SpGEMM solutions across both square and rectangular matrices, with reported speedups of $1.4{\times}$–$2.8{\times}$ on NVIDIA A100 and H100 platforms (2604.19004).

SpGEMM on GPU and Workflow Analysis

Efficient GPU SpGEMM uses Gustavson’s algorithm, leveraging the independence of output rows for parallelism and employing accumulators (hash-based or dense) to handle intermediate accumulation (Figure 1). The prevalent workflow consists of three phases: analysis, symbolic, and numeric, each entailing kernel launches and scratchpad memory allocations, guided by per-row structure analysis. Binning—mapping rows to kernels with suitable capacities—depends critically on precise per-row size predictions.

Figure 1: Computation of a single output row using Gustavson’s algorithm with both hash-based and dense accumulator implementations.

A critical constraint is that the numeric phase requires prior knowledge of the output row sizes for optimal memory allocation and minimizing post-compaction overhead due to CSR format demands. When the symbolic step's prediction cost outweighs benefits, the design’s efficiency is fundamentally compromised. Furthermore, existing accumulator designs are not uniformly efficient for both short (less than warp-sized) and very long (exceeding scratchpad limits) rows.

Figure 2: Common GPU SpGEMM workflow, highlighting the standard two-pass (symbolic and numeric) approach.

Estimation-Based Output Prediction Using HyperLogLog

Ocean substitutes the symbolic step with per-row estimation via HyperLogLog probabilistic sketches. HLL is appropriate for three reasons: (1) error can be tuned via the number of registers; (2) memory requirements are constant, enabling robust handling of long rows; (3) sketch updates are inherently parallelizable using atomic instructions.

Ocean constructs an HLL sketch for each row of input matrix $B$ and merges the relevant sketches according to the input pattern of $A$ rows (i.e., for row $i$ of $A$, the sketches corresponding to the nonzero columns in row $i$ are merged). This construct-and-merge pattern transforms the irregular and memory-bound symbolic pass into cache-friendly, regular operations (Figure 3), resulting in high estimation performance and parallel scalability.

Figure 3: Per-row output nonzero estimation using HLL construct-and-merge: sketches for each $B$ row are merged based on $A$’s nonzeros, giving an estimation for each $C$ row.

Ocean Workflow: Dynamic Analysis, Overflow Handling, and Workflow Selection

Ocean generalizes the workflow, orchestrating an initial analysis to dynamically decide between estimation-based, symbolic-based, or upper-bound-only workflows depending on matrix characteristics (Figure 4). Overflow—when HLL underestimates output size—is handled via fallback accumulation kernels operating on large dense configurations. The key statistics used for workflow selection are:

• Input Expansion Ratio (ER): Ratio of intermediate products to nonzeros in $2.8{\times}$0.
• Output Compression Ratio (CR): Ratio of intermediate products to output nonzeros in $2.8{\times}$1.

Sampling is used to estimate these ratios efficiently (Figure 5). If ER and sampled CR surpass specified thresholds, estimation-based execution is chosen; otherwise, symbolic or upper-bound approaches are employed.

Figure 4: Overview of Ocean’s SpGEMM pipeline; estimation-based stages (HLL estimation, overflow fallback, analysis-driven binning) are highlighted.

Figure 5: Calculation of Input Expansion Ratio (ER) and Output Compression Ratio (CR) for workflow selection.

Hybrid Accumulator Design and Kernel Selection

Ocean’s hybrid accumulator system improves both at the long- and short-row extremes:

• Enhanced Hash-Based Accumulator: For very long rows, Ocean uses a design where only indices remain in scratchpad, and values are stored in global memory, exploiting the fact that atomicAdd on global memory is not performance-prohibitive for these workloads due to fire-and-forget access patterns and the nature of update semantics.
• ESC Accumulator: For extremely short rows, "expand-sort-compact" (ESC) handling is employed for maximal memory access efficiency.

The accumulator type (hash, dense, or ESC) and kernel configuration are selected per row according to estimated statistics and binning outcomes.

Performance Results and Analysis

Ocean’s implementation was benchmarked on more than 400 matrices, encompassing both structure-rich and highly irregular patterns (Figure 6, Figure 7). The results are as follows:

• On the A100, across 337 square matrices, Ocean is best on 86% and second-best on an additional 11%. The average geometric mean performance is 63.7 GFLOPS, consistently outperforming state-of-the-art libraries (spECK, opSparse, TileSpGEMM, HSMU-SpGEMM) by factors up to $2.8{\times}$2.
• HLL-based estimation yields average per-row nonzero prediction errors below 13% for 32-register configurations, decreasing with more registers, and overflow handling is required in less than 1% of cases for 64-register settings.
• Ablation studies show the estimation-based workflow is responsible for the majority of the speedup, supported by the hybrid accumulation and adaptive analysis steps.

  Figure 6: Smoothed GFLOPS across square matrices, ordered by intermediate product count, showing Ocean’s robust scaling and absolute performance advantage.

  

Figure 7: Ocean’s speedup over spECK across (left) square matrices (H100) and (right) rectangular matrices (A100).

Estimation error and overflow distributions are visualized in Figure 8 and Figure 9.

Figure 8: CDF of average per-row HLL estimation error across numbers of registers.

Figure 9: CDF of overflow fraction across matrices for different HLL register sizes.

Theoretical and Practical Implications

Ocean introduces probabilistic estimation—previously unexplored at this scale and granularity—for SpGEMM on GPUs, fundamentally altering the cost trade-off for output structure prediction. It demonstrates that the uncertainty inherent in the HLL estimator is amortized by robust fallback mechanisms and strategically chosen expansion factors for accumulator sizing, resulting in rare and low-impact overflows. The hybrid accumulator design shows that judicious partitioning between on-chip and global memory can yield significant performance improvements, contrary to prior assumptions that value movement in global memory is prohibitive for SpGEMM.

On the theoretical side, this approach reframes output size prediction as a probabilistic streaming problem, making Ocean’s ideas applicable to a broader range of sparse algebraic primitives and analytics workloads—such as dynamic reordering and workload estimation. In practice, Ocean’s open-source implementation provides a performance-robust SpGEMM kernel for both single-GPU and as a local kernel in distributed SpGEMM systems.

Limitations and Future Work

While estimation drives strong average-case performance and reduces symbolic step overhead by up to $2.8{\times}$3, Ocean can incur higher memory overheads for some input configurations due to out-of-place compaction and peak memory reservation. Integrating finer-grained memory management or predictive switching to symbolic mode for memory-constrained cases could further improve robustness.

Key directions for further research include:

• Extending the estimation pipeline and hybrid kernel selection to other sparse primitives and chain products.
• Refining workflow selection with deeper matrix structural statistics and cost modeling for kernel configurations beyond scratchpad requirements.

Conclusion

Ocean establishes estimation-driven SpGEMM as a high-performance, scalable alternative on modern GPU architectures. By fusing probabilistic cardinality estimation with hybrid accumulation, it systematically sidesteps the heavy overhead of traditional symbolic passes. Rigorous evaluation and ablation on diverse real-world matrices validate both the conceptual and engineering advances of the framework.

Ocean’s contributions—estimation-driven workflow, dynamic adaptive analysis, and memory-aware hybrid accumulation—offer new paradigms for implementing irregular scientific kernels on massively parallel architectures.

Reference:

"Ocean: Fast Estimation-Based Sparse General Matrix-Matrix Multiplication on GPU" (2604.19004)