---
title: Estimation-Based SpGEMM Workflow
url: https://www.emergentmind.com/topics/estimation-based-spgemm-workflow
type: topic
---

# Estimation-Based SpGEMM Workflow

Searching arXiv for relevant papers on estimation-based SpGEMM workflows and related algorithm-selection work.
Estimation-based SpGEMM workflow denotes a family of sparse general matrix multiplication procedures in which structural properties of the product matrix are predicted or estimated before, or in lieu of, exact symbolic computation. In sparse matrix multiplication $C = A \cdot B$, the principal target of estimation is typically the output structure of $C$, especially the number of nonzeros per row or per column, because memory allocation, load balancing, and accumulator selection depend on it. Across CPU, GPU, vector, and accelerator settings, the workflow replaces or supplements exact symbolic phases with lightweight predictors derived from intermediate-product counts, sampled compression ratios, HyperLogLog sketches, or work metrics such as per-column operation counts. The approach is motivated by the cost of exact symbolic traversal and by the fact that SpGEMM performance is strongly conditioned on irregular sparsity, accumulator pressure, and communication patterns [2207.13848], [2604.19004], [1504.05022], [2303.02471], [1603.05627], [2406.10166].

## 1. Definition and motivating problem

Sparse general matrix multiplication (SpGEMM) computes $C = A \cdot B$ for sparse matrices and is a fundamental building block in scientific computing, graph analytics, and machine learning [2207.13848], [2303.02471], [2604.19004], [1504.05022]. A central difficulty is that the output sparsity pattern is unknown a priori. In CSR-based row-wise formulations, the critical quantity is the number of nonzeros per output row of $C$; in CSC-based column-wise vector formulations, the analogous quantity is the work or structural complexity per output column [2207.13848], [2303.02471].

The motivation for estimation is operational rather than merely analytical. Accurate prediction enables memory allocation for CSR row pointers and per-row buffers, load balancing across CPU threads or GPU warps and CTAs, and kernel selection among hash-based, merge-based, sorting-based, dense, or hybrid accumulators [2207.13848], [2604.19004]. In vector architectures, work estimation supports dynamic switching among SPA, SPARS, and HASH, while in communication-centric settings it guides the selection of 1D, 2D, or 3D algorithm families [2303.02471], [1603.05627]. In accelerator-oriented dataflow selection, estimation extends beyond output size to operational and traffic surrogates that distinguish inner-product, outer-product, and row-wise dataflows [2406.10166].

The immediate alternative to estimation is exact symbolic computation. Exact symbolic methods determine the output structure precisely but often perform a traversal similar in complexity to the numeric phase, thereby adding substantial overhead when used strictly for sizing or scheduling [2207.13848]. On GPUs, the symbolic pass alone accounts for roughly $28\%$ of total runtime on average in state-of-the-art baselines such as spECK, even though it produces only one number per row [2604.19004]. Upper bounds based on intermediate products are cheap, but can be extremely loose, particularly when the compression ratio is high [2207.13848]. Estimation-based workflows occupy the design space between these two extremes.

## 2. Structural quantities and estimator targets

The quantities estimated in SpGEMM workflows vary by execution model, but they all describe latent structure or work in the product.

For row-wise CSR SpGEMM, let $A \in \mathbb{R}^{m \times k}$, $B \in \mathbb{R}^{k \times n}$, and $C = A \cdot B$. The row-wise dataflow forms
$$
C_{i,:} = \sum_{k \in \mathrm{NZ}(A_{i,:})} A_{i,k} \cdot B_{k,:}.
$$
A fundamental quantity is the row-level intermediate-product count
$$
\mathrm{FLOP}_i = \sum_{k \in \mathrm{NZ}(A_{i,:})} \mathrm{nnz}(B_{k,:}),
$$
with total
$$
\mathrm{FLOP}_{\text{total}} = \sum_{i=1}^{m} \sum_{k \in \mathrm{NZ}(A_{i,:})} \mathrm{nnz}(B_{k,:}).
$$
The corresponding output compression ratio is
$$
\mathrm{CR} = \frac{\mathrm{FLOP}_{\text{total}}}{\mathrm{nnz}(C)}.
$$
These definitions underpin sampled compression-ratio estimation on CPUs and upper-bound workflows on GPUs [2207.13848], [1504.05022], [2604.19004].

In the GPU formulation of Ocean, the structural target for row $i$ is
$$
c_i = \left| \bigcup_{j \in \mathcal{N}_A(i)} \mathcal{N}_B(j) \right|,
$$
the number of distinct column indices reachable by two-hop expansion through the nonzeros of row $i$ of $A$ into rows of $B$ [2604.19004]. This is precisely the quantity needed for per-row capacity and bin assignment in Gustavson-style numeric kernels.

In the vector-architecture workflow using CSC storage, the principal estimator is per-column work:
$$
Op_j = \sum_{k \in \mathrm{supp}(B(:,j))} Z^A_k,
$$
where $Z^A_k$ is the number of nonzeros in column $k$ of $A$ [2303.02471]. Here the estimation target is not directly $\mathrm{nnz}(C(:,j))$ but the intermediate-product count governing whether SPA, SPARS, or HASH is preferable.

The 2015 GPU framework uses
$$
W_i = \sum_{p \in \mathrm{nz}(A_{i,:})} \mathrm{nnz}(B_{p,*}),
$$
which serves simultaneously as a work estimate and an upper bound on $\mathrm{nnz}(C_i)$ via $\mathrm{nnz}(C_i) \le \min(n, W_i)$ [1504.05022]. This bound is structurally loose but easy to compute in $O(\mathrm{nnz}(A))$ time.

A broader implication is that estimation-based workflows can be organized around three estimator families: direct size prediction, work prediction, and communication prediction. That taxonomy is not stated explicitly in the sources, but it is consistent with the sampled-$\mathrm{CR}$ method [2207.13848], the $Op_j$ vector work model [2303.02471], and the hypergraph communication framework [1603.05627].

## 3. Ratio-based estimation of output structure

A prominent estimation-based workflow is the sampled compression-ratio method for CSR SpGEMM proposed in "Predicting the Output Structure of Sparse Matrix Multiplication with Sampled Compression Ratio" [2207.13848]. Its central idea is to estimate $\mathrm{CR}$ from a sampled result computed on a small uniform random subset of rows of $A$, while keeping all rows of $B$ visible. For a sampled row set $\mathcal{S}$,
$$
\mathrm{FLOP}_s = \sum_{i \in \mathcal{S}} \mathrm{FLOP}_i,\qquad
\mathrm{nnz}(C_s) = \sum_{i \in \mathcal{S}} \mathrm{nnz}(C_{i,:}),
$$
and the sample-based compression-ratio estimator is
$$
\widehat{\mathrm{CR}} = \frac{\mathrm{FLOP}_s}{\mathrm{nnz}(C_s)}.
$$
Predicted row lengths then follow as
$$
\widehat{\mathrm{nnz}(C_{i,:})} = \frac{\mathrm{FLOP}_i}{\widehat{\mathrm{CR}}},
$$
with total prediction
$$
\widehat{\mathrm{nnz}(C)} = \frac{\mathrm{FLOP}_{\text{total}}}{\widehat{\mathrm{CR}}}.
$$
The workflow computes $\mathrm{FLOP}_i$ by precomputing the row degrees $d_k = \mathrm{nnz}(B_{k,:})$ and summing them over the nonzeros in each row of $A$ [2207.13848].

The distinguishing argument is error neutralization. If the sampled FLOP and sampled nonzero counts are perturbed by relative errors $\varepsilon_f$ and $\varepsilon_1$, then the induced relative error in the estimated total is
$$
\varepsilon_2 = \frac{\varepsilon_1 - \varepsilon_f}{1+\varepsilon_f}.
$$
When $\varepsilon_f \approx \varepsilon_1$, the estimator error becomes small. The paper reports a positive correlation between sampled FLOP and sampled nnz with correlation coefficient $97.01\%$ across $625$ test cases, which supports the ratio-based design [2207.13848].

The sampling strategy uses a small, uniform random subset of rows of $A$ with
$$
\text{sample\_num} = \min(0.003 \cdot m, 300),
$$
and computes $\mathrm{nnz}(C_s)$ exactly on the sample using a per-row hash-based accumulator with linear probing [2207.13848]. Because both matrices are typically in CSR, row sampling is preferred over column sampling, and row-wise dataflow is described as faster and simpler than inner-product dataflow for CSR inputs [2207.13848].

The empirical results are specific. On $625$ test cases derived from $25$ representative SuiteSparse matrices, the absolute relative error of the proposed method is $1.56\%$ on average and $25\%$ in the worst case, versus $8.12\%$ average and $156\%$ worst-case for a reference sampling design that scales sampled nnz only and ignores sampled FLOP [2207.13848]. The parallel implementation adds on average $0.78\%$ of execution time relative to BRMerge-Precise; on matrix-square benchmarks, computing $\mathrm{FLOP}_i$ costs $1.68\%$ on average, up to $4.12\%$, and the prediction step costs $0.72\%$ on average, up to $1.89\%$ [2207.13848].

The workflow then integrates prediction into downstream execution. It allocates CSR with
$$
\mathrm{row\_ptr}[i+1] = \mathrm{row\_ptr}[i] + \lceil \beta \cdot \widehat{\mathrm{nnz}(C_{i,:})} \rceil,
$$
where $\beta \ge 1$ is a small guard factor such as $\beta = 1.05$–$1.10$, and uses predicted row sizes to choose hash-based, merge-based, or sorting-based accumulation and to distribute rows proportionally to $\widehat{\mathrm{nnz}(C_{i,:})}$ or $\mathrm{FLOP}_i$ [2207.13848].

## 4. Upper-bound and hybrid pre-allocation workflows

An earlier and influential GPU framework estimates structure using per-row work bounds rather than direct output-size prediction [1504.05022]. In "A Framework for General Sparse Matrix-Matrix Multiplication on GPUs and Heterogeneous Processors" [1504.05022], all stages use CSR, and the symbolic estimation phase computes $U[i] := W_i$ on the GPU with one thread per output row. Because $W_i$ depends only on row lengths read from $B$’s row pointers, the estimation cost is $O(\mathrm{nnz}(A))$ [1504.05022].

The framework then performs host-side binning into $38$ bins organized into five groups:

| Bin group | Range of $U[i]$ | Allocation policy |
|---|---:|---|
| 1 | $U=0$ | Empty rows |
| 2 | $U=1$ | Single-entry rows |
| 3 | $U \in [2,32]$ | Exact-size upper bound per row |
| 4 | $U \in [33,64], [65,128], [129,256], [257,512]$ | Exact-size upper bound per row |
| 5 | $U > 512$ | Initial capacity $256$, then progressive growth |

For bin groups $1$–$4$, the temporary CSR matrix $\widetilde{C}$ allocates $nnz(\widetilde{c}_{i,*}) = U[i]$. For bin group $5$, capacity is initialized to $256$ and grown by doubling when needed [1504.05022]. This is described as a hybrid strategy interpolating between upper-bound allocation and progressive allocation. The rationale is explicit: full upper-bound allocation at the intermediate size $nnz(\widehat{C})$ can exceed device memory by orders of magnitude on irregular matrices, whereas precise pre-computation doubles the work and time [1504.05022].

The estimator is integral to algorithm dispatch. Rows with $U \le 32$ use a heap-based shared-memory accumulator. Rows with $33 \le U \le 512$ use bitonic ESC (Expansion–Sorting–Compression). Rows with $U > 512$ use merge-based insertion with GPU merge path, with checkpointing, dumping to global memory, re-allocation by $2\times$, and reload for long rows [1504.05022]. Load balancing follows the same estimate: short rows use one thread per row, moderate and long rows use one thread group per row, and per-bin counters avoid launching empty kernels [1504.05022].

This workflow does not estimate $\mathrm{nnz}(C_i)$ tightly. Instead, it treats arithmetic work as a safe upper bound and uses bin-specific accumulation strategies to absorb overestimation. A plausible implication is that it belongs to the same estimation-based lineage as later ratio-based and HLL-based methods, but emphasizes robust scheduling and memory control over prediction accuracy. That interpretation is consistent with the later critique that upper bounds can be “extremely loose” when compression ratios are large [2207.13848].

The framework also includes a heterogeneous-processor variant using re-allocatable shared virtual memory on APUs. For matrices that need re-allocation, this yields average $1.2\times$ and up to $1.8\times$ speedup on AMD A10-7850K by enabling in-place growth and shrink without copy [1504.05022].

## 5. Sketch-based estimation and workflow selection on GPUs

The GPU system Ocean introduces a more elaborate estimation-based symbolic workflow that replaces most exact symbolic accumulation with HyperLogLog (HLL) sketches [2604.19004]. Ocean constructs one HLL sketch per row of $B$, using $m = 2^p$ registers with $b=8$ bits per register. It uses $m \in \{32, 64\}$ adaptively: $p=5$ for $m=32$ when the input expansion ratio $ER < 48$, otherwise $p=6$ for $m=64$ [2604.19004].

For each row $i$ of $A$, Ocean merges the sketches of the rows of $B$ referenced by $\mathcal{N}_A(i)$ using elementwise maxima and estimates the union cardinality with the standard raw HLL estimator
$$
\hat{E} = \alpha_m \cdot m^2 \cdot \left(\sum_{r=1}^{m} 2^{-M_r}\right)^{-1},
$$
where $\alpha_m \approx 0.7213/(1 + 1.079/m)$ [2604.19004]. Linear counting is used when $\hat{E} < 2.5m$, and large-range corrections are enabled when $\hat{E} \gtrsim 2^{32}/30$ [2604.19004].

Ocean embeds HLL into a matrix-level workflow selector driven by two ratios. The input expansion ratio is
$$
ER = \frac{\sum_i \sum_{j \in \mathcal{N}_A(i)} |\mathcal{N}_B(j)|}{\mathrm{nnz}(A)},
$$
and the output compression ratio is
$$
CR = \frac{\sum_i \sum_{j \in \mathcal{N}_A(i)} |\mathcal{N}_B(j)|}{\mathrm{nnz}(C)}.
$$
Ocean estimates $CR$ cheaply via sampling: it builds HLL sketches for all rows of $B$ and randomly samples $3\%$ of rows of $A$, with minimum $600$ and maximum $10{,}000$, to merge their sketches and estimate approximate $c_i$ values [2604.19004].

The workflow then chooses among three modes:

| Workflow | Selection criterion |
|---|---|
| Upper-bound estimation | Average intermediate products per row $< 64$ |
| HLL estimation | Average intermediate products per row $\ge 64$, $ER \ge 8$, and sampled $CR \ge 8$ |
| Exact symbolic | Otherwise |

The rationale is architectural. Each HLL register is a byte, each CSR column index is $4$ bytes, and $ER \approx 8$ makes merging $m \approx 32$ registers per referenced row about as cheap as reading approximately $32$ symbolic indices [2604.19004]. Likewise, $CR \ge 8$ is used because estimation-based workflows delay the final CSR layout, and late compaction overhead grows when compression is low [2604.19004].

Ocean treats estimator uncertainty through allocation slack. Per-row capacity is
$$
cap_i = \mathrm{round\_up\_to\_bin}((1+S)\cdot \hat{c}_i),
$$
with default $S=0.5$, and $S$ may be raised to $1.0$ when $m=32$ [2604.19004]. The relative standard error of HLL is approximated by $1.04/\sqrt{m}$, yielding about $18.4\%$ for $m=32$ and $13.0\%$ for $m=64$ [2604.19004]. Measured average relative per-row errors are about $0.13$, $0.10$, and $0.07$ for $m \in \{32,64,128\}$, with overflow rates after binning with $1.5\times$ expansion and upward rounding of $1.2\%$, $0.3\%$, and less than $0.1\%$, respectively [2604.19004].

The numeric phase uses the estimator to guide accumulator choice. Ocean provides hash-based accumulators with hybrid memory, dense accumulators, and ESC accumulators, each with multiple bins [2604.19004]. Dense is preferred over hash when both fit the same bin, except in the largest bin where hybrid hash is favored; ESC is used only in the upper-bound workflow [2604.19004]. Overflow rows are handled by a fallback kernel using the largest dense accumulator configuration and a global buffer sized by the row’s intermediate product count [2604.19004].

The performance data are end-to-end. On $337$ square matrices and $64$ rectangular matrices from SuiteSparse, Ocean achieves on A100 a geometric mean of $63.7$ GFLOP/s for square cases and $36.3$ GFLOP/s for rectangular cases, and on H100 $108.0$ GFLOP/s for square cases [2604.19004]. It wins on $294/337$ square matrices on A100 and $312/337$ on H100, with average speedups of $1.4\times$–$2.8\times$ over leading GPU baselines across the study; more specifically, on A100 square matrices the average speedups over cuSPARSE, spECK, opSparse, TileSpGEMM, and HSMU are $18.8\times$, $1.4\times$, $2.6\times$, $3.5\times$, and $2.0\times$, respectively [2604.19004]. The estimation step costs approximately $4\%$ of runtime for selected matrices, compared with symbolic-step costs of about $28$–$30\%$ in spECK or the baseline Ocean variant [2604.19004].

## 6. Estimation for kernel, dataflow, and communication selection

Estimation-based SpGEMM workflows are not confined to output-size prediction. In several settings, the estimated quantity is used to select among algorithmic dataflows, accumulator families, or communication decompositions.

On RISC-V vector architectures, the workflow in "Optimization of SpGEMM with Risc-V vector instructions" [2303.02471] computes $Op_j$ for every output column and sorts columns by decreasing work. This work estimate drives dynamic switching between SPA, SPARS, and HASH, yielding hybrid algorithms H-SPA$(t)$ and H-HASH$(t)$ governed by a threshold $t$ [2303.02471]. If $Op_j < t$, the system uses SPARS or HASH inside blocks formed to homogenize column work; otherwise it uses SPA [2303.02471]. The threshold guidance comes from sensitivity experiments: $t$ in $[40,60]$ yields the best average, and $t=40$ is used in the main evaluation [2303.02471]. H-SPA$(t)$ and H-HASH$(t)$ obtain average speedups of $1.24\times$ and $1.57\times$ over SPA on $40$ SuiteSparse matrices, and $1.42\times$ and $1.99\times$ on the $22$ most sparse matrices [2303.02471].

In Misam, estimation is fused with machine learning for dataflow selection among inner-product, outer-product, and row-wise SpGEMM [2406.10166]. The system extracts light-weight structural features such as sparsityA, sparsityB, avg_row_lengthA, avg_row_lengthA_var, avg_col_lengthB, and blocks_accessed, then predicts the best dataflow per block using a decision tree or DQN [2406.10166]. The decision tree, limited to depth $9$, achieves $94\%$ selection accuracy; the DQN achieves $90\%$ [2406.10166]. Relative to any single fixed dataflow, the decision tree produces average speedups of $2.7\times$ over IP, $2.1\times$ over OP, and $2.64\times$ over RW, with $1.13\times$ improvement over the heuristic baseline [2406.10166]. The abstract reports gains of up to $28$ times [2406.10166]. Although Misam frames the problem as dataflow choice rather than symbolic estimation, it exemplifies the same principle: low-cost structural surrogates are used to avoid committing to a single static execution strategy.

At a larger scale, the hypergraph framework for SpGEMM estimates communication and memory-hierarchy traffic from sparsity patterns [1603.05627]. The fine-grained model associates computation vertices with nontrivial products and data vertices with nonzeros of $A$, $B$, and $C$, then encodes dependencies as hyperedges [1603.05627]. Estimation in this setting concerns communication volume along the critical path, total connectivity, and sequential I/O lower bounds. The workflow computes sparsity features, estimates flops $f = \sum_k \mathrm{nnz}(A(:,k)) \cdot \mathrm{nnz}(B(k,:))$, constructs coarse 1D or 2D hypergraph models, partitions or approximates them, and selects the algorithm family minimizing predicted communication subject to load balance [1603.05627]. The paper’s empirical conclusion is that at least one coarse-grained model is as communication efficient as the fine-grained model in each of the three application domains studied, but different applications favor different families [1603.05627].

Taken together, these works suggest that “estimation-based SpGEMM workflow” is best understood as a general control architecture: estimate a structural proxy, map it to a resource or algorithm decision, then execute the numeric phase under that guidance. The proxy may be $\widehat{\mathrm{nnz}}$, $\widehat{CR}$, $Op_j$, feature vectors for ML classifiers, or hypergraph cuts.

## 7. Accuracy limits, guardrails, and open directions

The principal limitation of estimation-based workflows is that their reliability depends on how well the estimator captures overlap and skew in the product structure. In the sampled compression-ratio method, when the correlation between sampled FLOP and sampled nnz weakens, error neutralization diminishes; the observed worst-case error reaches $25\%$ in extreme cases [2207.13848]. The paper identifies rows with extreme degree skew and duplicate aggregation patterns as contributors to residual error, and suggests adaptive sampling, importance sampling proportional to $\mathrm{FLOP}_i$, and per-cluster compression-ratio estimation as future directions [2207.13848].

Ocean incorporates explicit guardrails for such uncertainty [2604.19004]. It uses slack factors, upward bin rounding, and a robust overflow kernel. Under-estimation is therefore treated as a performance issue rather than a correctness issue: overflow rows are rerun with safe upper-bound capacity, while over-estimation increases compaction cost modestly but does not violate CSR correctness [2604.19004]. This suggests a design principle common to modern workflows: approximate sizing is acceptable when backed by bounded fallback mechanisms.

Upper-bound workflows have the opposite problem. They are safe but can over-allocate severely when the compression ratio is large, consuming memory and skewing binning or scheduling [2207.13848], [1504.05022]. Ocean formalizes part of this trade-off with $ER$ and $CR$ thresholds, choosing exact symbolic whenever late compaction would not amortize or when HLL merging would not be traffic-efficient [2604.19004].

Work-estimation workflows for vector architectures face different failure modes. The synthetic cross-overs reported for SPA versus SPARS and HASH depend on assumptions about uniform sparsity and sorted column degrees, and the paper notes that the general case deviates depending on column nnz distributions [2303.02471]. Hence, $Op_j$ is a useful but not exhaustive predictor. Sorting and blocking reduce imbalance but do not eliminate it [2303.02471].

ML-based selectors such as Misam face hardware dependence and distribution shift [2406.10166]. The learned efficiencies and feature-action mappings are tied to buffer sizes, interconnect, and PE scheduling, so retraining per accelerator configuration is advisable [2406.10166]. That caveat is consistent with the broader literature: estimated structure is platform-agnostic only at a coarse level, while thresholding, accumulator choice, and occupancy interactions are platform-specific.

Several future directions recur across the sources. These include adaptive or stratified sampling for nonuniform compression ratios [2207.13848]; alternative sketches with lower relative standard error at equal bandwidth, hierarchical sketches, and better memory-footprint modeling for multi-GPU settings [2604.19004]; predictive models for kernel efficiency beyond shared-memory capacity on vector hardware [2303.02471]; and tighter theoretical error bounds or uncertainty-aware fallbacks that bridge statistical estimation and resource provisioning [2207.13848], [2604.19004]. A plausible implication is that future estimation-based SpGEMM systems will combine multiple estimators—sampling, sketches, analytical bounds, and learned predictors—inside a unified controller rather than relying on a single proxy.

In aggregate, the literature defines an estimation-based SpGEMM workflow as a multistage procedure in which symbolic exactness is selectively replaced by predictive structure analysis, with the resulting estimates steering memory layout, load balancing, accumulator design, kernel dispatch, dataflow selection, and, in distributed settings, communication decomposition. The technical diversity of these workflows reflects the fact that the unknown output structure of sparse matrix products is not merely a combinatorial nuisance; it is the principal control variable of high-performance SpGEMM.

Source: https://www.emergentmind.com/topics/estimation-based-spgemm-workflow