A GPU Engine for Sparse Matrix-Vector Products on $k^2$-Tree-Compressed Graphs

A GPU Engine for Sparse Matrix-Vector Products on 𝑘 2-Tree-Compressed Graphs Francesco Tosoni

Gabriele Mencagli

Francesco.Tosoni@santannapisa.it Sant’Anna School of Advanced Studies, L’EMbeDS Pisa, Italy Department of Computer Science, University of Pisa Pisa, Italy

gabriele.mencagli@unipi.it Department of Computer Science, University of Pisa Pisa, Italy

arXiv:2611.03610v1 [cs.DS] 2 Nov 2026

Abstract

algebra (CLA) [26, 27], computation-friendly matrix formats [4– 7, 10, 11, 25, 29, 31, 41, 42, 56–58] run such products directly on compressed representations, though not always in time proportional to their size [1]; thus, Boldi–Vigna-style encodings of adjacency lists, by gaps between successors and references to similar lists [16, 30, 59], grammar-compressed matrices [29, 56] and wavelet trees [47, §6.2] have reached the GPU [32, 53, 58]; Navarro surveys this line of work [48]. Standard sparse formats pay for every non-zero: the coordinate format (COO) stores both coordinates, the compressed sparse row format (CSR) a column index plus one pointer per row, and library implementations of both also store a value, even when every value is one. A CSR serves products by rows, so the left product 𝑥 ⊤𝐴 needs a transposed kernel, which scatters its outputs, or a second, columnoriented copy (CSC). The 𝑘 2 -tree [18][47, §9.2.1] takes another route for Boolean matrices: each node splits its submatrix into 𝑘 2 blocks and keeps one bit per block, set iff the block holds a non-zero. Empty regions cost nothing and the position of a non-zero is implicit in its path from the root, so neither values nor coordinates are stored, and rows and columns are served alike. Introduced for web graphs more than fifteen years ago [18, 19], the 𝑘 2 -tree has since represented ternary relations, set operations and spatial data [17, 23, 33, 39], and supported sparse Boolean algebra: regular path queries (RPQs) in graph databases [6, 7], sums and products in cache-friendlier layouts [20, 28], and PageRank, where it balances time, space and energy between grammar compression and Zuckerli [59], a Boldi– Vigna-style compressor with entropy coding [57]. To the best of our knowledge, it has never been implemented on a GPU. The obstacle is the layout that makes the tree compact. Nodes are stored without pointers, level by level, and the children of a node are found by rank, the count of set bits up to a position in a bitvector [47, §4.2]. Rank is cheap in isolation, also on GPUs [32], but a top-down descent issues one lookup per visited node, each at an address known only after the previous one. This chain of dependent, irregular loads can be neither coalesced nor overlapped, and a GPU, which hides latency only behind many independent accesses, spends the descent waiting on memory.

The 𝑘 2 -tree is a compact representation of Boolean matrices: it stores the recursive hierarchy of their non-empty submatrices, with neither values nor coordinates, and serves rows and columns alike. Its traversal, however, locates the children of each node by rank queries, a chain of dependent irregular loads that leaves a GPU waiting on memory. We implement sparse matrix–vector multiplication (SpMV) on 𝑘 2 -trees for the GPU. Since the traversal does not depend on the vector, the positions of the children are static: we descend the tree one level at a time and obtain them from a parallel prefix sum or from a sampled directory of about 0.1 bits per node. Descending in tiles bounds the temporary memory by the depth of the tree, so the device memory of a product is known before it runs. The same traversal computes right and left products, for one vector or a batch, over the plus-times, Boolean and tropical semirings. We evaluate the engine on an NVIDIA GB10 and an H100 with fourteen Wikidata, Software Heritage, web and social graphs of up to 1.2 billion edges, against cuSPARSE, a grammar-compressed GPU engine and a Boolean CSR. Our results show that the tree is 2.3–27× smaller than the CSR of cuSPARSE and that, with its last internal level cached, a product takes 0.7–2.4× the time of cuSPARSE on the right and 0.27–2.6× on the left. On a 1.2-billion-edge graph, a product needs 1.1 GB on the device, 9× less than cuSPARSE. The engine is memory bounded: the 𝑘 2 -tree is the format of choice when device memory, not time, is the constraint.

CCS Concepts • Computing methodologies → Parallel algorithms; • Theory of computation → Data compression; • Mathematics of computing → Graph algorithms.

Keywords 𝑘 2 -tree, compressed sparse matrices, GPU, semiring SpMV, knowledge graphs, reachability

1

Introduction

Our approach. We traverse the tree level by level: the nodes of a level form a frontier that a kernel expands in parallel into the next. Our key observation is that this traversal is the same for every input vector, so what rank provides, where the children of each node lie, is static. On it we build variants that trade memory for time. The most direct, tree-only, keeps only the tree and recomputes the child offsets at every product by an exclusive prefix sum (scan) over the frontier, the data-parallel primitive set out by Blelloch [14]. Alternatively,

Iterative graph analytics are, in their inner loop, sparse matrix– vector products (SpMV): a step of breadth-first search is a Boolean product, a step of Bellman–Ford the same product over the tropical (min,+) semiring, and an iteration of PageRank the plus-times product [22][34, §7.3]. On GPUs this product is bound by memory: bandwidth sets its speed, capacity which graphs fit at all. Numerical linear algebra has thus long worked to move fewer bits per value or index [3, 36, 37, 52]. Following compressed linear 1

Tosoni et al.

a directory of offsets sampled every 256 nodes adds about 0.1 bits per node: the scan is faster for batches, the directory for single vectors. The frontier is the only memory a product needs beyond the matrix and the vectors, and descending in tiles, depth-first over contiguous ranges of each level, bounds it by the depth of the tree rather than its size. When memory is plentiful, the Morton codes (coordinates) of the last internal level can be stored with the matrix, at eight bytes per node: with this cached frontier, a single-vector product is one kernel over the leaves. The same traversal computes the left product by swapping two coordinates at the leaves, and serves batches and any semiring whose addition the GPU applies atomically. None of these choices depends on the vector, so the device memory of a product (tree, directory, cached codes and traversal buffers, sized by the host when it plans the tiles) is fixed by the tree and the tile size, besides the dense vectors: it is known after construction, before any kernel runs, and is the same on every GPU (Tables 2 and 3). Time and energy depend on the GPU, and we measure them (Sections 5 and 5.7). Being lossless, the tree also composes with lossy stages that produce a Boolean pattern, such as the pruning of a learned model, and leaves the accuracy of downstream tasks unchanged.

graphs on the right and on all of them in batch, so the 𝑘 2 -tree is the format for when memory, not time, is the constraint.

2 Background 2.1 The 𝑘 2 -tree Let 𝐴 ∈ {0, 1}𝑛×𝑛 with 𝑛 = 2ℓ (a rectangular matrix is padded with empty rows and columns, which cost nothing). The 𝑘 2 -tree with 𝑘 = 2 [18][47, §9.2.1] is the tree of the recursive 2 × 2 decomposition of 𝐴: the root has a child for each non-empty quadrant, and the recursion stops at single cells. Each internal node stores a signature 𝑠 of 𝑘 2 = 4 bits: with the quadrants numbered 𝑏 = 2𝑏 1 + 𝑏 0 , where 𝑏 1 selects the lower half of the rows and 𝑏 0 the right half of the columns, bit 𝑏 of 𝑠 is set iff quadrant 𝑏 is non-empty. In Figure 1, the root has 𝑠 = 1111 (bits from position 0), and its child 01, the topright quadrant, has 𝑠 = 1010, since only its blocks 00 and 10 hold non-zeros. Concatenating the signatures in level order yields the pointerless bitstring𝑇 . The children of the node whose bit is at offset 𝑝 of 𝑇 start at offset 𝑘 2 · rank1 (𝑇 , 𝑝), where rank1 (𝑇 , 𝑝) counts the set bits of 𝑇 [0..𝑝]: every step down the tree is a rank query [7][47, §9.2.1]. A matrix with 𝑚 non-zeros induces O (𝑚 log(𝑛 2 /𝑚)) nodes; our trees take from 2.4 bits per edge on a web crawl to 28 on a citation relation (Table 2).

Contributions. • A level-synchronous 𝑘 2 -tree SpMV for the GPU (Section 3), whose child offsets, obtained by rank in a descent, are static data, scanned at each product or read from a sampled directory, and whose memory beyond the matrix and the vectors is bounded by the depth of the tree. • Right and left products from one structure (Section 3.6): unlike CSR and Boldi–Vigna-style compressors, which serve one direction, the tree computes 𝐴𝑥 and 𝑥 ⊤𝐴 from the same data, with the same memory up to the dense 𝑥 and 𝑦, for one vector or a batch, over the plus-times, Boolean and tropical semirings. • An evaluation on fourteen graphs of up to 1.2×109 edges (Section 5) on an NVIDIA GB10 and an H100: Wikidata relations, a subgraph of the Software Heritage (SWH) archive and six web and social graphs, against the cuSPARSE CSR [49], which stores single-precision values, a recently proposed grammar-compressed GPU engine [58], bCSR, a Boolean CSR without values that we implement as a simple baseline, and a multicore CPU descent of the same tree; it includes a step-by-step ablation, a placement against measured hardware ceilings, the construction cost and an end-to-end multisource breadth-first search. The tree is 2.3–27× smaller than the cuSPARSE CSR, from the citation relation to the web crawls, as its compression follows the clustering of the non-zeros, and is built on one host core in 0.35– 312 s, about as fast as a RePair grammar [40]. With the frontier cached, a product takes 0.7–2.4× the time of cuSPARSE on the right, 0.27–2.6× on the left and 1.8–3.2× in batch, where it beats the grammar engine on every graph. Descended in tiles, the tree-only engine holds the fewest device bytes of any engine on eight of the nine largest graphs, 1.3–6.0× fewer than bCSR; on the SWH subgraph it needs 1.1 GB, against 10.3 GB for cuSPARSE. The tree trades time for these bytes: bCSR is faster on thirteen of fourteen

2.2

Semiring products

A semiring S = (𝐷, ⊕, ⊗, 0¯, 1¯ ) is a set 𝐷 with two operations such that (𝐷, ⊕, 0¯ ) is a commutative monoid with identity 0¯ , (𝐷, ⊗, 1¯ ) ¯ a monoid with identity 1¯ , ⊗ distributes over ⊕, Éand 0 annihilates under ⊗. It defines the product 𝑦 = 𝐴𝑥 by 𝑦𝑖 = 𝑗 𝐴𝑖 𝑗 ⊗ 𝑥 𝑗 . Our 𝐴 is the Boolean adjacency matrix of a graph without edge weights: ¯ an absent É edge is 0 and every edge carries the same weight 𝑤, so 𝑦𝑖 = 𝑤 ⊗ 𝑥 𝑗 . Contributions reach 𝑦𝑖 from many threads 𝑗: 𝐴𝑖 𝑗 =1 in no fixed order, which the associativity and commutativity of ⊕ make harmless (in floating point, up to rounding: Section 5.9); beyond the axioms, our kernel only requires an atomic ⊕ on the GPU. Three semirings are of direct use (Table 1). Plus-times, with 𝑤 = 1¯ = 1, is the kernel of PageRank and of any power iteration on an unweighted graph [34, §7.3]; Boolean (∨, ∧), also with 𝑤 = 1¯ , is a step of breadth-first search [21, §22.2] or reachability; tropical (min, +), whose 1¯ is 0, with unit edge length 𝑤 = 1, is a relaxation round of Bellman–Ford [21, §24.1] on an unweighted graph, a hopcount shortest-path step. Shortest paths as (min, +) matrix products are textbook material [21, §25.1], generalized to any semiring by GraphBLAS [22]. Only the matrix is Boolean: 𝑥 and 𝑦 hold arbitrary values of 𝐷 (single-precision reals here), so the tree serves any analytics on graphs without edge weights; weighted shortest paths would need edge values, which the 𝑘 2 -tree does not keep.

3

Level-Synchronous 𝑘 2 -tree SpMV

A GPU rewards independent work and coalesced accesses and punishes chains of dependent, irregular loads such as a top-down descent. We thus propose to descend the tree one level at a time (Algorithm 1): the nodes of a level form a frontier, held in level order, and one kernel expands all of them into the next level, with independent lookups issued by adjacent threads at adjacent positions. 2

A GPU Engine for Sparse Matrix-Vector Products on 𝑘 2 -Tree-Compressed Graphs

1

1

1

1

1

1

1

1

1

empty path f = 4, off = 0

0

1

2

3

code T , level 1 f off

00 1000 1 0

01 1010 2 1

10 1111 4 3

11 0001 1 7

(code ≪ 2) | b1 b0

0

1

2

3

4

5

6

7

code T , leaves (r, c)

0000 0001 (0,0)

0100 1001 (0,4)

0110 0001 (2,4)

1000 0001 (4,0)

1001 0010 (4,2)

1010 0100 (6,0)

1011 1000 (6,2)

1111 1111 (6,6)

leaf number off + j (code ≪ 2) | b1 b0

1

1

ε 1111

code T , level 0

1

popcount(T ) exclusive scan of f

odd bits → r, even bits → c path of cell (3, 5): code 01 10 11

Figure 1: A Boolean matrix (left) and its 𝑘 2 -tree (right) as the level-order bitstring 𝑇 , one level per block, each node above the range of its children. Each level shows what the traversal of Section 3 computes: the Morton code of a node (its parent’s code followed by the two bits 𝑏 1𝑏 0 of its quadrant, in color), its fan-out 𝑓 , and the offset of its first child, which a descent obtains by rank (boxed: the offsets a directory with 𝑆 = 2 stores). At the leaves, the code is deinterleaved into the origin (𝑟, 𝑐) of the 2 × 2 block. Highlighted: the path of cell (3, 5), through quadrant 01 of the root, block 10 of node 01 and cell 11 of leaf 2; its code 01 10 11 has row bits 011 = 3 and column bits 101 = 5. Table 1: The three semirings we implement: domain, operations, identities, and the contribution 𝑤 ⊗ 𝑥 𝑗 of a non-zero, applied by one GPU atomic (∞ is 109 ; non-negative floats compare as their bit patterns read as integers, so integer atomics serve max and min). All three are cross-checked (Section 5.9); timings use plus-times, the end-to-end search (Section 5.8) the Boolean semiring. semiring

𝐷

⊕

⊗ 0¯ 1¯

leaf

Algorithm 1 Level-synchronous 𝑘 2 -tree SpMV 𝑦 = 𝐴𝑥 over a semiring S. Node 𝑡 of the frontier has Morton code 𝐹𝑡 , signature 𝑠𝑡 , fan-out 𝑓𝑡 = popcount(𝑠𝑡 ), and first child off𝑡 in the next level; each parallel for gives one thread to each node. Require: tree 𝑇 , one signature (nibble) per node in level order, with levels 0 (root) to ℓ − 1 (leaves); vector 𝑥; semiring S with edge weight 𝑤; frontier 𝐹 , the Morton codes of the nodes of a cached level 𝑑 0 ¯𝑛 1: 𝑦 ← 0 ⊲ fill kernel 2: for 𝑑 = 𝑑 0 to ℓ − 2 do ⊲ one level per step Í 3: off𝑡 ← 𝑢<𝑡 𝑓𝑢 for every node 𝑡 of 𝐹 ⊲ ChildOffset 4: if 𝑑 < ℓ − 2 then ⊲ expand to the next frontier 5: parallel for each node 𝑡 of 𝐹 do 6: for the 𝑗-th set bit of 𝑠𝑡 , at position 𝑏 do ′ 7: 𝐹 off ← (𝐹𝑡 ≪ 2) | 𝑏 ⊲ code of child 𝑗 𝑡 +𝑗 8: end for 9: end for 10: 𝐹 ← 𝐹′ ⊲ two buffers, alternated 11: else ⊲ last internal level: leaves visited in place 12: parallel for each node 𝑡 of 𝐹 do 13: for the 𝑗-th set bit of 𝑠𝑡 , at position 𝑏 do 14: 𝜆 ← signature of leaf off𝑡 + 𝑗
15: (𝑟, 𝑐) ← DeinterleaveMorton (𝐹𝑡 ≪ 2) | 𝑏 16: for each set bit of 𝜆, at position
2 𝑑𝑟 + 𝑑𝑐 do 17: 𝑦𝑟 +𝑑𝑟 ← 𝑦𝑟 +𝑑𝑟 ⊕ 𝑤 ⊗ 𝑥𝑐+𝑑𝑐 ⊲ atomic 18: end for 19: end for 20: end for 21: end if 22: end for 23: return 𝑦

atomic

plus-times R + × 0 1 𝑥𝑗 atomicAdd Boolean {0, 1} ∨ ∧ 0 1 [𝑥 𝑗 ≠ 0] atomicMax tropical R ≥0 ∪ {∞} min + ∞ 0 1 + 𝑥 𝑗 atomicMin

Algorithm 1 uses the following notation. Levels are numbered from 0, the root, to ℓ − 1, the leaves, each a 2 × 2 block of cells, so level ℓ − 2 is the last internal one; within a level, nodes are numbered from 0 in level order. The frontier 𝐹 holds the nodes of the current level 𝑑: 𝐹𝑡 is the Morton code of node 𝑡, which interleaves the row and column bits of its block, two bits 𝑏 1𝑏 0 per level of its path (Section 3.1). Its signature 𝑠𝑡 is its nibble in 𝑇 (Section 2), and its fan-out 𝑓𝑡 = popcount(𝑠𝑡 ), the number of set bits of 𝑠𝑡 , is its number of children. Level order lists the children of consecutive nodes consecutively, so those of node 𝑡 occupy positions off𝑡 to Í off𝑡 + 𝑓𝑡 − 1 of level 𝑑 + 1, where off𝑡 = 𝑢<𝑡 𝑓𝑢 . Computing these offsets for a level is the step ChildOffset (line 3), by a scan or from a directory (Section 3.2); it replaces the rank queries of a descent. Child 𝑗 of node 𝑡 (0 ≤ 𝑗 < 𝑓𝑡 ) is the quadrant 𝑏 of the 𝑗-th set bit of 𝑠𝑡 , and its code is its parent’s followed by the two bits of 𝑏 (line 7). The leaves, children of the last internal level, are never written to a frontier: the thread of node 𝑡 reads the signature 𝜆 of each of its leaves from 𝑇 (line 14), decodes the leaf’s code into the top-left cell (𝑟, 𝑐) of its block (line 15), and adds the contribution of each non-zero cell (𝑟 + 𝑑𝑟, 𝑐 + 𝑑𝑐), 𝑑𝑟, 𝑑𝑐 ∈ {0, 1}, to 𝑦 by one atomic operation of the semiring (line 17, Table 1). The loop starts from a cached level 𝑑 0 , so the device never expands the top of the tree (Section 3.4); with the cached frontier of Section 3.5, 𝑑 0 = ℓ − 2 and only the else branch runs. Section 3.3 follows the algorithm on Figure 1.

3.1

Frontier, Morton codes and leaves

A frontier entry is a 64-bit Morton code, the quadrants on the path from the root, two bits 𝑏 1𝑏 0 per level. Descending to quadrant 𝑏 appends its bits, code′ = (code ≪ 2) | 𝑏 1𝑏 0 : one shift and one OR on one word, where a coordinate pair would need two updates, 3

Tosoni et al.

𝑟 ′ = 2𝑟 + 𝑏 1 and 𝑐 ′ = 2𝑐 + 𝑏 0 , and the frontier stays one array that adjacent threads read with coalesced loads. At the leaves, the code is decoded once: the row concatenates the 𝑏 1 bits of the path, the column the 𝑏 0 bits. In Figure 1, the leaf of code 01 10 is block (1, 2), with origin (2, 4), and its only set cell contributes 𝑥 5 to 𝑦3 . Signatures are stored one per nibble, so reading a node is a byte load, a shift and a mask; the signature of a leaf marks which of its four cells are non-zero. Since the tree stores neither values nor coordinates, the leaves read no values, which a valued format loads for every non-zero, possibly through a decompressing accessor [36], and no indices.

3.2

In general, the threads of the last kernel are the non-empty 4 × 4 blocks, at most min(𝑚, 𝑛 2 /16) and from 0.9 to 296 million on our graphs; the tree-only engine launches at most 220 of them at a time, still far more than a GPU keeps resident. Load balance. The work of a thread is bounded by the arity of the tree, not by the degrees of the graph: an internal node writes at most 𝑘 2 = 4 codes, and a node of the last internal level issues 1 to 𝑘 4 = 16 updates. A CSR kernel with one thread per row inherits the degree distribution, so a hub row stalls its warp [13]; the usual remedies are a warp per row [13], a merge-based split of rows and non-zeros into equal shares [43], layouts that balance non-zeros across threads [3], or rows sorted by length (bCSR takes the faster of a thread and a warp per row). The 𝑘 2 -tree splits a hub row across one thread per non-empty 4 × 4 block, so the imbalance within a warp is at most sixteen to one on any graph. Sorting the frontier by work would break the contiguity of children on which the offsets rely; with a bound of sixteen, we do not sort. Memory access. Adjacent threads read adjacent data: the codes 𝐹𝑡 (256 contiguous bytes per warp), the signatures (16 bytes per warp), the offsets from the scan buffer or from a directory sample shared by 𝑆 nodes, and the signatures of the leaves, which follow those of the previous thread. Only the accesses to 𝑥 and 𝑦 are irregular, as in any SpMV; but consecutive codes are consecutive in Z-order, so on a clustered matrix the rows and columns a warp touches fall in a narrow range and share cache lines. The loops over the four bits of a signature are unrolled, so the threads of a warp diverge only in how many updates they issue. Granularity. One thread per leaf would write and read back the widest level of the tree (Section 3.4); one thread per node of a higher level descends several levels on its own, which the directory permits. On the seven Wikidata relations and swh, descending two levels per thread instead of one costs 1.02–1.28× the time tree-only and 1.34–2.36× with the frontier cached, and each further level costs more: the work bound of a thread grows fourfold per level while the number of threads shrinks.

Child offsets: a scan or a sampled directory

A level-synchronous traversal visits the nodes of a level in left-toright order, which is exactly the order rank counts. In Figure 1, the four nodes of level 1 have signatures 1000, 1010, 1111, 0001 and fan-outs 1, 2, 4, 1. The exclusive prefix sum (scan) of the fan-outs, the data-parallel primitive set out by Blelloch [14], gives each node 𝑡 the offset off𝑡 of its first child (0, 1, 3, 7), and their total (8) is the size of the next frontier. One scan per level thus replaces every rank of the level with no structure beyond the tree; it is the data-parallel form of the sequential traversal that lets a CPU SpMV over the whole tree do without rank [57], and the way GPU breadth-first search places the output of a frontier [45]. We compute it with the single-pass scan with decoupled look-back of Merrill and Garland in the CUB library [44], rather than with Blelloch’s two-sweep algorithm. Because the offsets depend only on the tree, they can also be stored. One 32-bit offset per node would cost eight times the tree, so the directory keeps one every 𝑆 = 256 nodes: for each block of 𝑆 consecutive nodes, the offset of the first child of its first node. The offset of any node is the sample of its block plus the popcount of the signatures that precede it in the block, i.e., of at most 𝑆/16 64-bit words of the tree. In Figure 1, with 𝑆 = 2, the directory stores the boxed offsets 0 and 3, and node 01 finds its first child at 0 + popcount(1000) = 1. This one-level sampled rank directory [47, §4.2.1] is lighter than the two-level structures that answer arbitrary rank queries on GPUs [32]: it costs 0.06–0.11 bits per node on our graphs, and a level becomes one kernel, with no scan. It is faster for single vectors; in a batch, the scan is paid once while the directory lookup is repeated per vector, and the scan wins (Section 5.5).

3.3

3.4

Bounded scratch: leaves in place, a host-side top, and tiles

Besides the matrix and the vectors, a product needs the scratch: the Morton codes of the frontier and, when scanning, the offsets. Three properties of the level-order layout keep it small. First, the scan reads the fan-outs straight from the tree, so no fan-out array is stored. Second, levels never shrink with depth, since every node has a child, so the leaves are the widest level; but the leaves of a node of the last internal level are contiguous, so its kernel visits them in place (the else branch of Algorithm 1) and the widest frontier is never written. Third, a kernel launch costs more than the work of the first, tiny levels, so we expand them once on the host and start every product from the deepest level with at most 216 nodes, whose codes take at most 216 × 8 bytes = 219 bytes (0.52 MB). The remaining scratch is as wide as the last internal levels and grows with the matrix (Figure 2a); tiles remove this dependence (Figure 2b). In level order, the children of a contiguous range of nodes form a contiguous range, so a level can be cut into tiles, ranges of nodes with at most 𝑀 children, each descended to the leaves before the next is expanded: a depth-first traversal over ranges

Threads and work, on the running example

Every kernel gives one thread to each node of its frontier, or 𝐵 threads in a batch, one per vector. In Figure 1, ℓ = 3 and the root is expanded on the host (Section 3.4), so a product starts from the codes 00, 01, 10, 11 of level 1 and is one kernel of four threads, one per non-empty 4 × 4 block, running the else branch of Algorithm 1. The fan-outs 1, 2, 4, 1 scan to the offsets 0, 1, 3, 7 (line 3), so thread 10 owns leaves 3 to 6. For each, it reads the leaf’s signature, appends the quadrant to its code, deinterleaves the result into the origin of a 2 × 2 block and issues one atomic per set cell (lines 14–17): leaf 5, of code 10 10 and origin (6, 0), holds the single cell (6, 1) and adds 𝑥 1 to 𝑦6 . The four threads issue 1, 3, 4 and 4 updates, which collide: 𝑦1 receives cell (1, 1) from thread 00 and (1, 5) from thread 01, and 𝑦6 two cells each from threads 10 and 11; hence the semiring’s atomic. 4

A GPU Engine for Sparse Matrix-Vector Products on 𝑘 2 -Tree-Compressed Graphs

(b) tiles of M =2 children

node. The left product reads it by rows but accumulates by columns: a transposed kernel gives row 𝑖 to a thread, which adds 𝑥𝑖 to 𝑦 𝑗 for every column 𝑗 of the row, so threads whose rows share a column update the same 𝑦 𝑗 atomically and serialize when they collide; the alternative is a second, column-oriented copy. Both directions are needed together by direction-optimizing traversals [12] and by regular path queries, which follow edges and inverse edges [7, 35]. A batch of 𝐵 vectors is an 𝑛 × 𝐵 row-major array, so a leaf visit updates 𝐵 contiguous outputs from 𝐵 adjacent threads: the traversal, the scans and the reads of the tree are paid once per batch, and only the leaf updates grow with 𝐵. Batched products take their offsets from the scan.

launches

(a) whole levels

top

d+1

last

(a) 3 launches, but each buffer is as wide as its level. (b) 10 launches, one tile at a time, depth-first, but only M slots per level. Lines: buffer slots. Rows: kernel launches. Segments: a node writes a child. Circles: leaves visited in place. Bands: a code waiting in its slot.

top

d+1

M last

4

Platform and protocol. All experiments run on an NVIDIA GB10 (Grace–Blackwell, sm_121, 128 GB of memory shared by CPU and GPU), with CUDA 13.0 and nvcc -O3; Section 5.7 repeats them on an H100. Following [58], we compress once and multiply many times: a matrix is built once on the host and multiplied by many vectors not known in advance. We time construction separately (Section 5.1) and products on the device, between two CUDA events. Times are medians over ten runs of 50 products for our engine and over five runs of each baseline binary; the interquartile range stays within 9% of the median for us and 32% for the baselines. Batches have 𝐵 = 128 vectors, or 16 on swh and it-2004, whose dense blocks 𝑋, 𝑌 ∈ R𝑛×128 would take 47 and 42 GB, too much to run every engine next to its matrix on the GB10. Dense vectors are FP32 in every engine. Energy per product is read from the NVML energy counter around a sustained loop of products of at least 1.5 s, run after the timings, as in [58]. Following [38], the tables give absolute times. Unless stated otherwise, single-vector products use the directory (𝑆 = 256), batched ones the scan, and tiles hold 220 nodes. Sizes are decimal (1 MB = 106 bytes, 1 GB = 109 bytes). Our code, the raw logs and the scripts that regenerate every table and figure will be released at https://github.com/ftosoni/g-k2-tree (to appear).

Figure 2: The descent of a small tree as a data-flow network, drawn after the scan networks of Merrill and Garland [44]: whole levels against tiles.

of nodes. One buffer of at most 𝑀 codes per level suffices, so the scratch is O (𝑀ℓ), bounded by the depth of the tree rather than its size: at most 8𝑀 bytes of codes per level below the cached top (8.4 MB with 𝑀 = 220 ), plus, when scanning, 4𝑀 bytes of offsets and the temporary storage of the scan. Each tile costs one launch, so the tile size trades scratch for launches. The sequence of tiles depends only on the tree: the host plans it, and allocates the exact scratch, once, before any product. For single vectors, we also record the launches of a product once as a CUDA graph, which the host resubmits with a single call instead of one call per kernel. We use 𝑀 = 220 ; matrices whose levels are narrower than 𝑀 are traversed as without tiles.

3.5

Tree-only and cached frontier

How much of the traversal to keep between products is a choice between memory and time, and we evaluate both ends. In the treeonly configuration (labeled “𝑘 2 ” in the tables), a matrix keeps the tree, its directory and the cached top levels, and the scratch is transient: one scratch can serve every matrix on the device, as suits a knowledge graph whose thousands of relations are thousands of matrices. In the cached-frontier configuration (labeled “𝑘 2 +F”), the Morton codes of the last internal level are computed when the matrix is loaded and kept with it, at eight bytes per node of that level. A product then expands no internal level: each thread takes its code from the cached array and its leaf offset from the directory, so a single-vector product is one kernel over the leaves (after the initialization of 𝑦). Both configurations are exact.

3.6

Experimental Setup

Datasets. We use fourteen graphs in two groups (Table 2). The first holds Boolean relations, the data the 𝑘 2 -tree is meant for in knowledge-graph engines. Seven are relations of the Wikidata knowledge graph [60], in the version of [8] (9.6×108 edges over 5,419 relations), from the data package of [58] (Zenodo doi:10.5281/ zenodo.22677747): six relations of shapes from 1.7×105 × 1.4×105 to 1.0 × 107 × 1.7 × 103 , and the citation relation wd_cites_work (1.7×108 edges). The eighth, swh, from the same package, is a subgraph of the Merkle DAG [51] of Software Heritage (SWH), the universal archive of public source code [2, 24]: the forward artifact graph of its popular-3k-python export, with 4.6×107 nodes and 1.2×109 edges, distributed in BVGraph, the format of the WebGraph framework of Boldi and Vigna [16, 30]. The second group holds six graphs of the Laboratory for Web Algorithmics (LAW), as distributed (BVGraphs in layered-label-propagation order [15]): four web crawls (indochina-2004, uk-2002, arabic-2005, it-2004, up to 1.2×109 edges), the data the 𝑘 2 -tree was designed for, and two harder cases, a social network (ljournal-2008) and the links of the English Wikipedia (enwiki-2013). The package of [58] ships

Left products and batches

The 𝑘 2 -tree treats rows and columns alike, so the left product 𝑦 ⊤ = 𝑥 ⊤𝐴 needs no other structure: the kernel reads the same tree, frontier and scratch, and only swaps the decoded coordinates at the leaf, so that cell (𝑟, 𝑐) adds 𝑥𝑟 to 𝑦𝑐 . Its memory differs from that of a right product only in the dense 𝑥 and 𝑦, whose lengths swap for a rectangular matrix. A CSR, like Boldi–Vigna-style compressors [16, 30], stores for each row its non-zero columns, the successors of a 5

Tosoni et al.

RePair grammars of the Wikidata relations and swh; we compress the LAW graphs with the RePair builder of [29].

Representation size (MB)

104

Baselines. For cuSPARSE [49], the standard, we time SpMV, transposed SpMV for the left product, and, for batches, the fastest of SpMM with its default algorithm, CSR_ALG2 and CSR_ALG3; its released baseline in [58] allocates managed as well as device memory, so we also time the same calls on device memory in our harness and report the faster. The grammar engine [58]1 is a recently proposed level-synchronous GPU engine over a RePair [40] grammar of the matrix; we run its released binaries, and its grammar file, entropy-coded with asymmetric numeral systems (ANS) [29], is the compressed baseline for space. A grammar could serve left products too, but the released engine implements only right ones, hence no left column in Table 3. We implemented two more baselines. bCSR, a Boolean CSR with row pointers and column indices but no values, is a simple baseline, not part of our proposal, that isolates what a CSR saves without values; it multiplies with one thread or one warp per row, the scalar and vector kernels of Bell and Garland [13], whichever is faster, with atomic additions for the left product and one thread per row and vector for a batch. A CPU descent of the same tree and directory, by rank, on the 20 host cores, shows what the GPU contributes.

103

CSR (valued) CSR (Boolean) grammar k2 -tree (ours) BVGraph

102 web crawls

101

106

107

108

109

Number of edges (nnz)

Figure 3: Representation size vs. number of edges (log– log). The 𝑘 2 -tree stays below the cuSPARSE CSR and bCSR throughout, more than tenfold on the web crawls; the BVGraph, where available, is smaller still.

5.2

Right, left and batched products

Right product. 𝑘 2 +F takes 0.8–1.7× the time of cuSPARSE on the six relations, faster on two, and 0.7–2.4× over the fourteen graphs, fastest on the citation relation and slowest on swh (2.4×), where the leaves dominate (Table 3). Tree-only, re-deriving the frontier at every product costs 1.2–5.6× the time of cuSPARSE, and the gap to 𝑘 2 +F follows the nodes per non-zero: 1.1–1.3× on the web crawls, with fewer than one, and 4.5× on the citation relation, with seven. 𝑘 2 +F is 1.1–2.9× faster than the grammar engine on ten graphs, including all seven Wikidata relations, and 1.5–2.3× slower on swh and three web crawls, whose repeated rows a grammar evaluates once. The CPU descent of the tree on 20 cores takes 5.4–23× the time of the tree-only GPU engine. bCSR is the fastest engine on twelve graphs (0.31–1.02× the time of cuSPARSE, from half the traffic): 𝑘 2 +F beats it only on the citation relation, while the treeonly engine keeps 1.1–14× fewer bytes per matrix. Left product. From the same tree and scratch, 𝑘 2 +F computes ⊤ 𝑥 𝐴 in 0.27–1.35× the time of cuSPARSE’s transposed kernel on the Wikidata relations (3.7× faster on the citation relation) and in 1.1–2.6× on the other seven graphs. The left product turns the gather of a right product into a scatter: 𝑦 𝑗 collects the contributions of every row with a non-zero in column 𝑗. This is a property of the product, not of either format, and it matters most on relations with many more rows than columns, whose updates collide on few entries of 𝑦 and serialize. It slows cuSPARSE by 1.2–12× over its own right product, ours by up to 9.1×, and bCSR, which scatters with atomic additions, by 1.5–21×: bCSR then needs 0.95–1.6× the time of cuSPARSE, the fastest engine on the left on ten graphs, and 𝑘 2 +F beats it on five, by 3–4% on three, 1.3× on wd_subclass_of and 3.8× on the citation relation. A CSR recovers its right-product speed on the left only with a second, transposed copy, at twice the memory; the tree serves both from one structure an order of magnitude smaller. Batches. Batching amortizes the traversal, which is where the two configurations differ: at 𝐵 = 128 the tree-only engine is within 16% of 𝑘 2 +F, which takes 1.8–3.2× the time of cuSPARSE’s best SpMM and is 1.1–2.4× faster than the grammar engine on every graph. bCSR, with one thread per row and vector, is the fastest engine

5 Results 5.1 Space and construction The size of the tree follows the clustering of the non-zeros rather than their number (Table 2 and fig. 3): from 2.4 bits per edge on a web crawl to 28 on the citation relation. It is 2.3–27× smaller than the cuSPARSE CSR and 1.2–14× smaller than bCSR, most on the web crawls (21–27× and 11–14×). Narrower indices would not close the gap: packed into ⌈log2 𝑛⌉ bits, the column indices of a web crawl take 23–26 bits per edge, 8–10× the tree. Against compressors that exploit repetition across rows, the picture is mixed. The tree is smaller than the RePair grammar, by up to 1.75×, on the six LAW graphs (third block of Table 2) and on wd_subclass_of and wd_country; on the other six graphs of the first two blocks, whose rows repeat most, the grammar is 0.62–0.84× the tree. BVGraph, which encodes each adjacency list by gaps and by reference to a similar nearby list, is 0.60–0.80× the tree wherever available, but serves successors only. The tree-only engine adds 1.5–13% to the tree; the cached frontier adds eight bytes per node of the last internal level, reaching 0.17–1.31× the cuSPARSE CSR: a faster configuration, not a more compact one. The host builds the tree on one core in 0.35–2.1 s on the six relations and at most 312 s on swh, the time of 1,400–3,300 singlevector products; caching the frontier adds 0.7–8.1%. RePair takes 0.65–0.95× as long on the Wikidata relations and 1.3–1.5× on the other seven graphs. Since we compress once and multiply many times (Section 4), construction amortized over more than 1,400– 3,300 products costs less per product than a product, and neither construction time decides between the formats.

1 https://github.com/ftosoni/g-mm-repair

6

A GPU Engine for Sparse Matrix-Vector Products on 𝑘 2 -Tree-Compressed Graphs

Table 2: The fourteen graphs, in three blocks (Wikidata relations; SWH subgraph; LAW web and social graphs): non-zeros (nnz) and 𝑘 2 -tree structure; persistent size (MB) of the bare tree, of our engine (𝑘 2 : tree, directory and cached top levels; 𝑘 2 +F: plus the cached frontier), of the RePair grammar, of the BVGraph (BV) where the graph is distributed as one (the rectangular Wikidata relations come as triples), of the cuSPARSE CSR (cuSP., FP32 values) and of bCSR; construction time on one host core (s), from the parsed edge list to the tree-only engine on the device, and of the RePair grammar [29]. 𝑘 2 -tree nnz

nodes

levels

bits/edge

tree

1.14M 1.03M 3.06M 4.60M 2.02M 10.09M 166.68M

5.82M 5.82M 7.72M 17.29M 5.71M 16.93M 1.18G

19 18 22 22 21 24 24

20.5 22.5 10.1 15.0 11.3 6.7 28.3

2.91 2.91 3.86 8.65 2.85 8.46 590

3.25 3.16 4.25 9.06 3.22 9.09 606

11.5 10.5 19.5 39.2 12.1 43.4 1,784

graph wd_sports_team wd_cast_member wd_citizenship wd_occupation wd_subclass_of wd_country wd_cites_work swh ljournal-2008 enwiki-2013 indochina-2004 uk-2002 arabic-2005 it-2004

persistent size (MB) k2

construction (s)

k2 +F

gram.

BV

cuSP.

bCSR

k2

RePair

1.93 2.21 3.25 5.38 3.19 9.63 459

– – – – – – –

10.4 8.96 36.0 50.6 22.1 121 1,362

5.87 4.83 23.7 32.2 14.0 80.6 695

0.349 0.346 0.644 1.34 0.381 2.10 88.4

0.310 0.300 0.560 0.870 0.360 1.69 83.8

1.22G

1.27G

26

4.2

633

644

3,009

490

378

9,931

5,057

312

446

79.02M 101.36M 194.11M 298.11M 640.00M 1.15G

281.82M 413.42M 116.75M 232.62M 439.70M 795.23M

23 23 23 25 25 26

14.3 16.3 2.4 3.1 2.7 2.8

141 207 58.4 116 220 398

145 212 59.4 118 223 404

533 805 261 520 988 1,786

180 220 79.1 203 302 555

109 166 35.4 84.6 147 261

654 828 1,583 2,459 5,211 9,371

338 422 806 1,267 2,651 4,768

27.4 35.3 37.9 61.9 138 272

38.3 54.1 49.3 90.8 193 345

Table 3: Per-vector time and working set on the fourteen graphs: right product 𝐴𝑥, left product 𝑥 ⊤𝐴 (cuSPARSE: transposed SpMV; bCSR: transposed kernel; the released grammar engine has none) and batched right product (𝐵=128; † : 𝐵=16), for our engine (tree-only 𝑘 2 , cached frontier 𝑘 2 +F), the grammar engine (gram.), bCSR and cuSPARSE (cuSP.). Working set: matrix, dense vectors and scratch on the device during a single product. left product 𝑥 ⊤𝐴 (ms)

right product 𝐴𝑥 (ms) k2

working set (MB)

gram.

bCSR

cuSP.

k2

k2 +F

bCSR

cuSP.

k2 +F

gram.

bCSR

cuSP.

k2

k2 +F

gram.

bCSR

cuSP.

wd_sports_team 0.252 0.079 wd_cast_member 0.233 0.074 wd_citizenship 0.421 0.177 wd_occupation 0.909 0.356 wd_subclass_of 0.213 0.082 wd_country 1.16 0.776 wd_cites_work 64.2 14.3

0.134 0.129 0.482 0.706 0.239 1.57 42.1

0.014 0.035 0.154 0.204 0.046 0.518 17.9

0.047 0.045 0.187 0.309 0.076 0.948 20.0

0.260 0.230 1.67 2.02 0.869 6.29 65.3

0.090 0.072 1.39 1.27 0.746 5.67 14.9

0.080 0.053 1.26 1.32 0.975 5.84 57.2

0.084 0.053 1.14 1.08 0.940 5.16 54.3

0.030 0.027 0.161 0.216 0.091 0.555 6.02

0.058 0.054 0.279 0.362 0.166 1.03 7.11

0.009 0.011 0.069 0.082 0.037 0.239 2.61

0.016 0.014 0.074 0.089 0.051 0.263 2.83

21.1 19.1 39.9 64.1 24.0 81.3 787

12.9 11.8 31.0 53.1 18.4 83.6 1,861

9.33 11.1 34.6 43.4 19.7 121 1,478

7.32 6.10 35.3 46.1 20.3 121 772

11.9 10.2 47.5 64.5 28.4 161 1,439

swh

138

108

51.9

40.1

44.9

192

161

72.6

62.8

16.3†

18.7†

4.41†

8.25†

1,115

3,374

2,793

5,422

10,296

ljournal-2008 enwiki-2013 indochina-2004 uk-2002 arabic-2005 it-2004

17.2 23.6 17.0 19.9 46.9 83.1

6.70 8.33 15.0 14.8 37.4 65.9

10.7 11.0 6.56 16.8 24.2 44.2

3.66 3.76 6.51 9.67 20.1 36.4

3.60 4.20 7.08 11.9 23.8 42.5

17.2 31.7 27.5 33.1 101 175

6.81 14.6 25.6 27.7 91.7 158

7.03 10.0 14.4 26.1 53.3 92.7

6.07 6.78 10.9 16.4 42.6 73.0

1.61 2.06 1.57 2.84 5.39 12.8†

3.05 3.19 2.73 6.90 9.35 17.0†

0.561 0.502 0.436 0.891 1.33 3.58†

0.681 0.636 0.637 1.33 1.88 6.90†

283 576 338 839 174 321 366 668 488 1,170 850 2,116

626 715 430 1,050 1,471 2,595

380 456 865 1,415 2,833 5,098

697 861 1,642 2,607 5,393 9,701

in batch on every graph, at 0.52–0.93× the time of cuSPARSE and 2.3–4.1× faster than 𝑘 2 +F: a batch amortizes its matrix reads as it does ours, and it has no traversal to pay. Across batch widths on the Wikidata relations (Figure 4), 𝑘 2 +F, already a single kernel, gains less than cuSPARSE’s SpMM: it matches cuSPARSE up to 𝐵 = 4 (geometric mean 0.98–1.10×) and settles at about twice its time from 𝐵 = 16. The tree-only engine gains most, 4.4× per vector from 𝐵 = 1 to 𝐵 = 128, and from 𝐵 = 64 is within about 10% of 𝑘 2 +F on average.

5.3

batched 𝐴𝑋 (ms/vector)

k2 +F

graph

needs less than cuSPARSE on twelve graphs of fourteen: 0.09– 0.14× on the web crawls and swh, 0.39–0.41× on ljournal-2008 and enwiki-2013, 0.55× on the citation relation and 0.50–0.99× on four other Wikidata relations. On the nine largest graphs it holds fewer bytes than the cached configuration, the grammar engine (1.5–3.1× fewer) and cuSPARSE (7–11× fewer on the web crawls and swh); on swh, 1.1 GB (the 633 MB tree, the two vectors, and 116 MB of directory and scratch) against 10.3 GB for cuSPARSE and 2.8 GB for the grammar engine. The cached frontier needs 0.2–1.3× the working set of cuSPARSE. On the two smallest relations both configurations hold more than cuSPARSE (1.8–1.9× tree-only), by at most 9.3 MB: their levels barely exceed one tile, and eight bytes per frontier code equal a CSR’s eight bytes per non-zero. bCSR holds 0.32–0.88× the working set of the tree-only engine on five of the six smaller Wikidata relations, whose trees of a few megabytes

Working set

In the memory-bounded setting, what matters is what an engine holds on the device during a product: the matrix, the dense vectors and the scratch (for our engine, every byte allocated). Figure 5 plots it against time, both relative to cuSPARSE. The tree-only engine 7

Tosoni et al.

Per-vector time / cuSPARSE

k2 (tree-only)

grammar

k2 +F (cached frontier)

the internal levels add only 1.1–1.3× to the time of the leaves. On the two smallest relations and the citation relation, a block holds 1.1 non-zeros and the tree 5.1–7.1 nodes per non-zero, so every update pays for its own path: the tree is only 2.3–3.6× smaller than the cuSPARSE CSR, and tree-only takes 3.2–5.4× the time of cuSPARSE. The second property is the length of the rows, which sets the efficiency of a CSR: it reads one row pointer and writes 𝑦 once per row, costs that long rows amortize and short rows pay per non-zero. On swh and the LAW graphs, with 15–28 non-zeros per row, cuSPARSE streams its matrix at 83–99% of the copy bandwidth (Section 5.6), and 𝑘 2 +F takes 1.2–2.4× its time despite the clustering, as the tree updates 𝑦 once per non-zero. On wd_citizenship, wd_occupation, wd_subclass_of and wd_country, with 1.0–1.4 non-zeros per row, 𝑘 2 +F takes 0.82–1.15× the time of cuSPARSE, less on two. The third property is the locality of the columns within a row, which sets the cost of gathering 𝑥. On the citation relation, the 24 non-zeros of a row fall on columns spread over 12 million, and cuSPARSE reaches only 71 GB/s; the tree visits the non-zeros in Z-order, which keeps the entries of 𝑥 and 𝑦 a warp touches in a narrow range (Section 3.3), and 𝑘 2 +F takes 0.72× the time of cuSPARSE. The left product adds the distribution of the non-zeros over the columns. wd_citizenship, wd_occupation and wd_country hold 430–5,800 non-zeros per column on 1,747–10,610 columns: their updates collide on few entries of 𝑦, and both engines take 3.5– 7.9× the time of their own right product. Where the non-zeros spread over many columns (wd_sports_team, wd_cast_member, ljournal-2008, the citation relation), the tree, which treats rows and columns alike, takes 0.97–1.14× its right-product time, and the transposed kernel of cuSPARSE 1.2–2.7×; hence the 3.7× lead of the tree on the left on the citation relation.

4

2

1 cuSPARSE

1

4

16

64

128

256

Batch width B

Figure 4: Per-vector time against batch width, relative to cuSPARSE at the same width (its SpMV at 𝐵 = 1, its fastest SpMM above): geometric mean (line) and range (band) over the seven Wikidata relations (log–log).

Single-vector time, relative to cuSPARSE

8 enwiki-2013 ljournal-2008

4 swh

2 web crawls

1 smaller and faster

0.5

1

5

swh 6

2

4 3

wd cites work

Boolean CSR grammar k2 (tree-only)

0.25

k2 +F (cached frontier) 0.1

0.2

cuSPARSE

0.5

1

2

Working set, relative to cuSPARSE

5.5

Figure 5: Working set during a single right product (matrix + dense vectors + scratch) against its time, both relative to cuSPARSE (log scales); a line joins the two configurations of a graph. Points 1–6 are the six relations in the row order of Table 2; the larger graphs are named.

We switch the optimizations of Section 3 on one at a time on the six relations, each step against the previous one. Feeding a preallocated scan directly from the tree removes most of the cost of a naive implementation (2.6–7.4×), allocation and launch overhead; visiting the leaves in place and expanding the top levels on the host gain 1.2–1.5×; tiles of 220 nodes cost at most 6% (0.83–1.06× the time); the directory gains 1.3–1.6× and the CUDA graph 1.02– 1.11×, completing the 𝑘 2 configuration; the cached frontier gains 1.5–3.4×, from wd_country to the two smallest relations. In all, a single-vector product drops from 3–8 ms to 0.07–0.77 ms. On the seven Wikidata relations and swh, the directory beats the scan on every single-vector product, by 1.10–1.57× tree-only and 1.03–1.45× with the frontier cached, and loses in batch by up to 29%, hence our split; on a single vector, its block size hardly matters (under 6% between 𝑆 = 64 and 256). Shrinking the tile from whole levels to 216 nodes (Figure 6) cuts the scratch 13–457× and makes it independent of the matrix, and down to 220 nodes the time barely changes (0.82–1.04× the untiled time). Below that, the fixed cost of a tile takes over: about 1 𝜇s with the CUDA graph, 3 𝜇s with separate launches, and 10 𝜇s with a scan per tile. On these eight graphs, 220 nodes is the fastest tile or within 8% of it.

weigh less than the scratch, and is within 2% of it on the citation relation; on wd_country, swh and the six LAW graphs, the tree-only engine holds 1.3–6.0× fewer bytes. On these eight graphs no engine holds fewer: the tree-only engine is the low-memory end of the trade-off of Figure 5, at 1.2–5.6× the time of cuSPARSE, and the cached frontier spends memory to move along it.

5.4

Where the time goes

Which graphs favor the tree

Three properties of a graph explain where the tree wins and loses. The first is the clustering of the non-zeros, which sets the work of the tree per update; it shows in the nodes per non-zero (Table 2) and, at the finest scale, in the non-zeros per non-empty 4 × 4 block, i.e., per thread of the last kernel. On the web crawls and swh, a block holds 4.1–7.7 non-zeros and the tree 0.60–1.04 nodes per non-zero: the tree is 16–27× smaller than the cuSPARSE CSR, and 8

A GPU Engine for Sparse Matrix-Vector Products on 𝑘 2 -Tree-Compressed Graphs

Time, relative to untiled

Time, relative to untiled

wd occupation

Table 4: Where the time of a single right product goes on the eight largest graphs: the share of the tree-only time spent expanding internal levels, the rate of the leaves (non-zeros per second of the leaf kernel), and the bandwidth of the two CSR engines over the bytes they must move (the matrix, 𝑥 read once, 𝑦 read and written once). Measured ceilings: consecutive atomic updates per second, and the copy bandwidth.

wd country

32

32

16

16

8

8

4

4

2

2

1

1 100

101

102

100

101

wd cites work

102

levels

swh

𝑘2

𝑘 2 +F

CSR (GB/s)

32

32

graph

bCSR

cuSP.

16

16

wd_cites_work

71

8.73

11.2

44

71

8

8

swh

20

11.1

11.4

135

238

4

4

2

2

ljournal-2008 enwiki-2013 indochina-2004 uk-2002 arabic-2005 it-2004

55 60 12 24 18 19

10.3 10.7 13.0 19.9 16.7 17.4

11.6 12.0 13.0 20.6 17.3 17.9

109 121 132 153 142 140

200 209 237 230 235 233

1

1

100

101

102

103

100

Frontier scratch (MB) scan per piece directory, one kernel graph

101

102

103

104

Frontier scratch (MB) directory

measured ceiling

tile of 220 nodes

Figure 6: Tile sweep of the tree-only configuration on four graphs: frontier scratch against single-vector time relative to the untiled traversal, from whole levels (rightmost point) down to 214 nodes (leftmost), with offsets from a scan per tile, from the directory, and from the directory replayed as one graph of kernels (the default). The tile of 220 nodes is circled.

5.6

leaves (G nnz/s)

(% of 𝑘 2 )

19.0

241

bytes of vector traffic per non-zero, against a fraction of that for the matrix. Tree-only, expanding the internal levels takes a further 12–24% of a product on the web crawls and swh, and 55–71% on the social graphs and the citation relation, which have more nodes per non-zero. The compressed matrix has thus left the critical path: a faster tree needs fewer bytes per update, not per node, such as leaves larger than 2×2 whose row segments accumulate in registers, or narrower vector entries.

Hardware efficiency

5.7

We time every kernel of a product with Nsight Systems and compare each engine with two ceilings that we measure (Table 4): the copy bandwidth, 241 GB/s, and the rate of vector updates 𝑦𝑜 += 𝑥𝑖 , which a microbenchmark puts at 19.0 G/s for consecutive entries with atomics, 19.4 G/s with plain stores, and 0.98 G/s for scattered entries. We report the eight largest graphs, whose CSR exceeds the caches many times over; on smaller ones, the caches of the GB10 serve part of the matrix. The CSR engines are bound by the bytes of their matrix. Counting only what they must move (the matrix, 𝑥 read once, 𝑦 read and written once), cuSPARSE streams 200–238 GB/s on seven graphs, 83–99% of the copy bandwidth, and bCSR 109–153 GB/s; on the citation relation, whose scattered reads of 𝑥 this count leaves out, 71 and 44 GB/s. Their time falls only with the bytes of the matrix, eight per non-zero for cuSPARSE and four for bCSR. The tree is bound by its vectors instead. It reads 0.30–0.52 bytes of matrix per non-zero on the web crawls and swh (1.2–2.2 with the frontier cached, whose codes it reads instead of the internal levels), and its kernels keep the GPU busy for at least 99.8% of a product (96.5% on the smallest relations). Its leaves sustain 13.0–20.6 G non-zeros per second on the web crawls, the rate of consecutive updates (slightly above it on uk-2002, as the four cells of a leaf share cache lines), and 8.7–12.0 G/s on the other graphs, whose non-zeros cluster less. A CSR accumulates a row in a register and writes 𝑦 once per row; the tree enumerates the non-zeros of a block in Morton order, so it updates 𝑦 once per non-zero: up to twelve

Across GPU architectures

We repeat the measurements on an NVIDIA H100 of the GRICAD Kraken cluster2 (sm_90, discrete HBM, CUDA 13.1), with the same sources, inputs, protocol and batch widths (Table 5); every crosscheck passes again. The H100 copies 8.9× faster and launches kernels slightly faster, but replays a CUDA-graph node in 0.79 𝜇s against 0.41. Every engine gains, unequally: on a single-vector right product, cuSPARSE runs 6.8× faster than on the GB10 (geometric mean, 2.6–12.8×), bCSR 6.3×, the grammar engine 5.4×, and our engine 3.5× in both configurations (2.1–5.2×). Our factor to cuSPARSE on a single-vector right product thus about doubles: 𝑘 2 +F takes 2.68× its time (1.24–5.09×) against 1.40×; in batch it grows from 2.20× to 3.45×, and on the left from 1.34× to 1.71×. The grammar engine now beats 𝑘 2 +F on a single vector on eight graphs, against four, while in batch we stay faster on eleven, against fourteen. This matches what each engine reads: a CSR is streamed and follows bandwidth, whereas a tree product is a chain of short dependent kernels with gathered reads and atomic updates. The left product gains least, 2.1× for us and 2.7× for cuSPARSE, and nothing on the three relations whose 3.1–10.1 million non-zeros fall on at most 10,610 columns (0.96–1.18× their GB10 time for both), as colliding atomics serialize whatever the memory; on wd_subclass_of, ours even takes 1.28× its GB10 time, and cuSPARSE 1.02×. Of the launch-dependent settings, the source of offsets carries over: on 2 https://gricad.univ-grenoble-alpes.fr/en/news/2026/inauguration-kraken/

9

Tosoni et al.

Table 5: The engine on two GPUs: the GB10 (CUDA 13.0, driver 580.95.05), with memory shared by CPU and GPU, and an H100 with discrete memory (CUDA 13.1, driver 610.57.04). Device properties, measured copy bandwidth and launch costs; the time of each engine divided by that of cuSPARSE on the same GPU, i.e., the speedup of cuSPARSE over it (below 1: the engine is faster), as geometric mean over the fourteen graphs with the range in brackets; and the launch-dependent settings, chosen by geometric mean over the Wikidata relations and swh, tree-only.

GPU compute capability SMs memory (GB) L2 cache (MiB) copy bandwidth (GB/s) kernel launch (𝜇s) graph node (𝜇s)

GB10

H100

NVIDIA GB10 12.1 48 128 24 241 2.04 0.41

NVIDIA H100 9.0 132 100 60 2152 1.72 0.79

Table 6: Multi-source search to the fixed point, along the edges (fwd) and against them (rev): sources 𝐵, levels (∗ : stopped at 32 levels), time, and device memory (matrix and dense frontier tables). cuSPARSE uses one CSR per direction.

graph wd_subclass_of fwd rev wd_cites_work fwd rev swh fwd rev

5.9

device (GB) k2 k2 +F cuSP.

128 64 1.08 128 44 0.68 128 14 10.10 128 13 9.58 16 32∗ 8.18 16 32∗ 8.12

1.06 0.67 9.02 8.51 7.06 6.95

1.01 1.93 1.92 1.93 0.57 1.93 1.92 1.94 6.56 16.39 17.46 17.09 6.29 16.39 17.46 17.09 7.03 8.06 10.32 17.44 4.82 8.06 10.32 17.61

Correctness

Every engine is checked on a shared input: the grammar engine emits 𝑥 and its 𝑦, and cuSPARSE, a CPU CSR oracle and our engine recompute 𝑦 on the same 𝑥 (for the tropical semiring, which cuSPARSE lacks, against a CPU min-plus product, as in [58]). They agree on every relation and semiring, bit for bit on the Boolean and tropical ones and to single precision on plus-times (|Δ| < 10−6 , from the order in which atomics land). We also check both configurations, right and left, on all three semirings against a CPU product on all fourteen graphs, and both offset sources and every tile size on the Wikidata relations and swh. On all fourteen graphs, bCSR matches the CPU reference on its right, left and batched products, and so does the CPU descent on its right product.

directory scan 222

the Wikidata relations and swh, tree-only, the directory takes 0.67× the time of a scan per tile on a single vector (0.73× on the GB10), and in batch the scan stays faster, the directory taking 1.34× its time (1.09×). The tile size does not: with faster kernels and slower graph nodes, the fixed cost of a tile weighs more, and tiles of 220 nodes take 1.16× the time of whole levels tree-only (0.96× on the GB10), which tiles of 222 nodes bring to 1.02× at four times the scratch; the cached configuration stays within 1% at every tile on both GPUs. All H100 times above use tiles of 220 nodes. Every one of our matrices fits in the memory of the H100, so the case that motivates the tree on discrete GPUs, a matrix whose CSR does not fit, remains untested (Section 6).

5.8

time (s) k2 k2 +F cuSP.

uses right and left products of one tree; cuSPARSE SpMM receives a CSR oriented for each direction, so it always runs its fast path. Both produce identical distance tables, which a CPU search confirms on four of the sources. With the frontier cached, we take 1.00–1.44× the time of cuSPARSE (1.07–1.69× tree-only): at this batch width the dense tables, 42–99% of what cuSPARSE holds on the device, set the pace. The tree pays off in the matrix: on swh, the search holds 8.1 GB on the device against 17.4 GB.

time relative to cuSPARSE right, 𝑘 2 2.87 (1.22–5.62) 5.47 (2.17–14.3) right, 𝑘 2 +F 1.40 (0.713–2.41) 2.68 (1.24–5.09) right, grammar 1.87 (0.927–3.15) 2.35 (1.12–3.96) left, 𝑘 2 +F 1.34 (0.274–2.57) 1.71 (0.869–2.92) batched, 𝑘 2 +F 2.20 (1.81–3.23) 3.45 (1.82–7.87) settings that depend on the cost of a launch offsets, single directory offsets, batched scan tile (nodes) 220

𝐵 lev.

6

Limitations (1) Time against a CSR without values. bCSR is the fastest engine on the right on twelve graphs and in batch on all fourteen, 2.3–4.1× faster than 𝑘 2 +F; on the left, where it must scatter, the tree beats it on five. The tree buys memory with that time: on eight of the nine largest graphs it holds 1.3–6.0× fewer bytes than bCSR. (2) The tree needs clustering. On the citation relation, whose non-zeros show little locality, the tree has seven nodes per non-zero and is only 2.3× smaller than the cuSPARSE CSR, against 3.0× for the grammar; tree-only it is 3.2× slower than cuSPARSE, and the cached frontier recovers the speed only by holding more bytes than the CSR. (3) Compressors that exploit repetition are smaller. BVGraph is 0.60–0.80× the tree, the grammar is smaller on six graphs, and the tree needs 14–16 bits per edge on the social network and the English Wikipedia. In exchange, the tree offers right and left products from one structure and a data-parallel traversal.

End to end: multi-source search with right and left products

We close with an end-to-end application: hop distances from 𝐵 sources at once, a multi-source breadth-first search, by iterating the batched Boolean product to its fixed point, one level per product, along the edges (right products) and against them (left products). It yields the distances that Bellman–Ford computes on unit lengths by iterating the tropical product (Section 2), but each step propagates only the Boolean frontier of newly reached entries, which a small kernel marks in dense 𝑛 × 𝐵 tables of frontier and distances. We run it on the two Wikidata relations whose rows and columns are the same entities and on swh, the largest graph (Table 6). Our engine 10

A GPU Engine for Sparse Matrix-Vector Products on 𝑘 2 -Tree-Compressed Graphs

(4) Time on the largest graphs. There the leaves dominate: 𝑘 2 +F takes 1.2–2.4× the time of cuSPARSE on a single-vector right product and 1.5–2.7× its energy, and the grammar engine is faster on swh and three web crawls. The gain is in the working set. (5) Two GPUs, both larger than our matrices. All our matrices fit on both GPUs, so the reach beyond device memory remains to be shown, as does a parallel construction (now up to 312 s on one core). On the H100, cuSPARSE gains more from the faster memory, our factor to it on a single-vector right product about doubles, and the best tile is four times larger (Section 5.7).

7

at tens rather than thousands of GB/s. And the level-synchronous traversal could extend to products of two trees, whose CPU counterparts are the recursive 𝑘 2 -tree product [7] and its prejoining-based successor on qdags [6]: a frontier of node pairs would expand into matching quadrant sub-products, keeping Boolean 𝐴 · 𝐵, and by repeated squaring the transitive closure, compressed.

Acknowledgments All the H100 computations presented in this paper were performed using the GRICAD infrastructure (https://gricad.univ-grenoblealpes.fr), which is supported by Grenoble research communities. Francesco Tosoni was funded by the Alfred P. Sloan Foundation with the grant #G-2025-25193 (sloan.org). Francesco Tosoni has received travel reimbursement from SWH (Inria Paris), with which he collaborates on archival and retrieval technologies, and a travel scholarship from the Wikimedia Foundation; the SWH archive, Wikidata, and Wikipedia are among the data sources of this paper. These relationships are academic and did not influence the work reported here. Gabriele Mencagli declares no competing interests.

Related Work

Sums and products of 𝑘 2 -trees, and the regular path queries built on them, run on the CPU with one dependent rank per visited node [7], next to other compact query engines [8]; depth-first [20] and depthfirst unary degree sequence (DFUDS) [28] layouts improve their locality. Pointer-based quadtrees with CSR or dense leaves store valued sparse matrices, but their control overhead makes SpMV slower than CSR [61]; the 𝑘 2 -tree is their succinct, pointerless form. Regular path queries have been cast as sparse Boolean linear algebra [35], also over 𝑘 2 -trees [7] and over qdags, their generalization to higherarity relations [6], and run on GPUs over uncompressed graphs [50]. GPU breadth-first search [45, 46] places the output of each frontier element by a prefix sum [14]. GPU-based traversal on compressed graphs (GCGT) [53] decodes adjacency lists encoded Boldi–Vignastyle, by intervals and gaps under variable-length codes, to traverse graphs beyond device memory, and schedules the decoding to tame divergence; a 𝑘 2 -tree node is four bits read in place, and the tree serves right and left products. The grammar engine [58] has the level-synchronous, double-buffered structure we adopt, over a RePair grammar [40]. Numerical work on memory-bound kernels shortens the indices and values of sparse matrices [3, 52] or stores values, such as a Krylov basis, in fewer bits than they are computed in [36, 37]; we take the Boolean, lossless extreme, the memory corner of the read–update–memory trade-off [9]. GPU engines for iterative relational algebra and Datalog compute closures over uncompressed relations as fixed points of joins [54, 55].

8

References [1] Amir Abboud, Arturs Backurs, Karl Bringmann, and Marvin Künnemann. 2020. Impossibility Results for Grammar-Compressed Linear Algebra. In Advances in Neural Information Processing Systems, H. Larochelle, M. Ranzato, R. Hadsell, M.F. Balcan, and H. Lin (Eds.), Vol. 33. Curran Associates, Inc., Red Hook, NY, USA, 8810–8823. https://proceedings.neurips.cc/paper_files/paper/2020/file/ 645e6bfdd05d1a69c5e47b20f0a91d46-Paper.pdf [2] Jean-François Abramatic, Roberto Di Cosmo, and Stefano Zacchiroli. 2018. Building the universal archive of source code. Commun. ACM 61, 10 (2018), 29–31. doi:10.1145/3183558 [3] José Ignacio Aliaga, Hartwig Anzt, Enrique S. Quintana-Ortí, Andrés E. Tomás, and Yuhsiang M. Tsai. 2021. Balanced and Compressed Coordinate Layout for the Sparse Matrix-Vector Product on GPUs. In Euro-Par 2020: Parallel Processing Workshops, Bartosz Balis, Dora B. Heras, Laura Antonelli, Andrea Bracciali, Thomas Gruber, Jin Hyun-Wook, Michael Kuhn, Stephen L. Scott, Didem Unat, and Roman Wyrzykowski (Eds.). Springer International Publishing, Cham, 83– 95. [4] João Nuno Ferreira Alves, Samir Moustafa, Siegfried Benkner, Alexandre P. Francisco, Wilfried N. Gansterer, and Luís M. S. Russo. 2025. Accelerating Graph Neural Networks Using a Novel Computation-Friendly Matrix Compression Format. In IEEE International Parallel and Distributed Processing Symposium, IPDPS 2025, Milano, Italy, June 3-7, 2025. IEEE, Milano, Italy, 1091–1103. doi:10. 1109/IPDPS64566.2025.00100 [5] Nicolás Araya, José Zagal, Cecilia Hernández, and Gonzalo Navarro. 2026. Compressed Boolean Matrix Multiplication using Bicliques. arXiv:2610.32836 [cs.DS] https://arxiv.org/abs/2610.32836 [6] Diego Arroyuelo, José Cazorla, and Gonzalo Navarro. 2026. Boosting Graph Joins and Matrix Multiplications in Little Space. In Proceedings of the 9th Joint Workshop on Graph Data Management Experiences & Systems (GRADES) and Network Data Analytics (NDA) (GRADES-NDA 2026). ACM, Bengaluru, India. [7] Diego Arroyuelo, Adrián Gómez-Brandón, and Gonzalo Navarro. 2024. Evaluating regular path queries on compressed adjacency matrices. The VLDB Journal 34, 1 (26 Nov 2024), 2. doi:10.1007/s00778-024-00885-6 [8] Diego Arroyuelo, Aidan Hogan, Gonzalo Navarro, and Javiel Rojas-Ledesma. 2022. Time- and Space-Efficient Regular Path Queries. In 2022 IEEE 38th International Conference on Data Engineering (ICDE). IEEE, Kuala Lumpur, Malaysia, 3091–3105. doi:10.1109/ICDE53745.2022.00277 [9] Manos Athanassoulis, Michael S. Kester, Lukas M. Maas, Radu Stoica, Stratos Idreos, Anastasia Ailamaki, and Mark Callaghan. 2016. Designing Access Methods: The RUM Conjecture. In Proceedings of the 19th International Conference on Extending Database Technology (EDBT). OpenProceedings.org, Bordeaux, France, 461–466. doi:10.5441/002/edbt.2016.42 [10] Sebastian Baunsgaard and Matthias Boehm. 2023. AWARE: Workload-aware, Redundancy-exploiting Linear Algebra. Proc. ACM Manag. Data 1, 1, Article 2 (May 2023), 28 pages. doi:10.1145/3588682 [11] Sebastian Baunsgaard and Matthias Boehm. 2025. Morphing-Based Compression for Data-Centric ML Pipelines. Proc. VLDB Endow. 19, 3 (Nov. 2025), 440–454. doi:10.14778/3778092.3778104 [12] Scott Beamer, Krste Asanović, and David Patterson. 2013. Direction-Optimizing Breadth-First Search. Scientific Programming 21, 3-4 (2013), 702694. doi:10.3233/

Conclusions

A 𝑘 2 -tree can be multiplied efficiently on a GPU once its traversal is organized by levels: the traversal does not depend on the vector, so the child offsets that a descent obtains by dependent rank lookups become static data, scanned at each product or read from a sampled directory, and tiles bound the memory of a product by the depth of the tree. The resulting engine computes right and left products, single or batched, over three semirings, from a structure 2.3–27× smaller than the cuSPARSE CSR, within a small factor of cuSPARSE in time and faster than a grammar-compressed engine in batch; bCSR is faster on most graphs but holds up to 6× more bytes, so the 𝑘 2 -tree is the format for when memory, not time, is the constraint. Two directions remain for future work. On discrete GPUs, a bounded scratch over a tree an order of magnitude below the CSR could admit graphs whose CSR exceeds device memory and would otherwise cross the host interconnect at every product, 11

Tosoni et al.

[34] Gene H. Golub and Charles F. Van Loan. 2013. Matrix Computations - 4th Edition. Johns Hopkins University Press, Philadelphia, PA. doi:10.1137/1.9781421407944 [35] Semyon V. Grigorev, Georgiy Belyanin, and Rodion Suvorov. 2025. Single-Source Regular Path Querying in Terms of Linear Algebra. In Proceedings of Workshops at the 51st International Conference on Very Large Data Bases, VLDB 2025, London, UK, September 1-5, 2025, Norman W. Paton and John Paparrizos (Eds.). VLDB.org, London, UK. https://www.vldb.org/2025/Workshops/VLDB-Workshops-2025/ LSGDA/LSGDA25_12.pdf [36] Thomas Grützmacher, Hartwig Anzt, and Enrique S. Quintana-Ortí. 2023. Using Ginkgo’s memory accessor for improving the accuracy of memory-bound low precision BLAS. Software: Practice and Experience 53, 1 (2023), 81–98. doi:10. 1002/spe.3041 [37] Thomas Grützmacher, Robert Underwood, Sheng Di, Franck Cappello, and Hartwig Anzt. 2024. FRSZ2 for In-Register Block Compression Inside GMRES on GPUs. In SC24-W: Workshops of the International Conference for High Performance Computing, Networking, Storage and Analysis. IEEE, Atlanta, GA, USA, 240–249. doi:10.1109/SCW63240.2024.00038 [38] Torsten Hoefler and Roberto Belli. 2015. Scientific benchmarking of parallel computing systems: twelve ways to tell the masses when reporting performance results. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (Austin, Texas) (SC ’15). Association for Computing Machinery, New York, NY, USA, Article 73, 12 pages. doi:10.1145/2807591.2807644 [39] Susana Ladra, José R. Paramá, and Fernando Silva-Coira. 2017. Scalable and queryable compressed storage structure for raster data. Information Systems 72 (2017), 179–204. doi:10.1016/j.is.2017.10.007 [40] N. Jesper Larsson and Alistair Moffat. 1999. Offline Dictionary-Based Compression. In Proceedings of the Conference on Data Compression (DCC ’99). IEEE Computer Society, USA, 296. [41] Fengan Li, Lingjiao Chen, Yijing Zeng, Arun Kumar, Xi Wu, Jeffrey F. Naughton, and Jignesh M. Patel. 2019. Tuple-oriented Compression for Large-scale Mini-batch Stochastic Gradient Descent. In Proceedings of the 2019 International Conference on Management of Data (Amsterdam, Netherlands) (SIGMOD ’19). Association for Computing Machinery, New York, NY, USA, 1517–1534. doi:10.1145/3299869.3300070 [42] Giosuè Cataldo Marinò, Flavio Furia, Dario Malchiodi, and Marco Frasca. 2023. Efficient and Compact Representations of Deep Neural Networks via Entropy Coding. IEEE Access 11 (2023), 106103–106125. doi:10.1109/ACCESS.2023.3317293 [43] Duane Merrill and Michael Garland. 2016. Merge-based parallel sparse matrixvector multiplication. In Proceedings of the International Conference for High Performance Computing, Networking, Storage and Analysis (Salt Lake City, Utah) (SC ’16). IEEE Press, Article 58, 12 pages. doi:10.1109/SC.2016.57 [44] Duane Merrill and Michael Garland. 2016. Single-pass Parallel Prefix Scan with Decoupled Look-back. Technical Report NVR-2016-002. NVIDIA Corporation. https://research.nvidia.com/sites/default/files/pubs/2016-03_Single-passParallel-Prefix/nvr-2016-002.pdf [45] Duane Merrill, Michael Garland, and Andrew Grimshaw. 2012. Scalable GPU graph traversal. In Proceedings of the 17th ACM SIGPLAN Symposium on Principles and Practice of Parallel Programming (New Orleans, Louisiana, USA) (PPoPP ’12). Association for Computing Machinery, New York, NY, USA, 117–128. doi:10. 1145/2145816.2145832 [46] Duane Merrill, Michael Garland, and Andrew Grimshaw. 2015. HighPerformance and Scalable GPU Graph Traversal. ACM Trans. Parallel Comput. 1, 2, Article 14 (Feb. 2015), 30 pages. doi:10.1145/2717511 [47] Gonzalo Navarro. 2016. Compact Data Structures: A Practical Approach. Cambridge University Press, Cambridge, UK. [48] Gonzalo Navarro. 2026. Computation-Friendly Compression of Matrices: Speeding up Linear and Boolean Algebra. arXiv:2610.35362 [cs.DS] https://arxiv.org/ abs/2610.35362 [49] NVIDIA Corporation. 2026. cuSPARSE Library. https://docs.nvidia.com/cuda/ cusparse/index.html. https://developer.nvidia.com/cusparse Version 13.x (part of CUDA Toolkit). [50] Sungwoo Park, Seohyeon Kim, and Min-Soo Kim. 2026. cuRPQ: A HighPerformance GPU-Based Framework for Processing Regular and Conjunctive Regular Path Queries. Proc. ACM Manag. Data 4, 3, Article 156 (May 2026), 28 pages. doi:10.1145/3802033 [51] Antoine Pietri, Diomidis Spinellis, and Stefano Zacchiroli. 2020. The Software Heritage Graph Dataset: Large-scale Analysis of Public Software Development History. In Proceedings of the 17th International Conference on Mining Software Repositories (Seoul, Republic of Korea) (MSR ’20). Association for Computing Machinery, New York, NY, USA, 1–5. doi:10.1145/3379597.3387510 [52] Emil Schatzle, Tommaso Pegolotti, and Markus Puschel. 2026. Fast Entropy Decoding for Sparse MVM on GPUs . In 2026 IEEE International Parallel and Distributed Processing Symposium (IPDPS). IEEE Computer Society, Los Alamitos, CA, USA, 416–432. doi:10.1109/IPDPS65963.2026.00044

SPR-130370 [13] Nathan Bell and Michael Garland. 2009. Implementing sparse matrix-vector multiplication on throughput-oriented processors. In Proceedings of the Conference on High Performance Computing Networking, Storage and Analysis (Portland, Oregon) (SC ’09). Association for Computing Machinery, New York, NY, USA, Article 18, 11 pages. doi:10.1145/1654059.1654078 [14] Guy E. Blelloch. 1990. Prefix Sums and Their Applications. Technical Report CMU-CS-90-190. School of Computer Science, Carnegie Mellon University. [15] Paolo Boldi, Marco Rosa, Massimo Santini, and Sebastiano Vigna. 2011. Layered label propagation: a multiresolution coordinate-free ordering for compressing social networks. In Proceedings of the 20th International Conference on World Wide Web (Hyderabad, India) (WWW ’11). Association for Computing Machinery, New York, NY, USA, 587–596. doi:10.1145/1963405.1963488 [16] P. Boldi and S. Vigna. 2004. The webgraph framework I: compression techniques. In Proceedings of the 13th International Conference on World Wide Web (New York, NY, USA) (WWW ’04). Association for Computing Machinery, New York, NY, USA, 595–602. doi:10.1145/988672.988752 [17] Nieves R. Brisaboa, Guillermo De Bernardo, Gilberto Gutierrez, Susana Ladra, Miguel R. Penabad, and Brunny A. Troncoso. 2015. Efficient Set Operations over 𝑘 2 -Trees. In 2015 Data Compression Conference. IEEE, Snowbird, UT, USA, 373–382. doi:10.1109/DCC.2015.9 [18] Nieves R. Brisaboa, Susana Ladra, and Gonzalo Navarro. 2009. k2-Trees for Compact Web Graph Representation. In String Processing and Information Retrieval, Jussi Karlgren, Jorma Tarhio, and Heikki Hyyrö (Eds.). Springer Berlin Heidelberg, Berlin, Heidelberg, 18–30. [19] Nieves R. Brisaboa, Susana Ladra, and Gonzalo Navarro. 2014. Compact representation of Web graphs with extended functionality. Information Systems 39 (2014), 152–174. doi:10.1016/j.is.2013.08.003 [20] Gabriel Carmona and Giovanni Manzini. 2025. Depth First Representations of k2-trees. In String Processing and Information Retrieval: 32nd International Symposium, SPIRE 2025, London, UK, September 8–11, 2025, Proceedings (London, United Kingdom). Springer-Verlag, Berlin, Heidelberg, 28–44. doi:10.1007/978-3032-05228-5_4 [21] Thomas H. Cormen, Charles E. Leiserson, Ronald L. Rivest, and Clifford Stein. 2009. Introduction to Algorithms (3rd ed.). The MIT Press, Cambridge, MA. [22] Timothy A. Davis. 2019. Algorithm 1000: SuiteSparse:GraphBLAS: Graph Algorithms in the Language of Sparse Linear Algebra. ACM Trans. Math. Software 45, 4 (2019), 44:1–44:25. doi:10.1145/3322125 [23] Guillermo de Bernardo, Travis Gagie, Susana Ladra, Gonzalo Navarro, and Diego Seco. 2023. Faster compressed quadtrees. J. Comput. System Sci. 131 (2023), 86–104. doi:10.1016/j.jcss.2022.09.001 [24] Roberto Di Cosmo. 2020. Archiving and Referencing Source Code with Software Heritage. In Mathematical Software – ICMS 2020, Anna Maria Bigatti, Jacques Carette, James H. Davenport, Michael Joswig, and Timo de Wolff (Eds.). Springer International Publishing, Cham, 362–373. [25] Daniel Donenfeld, Stephen Chou, and Saman Amarasinghe. 2022. Unified Compilation for Lossless Compression and Sparse Computing. In 2022 IEEE/ACM International Symposium on Code Generation and Optimization (CGO). IEEE, Seoul, Republic of Korea, 205–216. doi:10.1109/CGO53902.2022.9741282 [26] Ahmed Elgohary, Matthias Boehm, Peter J. Haas, Frederick R. Reiss, and Berthold Reinwald. 2016. Compressed linear algebra for large-scale machine learning. Proc. VLDB Endow. 9, 12 (Aug. 2016), 960–971. doi:10.14778/2994509.2994515 [27] Ahmed Elgohary, Matthias Boehm, Peter J. Haas, Frederick R. Reiss, and Berthold Reinwald. 2019. Compressed linear algebra for declarative large-scale machine learning. Commun. ACM 62, 5 (April 2019), 83–91. doi:10.1145/3318221 [28] Antonio Fariña, Adrián Gómez-Brandón, Asunción Gómez-Colomer, and Gonzalo Navarro. 2026. Cache-Friendly Compressed Boolean Matrices. In String Processing and Information Retrieval, Golnaz Badkobeh, Jakub Radoszewski, Nicola Tonellotto, and Ricardo Baeza-Yates (Eds.). Springer Nature Switzerland, Cham, 95–108. doi:10.1007/978-3-032-05228-5_9 [29] P. Ferragina, G. Manzini, T. Gagie, D. Köppl, G. Navarro, M. Striani, and F. Tosoni. 2022. Improving matrix-vector multiplication via lossless grammar-compressed matrices. Proc. VLDB Endowment 15, 10 (2022), 2175–2187. [30] Tommaso Fontana, Sebastiano Vigna, and Stefano Zacchiroli. 2024. WebGraph: The Next Generation (Is in Rust). In Companion Proceedings of the ACM Web Conference 2024 (Singapore, Singapore) (WWW ’24). Association for Computing Machinery, New York, NY, USA, 686–689. doi:10.1145/3589335.3651581 [31] Alexandre P. Francisco, Travis Gagie, Dominik Köppl, Susana Ladra, and Gonzalo Navarro. 2022. Graph Compression for Adjacency-Matrix Multiplication. SN Computer Science 3, 3 (21 Mar 2022), 193. doi:10.1007/s42979-022-01084-2 [32] Marco Franzreb, Martin Burtscher, and Stephan Rudolph. 2025. GPU Implementation of the Wavelet Tree. arXiv:2505.03372 [cs.DS] https://arxiv.org/abs/2505. 03372 [33] Sandra Álvarez García, Nieves R. Brisaboa, Guillermo de Bernardo, and Gonzalo Navarro. 2014. Interleaved 𝑘 2 -Tree: Indexing and Navigating Ternary Relations. In 2014 Data Compression Conference. IEEE, Snowbird, UT, USA, 342–351. doi:10. 1109/DCC.2014.56 12

A GPU Engine for Sparse Matrix-Vector Products on 𝑘 2 -Tree-Compressed Graphs

’16). Association for Computing Machinery, New York, NY, USA, 1875–1884. doi:10.1145/2939672.2939864 [57] F. Tosoni, P. Bille, V. Brunacci, A. De Angelis, P. Ferragina, and G. Manzini. 2025. Toward greener matrix operations by lossless compressed formats. IEEE Access 13 (2025), 56756–56773. [58] Francesco Tosoni and Gabriele Mencagli. 2027. Streaming Right Multiplication over Grammar-Compressed Matrices: A Memory-Bounded GPU Engine for Genotype and Graph Data. In Proceedings of the SIAM Symposium on Algorithm Engineering and Experiments (ALENEX). SIAM. To appear. [59] Luca Versari, Iulia-Maria Comsa, Alessio Conte, and Roberto Grossi. 2020. Zuckerli: A New Compressed Representation for Graphs. IEEE Access 8 (2020), 219233– 219243. doi:10.1109/ACCESS.2020.3040673 [60] Denny Vrandečić and Markus Krötzsch. 2014. Wikidata: a free collaborative knowledgebase. Commun. ACM 57, 10 (Sept. 2014), 78–85. doi:10.1145/2629489 [61] Ivan Šimeček. 2009. Sparse Matrix Computations Using the Quadtree Storage Format. In 2009 11th International Symposium on Symbolic and Numeric Algorithms for Scientific Computing (SYNASC). IEEE, Timisoara, Romania, 168–173. doi:10.1109/SYNASC.2009.55

[53] Mo Sha, Yuchen Li, and Kian-Lee Tan. 2019. GPU-based Graph Traversal on Compressed Graphs. In Proceedings of the 2019 International Conference on Management of Data (Amsterdam, Netherlands) (SIGMOD ’19). Association for Computing Machinery, New York, NY, USA, 775–792. doi:10.1145/3299869.3319871 [54] Ahmedur Rahman Shovon, Thomas Gilray, Kristopher Micinski, and Sidharth Kumar. 2023. Towards Iterative Relational Algebra on the GPU. In 2023 USENIX Annual Technical Conference (USENIX ATC 23). USENIX Association, Boston, MA, 1009–1016. https://www.usenix.org/conference/atc23/presentation/shovon [55] Yihao Sun, Ahmedur Rahman Shovon, Thomas Gilray, Sidharth Kumar, and Kristopher Micinski. 2025. Optimizing Datalog for the GPU. In Proceedings of the 30th ACM International Conference on Architectural Support for Programming Languages and Operating Systems, Volume 1 (Rotterdam, Netherlands) (ASPLOS ’25). Association for Computing Machinery, New York, NY, USA, 762–776. doi:10. 1145/3669940.3707274 [56] Yasuo Tabei, Hiroto Saigo, Yoshihiro Yamanishi, and Simon J. Puglisi. 2016. Scalable Partial Least Squares Regression on Grammar-Compressed Data Matrices. In Proceedings of the 22nd ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (San Francisco, California, USA) (KDD

13

Related documents

• Computation-Friendly Compression of Matrices: Speeding up Linear and Boolean Algebra#1458832
• UCCL-Zip: Lossless Compression Supercharged GPU Communication#124033
• Simultaneous Matrix-Free HSS Compression and Factorization with Adaptive Ranks#1341795
• Mathematical foundations of the GraphBLAS#1366569
• Resource allocation in sparse graphs#306912