arXiv CSopen access

Right Multiplication on Grammar-Compressed Matrices: A Streaming, Memory-Bounded GPU Engine

Right Multiplication on Grammar-Compressed Matrices: A Streaming, Memory-Bounded GPU Engine

arXiv:2607.24971v1 [cs.MS] 27 Jul 2026

Francesco Tosoni∗

Gabriele Mencagli†

1 Introduction Following [21, 22], research into computation-friendly matrix formats [1, 4, 5, 7, 10, 19, 23, 25, 33, 42] enables linear algebra directly on compressed representations. Similar principles drive compressed-text engines on GPUs and FPGAs [46, 47] and indexing repetitive collections [37, 38], where space scales with repetitiveness rather than raw size, as for the r measure of the r-index [26, 27] or run-length-BWT machinery [16] adapted for genotype data. The mm-repair scheme [23] encodes the structure of a matrix M ∈ Rn×m as a string, grammar-compressed via RePair [31] into (C, R, V ) (Equations (2.1), (2.3), and (2.4)). This supports right products y = M x and left products xT = y T M in O(|C| + |R|) time and O(|R|) space [42]. Prior works [23, 42] parallelize this by partitioning M into row blocks, compressing them separately, and assigning each to a thread. We propose an orthogonal, GPU-specific parallelization within a single grammar. We focus on memory-constrained engineering: maximizing in-place computation near physical memory limits. Since the grammar is a DAG of out-degree 2 (sections 2 and 3), the right product is a bottom-up, conflict-free gather. The dual left product is a scatter with write contention on highly shared nodes: an obstacle on GPUs. This up- and down-sweep pattern recalls Blelloch’s parallel prefix sum [11]: our right-product gather mirrors his reduction, while his distribution mirrors the leftproduct scatter. Yet, though balanced trees allow local rewrites, our DAG’s high-in-degree nodes force accumulation under contention. We thus restrict this paper to the right product: it is essential for power iteration (e.g., PageRank) [23, 42, 25] and avoids scatter bottlenecks. Streaming on a properly layered grammar. The sweep is embarrassingly parallel within each level (Lemma 3.1), but the target regime needs more than node independence: a properly layered grammar, where every nonterminal child sits one level below its parents (Definition 3.2). Then each level reads only the level below and writes the next, so traversing bottom-up every frontier is freed once consumed and two buffers alternate ∗ Sant’Anna School of Advanced Studies, L’EMbeDS, p.zza Martiri della Libertà 33, 56127 Pisa PI, Italy read- and write-only (Lemma 3.3). We retain RePair’s compression capabilities and impose proper layering via ([email protected]). † Department of Computer Science, University of Pisa, L.go B. pass-through completion (section 4), inserting identity Abstract. Grammar-compressed matrices (the mm-repair family [23]) store a matrix’s non-zero structure as a RePair straight-line program (SLP) [31], supporting matrix–vector products in time and space proportional to the compressed size. We target the regime where this is decisive on a GPU: when the uncompressed matrix exceeds device memory, so footprint (not floating-point throughput) is the binding constraint. Our SLP is a directed acyclic graph (DAG) of out-degree 2, and the right product y = M x is a single bottom-up sweep (leaves → roots): a conflictfree gather. The dual left product is a scatter with write contention, which we deliberately avoid. We make the grammar properly layered (every nonterminal child one level below its parent) via pass-through completion, which inserts identity nodes to carry values upward until consumed. This yields a streaming evaluation in which each level reads only the level below and writes the next, so the live set fits in two alternating read-only/write-only buffers instead of scaling with the whole grammar; the per-level width equals the live set. On genotype matrices, where a polygenic score is exactly the right product y = Gβ [14], a CUDA implementation shows a clear space advantage: a device footprint 4 to 8 times smaller than a materialized cuSPARSE [39] CSR baseline, single-vector times within a small factor of cuSPARSE, and consistently lower energy. Because the sweep needs only an associative combine, the same engine and schedule evaluate any monoid homomorphism over the grammar by swapping a small leaf/combine/emit policy; the same reachability sweep then scales to the billion-edge Software Heritage graph (261 TB dense and unmaterializable, 21× smaller serialized than CSR), where the memory argument holds. We frame this as an algorithm-engineering case study: structural metrics (depth, live-set width, completion cost) are measured, architecture-independent grammar properties, whereas time and energy are profiled on a single board.

Pontecorvo 3, 56127 Pisa PI, Italy ([email protected]).

nodes that propagate values upward until they are consumed; the resulting per-level width precisely matches the live set (i.e., the number of active nonterminals). Scope and non-goals. We pursue a “compress-once, multiply-many-times” scenario in which the input vectors are not known a priori. The GPU serves as our sole target architecture; a host-constructed grammar feeds a level-synchronous GPU engine, which subsequently collects all throughput metrics. We evaluate performance against the 16-thread CPU mm-repair baseline from [23] and the NVIDIA cuSPARSE CSR kernel [39] as references. The depth, live-set width, and structural sizing metrics are architecture-independent grammar properties computed on the host during construction; these a priori grammar features enable prediction of GPU execution costs. We exclude left multiplication and FPGA-based deployments, deferring them to future work (section 8). We delineate in sections 6 and 7 board-dependent versus invariant metrics, and validate on two domains stressing different axes: genotype matrices, whose dense form is prohibitive at biobank scale so y = Gβ wins on space at time parity, and repetitive graphs, where the same monoid sweep reaches a billionedge software graph. 2 Background: grammar-compressed matrices CSRV. Following the methodology in [23], let V [1, k] store the distinct non-zero values of M . Scanning M row by row, each non-zero entry M [i][j] = V [ℓ] is emitted as a pair ⟨ℓ, j⟩. A delimiter $ terminates each row, thereby producing a string S composed of value-column pairs and $. Equations (2.1) and (2.2) give the small matrix used as our running example (an 8 × 7 matrix with |V | = 4 distinct non-zero values and three distinct rows, repeated to expose shared structure) and its CSRV string S; the value array is V = [ 1.2, 2.3, 3.4, 4.5 ]. In (2.2) we write out only the three distinct rows (1, 2, 6); the · · · mark where the repeated row blocks recur: rows 3, 5 repeat row 2, rows 4, 8 repeat row 1, and row 7 repeats row 6. The pair ⟨ℓ, j⟩ encodes value V [ℓ] in column j; for example, ⟨2, 1⟩ is V [2] = 2.3 in column 1, while the same value in column 6 is ⟨2, 6⟩. Only equal values in the same column share a pair. A row delimiter $ closes each row. 

(2.1)

 2.3 0 1.2 1.2 1.2 2.3 1.2 1.2 4.5 3.4 1.2 1.2 2.3 4.5   1.2 4.5 3.4 1.2 1.2 2.3 4.5   2.3 0 1.2 1.2 1.2 2.3 1.2   M =  1.2 4.5 3.4 1.2 1.2 2.3 4.5  0 4.5 3.4 2.3 2.3 0 4.5    0 4.5 3.4 2.3 2.3 0 4.5 2.3 0 1.2 1.2 1.2 2.3 1.2

S = ⟨2, 1⟩⟨1, 3⟩⟨1, 4⟩⟨1, 5⟩⟨2, 6⟩⟨1, 7⟩ $ (2.2)

⟨1, 1⟩⟨4, 2⟩⟨3, 3⟩⟨1, 4⟩⟨1, 5⟩⟨2, 6⟩⟨4, 7⟩ $ · · · ⟨4, 2⟩⟨3, 3⟩⟨2, 4⟩⟨2, 5⟩⟨4, 7⟩ $ · · ·

Grammar. RePair [31] compresses S by repeatedly replacing the most frequent pair of adjacent symbols with a new nonterminal, modified to ensure it never pairs elements across a $ boundary [23]. This process yields a set of production rules Ni → Ai Bi (where each Ai , Bi is either a terminal ⟨ℓ, j⟩ or a preceding nonterminal N<i ) and a final sequence C whose constituent symbols expand directly to the matrix rows. The indexing follows a strict topological ordering: i < j whenever Ni appears on the right-hand side of the rule for Nj . Executing RePair on the string S of (2.2) yields:  N2 → N1 ⟨2, 6⟩      N4 → ⟨1, 1⟩N3     N6 → ⟨2, 1⟩N5  N8 → N4 N7     N10 → ⟨2, 4⟩⟨2, 5⟩     N12 → N11 ⟨4, 7⟩

(2.3)

 N1 → ⟨1, 4⟩⟨1, 5⟩      N3 → ⟨4, 2⟩⟨3, 3⟩     N → ⟨1, 3⟩N 5 2 R=  N7 → N2 ⟨4, 7⟩      N9 → N6 ⟨1, 7⟩    N11 → N3 N10

(2.4)

C = N9 $ N8 $ N8 $ N9 $ N8 $ N12 $ N12 $ N9 $.

In this specific instance, each row within C reduces to a single nonterminal symbol; thus, C tracks the eight distinct row roots drawn from the set {N9 , N8 , N12 }. In practice, C can be longer and may incorporate bare terminals. Indeed, RePair stops when no more pairs of consecutive symbols appear more than once [23]. The deepest root, N9 , encapsulates a chain of five nested rules (N9 → N6 → N5 → N2 → N1 ), indicating that the grammar spans five logical levels. The right-product sweep. Let evalx (⟨ℓ, j⟩) = V [ℓ] · x[j] denote the scalar contribution of an individual matrix entry to the output product. The underlying grammar transforms the multiplication y = M x into a single bottom-up evaluation pass governed by three key properties established in [23]. Lemma 2.1 (Additivity of evalx [23]). For any grammar rule Ni → Ai Bi , the evaluation distributes additively: evalx (Ni ) = evalx (Ai ) + evalx (Bi ). Lemma 2.2 (Rows as roots [23]). Given a compressed sequence C = Ni1 $ · · · Nin $ and a product y = M x, the output vector matches the root evaluations: y[r] = evalx (Nir ) for every row index r = 1, . . . , n. Theorem 2.3 (Compressed-time right product [23]). Given the grammar components (C, R, V ) representing a matrix M ∈ Rn×m and an input vector

x ∈ Rm , the product y = M x can be evaluated in O(|C| + |R|) time using O(|R|) words of auxiliary storage.

a value two frontiers below it. Pass-through completion (subsection 4.2, Figure 3.1) splits exactly such edges so the streaming sweep of Lemma 3.3 applies.

In practice, the CPU implementation of [23] populates an auxiliary array W [1, |R|] such that W [i] = evalx (Ni ) via a single forward traversal. Because the topological numbering ensures child nodes precede their respective parents, Lemma 2.1 resolves sequentially as W [i] = cv(Ai ) + cv(Bi ), where cv(Nj ) = W [j] and cv(⟨ℓ, j⟩) = V [ℓ]x[j]. The final answers are then extracted via y[r] = W [rootr ], as in Lemma 2.2. Preserving the explicit $ boundaries enforces a clean, row-separated representation that underpins the rightproduct execution; all our structural GPU-specific modifications respect this property. Theorem 2.3 provides the operational bound maintained by our engine, and the streaming framework covered in subsection 4.3 represents its level-synchronous, parallel GPU realization.

Lemma 3.1 (Levels are antichains). If lvl(u) = lvl(v) for distinct nodes u ̸= v, then G contains no directed edge connecting u and v.

3 The grammar as a DAG, and its levels We interpret the rules in (2.3) as a directed acyclic graph G: each nonterminal maps to an internal node, each distinct terminal forms a leaf node, and directed edges span from every parent Ni to its two children. Consequently, nonterminals maintain an out-degree of 2, terminals remain childless leaves, and the set of roots corresponds to the symbols active in C. Crucially, G forms a DAG rather than a tree, as reused subexpressions exhibit an in-degree > 1. For instance, N2 drives both N5 and N7 , N3 feeds both N4 and N11 , and individual roots like N9 , N8 , N12 serve multiple rows within C. From the perspective of the right product, this multi-parent sharing is benign because it represents pure read-sharing. Levels. We assign every node an integer level corresponding to its maximal path distance from a leaf:

Proof. Every valid edge in G terminates at a child node possessing a strictly smaller level according to (3.1); hence, the source and destination endpoints of an edge can never occupy the same level. Lemma 3.1 serves as the core parallelization driver for our architecture: all nodes residing within a given level can be evaluated simultaneously. The total number of levels, meaning the overall grammar depth L = maxi lvl(Ni ), defines the critical path of the computation and determines the number of sequential GPU kernel launches. Yet, node independence alone does not fully satisfy our design objectives; we require a more rigid structural configuration. Definition 3.2 (Proper layering). A grammar is properly layered if every nonterminal child node is positioned exactly one level beneath its parent node: for any rule N → AB where lvl(N ) = k, any child that is a nonterminal must satisfy lvl = k − 1. Terminal children face no such constraints; they interface directly with the persistent inputs x and V , which remain globally resident in device memory and are never deallocated.

Lemma 3.3 (Liveness / double buffering). When executing over a properly layered grammar, the rightproduct sweep requires only two active value buffers in working memory. While evaluating level k, the kernel reads exclusively from level k − 1 and writes to level k. Upon completing level k, the frontier at level k − 1 has no remaining consumers and can be safely deallocated; any root node finalized at level k is immediately emitted to the output vector y. Consequently, the transient lvl(terminal) = 0,
(3.1) memory footprint is bounded by O(maxk wk ), where wk lvl(Ni → Ai Bi ) = 1 + max lvl(Ai ), lvl(Bi ) . is the width of level k, rather than scaling as O(|R|), and Applying this definition to the rules in (2.3) gener- the execution alternates between two dedicated read-only ates five distinct strata: {N1 , N3 , N10 }, {N2 , N4 , N11 }, and write-only buffers. {N5 , N7 , N12 }, {N6 , N8 }, and {N9 }, as illustrated in Proof. By Definition 3.2, all parent nodes consumFigure 3.1. The left spine N9 → N6 → N5 → N2 → N1 reing a nonterminal from level k − 1 must reside exactly alizes the five levels. Nodes are shared (N2 feeds N5 and at level k. Once level k is fully processed, all potential N7 ; N3 feeds N4 and N11 ; the terminal ⟨4, 7⟩ feeds N7 consumers of the level-(k − 1) buffer have finished exand N12 ). Many edges skip a level by reaching a termiecution, allowing it to be cleared. Terminal nodes are nal (e.g., N9 , N6 , N5 each drop to a level-0 leaf): these read from persistent global memory and do not depend are harmless, as terminals read the persistent input and on transient buffers. are never freed (Definition 3.2). The critical edges are We use pass-through completion (section 4) to the nonterminal level-skipping edges, such as N8 → N4 (red) from level 4 to level 2. This skips a level and vio- ensure proper layering, enabling streaming, doubleMinimizing L reduces kernel lates proper layering, since evaluating N8 would require buffered execution.

level 5

N9

level 4

N6

level 3

N5

level 2

N2

N7

P

N12

N4

N11

N1

level 1

level 0

N8

⟨1, 4⟩

⟨1, 5⟩

N3

⟨2, 6⟩

⟨1, 3⟩

⟨2, 1⟩

⟨1, 7⟩

⟨1, 1⟩

⟨4, 2⟩

N10

⟨3, 3⟩

⟨4, 7⟩

⟨2, 4⟩

⟨2, 5⟩

Figure 3.1: The grammar DAG of (2.3) (L = 5 levels, background bands; (3.1)). Circles represent non-terminals, while boxes represent terminal leaves. Dashed gray edges are harmless terminal level-skips (Definition 3.2). The lone nonterminal level-skipping edge N8 → N4 (red) violates proper layering. Passthrough completion (subsection 4.2) splits it through the synthesized node P (orange path N8 → P → N4 ), making every nonterminal edge level-adjacent. launches, while bounding the maximum level width wk ensures live sets fit within device memory. This approach mirrors sparse dynamic programming on smallwidth DAGs [32], where costs scale with the width of a minimum path cover: here, the per-level antichains of Lemma 3.1. 4 A proper-layered streaming engine Our central thesis is that we can fully preserve the compression advantages of the mm-repair format while enabling the streaming, double-buffered execution model of Lemma 3.3. We achieve this by transforming the raw RePair output into a properly layered grammar during a post-processing phase. The overhead introduced consists of the synthesized pass-through nodes; yet, our experiments demonstrate that this cost remains low for highly compressible matrices and scales up only on dense, poorly compressible datasets (subsection 4.4). 4.1 Step 1: RePair (unchanged) We execute RePair on the input sequence S exactly as described in [23], processing the matrix row by row and strictly enforcing the $ boundaries. We treat it here as an immutable black box. Its primary drawback for highthroughput parallelization is structural: its greedy, global selection of the most frequent digram yields deep graph topologies containing level-skipping edges (Figure 3.1). We report the overall depth L and the exact distribution of edge spans for our datasets in subsection 4.4.

ering guarantees, we explicitly add nodes to the grammar graph. For every nonterminal directed edge N → . . . M . . . where the level differential d = lvl(N )−lvl(M ) is strictly greater than 1, we intercept the edge by inserting a sequential chain of d − 1 pass-through nodes denoted M = M0 , M1 , . . . , Md−1 . These artificial nodes implement identity rules of the form Mi → Mi−1 , ensuring that evalx (Mi ) = evalx (Mi−1 ). The parent node N is then modified to read from Md−1 at level lvl(N ) − 1; terminal edges are left unaltered (Definition 3.2). Completion runs independently per row, preserving the row roots in C and the alignment y[r] = W [rootr ]; it is an offline, amortized preprocessing cost, not an online penalty (subsection 4.4). Figure 3.1 shows the result: the single level-skipping edge N8 → N4 (red) is split by one pass-through P at level 3 (P → N4 ), so N8 reads P one level below and every nonterminal edge becomes level-adjacent (Lemma 3.3). The passthrough is not overhead but bookkeeping: it is N4 held live across level 3 for its higher consumer N8 , so each band’s width is exactly the live set. Figure 4.2 replays this as a numeric, round-by-round computation. Such insertions are not wasteful: a pass-through crossing level k is a value alive but unconsumed there, so the completed width wk tracks the live set and the buffers hold exactly the values in flight (on GPUs, forwarded frontier slots; on an FPGA pipeline they would correspond to pipeline registers). Worst-case completion adds O(|R| L) nodes, but the actual inflation follows the edge-span distribution, which we measure empirically (subsection 4.4). Depth still matters, but is not enough. Grammar height governs both the sequential launch count L and the length of individual pass-through chains. Thus, pre-balancing to a logarithmic height bound of O(log n) [28, 35] is an attractive complementary optimization. We do not apply it in this work (our pipeline consumes the raw RePair grammar unchanged), and we stress that it would not replace completion. Structural balancing bounds only the depth, not the active live set, meaning a balanced grammar can still host nonterminal edges that skip multiple levels. Because only explicit completion guarantees compliance with the streaming model, balancing is orthogonal to our contribution; we leave its integration to future work (section 8).

4.3 Step 3: the streaming, double-buffered sweep During the offline phase, we evaluate (3.1) on the completed grammar. Nodes are laid out level by level, with rules concatenated into a flat array and partitioned by offsets so each level k occupies a contiguous block. Node children are encoded as 4.2 Step 2: pass-through completion (layer- absolute offsets relative to the preceding frontier array, ing) Because standard RePair offers no structural lay- with tags identifying terminal, nonterminal, or pass-

kernel round(level k, prev, cur): # one thread t per node of level k r = rules[k][t] cur[t] = r.coeff * cv(r.left) + cv(r.right) # one fused MAD, no branch kernel emit(level k, cur): # one thread per C-occurrence at level k atomicAdd(&y[row[e]], cur[pos[e]]) # emit on the spot, then free # host: zero y; for k = 1..L { round(k,prev,cur); emit(k,cur); swap(prev,cur) }

tion keeps every lookup on the adjacent level rather than scattering across the DAG. The sequential launch count is bounded by L, which height-balancing would cap logarithmically. Per-round width tapers toward the roots; occupancy peaks where most of the reduction happens. k=1

k=2

k=3

read

read

read

buf α

WRITE

READ

WRITE

buf β

idle

WRITE

READ

y: output

emit

emit

emit

WRITE (cur)

persistent

2

T : input

Listing 1: Simplified streaming round: read prev (level k − 1), write cur (level k); the host swaps the buffers after each launch. cv(c) = prev [c] if c ≥ 0, else T [−c−1]. A node of C is emitted into y the moment its level is computed (atomic add), then freed.

through types. Maintaining nodes in RePair index order ensures that child reads remain monotone and largely contiguous across warps. Execution of y = M x follows Listing 1 using alternating buffers, prev and cur . For each level k = 1, . . . , L, a single kernel launch reads level k − 1 and writes level k, with the host swapping buffer pointers to alternate modes (Lemma 3.3). The maximum transient memory footprint is O(maxk wk ). To prevent warp divergence, we handle all node variants using a branch-free, fused multiply–add sequence. We allocate a dedicated zero terminal at T [0] = 0 and encode all child references as signed array indices: a non-negative index c indicates a lookup in the previous frontier buffer prev [c], whereas a negative index c routes the lookup to the persistent terminal array T [−c − 1], where T [p] = V [ℓp ] x[jp ]. Every node can then be evaluated using the unified expression cur = coeff · cv(left) + cv(right), where the tuple parameters are assigned as follows: (coeff, right) = (1, B) for standard binary rules N → AB; (1, zero) for passthrough operations N → A; and (t, zero) for run-length rules N → At , giving evalx (N ) = t · evalx (B) from a single child read. The last case never arises for us: within a row, the column index is distinct at every position, so no pair (hence no digram) can repeat, and RePair emits only binary productions; the engine nonetheless supports it for run-length-enriched grammar constructions too (section 8). Figure 4.1 illustrates the execution on three consecutive levels of the running example (L = 5; Figure 4.2 runs it in full). The round kernel is a conflict-free gather: threads write distinct slots of cur , and concurrent reads of shared children in prev are harmless. A root of C is emitted into y by atomic addition the instant its level is finalized, then discarded; this keeps roots from propagating upward, so wk stays minimal rather than inflating to Θ(rows). Pass-through comple-

1

α reused

3

4

READ (prev)

Figure 4.1: The double-buffered, level-synchronous sweep (Listing 1), shown for the first three levels (k = 1 double 1, 2, 3). Buffers α, β swap roles every level. ○ 2 persistent terminals; ○ 3 branch-free fused buffering; ○ 4 emit-on-the-spot. MAD; ○ 4.4 Experimental results We evaluate the engine on large genotype datasets (GB10 platform, section 6): computing a polygenic risk score is the right product y = Gβ [14]. Biobank-scale panels (hundreds of thousands of samples, millions of variants) exceed physical memory, which is why genomics relies on specialized run-length/BWT-based haplotype formats [20, 16]. Our matrices G come from human chromosomes 22, 21, and 20 (from the 1000 Genomes Project [15]): individuals as rows, SNPs as columns in genomic order, each cell in {0, 1, 2}; adjacent columns are in heavy linkage disequilibrium (LD), with repetitiveness that RePair exploits. We use two scales (a 105 -variant subset and the full per-chromosome sequence), plus synthetic panels from the coalescent with recombination via msprime [30, 9]. The recombination-to-mutation ratio of these synthetic panels tunes LD block structure (reproducible via fixed seeds), with low ratios creating highly compressible, long-LD blocks and high ratios the inverse. We avoid full biobank cohorts (section 7); our largest synthetic panel (crossover_synth, 10,000 individuals, 700,000 variants, 1.00 G non-zeros) serves as a scale benchmark. Table 4.1 gives the grammars’ structural parameters (|R|, depth L, max frontier width w∗ , pass-through inflation +pt). Table 4.2 gives the average time per right product (over 100 vectors) for our GPU streaming engine (GPU); a parallel OpenMP CPU sweep (20 threads) and a sequential CPU one; the original

read

k=1

read

k=2

read

k=3

⟨1,4⟩

⟨1,5⟩

⟨2,6⟩

⟨4,2⟩

⟨3,3⟩

⟨2,4⟩

⟨2,5⟩

2.4

1.2

2.3

4.5

3.4

4.6

2.3

3.6

7.9

6.9

N1

N3

N10

3.6

7.9

6.9

N1

N3

N10

5.9

10.3

14.8

N2

N4

N11

5.9

10.3

14.8

N2

N4

N11

⟨2,6⟩

⟨1,1⟩

2.3

2.4

matrix

7.1

10.4

10.3

19.3

N5

N7

P (pt)

N12

⟨1,3⟩

⟨4,7⟩

1.2

4.5

⟨2,1⟩

read

7.1

10.4

10.3

N5

N7

P

Table 4.1: Structural figures, measured on the GB10 node. |R|: RePair non-terminals; L: grammar depth; w∗ : maximum live width after completion (the streaming-buffer size); +pt: pass-through nodes added by completion. K = 103 , M = 106 .

4.6

y = Mx 12.9

= N9

20.7

= N8

. . . 19.3

k=4

11.7

20.7

N6

. . .

N8 ⟨1,7⟩

read

11.7 N6

k=5

12.9

= N12

1.2

persistent T read (prev)

emitted root

write (cur)

pass-through

N9

Figure 4.2: The DAG grammar execution for input x = [ 2, 1, 1, 2, 1, 1, 1 ]⊤ and V = [ 1.2, 2.3, 3.4, 4.5 ]. Terminals evaluate as ⟨ℓ, j⟩ 7→ V [ℓ] · x[j]. Distinct row roots N12 , N8 , N9 evaluate to 19.3, 20.7, 12.9. Ellipses (. . .) elide repeated rows.

mm-repair (mmr) at 1 and 16 threads; and vendor cuSPARSE SpMV, with ×seq representing the GPU speed-up over sequential mm-repair. Table 4.3 adds device footprint and energy versus cuSPARSE. Depth and completion. The critical depth L remains low across all test configurations: 20 to 30 tiers for the grammars of Table 4.1, and L = 32 for the larger crossover_synth. Crucially, the maximum streaming frontier width w∗ stays well below the base rule count |R| (about one third on the real chromosomes, up to one half on the synthetics), enabling the dual buffers to fit within device memory. Pass-through nodes are added in proportion to the underlying edge-span profiles (subsection 4.2). As the genotype alphabet is discrete and small ({0, 1, 2}), terminal evaluation is cheap and does not bound performance, leaving traversal and emission atomics as the likely dominant kernel cost. Utilizing CUDA Unified Memory (cudaMallocManaged) prevents out-of-memory errors on the host during graph construction. Time vs. mm-repair and cuSPARSE. Our GPU streaming engine outperforms the single-threaded CPU mm-repair reference across all test cases, achieving speedups ranging from 5.4× to 15.5×. When compared

|R|

L

w∗

Real Genotypes (1000 Genomes) Chr22 437.6K 20 146.0K Chr22 full 4.29M 24 1.41M Chr21 404.6K 24 133.1K Chr21 full 4.30M 26 1.39M Chr20 438.3K 22 144.7K Chr20 full 6.65M 24 2.18M Synthetic Genotypes (Haplotypes) synth_small 625.3K 25 225.0K synth_large 3.94M 29 1.39M synth_ld_high 1.57M 30 498.4K synth_ld_low 3.47M 21 1.54M synth_ind_large 1.40M 30 489.8K

+pt 668.1K 7.02M 568.6K 7.29M 718.9K 11.12M 1.23M 8.29M 3.55M 4.85M 3.08M

against the heavily optimized, vendor cuSPARSE CSR kernel, our single-vector engine maintains competitive parity on the real chromosomes, from 1.41× faster to 1.20× slower than cuSPARSE. Further, our engine outperforms cuSPARSE on synthetic configurations characterized by high linkage disequilibrium or expanded sample sizes (achieving up to a 4.2× performance improvement on synth_ind_large), where the underlying compression factor is at its maximum. Space and energy vs. cuSPARSE. Table 4.3 reports analytic device footprint, time, and energy of the engine versus cuSPARSE CSR SpMV. The engine is consistently smaller and more energy-efficient, especially on highly repetitive or high-LD matrices. We observe a 4.0× to 4.7× footprint reduction across the real chromosomes, which correlates with lower energy usage (up to 35% on the full chromosomes, and roughly half on the 105 -variant subsets). This efficiency stems from arithmetic operations running directly on the compressed grammar, without materializing uncompressed vectors. On synthetics, the LD level sets the scaling: high LD (recombination rate 10−9 ) gives a 7.3× smaller space and beats cuSPARSE by 3.6× in time and 4.6× in energy. Conversely, low LD (recombination rate 10−7 ) drops to 3.5× space and runs at time parity. Widening the cohort (synth_ind_large, 10,000 individuals) reaches 7.9× smaller space and 4.7× less energy, at the speed-up already reported above. For biobank-scale problems, this memory advantage decides whether a matrix fits in high-speed GPU memory (Figure 4.3). Pushing to billion-nonzero scale. To probe the largest-scale regime, we generated a large synthetic matrix (crossover_synth, 10,000 × 700,000 genotypes,

Table 4.2: Genotype average time per right matrix–vector product (ms/vector) on the GB10 node. matrix

GPU

Real Genotypes (1000 Genomes) Chr22 (2504 × 0.10M) 0.50 Chr22 (2504 × 1.06M) 6.69 Chr21 (2504 × 0.10M) 0.50 Chr21 (2504 × 1.05M) 6.38 Chr20 (2504 × 0.10M) 0.48 Chr20 (2504 × 1.74M) 11.73 Synthetic Genotypes (Haplotypes) synth_small (2K × 50K) 0.51 synth_large (5K × 200K) 2.97 synth_ld_high (5K × 100K) 0.98 synth_ld_low (5K × 100K) 3.78 synth_ind_large (10K × 50K) 0.83

OpenMP

seq.

mmr (seq)

mmr (16th)

cuSPARSE

×seq

11.01 32.03 14.07 30.05 18.48 40.67

9.37 93.22 10.08 87.19 8.92 139.31

2.80 68.60 2.70 66.60 3.00 139.70

3.00 27.80 0.60 28.00 2.70 42.90

0.64 6.21 0.71 6.81 0.68 9.76

5.6 10.3 5.4 10.4 6.2 11.9

12.22 31.33 23.46 23.66 22.83

7.95 50.57 18.12 66.53 16.64

2.80 46.00 7.90 47.30 9.20

2.50 26.00 11.10 21.10 6.30

0.82 7.41 3.58 3.78 3.49

5.4 15.5 8.1 12.5 11.1

Batched evaluation (SpMM). Batching B righthand vectors (Y = M X, the kernel of block power/Krylov iteration [29, §7.3, §10.1, §10.3.6]) amortizes 103 grammar traversal. With both sides given their best configuration (the engine’s best B, and cuSPARSE’s best of ALG_DEFAULT/CSR_ALG2/CSR_ALG3), cuSPARSE 102 leads by 1.7× to 10.9×. This gap is narrowest on the full real chromosomes (1.7 to 2.4×) and widest on the least compressible synthetics (7.2 to 10.9×). Unlike 107 108 109 the single-vector case, the footprint advantage does not Number of Non-Zeros (NNZ) persist at large B, where the engine’s B-wide terminal array and streaming buffers overtake cuSPARSE’s Figure 4.3: Analytic device footprint, grammar en- B-invariant CSR (section 7). Full per-dataset batched gine vs. cuSPARSE CSR, on log-log scale; both se- times and the algorithm-sweep that justifies CSR_ALG3 ries are analytic, as in Table 4.3. Even at the as the vendor baseline are in section A. largest scale (crossover_synth, 1.00 G non-zeros) cuS5 Beyond (+, ×): a monoid-homomorphism PARSE’s CSR (8.03 GB) still fits, but the Grammar engine Our engine does not rely on full-ring axioms, Engine’s 1.01 GB is 8.0× smaller. so the same schedule that computes the genotype product also evaluates graph reachability, with no change 1.00 billion non-zeros). In standard CSR representa- to the construction pipeline in section 4. Correctness tion, this matrix requires 8.03 GB of GPU memory; on over the shared DAG needs only that the combine ⊕ the GB10’s unified 119 GiB pool, cuSPARSE runs, but be associative (a monoid): then a reused nonterminal’s at 48.38 ms per vector and 1976 mJ/vec: the highest value is invariant to how its sub-expansion is parenthesized, so the intra-level parallelism of Lemma 3.1 and single-vector cost of any genotype matrix we test. In contrast, our Grammar Engine processes the ma- the liveness of Lemma 3.3 still hold. Commutativity is trix entirely in its compressed representation, requiring not required, since RePair preserves symbol order and an analytic device footprint of only 1.01 GB (measured never pairs across row boundaries. The leaf map (nonpeak: 0.98 GB), which is 8.0× smaller than CSR. Once zero terminals) is the multiplicative ⊗; an unmapped constructed, the engine evaluated the right product in entry takes the ⊕-identity (the semiring’s zero: 0 for 20.09 ms per vector (14.5× faster than the sequential (+, ×) and Boolean, +∞ for tropical). CPU reference and 2.41× faster than cuSPARSE), with Our engine evaluates any arbitrary monan energy consumption of 628 mJ/vec (3.15× less; see oid homomorphism over the grammar without Users swap low-level decrossover point in Figure 4.3). At this scale, the ad- schedule adjustments. vantage is no longer confined to space but also to time vice primitives via a unified configuration struct Two and energy; at a true biobank scale (where the uncom- (SEMIRING={plustimes,boolean,tropical}). pressed matrix exceeds device memory), it is what keeps semirings are of immediate interest: the Boolean the matrix resident at all. semiring (⊕ = or, ⊗ = and), where one matrix–vector Memory Footprint (MB)

104

cuSPARSE CSR Grammar Engine (Ours)

Table 4.3: Genotype right product. Single-vector engine vs. cuSPARSE CSR SpMV in unified memory on the GB10 node. nnz: non-zeros; M = 106 ; K = 103 . matrix (rows × cols)

nnz

Real Genotypes (1000 Genomes) Chr22 (2504 × 0.10M) 12.58M Chr22 (2504 × 1.06M) 127.58M Chr21 (2504 × 0.10M) 14.04M Chr21 (2504 × 1.05M) 140.20M Chr20 (2504 × 0.10M) 13.52M Chr20 (2504 × 1.74M) 201.04M Synthetic Genotypes (Haplotypes) synth_small (2K × 50K) 16.62M synth_large (5K × 200K) 152.44M synth_ld_high (5K × 100K) 73.47M synth_ld_low (5K × 100K) 77.01M synth_ind_large (10K × 50K) 71.60M crossover_synth (10K × 700K) 1002.97M

pass expands a single BFS frontier; and the Tropical semiring (⊕ = min, ⊗ = +), where a pass performs one Bellman–Ford edge relaxation. We measure a single pass: the natural unit of work; iterating to the BFS/SSSP fixpoint requires a square, endomorphic node remap, which we leave to future work (section 7). This recasts adjacency compression as an enabler for matrixalgebraic graph queries, such as WebGraph [12, 24], Zuckerli’s [43] intervalisation for adjacency–vector products [25, 42], k 2 -tree representations [13][36, §9.2.1], and Boolean-matrix algebras for regular path queries over compressed structures [7, 6, 8, 5]. This generalization reinforces our primary thesis: the same memory optimization that justifies grammar compression for genotypes (storing and operating on a highly compressed representation when uncompressed data structures exceed device limits) applies to any highly-repetitive graph topology. To demonstrate this capability, we evaluate the Boolean and Tropical semirings on five adjacency relation matrices built from Wikidata (the same Zenodo collection [8]), chosen to span scales, edge densities, and topologies: wd_sports_team (P54, member of sports team, 332,121 × 29,854), wd_cast_member (P161, cast member, 173,977 × 144,095), wd_citizenship (P27, country of citizenship, 2,874,250×2,556), wd_occupation (P106, occupation, 3,459,933 × 10,610), and wd_subclass_of (P279, subclass of, 1,487,709 × 73,417). The first three are bipartite, entity-to-entity relations with varying widths and degrees of repetition. The last two add the largest edge set and an ontologyhierarchy (rather than bipartite) topology. Because Wikidata subjects [41, 44, 45] number in the millions, these matrices are far too large to materialize densely;

space (MB)

time (ms)

energy (mJ)

eng.

cuS.

eng.

cuS.

eng.

cuS.

25.2 246.1 24.1 246.6 25.2 381.6

101.1 1025.4 112.8 1126.4 108.6 1616.1

0.50 6.70 0.50 6.40 0.47 11.72

0.64 6.21 0.71 6.81 0.68 9.76

14 199 14 183 14 328

26 257 31 280 28 439

31.6 201.4 81.0 173.8 72.4 1005.6

133.2 1221.0 588.5 616.8 573.4 8030.5

0.49 2.99 0.99 3.74 0.85 20.09

0.82 7.41 3.58 3.78 3.49 48.38

15 96 32 123 29 628

31 279 146 144 135 1976

we thus build them directly from the sparse per-source edge list into the CSRV/RePair grammar without ever forming the uncompressed matrix. Dataset characteristics and grammar structure (averaged over 50 iterations, CPU references verified bit-for-bit) are centralized in Table B.2, and runtime/space results in Table 5.2. We then push the space argument to the graph scale (Table 5.1) for the two largest Wikidata relations and for a billion-edge Software Heritage graph [40]. Compressibility tracks per-source repetition, as with the genotype matrices of subsection 4.4. RePair compresses wd_citizenship’s 3.06M edges to just 3,295 rules and the ontology hierarchy wd_subclass_of’s 2.02M edges to 6,825 rules. Both are relations in which each source shares its target structure with many others. wd_sports_team and wd_occupation compress moderately well (70,158 and 41,544 rules from 1.14M and 4.60M edges respectively), while wd_cast_member (where each film has a largely distinct cast) compresses least (51,677 rules from 1.03M edges, the only relation whose device footprint fails to beat CSR). Grammar depth ranges L = 5 to 8 and pass-through completion overhead ranges 11% to 21% across all five relations (Table B.2). None collapses to a flat, near-uncompressed structure that a degenerate, extremely narrow relation can produce, since all five expose genuine per-source repetition at this scale. Operating directly on these compressed DAG topologies, the engine’s analytic device footprint undercuts cuSPARSE CSR on four of the five relations (21% to 33% smaller on wd_sports_team, wd_citizenship, wd_occupation, and wd_subclass_of). Yet, it is larger on the poorly repetitive wd_cast_member (11.12 MB vs. 10.24 MB). This matches the failure mode seen on

badly-compressing genotype instances (subsection 4.4): once a relation lacks high repetitiveness, the per-rule bookkeeping overhead of the layered grammar outweighs the modest structural savings. Single-vector times remain within a small factor of cuSPARSE (Boolean 0.09 to 0.73 ms vs. 0.06 to 0.50 ms) while outrunning the sequential CPU reference by 17× to 34×. The batched B= 16 sweeps amortize the traversal further (e.g., Boolean wd_sports_team drops from 0.092 to 0.051 ms/vector; Table 5.2). These results confirm that the space advantage is real but conditional on genuine per-source repetitiveness: a measurable structural fact about each relation. Semiring-native baseline (GraphBLAS). For these semirings, we compare to SuiteSparse:GraphBLAS [17] running the same operation (lor_land/min_plus mxv) on the uncompressed matrix (20 threads, Table 5.2). Working on the compressed grammar, our GPU engine outperforms GraphBLAS on all five Boolean sweeps. For the Tropical semiring, it is 3× to 7× faster across all five relations (e.g., wd_occupation: 0.73 vs. 3.57 ms). A GPU-native cuGraph BFS/SSSP reference is discussed in section 7; it performs a full traversal rather than a single mat-vec, serving as an end-to-end comparison rather than a per-mat-vec baseline. Scaling to 10M–166M edges. To probe the space argument at scale, we build the two largest single Wikidata relations, bracketing the compressibility spectrum: the narrow categorical wd_country (10.1M edges) and the citation network wd_cites_work (166.7M edges); see Table 5.1. Their dense forms (2.2 GB, 10.8 TB) are unusable, making the sparse-to-grammar construction (section 5) essential to make them addressable. Two effects stand out. Compressibility tracks persource repetition: wd_country collapses to 2,002 rules (serializing 16.8× below CSR’s footprint), whereas wd_cites_work barely compresses (3.1×) due to its largely unique citation lists. More subtly, the deviceresident advantage is far smaller than the serialized one. The sweep still materializes one root per row plus the persistent terminal array, a Θ(rows) cost. Consequently, even the 16.8×-compressible wd_country is only ∼ 25% smaller than CSR in device bytes, and wd_cites_work is essentially tied. At graph scale, the grammar’s decisive advantage lies in its serialized/transfer size and in avoiding the dense matrix construction, rather than device-resident bytes, which are lower-bounded by the per-row root/terminal structure. This marks an honest boundary complementary to the genotype scale regime (subsection 4.4), where the engine’s device footprint stays several-fold below CSR even as both fit. A billion-edge software graph (Software Heritage). Our largest instance, and the only one from the soft-

ware domain, is drawn from the Software Heritage (SWH) archive [40], the global Merkle-DAG of publicly archived source code [34, §2.1][2, 18]. It is repetitive by construction. Source code evolves by small edits over many revisions and forks that reuse the same files and directory subtrees across the archive: exactly the redundancy RePair exploits. From the 2021-03-23-popular-3k-python export, we build a single boolean adjacency of 45.7M nodes and 1.22 billion edges (Table 5.1). Its dense form (≈ 261 TB) cannot be materialized, and even a boolean CSR needs ≈ 10 GB. RePair compresses it to a 30.3M-rule grammar serialized to 490 MB, representing a 21× reduction under CSR. Crucially, and unlike wd_cites_work, the advantage here also survives device-resident (2.79 GB, 3.7× under CSR), because SWH is simultaneously highly repetitive and dense enough (≈ 27 edges/row) to amortize the Θ(rows) root/terminal cost. Structurally, its grammar is deep and wide (L = 70, w∗ = 8.98M), presenting a distinct stress test compared to the shallow bipartite relations (L = 5 to 8). Correctness is verified bit-for-bit against the CPU reference on the full 45.7M-node, 1.22 billion-edge graph (Boolean semiring; the GPU output equals both the sequential and OpenMP CPU sweeps, with maximum absolute difference 0 over all 45.7M outputs). Our grammar serializes to ≈ 3.2 bits/edge; WebGraph [12, 24] stores it slightly smaller (2.485 bits/edge) but does not evaluate a semiring mat-vec on the compressed form on a GPU. 6 Experimental setup All metrics, including host-side grammar decomposition, structural analysis, and GPU execution throughput profiling, are collected on a device exclusively reserved for our experiments and equipped with an NVIDIA GB10 Grace–Blackwell superchip architecture. This hardware platform runs Ubuntu 24.04 LTS (Linux kernel 6.17) and features a 20-core Arm CPU (aarch64 architecture) tightly coupled to 119 GiB of unified coherent memory shared over a high-bandwidth NVLink-C2C interconnect. The onboard GPU hosts 48 streaming multiprocessors (SMs) targeting compute capability 12.1 (sm_121). The software stack comprises the CUDA 13.0 toolkit, g++ 13.3, and CMake 3.28. Energy is measured from the on-die power telemetry via NVML: we sample nvmlDeviceGetTotalEnergyConsumption on the Blackwell device before and after a sustained evaluation loop (repeated until at least 1.5 s of sustained execution) and divide the energy delta by the number of products, reporting GPU energy in mJ per vector; these figures are subject to the profiling variance noted in section 7. The most consequential architectural trait of the GB10 platform for our implementation is its fully unified and coherent physical memory layout: the host

Table 5.1: Largest graph relations that bracket the compressibility spectrum: two Wikidata relations and the billion-edge Software Heritage software graph. REANS: serialized grammar (the engine’s on-disk/transfer size); CSR / eng. dev.: analytic device footprint of cuSPARSE CSR and the engine (for these Boolean relations CSR’s serialized size and device footprint coincide up to the negligible x/y vectors, so the CSR column doubles as both baselines); Bool: single-vector Boolean time (engine / cuSPARSE), GB10 node. Dense forms (2.2 GB / 10.8 TB / 261 TB) are unusable. nnz

|R|

L

REANS (MB)

CSR (MB)

eng. dev. (MB)

Bool eng/cuS (ms)

10,089,284 166,682,725 1,218,488,928

2,002 9,919,927 30,292,590

5 15 70

9.6 459.0 490.1

161.2 1439.7 10301.0

120.8 1478.4 2792.9

1.63 / 1.41 42.71 / 34.56 52.79 / 65.42

relation wd_country (10.1M) wd_cites_work (166.7M) swh (1.22G)

Table 5.2: Graph right product under Boolean and Tropical semirings across five Wikidata relations on the GB10 node (single-vector time, slash-separated per method: Boolean CSR/eng/GB, Tropical eng/GB; batched engineonly at B=16; footprint columns analytic CSR/eng). CSR: uncompressed cuSPARSE (Tropical: no vendor kernel, N/A). GB: SuiteSparse:GraphBLAS, 20 threads (lor_land/min_plus mxv). Structure in Table B.2. relation wd_sports_team wd_cast_member wd_citizenship wd_occupation wd_subclass_of

analytic MB

Boolean

nnz

|R|

(CSR/eng)

single (CSR/eng/GB)

B=16 (eng)

single (eng/GB)

B=16 (eng)

1,136,249 1,033,124 3,063,058 4,596,658 2,024,347

70,158 51,677 3,295 41,544 6,825

11.87/9.33 10.24/11.12 47.52/34.63 64.51/43.35 28.40/19.65

0.064/0.092/0.267 0.056/0.116/0.144 0.374/0.501/0.684 0.501/0.729/0.818 0.173/0.222/0.550

0.051 0.046 0.310 0.390 0.156

0.092/0.507 0.125/0.383 0.438/3.021 0.735/3.570 0.247/1.778

0.051 0.059 0.280 0.395 0.179

CPU and destination GPU share the 119 GiB pool over the NVLink-C2C bus. This eliminates explicit host-to-device memory copies (cudaMemcpy) for both the base grammar and the active streaming frontiers (Lemma 3.3), validating the use of compressed formats to save memory footprint and bandwidth. We allocate the frontier arrays via cudaMallocManaged. Correctness. Every algorithm and configuration we report is checked on a shared input vector against independent implementations, which must agree: the host CPU reference sweep (sequential and 20-thread OpenMP), the mm-repair CPU grammar mat-vec, the materialized-CSR path (cuSPARSE for (+, ×) and Boolean, plus a CPU SpMV / min-plus over the same grammar-reconstructed matrix), the batched (SpMM) engine and cuSPARSE kernels, and, for the graph semirings, SuiteSparse:GraphBLAS (∨∧ and min, +). Across all genotype and Wikidata relations, these agree within float precision for (+, ×) and bit-for-bit for Boolean and Tropical; the full Software Heritage graph (45.7M nodes, 1.22 billion edges) is verified bit-for-bit against both the CPU reference sweep and the materialized CSR. 7 Limitations We delineate the main boundaries of the current evaluation. 1. Evaluation scale. Full biobank-scale matrices (hundreds of thousands of individuals) are our target

Tropical

but are modeled here by synthetics up to 10K samples. Whether pass-through inflation and buffer sizing remain stable at full scale without hurting occupancy remains to be verified. Mitigating this, our synthetic panels are not ad-hoc: they are simulated with msprime [30, 9], the communitystandard coalescent-with-recombination simulator, whose panels reproduce the haplotype-block and LD structure of real cohorts. What our synthetics leave untested is thus the scale of that structure, not its realism. 2. Host-side construction cost. The compressed format carries a one-time offline cost (Table B.1: the RePair build, plus milliseconds of completion and level-bucketing), amortized after a few hundred products and thus negligible in the many-vector regimes we target. 3. Batched throughput. cuSPARSE’s tuned SpMM outperforms our engine across all configurations (1.7× to 10.9×). The single-vector memory advantage also erodes at large batch: the engine holds its terminal array (one slot per each of the α distinct leaf ⟨ℓ, j⟩) and both streaming buffers B-wide, a cost proportional to (α + 2w∗ )B absent from cuSPARSE, whose CSR is B-invariant and grows only the dense X, Y . On Chr22 full the footprints

Per-vector throughput (ms/vector)

cross near B ≈ 50, past the engine’s throughputoptimal B=64, so there cuSPARSE is smaller; evaluating the leaves in the sweep rather than materializing them B-wide would drop the Θ(αB) term, which we leave to future work. What survives is structural: we never materialize the uncompressed matrix.

3 2.5

Grammar Engine (Ours) cuSPARSE (ALG DEFAULT / ALG2) cuSPARSE (CSR ALG3)

2 1.5 1

4. Graph analysis. The space advantage over CSR is conditional on the relation RePair-compressing 0.5 (Tables B.2 and 5.2) and reverses on poorly repeti16 32 64 128 256 tive ones (wd_cast_member); our GraphBLAS comBatch size B parison is a single mxv, so an end-to-end evaluation (fixpoint-iterated reachability/SSSP and batched Figure A.1: Batched right product throughput (ms/vecmulti-source, against the cuGraph BFS/SSSP ref- tor) vs. batch size B on Chr22 full. erence) is left to future work. 5. Profiling variance. Because the GB10 shares unified memory with the host (section 6), timings and energy vary with OS scheduling; the structural metrics (L, w∗ , pass-through counts; Tables B.2 and 4.1) are instead invariant, architectureindependent grammar traits.

CSR_ALG3 (abbreviated a2/a3) because the first two degrade severely as B scales up on short, ultra-wide matrix shapes, whereas CSR_ALG3 stays stable; Figure A.1 plots this degradation on Chr22 full, justifying CSR_ALG3 as our vendor baseline. With both environments optimized, cuSPARSE leads by 1.7× to 10.9× (×cuS in the table), the gap narrowing on the full real chromosomes 8 Conclusion and future work We showed and widening on the least compressible synthetics. that a right SpMV product over a grammar-compressed B Additional tables Table B.1 gives the twomatrix is a conflict-free, level-synchronous DAG gather, and that the efficient parallel model is properly lay- part host construction cost per genotype matrix, and ered : the sweep streams through two alternating buf- Table B.2 the structural figures of the five Wikidata fers, freeing frontiers immediately (Lemma 3.3). We relations. maintain RePair’s compression and ensure proper layReferences ering via post-processing pass-through completion. On large genotype matrices against cuSPARSE CSR, our [1] Amir Abboud, Arturs Backurs, Karl Bringengine’s primary advantage is space. On real chromomann, and Marvin Künnemann. Impossibility resomes, it occupies 4.0× to 4.7× less device memory at sults for grammar-compressed linear algebra. In single-vector parity (1.41× faster to 1.20× slower) and H. Larochelle, M. Ranzato, R. Hadsell, M.F. Ballower energy, scaling to 7.9× smaller and 4.2× faster on can, and H. Lin, editors, Advances in Neural Inforhighly compressible synthetics, while matching largemation Processing Systems, volume 33, pages 8810– scale polygenic scoring (y = Gβ). Requiring only an 8823. Curran Associates, Inc., 2020. associative combine, our engine extends this space advantage to repetitive graphs with up to a billion edges. [2] Jean-François Abramatic, Roberto Di Cosmo, and Stefano Zacchiroli. Building the universal archive of Extensions include: scaling to full biobank cohorts source code. Commun. ACM, 61(10):29–31, 2018. where CSR exceeds GPU memory; broadening the semiring framework (section 5) via GraphBLAS/cuGraph benchmarks; left multiplication (xT = y T M ) via pull- [3] Ankith Reddy Adudodla and Dominik Kempa. Engineering Fast and Space-Efficient Recompression based reduction; integrating height-balancing [28, 35] from SLP-Compressed Text, pages 222–232. to shorten L; exploring recompression-based constructions [3]; and targeting FPGA pipelines. [4] João Nuno Ferreira Alves, Samir Moustafa, Siegfried Benkner, Alexandre P. Francisco, WilA Batched evaluation (SpMM) Table A.1 refried N. Gansterer, and Luís M. S. Russo. Acports the engine’s optimal batch configurations (best celerating graph neural networks using a novel batch B, averaged over 40–50 iterations, verified computation-friendly matrix compression format. against the single-vector reference) against the vendor In IEEE International Parallel and Distributed cusparseSpMM kernel under its own most competitive Processing Symposium, IPDPS 2025, Milano, parameters. We sweep ALG_DEFAULT, CSR_ALG2, and

Table A.1: Batched right product (SpMM, Y = M X): best per-vector time (ms/vector) on the GB10 node. matrix

engine (B)

cuSPARSE (alg,B)

×cuS

0.024 (a3,256) 0.550 (a3,32) 0.020 (a3,256) 0.579 (a3,32) 0.025 (a3,256) 1.25 (a3,32)

4.2 2.4 4.9 2.2 4.2 1.7

0.019 (a3,256) 0.242 (a3,64) 0.121 (a3,128) 0.087 (a3,128) 0.079 (a3,256)

7.2 4.3 3.1 10.9 4.2

Real Genotypes (1000 Genomes) Chr22 (2504 × 0.10M) 0.101 (32) Chr22 (2504 × 1.06M) 1.33 (64) Chr21 (2504 × 0.10M) 0.096 (16) Chr21 (2504 × 1.05M) 1.30 (128) Chr20 (2504 × 0.10M) 0.104 (16) Chr20 (2504 × 1.74M) 2.13 (256) Synthetic Genotypes (Haplotypes) synth_small (2K × 50K) 0.141 (16) synth_large (5K × 200K) 1.04 (32) synth_ld_high (5K × 100K) 0.380 (32) synth_ld_low (5K × 100K) 0.945 (256) synth_ind_large (10K × 50K) 0.334 (16)

Table B.1: Two-part host construction cost per genotype matrix on the GB10 node. grammar RePair : one-time offline mm-repair build (shared with the CPU baseline); +pt build : the engine’s marginal completion step (passthrough, level-bucketing, terminal compaction). ≈ mat-vecs: +pt build over one right-product time (Table 4.3), i.e. the products after which it amortizes. matrix

+pt build (ms)

≈ mat-vecs

13.02 161.26 13.28 157.19 12.69 244.34

82.89 945.34 88.01 895.84 80.16 1435.69

167 141 176 140 169 122

9.72 94.36 46.50 56.04 43.90

111.85 820.39 325.59 668.65 284.71

228 274 328 179 337

grammar RePair (s)

Real Genotypes (1000 Genomes) Chr22 (2504 × 0.10M) Chr22 (2504 × 1.06M) Chr21 (2504 × 0.10M) Chr21 (2504 × 1.05M) Chr20 (2504 × 0.10M) Chr20 (2504 × 1.74M) Synthetic Genotypes (Haplotypes) synth_small (2K × 50K) synth_large (5K × 200K) synth_ld_high (5K × 100K) synth_ld_low (5K × 100K) synth_ind_large (10K × 50K)

Italy, June 3-7, 2025, pages 1091–1103. IEEE, 2025. [5] Diego Arroyuelo, José Cazorla, and Gonzalo Navarro. 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, June 2026. [6] Diego Arroyuelo, Adrián Gómez-Brandón, Aidan Hogan, Gonzalo Navarro, and Javiel RojasLedesma. Optimizing rpqs over a compact graph representation. The VLDB Journal, 33(2):349–374, Mar 2024.

[7] Diego Arroyuelo, Adrián Gómez-Brandón, and Gonzalo Navarro. Evaluating regular path queries on compressed adjacency matrices. The VLDB Journal, 34(1):2, Nov 2024. [8] Diego Arroyuelo, Aidan Hogan, Gonzalo Navarro, and Javiel Rojas-Ledesma. Time- and spaceefficient regular path queries. In 2022 IEEE 38th International Conference on Data Engineering (ICDE), pages 3091–3105, 2022. [9] Franz Baumdicker, Gertjan Bisschop, Daniel Goldstein, Graham Gower, Aaron P. Ragsdale, Georgia Tsambos, Sha Zhu, Bjarki Eldon, E. Castedo Ellerman, Jared G. Galloway, Ariella L. Gladstein, Gregor Gorjanc, Bing Guo, Ben Jeffery, Warren W. Kretzschumar, Konrad Lohse, Michael

Table B.2: Structural figures for the five Wikidata relation matrices, measured on the GB10 node (columns as in Table 4.1; total: layered rules after pass-through completion; dimensions in section 5). relation

nnz

|R|

total

wd_sports_team 1.14M 70.2K 80.9K wd_cast_member 1.03M 51.7K 62.3K wd_citizenship 3.06M 3.3K 3.7K wd_occupation 4.60M 41.5K 46.7K wd_subclass_of 2.02M 6.8K 7.6K

w∗

+pt

6 49.3K 8 42.4K 5 2.4K 8 19.5K 5 5.8K

10.7K 10.7K 432 5.1K 776

L

Matschiner, Dominic Nelson, Nathaniel S. Pope, and querying uk biobank data. Bioinformatics, Consuelo D. Quinto-Cortés, Murillo F. Rodrigues, 39(9):btad552, 09 2023. Kumar Saunack, Thibaut Sellinger, Kevin Thornton, Hugo van Kemenade, Anthony W. Wohns, [17] Timothy A. Davis. Algorithm 1000: SuiteSparse:GraphBLAS: Graph algorithms in the lanYan Wong, Simon Gravel, Andrew D. Kern, Jere guage of sparse linear algebra. ACM Transactions Koskela, Peter L. Ralph, and Jerome Kelleher. on Mathematical Software, 45(4):44:1–44:25, 2019. Efficient ancestry and mutation simulation with msprime 1.0. Genetics, 220(3):iyab229, 2022. [18] Roberto Di Cosmo. Archiving and referencing source code with software heritage. In Anna Maria [10] Sebastian Baunsgaard and Matthias Boehm. Bigatti, Jacques Carette, James H. Davenport, Aware: Workload-aware, redundancy-exploiting Michael Joswig, and Timo de Wolff, editors, Mathlinear algebra. Proc. ACM Manag. Data, 1(1), May ematical Software – ICMS 2020, pages 362–373, 2023. Cham, 2020. Springer International Publishing. [11] Guy E. Blelloch. Prefix sums and their applica[19] Daniel Donenfeld, Stephen Chou, and Saman Amations. Technical Report CMU-CS-90-190, School rasinghe. Unified compilation for lossless compresof Computer Science, Carnegie Mellon University, sion and sparse computing. In 2022 IEEE/ACM November 1990. International Symposium on Code Generation and Optimization (CGO), pages 205–216, 2022. [12] P. Boldi and S. Vigna. The webgraph framework i: compression techniques. In Proceedings of the [20] Richard Durbin. Efficient haplotype matching 13th International Conference on World Wide Web, and storage using the positional burrows–wheeler WWW ’04, page 595–602, New York, NY, USA, transform (pbwt). Bioinformatics, 30(9):1266– 2004. Association for Computing Machinery. 1272, 05 2014.

[13] Nieves R. Brisaboa, Susana Ladra, and Gonzalo [21] Ahmed Elgohary, Matthias Boehm, Peter J. Haas, Navarro. k2-trees for compact web graph represenFrederick R. Reiss, and Berthold Reinwald. Comtation. In Jussi Karlgren, Jorma Tarhio, and Heikki pressed linear algebra for large-scale machine learnHyyrö, editors, String Processing and Information ing. Proc. VLDB Endow., 9(12):960–971, August Retrieval, pages 18–30, Berlin, Heidelberg, 2009. 2016. Springer Berlin Heidelberg. [22] Ahmed Elgohary, Matthias Boehm, Peter J. Haas, [14] Shing Wan Choi, Timothy Shin Heng Mak, and Frederick R. Reiss, and Berthold Reinwald. ComPaul F. O’Reilly. Tutorial: a guide to performing pressed linear algebra for declarative large-scale polygenic risk score analyses. Nature Protocols, machine learning. Commun. ACM, 62(5):83–91, 15(9):2759–2772, sep 2020. April 2019. [15] The 1000 Genomes Project Consortium. A global [23] P. Ferragina, G. Manzini, T. Gagie, D. Köppl, reference for human genetic variation. Nature, G. Navarro, M. Striani, and F. Tosoni. Improving 526(7571):68–74, Oct 2015. matrix-vector multiplication via lossless grammarcompressed matrices. Proc. VLDB Endowment, [16] Davide Cozzi, Massimiliano Rossi, Simone Rubi15(10):2175–2187, 2022. nacci, Travis Gagie, Dominik Köppl, Christina Boucher, and Paola Bonizzoni. µ-pbwt: a [24] Tommaso Fontana, Sebastiano Vigna, and Stefano lightweight r-indexing of the pbwt for storing Zacchiroli. Webgraph: The next generation (is in

rust). In Companion Proceedings of the ACM Web Conference 2024, WWW ’24, page 686–689, New York, NY, USA, 2024. Association for Computing Machinery.

[35] G. Navarro, F. Olivares, and C. Urbina. Balancing run-length straight-line programs. In SPIRE, 2022. Also arXiv:2206.13027.

[36] Gonzalo Navarro. Compact Data Structures: A Practical Approach. Cambridge University Press, [25] Alexandre P. Francisco, Travis Gagie, Dominik 2016. Köppl, Susana Ladra, and Gonzalo Navarro. Graph compression for adjacency-matrix multiplication. [37] Gonzalo Navarro. Indexing highly repetitive string SN Computer Science, 3(3):193, Mar 2022. collections, part i: Repetitiveness measures. ACM Comput. Surv., 54(2), March 2021. [26] Travis Gagie, Gonzalo Navarro, and Nicola Prezza.

Optimal-time text indexing in bwt-runs bounded [38] Gonzalo Navarro. Indexing highly repetitive string space. In Proceedings of the Twenty-Ninth Annual collections, part ii: Compressed indexes. ACM ACM-SIAM Symposium on Discrete Algorithms, Comput. Surv., 54(2), February 2021. SODA ’18, page 1459–1477, USA, 2018. Society for Industrial and Applied Mathematics. [39] NVIDIA Corporation. cuSPARSE library. https:// docs.nvidia.com/cuda/cusparse/index.html, 2026. [27] Travis Gagie, Gonzalo Navarro, and Nicola Prezza. Version 13.x (part of CUDA Toolkit). Fully functional suffix trees and optimal text searching in bwt-runs bounded space. J. ACM, [40] Antoine Pietri, Diomidis Spinellis, and Stefano Zacchiroli. The software heritage graph dataset: 67(1), January 2020. Large-scale analysis of public software development [28] M. Ganardi, A. Jeż, and M. Lohrey. Balanchistory. In Proceedings of the 17th International ing straight-line programs. Journal of the ACM, Conference on Mining Software Repositories, MSR 68(4):27:1–27:40, 2021. ’20, page 1–5, New York, NY, USA, 2020. Association for Computing Machinery. [29] Gene H. Golub and Charles F. Van Loan. Matrix Computations - 4th Edition. Johns Hopkins Uni- [41] Philipp Scharpf, Corinna Breitinger, Andreas Spitz, Norman Meuschke, André Greiner-Petter, versity Press, Philadelphia, PA, 2013. Moritz Schubotz, and Bela Gipp. Entity linking with wikidata: A systematic literature review. [30] Jerome Kelleher, Alison M. Etheridge, and Gilean ACM Comput. Surv., 58(9), February 2026. McVean. Efficient coalescent simulation and genealogical analysis for large sample sizes. PLOS [42] F. Tosoni, P. Bille, V. Brunacci, A. De Angelis, Computational Biology, 12(5):e1004842, 2016. P. Ferragina, and G. Manzini. Toward greener matrix operations by lossless compressed formats. [31] N. Jesper Larsson and Alistair Moffat. Offline IEEE Access, 13:56756–56773, 2025. dictionary-based compression. In Proceedings of the Conference on Data Compression, DCC ’99, page [43] Luca Versari, Iulia-Maria Comsa, Alessio Conte, 296, USA, 1999. IEEE Computer Society. and Roberto Grossi. Zuckerli: A new compressed representation for graphs. IEEE Access, 8:219233– [32] Veli Mäkinen, Alexandru I. Tomescu, Anna Kuos219243, 2020. manen, Topi Paavilainen, Travis Gagie, and Rayan Chikhi. Sparse dynamic programming on dags with [44] Denny Vrandečić and Markus Krötzsch. Wikidata: small width. ACM Trans. Algorithms, 15(2), Feba free collaborative knowledgebase. Commun. ruary 2019. ACM, 57(10):78–85, September 2014.

[33] Giosuè Cataldo Marinò, Flavio Furia, Dario Malchiodi, and Marco Frasca. Efficient and compact representations of deep neural networks via entropy coding. IEEE Access, 11:106103–106125, 2023.

[45] Denny Vrandečić, Lydia Pintscher, and Markus Krötzsch. Wikidata: The making of. In Companion Proceedings of the ACM Web Conference 2023, WWW ’23 Companion, page 615–624, New York, NY, USA, 2023. Association for Computing Machinery.

[34] Tom Mens, Coen De Roover, and Anthony Cleve, editors. Software Ecosystems: Tooling and Analytics. Springer Cham, 1 edition, 2023. XXII, 314 [46] Feng Zhang, Zaifeng Pan, Yanliang Zhou, Jidong pages. Zhai, Xipeng Shen, Onur Mutlu, and Xiaoyong

Du. G-tadoc: Enabling efficient gpu-based text analytics without decompression. In 2021 IEEE 37th International Conference on Data Engineering (ICDE), pages 1679–1690, 2021. [47] Yanliang Zhou, Feng Zhang, Tuo Lin, Yuanjie Huang, Saiqin Long, Jidong Zhai, and Xiaoyong Du. F-tadoc: Fpga-based text analytics directly on compression with hls. In 2024 IEEE 40th International Conference on Data Engineering (ICDE), pages 3739–3752, 2024.