SpSYRK: Half the Work in Distributed Sparse Matrix Multiplication

arXiv:2608.09713v1 [cs.DC] 10 Aug 2026

SpSYRK: Half the Work in Distributed Sparse Matrix Multiplication Thomas McFarland

Julian Bellavita

Giulia Guidi

Cornell University Ithaca, NY, USA tfm62@cornell.edu

Cornell University Ithaca, NY, USA jb2695@cornell.edu

Cornell University Ithaca, NY, USA gg434@cornell.edu

Abstract—The symmetric rank-k update (SYRK), C = AA⊤ , computes the dot product between each pair of rows of A, producing the Gram matrix C. Its sparse variant underpins similarity search in machine learning, graph analytics, and genomics, including Jaccard similarity on datasets too large for a single node. Yet, despite the symmetry in its inputs and outputs, existing distributed sparse matrix multiplication algorithms, such as Sparse SUMMA, treat sparse SYRK as generic multiplication, leaving performance untapped. In this paper, we present distributed sparse SYRK approaches that leverage symmetry. The approach partitions the off-diagonal blocks of the output between the upper and lower triangular portions of the process grid and computes only the lowertriangular part of each diagonal block, reducing per-process communication and computation compared with state-of-the-art distributed SpGEMM. A second variant reorders communication to avoid materializing A⊤ . On 32 nodes of the Perlmutter supercomputer, the algorithm achieves a 2× speedup over an optimized Sparse SUMMA on matrices where local multiplication dominates the runtime; the advantage narrows on communication-bound inputs, a dependence the cost model predicts from the arithmetic intensity. Our variant rectifies this and consistently achieves superior scaling at high process counts. The approach is a drop-in replacement for any application computing C = AA⊤ via a distributed SpGEMM routine, and its triangular output can be consumed directly by subsequent operations, reducing both compute and memory footprint.

I. I NTRODUCTION Computing the product of a large sparse matrix with its transpose, C = AAT for A ∈ Rm×n , is a core operation in many scientific and engineering applications. This operation is known as the sparse symmetric rank-k update, or sparse SYRK. It shows up in computing Jaccard similarity [1], hypergraph partitioning [2], triangle counting [3], the normal equations AT A in least-squares problems, Gram matrix construction in statistics and machine learning, genome assembly [4], [5], and protein sequence alignment and similarity search [6]–[8]. In practice, A often has hundreds of millions of nonzeros, and the output C typically has 1–2 orders of magnitude more nonzeros than A. In these settings, distributedmemory parallelism is needed, and applications typically use a distributed-memory sparse general matrix-matrix multiply (SpGEMM) routine to perform sparse SYRK. This strategy fails to exploit the fact that C is symmetric. In particular, off-the-shelf distributed SpGEMM routines, such

as CombBLAS [9], compute the entire output matrix C, even though only the lower or upper triangular portion of C is strictly necessary. They also materialize the explicit transpose AT , incurring additional communication and memory costs that, as we show, are avoidable. In this paper, we develop and implement SpSYRK, the first distributed-memory parallel algorithm for sparse SYRK. SpSYRK uses a 2D partitioning scheme similar to the Sparse SUMMA algorithm [10] for distributed SpGEMM, but it exploits output symmetry by splitting the computation of offdiagonal blocks between processes in the upper and lower triangular regions of the process grid, while diagonal processes compute only the lower triangular portion of their block, theoretically reducing computation by 1/2. A key component is CommSpSYRK, a communicationreordered variant that eliminates the explicit materialization of A⊤ . Rather than transposing and redistributing the input, CommSpSYRK reorders the communication so that each process receives only the submatrices it needs. This variant has the lowest communication volume of the schemes presented and delivers the best performance at scale. The regions of the input matrices that must be communicated to each process are smaller, leading to reduced per-process communication by as much as a factor of 3/4. The algorithm can also materialize the full output matrix from the explicitly computed half, which is necessary for downstream tasks such as the overlap detection component of ELBA [5]. This materialization is often wasteful, though; many applications work directly with the triangular output, which always offers a memory advantage and, in specific pipelines, a computational one. The evaluation compares strong scaling against the distributed SpGEMM routine in CombBLAS [9] and an optimized Sparse SUMMA implementation that is usually faster than the CombBLAS version. Both were included to provide an independent algorithm and a SUMMA algorithm using the same kernels. On 32 nodes of NERSC Perlmutter, SpSYRK achieves twice the performance of optimized Sparse SUMMA on compute-bound inputs, and both SpSYRK and CommSpSYRK are up to 10× faster than CombBLAS. On bandwidthbound inputs, CommSpSYRK maintains a 2× speedup and, in most cases, achieves a 3× speedup. The main contributions are: 1) SpSYRK, the first distributed-memory algorithm for

Algorithm 1 Sparse SUMMA from the view of process (i, j) Require: A and B distributed on a pr × pc process grid Ensure: Matrix C = AB 1: (i, j) ← grid indices of this process 2: Initialize Clocal = 0 3: for stage = 0 → pc − 1 do 4: if j = stage then 5: Broadcast Alocal within process row 6: else 7: Receive Arecv from process (i, stage) 8: end if 9: if i = stage then 10: Broadcast Blocal within process column 11: else 12: Receive Brecv from process (stage, j) 13: end if 14: Clocal += Arecv Brecv 15: end for 16: return Clocal

B. Distributed SpGEMM

sparse SYRK that takes advantage of the symmetry in the output; 2) CommSpSYRK is a communication-reordered variant that eliminates the explicit materialization of A⊤ and achieves the best performance at scale; 3) A triangular output format that offers memory advantages for symmetric outputs and computational benefits in specific pipelines; 4) Up to 2× faster than Sparse SUMMA and up to 10× faster than CombBLAS on 32 nodes of Perlmutter.

Prior work has developed several distributed-memory SpGEMM algorithms. Sparse SUMMA [10] and Sparse Cannon [18] adapt the standard SUMMA and Cannon algorithms to sparse inputs, partitioning the inputs across a 2D process grid and casting SpGEMM as a sequence of block outer products. Both have scalable communication costs and rely on random permutation for load balancing, a strategy also used here. Sparse SUMMA is the basis for the SpGEMM implementation in CombBLAS [9]. Brock et al. [19] introduce a GPU-accelerated Sparse SUMMA that uses one-sided communication and work stealing for load balancing without random permutation. Bellavita et al. [13] present a distributed SpGEMM for GPU machines with hierarchical interconnects, using a hybrid 1D/2D trident partitioning and a communication schedule that reduces offnode traffic. Both target GPU platforms, which makes direct comparison difficult, and are largely orthogonal to this work. Hussain et al. [2] present a Sparse SUMMA variant for outputs too large for memory, generating a subset of output columns at a time to support applications that never require the full output, such as Markov clustering. Hong et al. [20] use 1D partitioning with sparsity-aware communication via onesided MPI, scaling better than Sparse SUMMA only when the input is highly structured. PETSc and Trilinos implement a similar 1D, point-to-point, sparsity-aware approach [21], [22]. Ranawaka et al. [23] develop a distributed SpGEMM algorithm for C = AB, where B is tall and skinny. None of these works exploit symmetry in the case of C = AAT to reduce communication and computation, which is our main contribution.

II. BACKGROUND AND R ELATED W ORK

III. A LGORITHM

In this section, we define sparse SYRK, review the SUMMA algorithm and prior work on distributed SpGEMM.

Our goal is to compute the sparse SYRK operation C = AAT for a sparse matrix A ∈ Rm×n , distributed across p √ processes on in a q × q grid, where q = p. For simplicity, we assume a square process grid and square blocks; results for rectangular inputs are discussed in Section V-D. The square-grid requirement is inherited from the CombBLAS backend [9]; individual processes hold rectangular submatrices when A is non-square. Process (i, j) owns an m/q × m/q block Cij of the output and an m/q × n/q block Aij of the input. Our starting point is Sparse SUMMA (Alg. 1), which computes C = AB with B = AT materialized explicitly: at stage s, process (i, s) broadcasts Ais along its process row and process (s, j) broadcasts Bsj along P its process column, and each process accumulates Cij = s Ais Bsj . If we apply this to sparse SYRK, it is wasteful because (i) C is symmetric, so computing all of it does twice the necessary work, and (ii) B = A⊤ contains no information not already in A, yet the baseline materializes it via a transpose (a local operation plus a communication round), which doubles the memory overhead of the input. SpSYRK eliminates the first redundancy through symmetry-aware partitioning of computation, and CommSpSYRK additionally eliminates the second through a reordered communication schedule that never

A. SpGEMM and Sparse SYRK SpGEMM computes C = AB where both inputs are sparse, and is a widely used primitive [1], [4], [7], [11]–[14]. This paper focuses on the sparse symmetric rank-k update, or sparse SYRK: given a sparse matrix A ∈ Rm×n , compute C = AAT ∈ Rm×m . The output C is symmetric, so only its lower (or upper) triangular part contains unique information. The sparse SYRK is a special case of SpGEMM in which the second operand is the transpose of the first. In the dense setting, exploiting symmetry in rank-k updates is standard. SYRK is a level-3 BLAS routine [15], [16], and its distributed counterpart PxSYRK in PBLAS, used throughout ScaLAPACK [17], computes only one triangle and halves the arithmetic relative to a general matrix multiply—no analogous routine exists in distributed sparse memory: sparse libraries expose only general SpGEMM, forcing applications to compute the full output and discard half. This work closes that gap, addressing challenges absent in the dense case, such as unpredictable output sparsity and load imbalance caused by the irregular nonzero distribution.

Fig. 1: Computation of each process in Sparse SUMMA (left) versus CommSpSYRK (right) on a 3 × 3 grid. Colors indicate communicators for the broadcast of AT . Diagonal processes compute only the lower triangle of their block; off-diagonal processes compute the complementary halves. ⊤

forms A . Both variants are agnostic to the local kernel: any local sparse SYRK or SpGEMM routine can be used, so the algorithm ports across hardware and library ecosystems. A. Symmetry-Aware Computation The symmetry Cij = C⊤ ji lets each process compute only half of its assigned block. For a diagonal process (i, i), the local accumulation is X X X Cii = Ais (A⊤ )si = Ais A⊤ SYRK(Ais ), is = s

s

s

a sum of local rank-k updates, each with a symmetric result. Diagonal processes therefore compute only the lower (or equivalently, upper) triangle of Cii . For an off-diagonal pair, Cij on process (i, j) is the transpose of Cji on (j, i), so the work is split evenly. The row/column split is chosen to reduce computation: the lowerm rows triangular process (i, j) with i > j computes the first 2q m of Cij , and its partner (j, i) computes the last 2q columns of Cji . Given that row r of Cij is column r of Cji , the two half-blocks are disjoint and together cover the full block pair (Figure 1). In a load balanced matrix, every process computes half of its Sparse SUMMA workload, with the remainder recovered by symmetry for diagonal processes or from the partner for off-diagonal processes. Load balanced computation can be induced via a symmetric random permutation, which we support in our implementation. On completion, each process holds half of its output block. The lower triangle on the diagonal and the complementary half-blocks off it. If an application requires the full matrix, SpSYRK materializes it by transposing. In many applications, however, the triangular output is used directly. Thus, skipping materialization halves the memory of C, which is typically much larger than A, and removes the output transpose from the performance-critical path. B. CommSpSYRK: Reduced-Volume Communication The computation split reduces the inputs required by each P m process. Given that the first 2q rows of Cij = s Ais (A⊤ )sj m depend only on the first 2q rows of each Ais , and the last

Fig. 2: Communication at stage 0 in Sparse SUMMA (left) versus CommSpSYRK (right) on a 3×3 grid. In CommSpSYRK, the column-direction data originates from the transposed grid positions and is trimmed to each receiver’s needs. n columns of Cji depend only on the last 2q columns of each ⊤ (A )sj , lower-triangular processes need only the top half of each received A block, and upper-triangular processes need only the right half of each received A⊤ block. Diagonal processes need no A⊤ at all, since their accumulation involves only products of the form Ais A⊤ is . One iteration of the resulting pattern is shown on the right in Figure 2. The overall result is that less data needs to be communicated at each stage. Notably, implementing this system cannot be done using MPI collectives. MPI Bcast sends identical data to all receivers, MPI Scatterv splits the root’s buffer into disjoint pieces, and MPI Alltoallv requires every member of a communicator to participate. Thus, we use MPI ISend and MPI IRecv to implement CommSpSYRK. CommSpSYRK also can eliminate AT . Observe that in the baseline, the block (AT )sj = ATjs broadcast at stage s is simply the transpose of a block of A that already resides on process (j, s). CommSpSYRK reorders the send pattern: at stage s, instead of process (s, j) sending Bsj down column j, process (j, s) sends Ajs down column j. The set of receivers is unchanged; only the sender is reflected across the diagonal. The local accumulation becomes: X Cij = Ais ATjs , s

pushing the transpose into the local kernel, where it can be fused with the multiply at negligible cost. This assumption is validated empirically in Section IV-B. Thus, CommSpSYRK never materializes A⊤ . This increases collective communication costs, since each receiver gets a different message and the sender is outside a fixed broadcast group, so CommSpSYRK uses q − 1 point-to-point messages per stage instead of an O(log q)-depth broadcast. C. Theoretical Analysis The analysis uses the α–β model [24], extended with local computation rates. A message of w words costs αcomm + βcomm w, where αcomm is the per-message latency and βcomm is the per-word transfer time. A pipelined-tree broadcast [25] is assumed, which, for a broadcast of size n among p processors, sends log(p) messages and O(n) words. Each local kernel

c τin τout σmul

Sparse SUMMA 2 log q 1 0 1

SpSYRK 2 log q 1 0 or 1† 1/2

CommSpSYRK q−1 0 0 or 1† 1/2

σA σB Received volume

TABLE I: Cost coefficients by algorithm. † Output transpose applies only when the application requires full materialization of C; the triangular-output mode sets τout = 0. x ∈ {mul, tr} costs αloc +βx ux , where αloc is a per-invocation overhead shared across local kernels and βx is the marginal cost per unit of work ux : multiply-adds for the local sparse multiply and nonzeros processed for the sparse addition and local transpose. Let η = nnz(A)/p be the per-process nonzero count, approximately uniform after random permutation; let m′ = m/q be the local row count; let γ = nnz(C)/nnz(A) be the output expansion factor; and let: n

φ=

2 1X nnz A(:, k) p k=1

be the per-process multiply-add count of AA⊤ under aP random row/column permutation, which distributes the 2 global multiply-add approximately unik nnz(A(:, k)) formly across p processes. φ is not proportional to η since it depends on the column-nonzero distribution of A; matrices with the same nnz(A) can differ by orders of magnitude in local multiplication cost. For simplicity, A is assumed square (m = n). The three components are then: Ttr (ν) = αloc + βtr ν + αcomm + βcomm (2ν + m′ ), (1) | {z } | {z } local transpose

pairwise exchange

Tcomp (f ) = αloc + βmul f,

(2) ′

Tcomm (ν) = c αcomm + βcomm (2ν + m ),

(3) ′

where a CSR block with ν nonzeros occupies 2ν + m words (index and value types of equal width), and c is the message multiplier of the communication pattern: c = 2 log q for a broadcast and c = q − 1 for a point-to-point round. The local computation is modeled with a single multiplyadd operation. The sparse multiply-accumulate performs one addition per partial product, so the number of additions equals the multiply-add count, φ. Every variant runs exactly q stages; c belongs inside Tcomm , multiplying αcomm only, as in Eq. 3. Sparse SUMMA, SpSYRK, and CommSpSYRK share the same high-level structure: an optional input transpose to form A⊤ , followed by q stages of communicating two blocks and performing a local multiply-and-accumulate, and then an optional output transpose. The algorithms differ in where and how much data moves, captured by: T = τin Ttr (η) + Tcomp (σmul φ)   + q Tcomm (σA η) + Tcomm (σB η) + τout Ttr (σmul γη). (4) The coefficients τin , τout ∈ {0, 1} toggle the input and output transposes, σmul scales the local work, and σA , σB

Upper 1 1/2 3/4

Lower 1/2 1 3/4

Diagonal 1 0 1/2

TABLE II: Communication scaling factors for CommSpSYRK by grid position, relative to Sparse SUMMA. For SpSYRK, σA = σB = 1 everywhere. scale the communicated volume of the row- and columndirection operand. Table I gives c, τin , τout , and σmul for each algorithm, and Table II gives σA and σB by grid position for CommSpSYRK. Because σmul = 1/2, both SpSYRK and CommSpSYRK halve the local work relative to Sparse SUMMA. CommSpSYRK further reduces communication volume: off-diagonal processes receive 3/4 of the baseline volume and diagonal processes receive 1/2. Given that diagonal processes make up only a 1/q fraction of the grid, the asymptotic per-process volume is 3/4 that of Sparse SUMMA, and the input transpose is eliminated. As q grows, the three terms in Eq. 4 scale differently: local work as 1/q 2 , communicated volume as nnz(A)/q, and message count as 2q log q under broadcasts or q(q−1) under pointto-point communication. Three regimes follow. For small q, the βmul φ term dominates and both variants are computebound, with CommSpSYRK also avoiding the input transpose. For intermediate q, CommSpSYRK’s reduced volume wins in the bandwidth-bound regime, while SpSYRK converges toward Sparse SUMMA, matching its communication. For large q, latency dominates and CommSpSYRK’s q −1 messages per stage become a liability compared to a 2 log q broadcast. The compute-to-bandwidth crossover, βcomm η ≈ αcomm q, shifts to higher process counts as nnz(A) increases. Empirical results demonstrate that CommSpSYRK scales the best despite this, thus the crossover point is presumed to be at a very high process count for most inputs. IV. I MPLEMENTATION AND O PTIMIZATION Our implementation relies on CombBLAS for I/O and baselines [9], MPI for inter-process communication, and Intel MKL for local computation [26]. Our implementation spans four configurations across two independent axes. The first axis selects the communication pattern, either the Sparse SUMMA broadcast or the point-to-point CommSpSYRK schedule, isolating the communication savings of the new schedule from the computational savings of the SYRK-aware kernel. The second axis selects whether C is fully materialized or stored halved, isolating the memory savings from the triangular output. A. Computation Algorithm 2 describes the local computation, with r and c denoting the row and column counts of the received blocks. Halving a CSR matrix along rows and halving it along columns differ in cost, and this asymmetry drives most of the off-diagonal process behavior. Halving along rows is free, since the first half of the row-pointer array already describes

Algorithm 2 Local SpSYRK computation on process (i, j)

Label

Require: Left block Arecv ; right block Rrecv when i ̸= j. In SpSYRK, Rrecv is a received block of A⊤ ; In CommSpSYRK, Rrecv = A′⊤ recv , transposed in-kernel. Ensure: Accumulates this process’s half of Cij 1: if i = j then 2: Clocal += SYRK(Arecv ) ▷ lower triangle only 3: else if i > j then 4: Clocal += Arecv [ : r/2, :] × Rrecv ▷ Row top half 5: else 6: Clocal += Arecv × Rrecv [ :, c/2 : ] ▷ Col right half 7: end if

vas stokes 2M stokes kron g500-logn21 rgg n 2 24 s0 relat9 rel9 spal 004

Algorithm 3 CommSpSYRK send at stage s, process (i, j) Require: A distributed on a q × q grid Ensure: Each process receives the input blocks required for its local computation Row send (supplies Ais along process row i): 1: if j = s then 2: for each receiver (i, j ′ ) in row i do 3: if i ≥ j ′ then ▷ diagonal or upper receiver: full block 4: IS END Alocal to (i, j ′ ) 5: else ▷ lower receiver: top half of rows 6: IS END Alocal [ : r/2, :] to (i, j ′ ) 7: end if 8: end for 9: end if Column send (supplies Ajs , standing in for Bsj = A⊤ js ): 10: if i = s then 11: for each off-diagonal receiver (i′ , j) in column j, i′ ̸= j do 12: if i′ < j then ▷ upper receiver: bottom half of rows 13: IS END Alocal [ r/2 :, :] to (i′ , j) 14: else ▷ lower receiver: full block 15: IS END Alocal to (i′ , j) 16: end if 17: end for 18: end if

a valid submatrix. Halving along columns is not free, because the first half of the columns is contiguous in neither the column nor the values array, which forces a full traversal and rebuild. A process computing the top rows of Cij obtains its input immediately, whereas one computing the right columns must either prune columns directly from the CSR structure or transpose so that columns become rows and pruning is cheap again. In our experiments, we found that local MKL transposes are at least an order of magnitude faster than MKL SpGEMM kernels, so the implementation transposes B, builds a handle from the latter half of its row pointers, and then transposes back to obtain the last columns of B. Because MKL cannot handle

Rows (K) Cols (K) NNZ(A) (K) NNZ(C) (K) Density (%) 2,146 11,449 2,097 16,777 12,360 9,888 10

2,146 11,449 2,097 16,777 549 274 321

65,129 349,321 182,082 265,114 38,955 23,667 46,168

1,146,604 5,937,244 5,772,532 845,416 5,516,977 2,038,738 51,899

γ Size (GB)

0.0014 18 0.00027 17 0.0041 32 0.000094 3.2 0.0057 142 0.00087 86 1.44 1.1

1.5 8.4 1.5 2.1 0.6 0.4 1.1

TABLE III: Test matrices. NNZ(C) counts nonzeros in the upper triangle of C = AAT . Density is NNZ(A)/mn. The compression factor γ = FLOP/2 NNZ(C) measures accumulation per output nonzero; larger γ means more compute bound.

the off-by-one inner dimensions that occur when blocks divide unevenly, the inner dimensions of both operands and the output are padded and then trimmed at the end. Diagonal processes (i = j) call mkl_sparse_syrk directly, while off-diagonal processes call the general routine mkl_sparse_sp2m. B. Communication and Transposes In Sparse SUMMA, the process broadcasting a block down a column is itself a member of that column, so MPI Bcast suffices. CommSpSYRK instead sends the underlying block of A: at stage s, process (j, s) sends Ajs down column j, playing the role that Bsj = A⊤ js plays in Sparse SUMMA, with the transpose deferred to the local kernel (Algorithm 3, column-send). Because the sender is not a member of the receiving column, a collective cannot be used, so individual sends are issued instead. Given that a sending process may simultaneously be a receiver, MPI ISend prevents deadlock. MKL provides no transpose routine, so two alternatives were implemented and evaluated: mkl_sparse_add with a zero matrix and mkl_convert_csr. We found an order of magnitude gap from SpGEMM kernels to mkl_convert_csr, and this justifies the mkl_convert_csr approach for both the compute-only routine and the full materialization of C. If a transpose across the process grid is unavoidable, as in SpSYRK’s input transpose and output materialization, a local transpose is followed by a one-to-one exchange. V. R ESULTS AND D ISCUSSION This section reports strong scaling and runtime breakdown for each algorithm against appropriate baselines. Throughout, “processes” refers to MPI processes. A. Experimental Setup Experiments were run on the CPU nodes of NERSC Perlmutter. Each node has 128 CPU cores and 512 GB of memory, connected through a Dragonfly topology. Runs used 32 nodes with 8 threads or 4 physical cores per process and 32 processes per node, varying the process count to study strong scaling. We chose this distribution as we found it to be a sweet spot on per process compute times. Two matrix properties govern the cost of sparse SYRK. The first is output expansion, nnz(C)/nnz(A), which determines how much runtime is spent in local kernels and how much

memory the result uses. The second is the density of the nonzeros, which determines how concentrated the communication is and therefore how much a symmetry-aware schedule can save. vas_stokes_2M and stokes are matrices from semiconductor process problems with localized patterns and modest expansion; kron_g500-logn21 is a Kronecker power-law graph whose hub vertices produce poor locality and large expansion; and rgg_n_2_24_s0 is a random geometric graph that is spatially localized but has low expansion. These four cover the two axes independently rather than varying them together. The remaining three matrices, relat9, rel9, and spal_004, are rectangular and are used in Section V-D to test the algorithm beyond the square cases, but also span expansion factors and locality, with spal_004 has low expansion but good locality, and relat9 and rel9 have similar expansion but different locality. For the baselines, we used CombBLAS’ Sparse SUMMA algorithm, except on kron due to out of memory errors, and an optimized Sparse SUMMA implementation built on top of CombBLAS but using MKL kernels. CombBLAS divides A into two halves and B into p row-wise partitions. In all implementations, matrix values are stored in double precision and the indices in int64. B. Load Balance Under load imbalance, CommSpSYRK continued to improve as the number of processes increased, whereas the reference Sparse SUMMA and SpSYRK stalled. The synchronous broadcasts in Sparse SUMMA force every process to wait for the most heavily loaded one, so even small nonzero imbalances limit scaling, while the point-to-point rounds of CommSpSYRK tolerate imbalance better. CommSpSYRK, in its smaller communication rounds, is “sparsity-aware,” and therefore improves more under imbalanced conditions; applying a random permutation to the rows and columns of A restores balance across all approaches. Therefore, we report load balance wherever possible, except in the case of a rectangular input matrix, where symmetric load balance is not possible. A nonsymmetric permutation is possible but is not offered by CombBLAS, so we do not report those results here. Our results show that load balancing generally improves on standard algorithms while not significantly improving CommSpSYRK, our best algorithm. Therefore we note that CommSpSYRK avoids random permutation as a pre-processing step. C. Comparison to Baselines Figure 3 shows strong-scaling results for all test matrices, with markers indicating mean runtime and error bars representing standard deviation. The gap between Sparse SUMMA and SpSYRK varied by matrix: rgg showed parity, spal_004 some regression, and the rest saw roughly a 2× speedup. Under CommSpSYRK, spal_004 and rgg improved significantly, while the rest showed parity or, in the case of kron, a slight regression. This variation follows arithmetic intensity. Halving local work pays off only when local multiplication dominates the runtime, a condition the cost

Fig. 3: Runtime scaling on the four test matrices. Error bars show variation across different runs.

model captures through φ, the multiply-add count of AA⊤ (Section III-C). The matrices with low φ relative to nnz(A) spend proportionally more time on communication and CSR row and column metadata, neither of which SpSYRK changes, so the halved computation is hidden. For low process counts, where local computation dominates, SpSYRK and CommSpSYRK are comparable and often outperform Sparse SUMMA. As the process count increases and communication begins to dominate, SpSYRK converges toward Sparse SUMMA, matching its communication volume, while CommSpSYRK pulls ahead and approaches the predicted 3/4 ratio. Our measurements therefore lie in the compute-dominated regime at lower process counts and approach the bandwidth-dominated regime in higher process counts. A further breakdown is shown in Section V-F. For some matrices, notably rgg and spal_004, the lowprocess SpSYRK gains fall short of the prediction. The cause is an inefficiency in the MKL kernel. MKL’s local sparse SYRK kernel is less efficient than its SpGEMM kernel on very sparse matrices, so halving the input does not halve the runtime. Our local benchmarks on random matrices of varying densities confirm this, showing that the MKL sparse SYRK kernel is slower than performing the same operation via mkl_sparse_sp2m for matrices with densities less than or equal to ≈ 0.005%. The slowdown scales with sparsity, with matrices of density lower than 0.0001% consistently showing over 2× slowdown. D. Rectangular Matrices In our analysis in Section III-C, the row count enters explicitly into Tcomm and Ttr through the m′ term, whereas the multiply-add count φ is determined by the column-nonzero distribution rather than by the row count. Because the CommSpSYRK bisects matrices row-wise, the row count is reduced at the same rate as the number of nonzeros when they are distributed uniformly across rows; in this regime, the per-process scaling factor for communication and transpose costs is the same for both quantities. Therefore, a rectangular matrix in our model would experience similar performance scaling between Sparse SUMMA and CommSpSYRK. Figure 3 confirms this on relat9, rel9, and spal_004 At high process counts CombBLAS outperforms on relat9 and rel9, which we suspect arises from the CSC format of CombBLAS giving a scaling advantage. E. Communication and Computation Breakdown Figure 5 shows the runtime breakdown: computation in MKL kernels, communication in transfer rounds and staging, and transpose for forming A⊤ and materializing the output. The breakdown is for a single matrix, but the relative trends hold across the entire suite. Sparse SUMMA and SpSYRK show similar trends with increasing process counts and similar scaling, which is expected since Section III-C implies that both achieve constant fractions of the Sparse SUMMA runtime in the high-process limit.

Figure 4 breaks communication time down by process type: diagonal, upper triangular, and lower triangular, which makes CommSpSYRK’s variance clearer. Diagonal processes have lower communication times, receiving only about 23 of the data that off-diagonal processes do. Grouping by type removes much of the variance; the remainder reflects CommSpSYRK’s asynchronous communication, which approaches a synchronous broadcast round as the process count increases. The stratification between the upper and lower off-diagonals reflects synchronization effects, since the communication system is identical to Sparse SUMMA and the split comes from processes idling while waiting for others to enter the round. F. Memory Footprint and the Materialization Tradeoff Figure 5 separates materialization from the rest of the runtime. The absolute cost of materialization falls with process count, as Section III-C predicts, since the output transpose scales with the per-process output ζ. As a fraction of total runtime, it stays roughly stable because the remaining terms shrink at a similar rate. For the single-process case, the breakdown highlights the relative inefficiency of the MKL sparse SYRK kernel noted in Section V-C, since computation there reduces to local multiplication. For Sparse SUMMA and SpSYRK, communication grows until it dominates, and the total runtime increases beyond about 256 processes, whereas CommSpSYRK’s communication stays small and its runtime keeps decreasing over the tested range. These observations argue against materializing the full output for two reasons. The first is memory. The output is far denser than the input: across our matrices, nnz(C) exceeded nnz(A) by roughly an order of magnitude (Table III), so C, not A, sets the peak memory footprint and thus determines the largest problem that fits in a given allocation; storing only the triangular half cuts this in half. The second issue concerns downstream use: consumers are often assumed to need the full matrix, but the motivating applications frequently do not. For instance, similarity search enumerates each candidate pair once [4], [5], [7], using only half of the overlap matrix, and the Cholesky factorization of the normal equations reads only a single triangular part. In addition, we show that it is possible to perform downstream sparse primitives using only the non-materialized C. Consider the product Cx, and let CL denote the stored lowertriangular output, including the diagonal, with D = diag(CL ). Given that C = CL + C⊤L − D,
⊤ Cx = CL x + x⊤ CL − Dx, (5) where the second term uses the transpose of the vector rather than of the matrix, so C⊤L is never formed, and Dx removes the double-counted diagonal through a local scaling on diagonal processes. The arithmetic is unchanged; the only extra cost is that x moves twice, once along process columns and once along process rows, instead of once. The same decomposition applies to a matrix operand X,
⊤ CX = CL X + X⊤ CL − DX, (6)

Fig. 4: Communication time, broken down by process location. The communication cost of each algorithm is shown, divided by whether a process is on the diagonal of the process grid, the upper triangular, or the lower triangular.

Fig. 5: Runtime split by operation for each algorithm. The runtime is divided into computation time, communication time, time spent computing transposes, and time spent materializing the full output. Each bar’s components are determined by the runtime of that component among the processes with the largest total runtime in the grid on RGG.

and here the second term reuses CommSpSYRK’s transposedsend pattern, so the blocks of CL it needs are routed by the schedule without additional machinery. The aggregate data volume of CL moved in Eq. 6 is 2 nnz(CL ) ≈ nnz(C), matching a SUMMA multiply on the materialized matrix, with X still moved once. Relative to materializing C and then multiplying, the savings are the output transpose Ttr and half the resident footprint of C, with equal product-phase communication and arithmetic cost. For Cx, the product can instead be computed as A(A⊤ x), avoiding C entirely, so there is no advantage there. For CX the situation is different, since A(A⊤ X) can require more arithmetic than (AA⊤ )X, making the triangular form worthwhile. Computing the full symmetric completion is thus unnecessary when the output feeds a subsequent kernel, which both lets us evaluate our algorithms without materialization and provides a general method for consuming symmetric output. VI. C ONCLUSION This paper presented SpSYRK and CommSpSYRK, the first distributed-memory algorithms for sparse SYRK that exploit the symmetry of the output matrix. SpSYRK divides each offdiagonal block between upper- and lower-triangular process pairs and computes only the lower triangle of each diagonal block, halving the computation compared to Sparse SUMMA. CommSpSYRK also reorganizes communication so that A⊤ is

not formed, reducing communication volume and eliminating the need to transpose the input. The triangular output can be materialized when a downstream component requires the complete matrix. On the Perlmutter supercomputer, the algorithms achieved a 2× speedup over an optimized Sparse SUMMA on 32 nodes for compute-bound input, with the advantage determined by arithmetic intensity. Because the partitioning and communication schemes are agnostic to the local kernel, future work will study GPU acceleration and asynchronous execution that overlaps communication with computation. The observed communication disparity suggests this could further improve CommSpSYRK. 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-AC0205CH11231, using NERSC award ASCR-ERCAP0030076. This material is based upon work supported by the U.S. Department of Energy, Office of Science, Office of Advanced Scientific Computing Research, Department of Energy Computational Science Graduate Fellowship under Award No. DESC0025528.

R EFERENCES [1] M. Besta, R. Kanakagiri, H. Mustafa, M. Karasikov, G. Rätsch, T. Hoefler, and E. Solomonik, “Communication-efficient jaccard similarity for high-performance distributed genome comparisons,” in 2020 IEEE International Parallel and Distributed Processing Symposium (IPDPS). IEEE, 2020, pp. 1122–1132. [2] M. T. Hussain, O. Selvitopi, A. Buluc, and A. Azad, “Communication-Avoiding and Memory-Constrained Sparse MatrixMatrix Multiplication at Extreme Scale,” in 2021 IEEE International Parallel and Distributed Processing Symposium (IPDPS). Portland, OR, USA: IEEE, May 2021, pp. 90–100. [Online]. Available: https://ieeexplore.ieee.org/document/9460450/ [3] A. Azad, A. Buluç, and J. Gilbert, “Parallel triangle counting and enumeration using matrix algebra,” in 2015 IEEE International Parallel and Distributed Processing Symposium Workshop. IEEE, 2015, pp. 804–811. [4] G. Guidi, M. Ellis, D. Rokhsar, K. Yelick, and A. Buluç, “Bella: Berkeley efficient long-read to long-read aligner and overlapper,” in SIAM Conference on Applied and Computational Discrete Algorithms (ACDA21). SIAM, 2021, pp. 123–134. [5] G. Guidi, O. Selvitopi, M. Ellis, L. Oliker, K. Yelick, and A. Buluç, “Parallel string graph construction and transitive reduction for de novo genome assembly,” in 2021 IEEE International Parallel and Distributed Processing Symposium (IPDPS). IEEE, 2021, pp. 517–526. [6] A. Buluç, “The ubiquitous sparse matrix-matrix products,” 2025. [Online]. Available: https://arxiv.org/abs/2508.04077 [7] O. Selvitopi, S. Ekanayake, G. Guidi, G. A. Pavlopoulos, A. Azad, and A. Buluç, “Distributed many-to-many protein sequence alignment using sparse matrices,” in SC20: International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE, 2020, pp. 1–14. [8] O. Selvitopi, S. Ekanayake, G. Guidi, M. G. Awan, G. A. Pavlopoulos, A. Azad, N. Kyrpides, L. Oliker, K. Yelick, and A. Buluç, “Extremescale many-against-many protein similarity search,” in SC22: International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE, 2022, pp. 1–12. [9] A. Buluç and J. R. Gilbert, “The Combinatorial BLAS: design, implementation, and applications,” The International Journal of High Performance Computing Applications, vol. 25, no. 4, pp. 496–509, Nov. 2011. [Online]. Available: https://journals.sagepub.com/doi/10.1177/1094342011403516 [10] ——, “Parallel sparse matrix-matrix multiplication and indexing: Implementation and experiments,” SIAM Journal on Scientific Computing, vol. 34, no. 4, pp. C170–C191, 2012. [11] J. R. Gilbert, S. Reinhardt, and V. B. Shah, “High-performance graph algorithms from parallel sparse matrices,” in Applied Parallel Computing. State of the Art in Scientific Computing, B. Kågström, E. Elmroth, J. Dongarra, and J. Waśniewski, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg, 2007, pp. 260–269. [12] T. A. Davis, “Graph algorithms via suitesparse: Graphblas: triangle counting and k-truss,” in 2018 IEEE High Performance extreme Computing Conference (HPEC), 2018, pp. 1–6. [13] J. Bellavita, L. Pichetti, T. Pasquali, F. Vella, and G. Guidi, “Communication-avoiding spgemm via trident partitioning on hierarchical gpu interconnects,” arXiv preprint arXiv:2603.21444, 2026. [14] Y. Li and G. Guidi, “Ocean: Fast estimation-based sparse general matrixmatrix multiplication on gpu,” arXiv preprint arXiv:2604.19004, 2026. [15] L. S. Blackford, A. Petitet, R. Pozo, K. Remington, R. C. Whaley, J. Demmel, J. Dongarra, I. Duff, S. Hammarling, G. Henry et al., “An updated set of basic linear algebra subprograms (blas),” ACM Transactions on Mathematical Software, vol. 28, no. 2, pp. 135–151, 2002. [16] J. J. Dongarra, J. Du Croz, S. Hammarling, and I. S. Duff, “A set of level 3 basic linear algebra subprograms,” ACM Transactions on Mathematical Software (TOMS), vol. 16, no. 1, pp. 1–17, 1990. [17] L. S. Blackford, J. Choi, A. Cleary, E. D’Azevedo, J. Demmel, I. Dhillon, J. Dongarra, S. Hammarling, G. Henry, A. Petitet et al., ScaLAPACK users’ guide. SIAM, 1997. [18] A. Buluc and J. R. Gilbert, “Challenges and advances in parallel sparse matrix-matrix multiplication,” in 2008 37th International Conference on Parallel Processing, 2008, pp. 503–510.

[19] B. Brock, A. Buluç, and K. Yelick, “Rdma-based algorithms for sparse matrix multiplication on gpus,” in Proceedings of the 38th ACM International Conference on Supercomputing, 2024, pp. 225–235. [20] Y. Hong and A. Buluç, “A sparsity-aware distributed-memory algorithm for sparse-sparse matrix multiplication,” in SC24: International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE, 2024, pp. 1–14. [21] S. Balay, S. Abhyankar, M. F. Adams, S. Benson, J. Brown, P. Brune, K. Buschelman, E. M. Constantinescu, L. Dalcin, A. Dener, V. Eijkhout, J. Faibussowitsch, W. D. Gropp, V. Hapla, T. Isaac, P. Jolivet, D. Karpeev, D. Kaushik, M. G. Knepley, F. Kong, S. Kruger, D. A. May, L. C. McInnes, R. T. Mills, L. Mitchell, T. Munson, J. E. Roman, K. Rupp, P. Sanan, J. Sarich, B. F. Smith, S. Zampini, H. Zhang, H. Zhang, and J. Zhang, “PETSc Web page,” https://petsc.org/, 2025. [Online]. Available: https://petsc.org/ [22] M. A. Heroux, R. A. Bartlett, V. E. Howle, R. J. Hoekstra, J. J. Hu, T. G. Kolda, R. B. Lehoucq, K. R. Long, R. P. Pawlowski, E. T. Phipps, A. G. Salinger, H. K. Thornquist, R. S. Tuminaro, J. M. Willenbring, A. Williams, and K. S. Stanley, “An overview of the trilinos project,” ACM Trans. Math. Softw., vol. 31, no. 3, p. 397–423, Sep. 2005. [Online]. Available: https://doi.org/10.1145/1089014.1089021 [23] I. Ranawaka, M. T. Hussain, C. Block, G. Gerogiannis, J. Torrellas, and A. Azad, “Distributed-memory parallel algorithms for sparse matrix and sparse tall-and-skinny matrix multiplication,” in SC24: International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE, 2024, pp. 1–17. [24] R. W. Hockney, “The communication challenge for mpp: Intel paragon and meiko cs-2,” Parallel Comput., vol. 20, no. 3, p. 389–398, Mar. 1994. [Online]. Available: https://doi.org/10.1016/S0167-8191(06)80021-9 [25] J. Bruck, R. Cypher, and C.-T. Ho, “Multiple message broadcasting with generalized fibonacci trees,” 1992. [26] E. Wang, Q. Zhang, B. Shen, G. Zhang, X. Lu, Q. Wu, and Y. Wang, “Intel math kernel library,” in High-Performance Computing on the Intel® Xeon Phi™: How to Fully Exploit MIC Architectures. Springer, 2014, pp. 167–188.

Related documents

• Distributed Algorithms#3468
• Butterfly Factorization Via Randomized Matrix-Vector Multiplications#1119895
• DB-SpMSpV: Dual-View Blocked Sparse Matrix-Sparse Vector Multiplication for Dynamic GPU Workloads#620708
• DB-SpMSpV: Dual-View Blocked Sparse Matrix-Sparse Vector Multiplication for Dynamic GPU Workloads#619029
• DB-SpMSpV: Dual-View Blocked Sparse Matrix-Sparse Vector Multiplication for Dynamic GPU Workloads#459612