Taurus∗: Accelerating Out-of-Core Graph Neural Network Inference on Billion-Scale Graphs † Pranjal Naman
and Yogesh Simmhan1
Department of Computational and Data Sciences (CDS), Indian Institute of Science (IISc), Bangalore 560012 India
arXiv:2607.17374v1 [cs.DC] 19 Jul 2026
Email:{pranjalnaman, simmhan}@iisc.ac.in
Abstract Graph Neural Network (GNN) inference on billion-scale graphs is challenging due to the large memory footprint of features and embeddings and high disk I/O costs in outof-core settings. Existing distributed GNN systems incur high communication times and infrastructure costs while disk-based GNN systems are primarily tailored to training and experience massive wasted reads during inference on the entire graph. We present Taurus, a single-machine system for GNN inference on graphs that do not fit in RAM, supporting both exact full-graph inference and fanout-sampled inference. To avoid random and repeated feature gathers, Taurus reformulates layer-wise inference as source-centric broadcasts over sequential SSD scans, backed by a pipelined GPU–CPU–SSD hierarchy, topology-aware reordering, pending-message eviction and a GPU-resident store for highdegree vertices. It further uses non-buffered sequential reads and GPU-backed writes to reduce page-cache pollution, host-memory pressure and write overheads. On out-of-core graphs with up to 269M vertices, 4B edges, and 514 GiB of features, Taurus outperforms the strongest layer-wise baseline, DGI, by 7–25×, and vertex-wise baselines by 40–140×.
1
Introduction
Graph Neural Networks (GNNs) are an effective tool for learning representations from graph data, capturing both topology and associated features [2, 3]. This makes them adept at performing a variety of tasks like detecting fraud in transaction networks [4, 5], predicting traffic flows and signaling in Intelligent Transportation Systems (ITS) [6, 7], and making e-commerce recommendations [8]. Most GNNs follow a message-passing paradigm, where each vertex iteratively gathers and aggregates information from its neighbors while applying a neural-network transformation at every layer. This iterative neighborhood aggregation incurs significant computational and memory overhead due to irregular graph accesses and repeated neural network computations. ∗
Taurus is the constellation that houses the star, Atlas, reflecting its evolution from our prior out-of-core GNN inference framework, ATLAS [1]. † Extended full-length version of paper that appeared at HPDC 2026: “ATLAS: Efficient Out-of-Core Inference for Billion-Scale Graph Neural Networks”, Pranjal Naman and Yogesh Simmhan, in the 35th ACM International Symposium on High-Performance Parallel and Distributed Computing (HPDC), 2026. DOI: https://doi.org/10.1145/3806645.3807597
1
Motivation Given these memory and computational costs, optimizing inference is critical for real-world deployments. GNNs are often deployed on evolving graphs requiring predictions to be periodically refreshed [9, 10]. While this refresh is typically not latency-sensitive, it must complete within practical timeframes, i.e., hours rather than days, using either exact or sampled inference depending on application needs. This is further complicated by the scale of realworld graphs, often containing millions–billions of vertices and edges, e.g., fintech transaction networks for fraud detection and social or e-commerce graphs for recommendation tasks. Although graph topology may fit in the RAM of a single machine, vertex features/embeddings often dominate the memory footprint. E.g., the IGB-Full citation graph [11] (269M vertices, 4B edges, 1024 FP16 features) requires 514 GiB of RAM (Table 1). To circumvent this, distributed GNN systems partition graphs across multiple servers and perform training or inference collaboratively on distributed subgraphs, incurring high infrastructure and network costs. Recent inference systems further optimize performance through probabilistic caching of remote features and communication-aware graph and feature partitioning [12, 13]. Disk-based GNN training has emerged as a promising alternative [14–18], along the lines of prior Out-of-Core (OOC) disk-based parallel graph processing [19, 20]. These focus on efficient data layouts and intelligent caching strategies to fully utilize memory, disk capacity, and bandwidth on a single machine. This is all the more relevant given the lower latency and higher bandwidth of Solid State Disks (SSDs) with capacities of 2 TiB+ common even for prosumer disks, at a much lower price point than RAM. While disk-based OOC GNN training has been extensively studied, inference presents distinct challenges that remain largely unexplored. (i) Working-Set Amplification During Inference: Existing OOC training systems optimize data transfer [21, 22], data organization through reordered movement and multi-tier caching [14, 16–18], and storage layouts tailored for efficient training [14]. While inference, at first glance, appears to be just the forwardpass phase of training, these are fundamentally different workloads. Training operates on a small labeled subset of vertices (e.g., ≈ 1% in OGBN-Papers100M [23]), whereas inference computes embeddings for the entire graph. As a result, OOC training optimizations that exploit a limited working set, such as computation-graph precomputation [14, 16], do not readily apply because materializing computation graphs for all vertices is prohibitively expensive. (ii) Sampling Effects on Inference: Training typically employs neighborhood sampling [3] to bound memory usage by aggregating over only a subset of neighbors. Inference may require exact full-neighborhood aggregation for accuracy-critical domains such as fintech and ITS, where sampling can introduce non-deterministic predictions [10,24,25]; other applications may tolerate fanout-sampled inference. In both modes, disk-resident inference still streams vertex features or intermediate embeddings across layers, and sampling mainly reduces message propagation. E.g., training a 2-layer GNN on IGB-Large with 1% labeled vertices touches only ≈ 8% of all vertices per epoch [26]. Consequently, OOC GNN inference cannot directly inherit training optimizations. It must control read amplification, active-state residency and output materialization. Challenges Existing out-of-core (OOC) GNN inference typically uses either vertex-wise [14, 18, 27] or layer-wise [28] gather-based execution. Vertex-wise methods (e.g., DGL [27]) recursively aggregate k-hop neighborhoods for each batch, while layer-wise methods (e.g., DGI [28]) compute embeddings for all vertices one layer at a time to eliminate redundant neural-network computations. Despite these differences, both suffer severe I/O bottlenecks on large diskresident graphs. We demonstrate these inefficiencies using DGI (our strongest layer-wise 2
1
0.0 DGI DGI 0 L1 L2
(a) Random and repeated accesses
Verts. Req. GiB
Bytes Req./Load (GiB, )
0.2
2
Verts./Blks. Req./Load (M)
0.4
3
Redundancy Factor
Acc. Density / Batch
0.6
500 Blks. Req. PA Size 500 Blks. Load. 400 400 300 300 200 200 100 100 0 DGI DGI 0 L1
L2
(b) Read amplification
Figure 1: Challenges limiting out-of-core GNN inference on Papers [23]. (a) Random and repeated accesses quantified by access density (left Y axis) and redundancy factor (stars, right Y axis). (b) Read amplification: requested vertices/blocks and loaded blocks (bars, left Y axis), requested/loaded bytes (markers, right Y axis), and Papers feature size (dashed line). baseline) for 2-layer GraphConv (fanout=10 ) [2] sampled inference on the RCMK-reordered PApers graph (Table 1) on a 32 GiB RAM machine: (1) Random Access: Neighbors required for inference are scattered across the feature store (Fig. 3b, top). Access density – the ratio of neighbors fetched per batch to the ID range they occupy – for DGI stays ≈ 0.2–0.3 across both layers (Fig. 1a, box plots, left Y axis), well below 1 even with RCMK reordering (1 = contiguous accesses). (2) Repeated Access: Neighborhoods fetched across batches have a high overlap leading to repeated fetches (Fig. 3b, center ). Redundancy factor – the total vertex fetches divided by unique vertices processed overall – is ≈ 3× for both layers using DGI (Fig. 1a, stars, right Y axis). (3) Read Amplification: Random and repeated accesses translate directly into excess disk I/O; since storage is accessed at block granularity (e.g., 4 KiB), scattered reads fetch data that remains unused (Fig. 3b, bottom). Across both layers, DGI requests ≈ 750 GiB of feature blocks (Fig. 1b, teal circles) against only ≈ 300 GiB of features actually needed by the requested vertices (orange circles). Although the OS page cache absorbs a portion of this, physical reads still total 568 GiB (green circles) – nearly 1.9× the useful bytes requested, and over 8.35× of PA’s 68 GiB total size. The common root cause is destination-centric gathering. Since each destination independently pulls neighbor embeddings, systems repeatedly move the same source data, lose sequentiality and amplify disk traffic even with graph reordering. We also empirically demonstrate these challenges in Fig. 2 by comparing our proposed Taurus with three baselines: DGI [28] (layer-wise inference), and Ginex [18] and DGL [27] (training frameworks adapted for vertex-wise inference), using sampled 2-layer GraphConv inference [2] over PApers and MAG-Cites (Table 1), on a GPU workstation with 128 GiB RAM, RTX 5090 GPU and 2 TiB SSD (§ 4.1). Both DGL [27] and Ginex [18] are unable to complete the inference even within a 4 h time budget on MA, taking an extrapolated ≈ 16 h (Fig. 2b, left Y axis), while DGI [28] took ≈ 2.5 h. In contrast, our Taurus framework completes this in < 0.5 h for both layers.
3
101
2
100
L1 L2 L1 L2 0 GX DG DI TA
105
60
104
45
103
30
102
15
101
L1 L2 L1 L2 0 GX DG DI TA
(a) PA/GCN2
Relative Speedup
4
Relative Speedup
102
Execution Time (s) [log]
103
Extrapolated 8 6
Execution Time (s) [log]
104
(b) MA/GCN2
Figure 2: Inference time (left Y axis, bars) and speedup relative to Taurus (right Y axis, markers) for 2-layer GraphConv inference with fanout=10 on disk-resident graphs (topology+features) reported for GineX, DGL, DGI, and TAurus (ours) on PApers and MAG-Cites using a 5090 GPU. Proposal To address these challenges, we leverage a key insight: layer-wise GNN inference can be reformulated as source-centric broadcasts, enabling sequential disk access instead of repeated random gathers. This reduces read amplification, but naïve broadcasts merely shift the bottleneck to partial-state memory pressure and random output writes. Thus, Taurus must preserve sequential reads while bounding active states and avoiding random-write amplification. Contributions We present Taurus, a disk-based GNN inference framework for billion-scale graphs on a single workstation. Taurus supports exact full-graph and fanout-sampled inference by replacing destination-centric gathers with source-centric broadcasts, combined with tiered GPU–RAM– SSD aggregation, topology-aware reordering and eviction, and pipelined I/O, aggregation, GPU compute and output. Specifically, we make the following contributions: 1. Broadcast-based inference model. Taurus introduces broadcast-based layer-wise inference that reads features and embeddings sequentially, reducing read amplification relative to gather-based approaches. It supports both exact full-graph and fanout-sampled inference while preserving their respective GNN semantics. 2. GPU-enhanced tiered runtime. Taurus employs a pipelined GPU–RAM–SSD hierarchy that overlaps I/O, aggregation, transformation, and output, reducing memory pressure and data movement overheads while maintaining high throughput under constrained memory. 3. Architecture support and topology-aware execution. We extend broadcast inference to Graph Attention Networks and memory-efficient GraphSAGE, and introduce a topologyaware graph reordering strategy that reduces partial-state residency and evict-reload cycles during execution. 4. Comprehensive evaluation. We evaluate Taurus on billion-scale citation and social-network graphs with feature sizes above 500 GiB under exact and sampled inference. We show substantial reductions in disk traffic and runtime over DGL GraphBolt-backed and outof-core baselines, with ablations of key runtime components.
4
2 1
ℎ
ℎ
(a) Graph ℎ
𝐵 ℎ
ℎ
𝐵 ℎ
ℎ
Legend
36
14
0 2 4
0
4
Scattered random reads per vertex
Repeated 𝐵 ℎ ℎ 14 36 Reads 0 4 0 4 2 𝐵 ℎ Amplified 3 𝐵 ℎ ℎ 0 Reads 4 On Disk 𝐵 ℎ ℎ Feature Store 2𝐵 ℎ ℎ
4
𝐵 ℎ
Random Reads
Feats. of 0 and 4 read twice
…
5
ℎ
𝐵
3
0
𝐵 ℎ
(b) Gather-based GNN Methods 21
0
1
Wasted reads (ℎ , ℎ , ℎ 𝑤𝑎𝑠𝑡𝑒𝑑)
3
0 0 4
0
𝐵 ℎ
ℎ
𝐵 ℎ
ℎ
𝐵 ℎ
ℎ
21 0
0 0
3
2
44
1
Completed
1
Aggregated Value (c) TAURUS Broadcast
𝐵 ℎ
ℎ
𝐵 ℎ
ℎ
𝐵 ℎ
ℎ
21 4
1
0 0 4
36
9
Feature Read
Figure 3: Gather-based versus broadcast-based execution for one GNN layer. This article extends our previous conference work ATLAS [1], which introduced broadcastbased layer-wise GNN inference using a RAM–SSD hierarchy and pipelined execution. Taurus extends it with: (1) GPU-resident aggregation and GPU-backed output materialization to reduce host-memory pressure and write overheads; (2) A topology-aware reordering objective that targets partial-state residency and eviction–reload cycles; (3) Support for fanout-sampled inference, memory-efficient GraphSAGE and Graph Attention Networks; and (4) Additional billion-scale datasets, a DGL GraphBolt baseline and expanded ablations.
2
Background
2.1
GNN Training and Inference
A GNN layer gathers and aggregates the embeddings of a vertex u’s in-neighbors (N − (u)) to produce an intermediate representation xlu (Eqn. 1), which is transformed by a learnable Update function and non-linear activation σ(·) to generate hlu (Eqn. 2). Training repeats this for L layers in the forward pass and then updates parameters through backpropagation. − xlu = Aggregatel ({hl−1 v , v ∈ N (u)})
(1)
l hlu = σ(Updatel (hl−1 u , xu ))
(2)
Aggregating all in-neighbors recursively leads to neighborhood explosion and out-of-memory (OOM) errors [3]. To mitigate this, training typically employs neighborhood sampling, where only a subset of neighbors is aggregated at each hop [3]. In contrast, GNN inference requires only the forward pass and may use either exact fullneighborhood aggregation for deterministic embeddings [10, 24, 25] or fanout sampling when approximation is acceptable. Full-neighborhood vertex-wise inference suffers from neighborhood explosion and redundant computation, motivating layer-wise inference [28], which materializes embeddings for all vertices one layer at a time. Taurus adopts and optimizes layer-wise inference for out-of-core execution.
2.2
Gather-based Execution Model
Most GNN systems implement message passing as gather-based execution: a destination vertex u retrieves in-neighbor embeddings {hl−1 | v ∈ N − (u)} to compute hlu . When embeddings v 5
Topology
Features/Embeddings 1
2 0
1
Partition 0
Read Queue
…
…
0
2
1
Partition 1
Sparse CSR Representation
0
Partition N O_DIRECT Reads
Chunk
Chunk
Ingest & Assemble
Orchestrator Transform Send Msgs. OG NS
HS
CP
EV
Vertex States
Memory Manager GPU Store
Pending Messages
taurus I/O
Graph Reader
Update Layer Transform
Writer Eviction Policy Min Heap
Hot Store
Cold Store Spilled Buffers
Local Control
Graduation Processor GPU Queue
Grad. Buffer
MMAP I/O
2
Write Queue
Write Buffers & Run Files
…
GDS
Figure 4: Taurus Architecture: Sequential graph/embedding reads, tiered aggregation, GPU transformation & run-file output. reside on disk, these destination-centric gathers cause irregular and repeated accesses (Fig. 3). Message-passing execution consists of sample, gather, transfer and compute. Sampling selects a subset of in-neighbors to limit neighborhood explosion [3]; exact inference instead uses all inneighbors. Gather retrieves selected or full-neighborhood embeddings, transfer moves them to the GPU, and compute applies aggregation followed by neural-network transformations. Popular architectures such as GraphConv [2], GraphSAGE [3], GIN [29] and GAT [30] follow this paradigm. Taurus targets their out-of-core inference.
2.3
Layer-wise Inference
A common inference approach reuses training code with the backward pass disabled (vertexwise inference) [27]. Under full-neighborhood aggregation, overlapping neighborhoods cause memory blowup and redundant computation. Layer-wise inference instead materializes embeddings for all vertices once per layer and reuses them in later layers [9, 10, 28]. However, layer-wise execution does not eliminate redundant data movement. In OOC settings, gatherbased execution still makes each destination independently fetch in-neighbor embeddings, so read volume grows with propagated messages (≈(sampled) edges) rather than unique vertices. Graph reordering [28, 31] improves locality but cannot remove these repeated gathers. Since inference is often lightweight, OOC execution becomes I/O-bound, and repeated layer-wise access patterns amplify cumulative read traffic.
3
System Design
3.1
Broadcast Execution Model and Challenges
Taurus replaces destination-centric gathers with source-centric broadcasts, streaming each source feature/embedding once per pass in vertex order (Fig. 3b), for both full-neighborhood and sampled message sets. In contrast, gather-based inference performs scattered reads, repeatedly fetches shared source features, and loads unused records from block-granularity storage. Broadcast is semantically equivalent to gather for the same message set, Aggregate, and Update functions, but realizing it for OOC graphs introduces several challenges. 1. Bounded memory. Broadcast execution eliminates repeated feature reads but creates many partially aggregated destination states per layer. For large graphs, these states cannot 6
be fully materialized in memory and must be managed under a fixed budget. 2. Vertex ordering. Source-vertex order determines when destinations become active and complete. Poor ordering lengthens partial-state lifetimes, increasing RAM pressure and SSD spills. 3. Output materialization. Vertices complete only after receiving all messages, so completion order is not sequential. Writing outputs immediately can therefore replace random-read amplification with random-write amplification. 4. I/O efficiency and overlap. Topology, embeddings, aggregation, transformation and output must be organized and pipelined so that sequential I/O, compute and writes do not stall each other.
3.2
Data Layout
The on-disk layout is central to OOC performance. Since Taurus broadcasts along vertex out-edges, it stores topology in Compressed Sparse Row (CSR) format (Fig. 4, top), requiring O(|V | + |E|) space and enabling sequential source-vertex scans via file offsets by the graph reader. For features and intermediate embeddings, vertices do not complete in vertex-ID order under broadcast execution. A dense in-memory embedding array is infeasible, while dense on-disk placement causes costly random writes. Externally sorting completed embeddings would add large temporary storage and I/O. To address these challenges, Taurus partitions the vertex-ID space into fixed ranges and maintains embeddings in sorted order within each range (Fig. 4, top, Partition 0–N ). Since an individual range may still exceed available memory, each partition is materialized as multiple sorted run files, generated by sorting buffered embeddings in-memory and flushing them sequentially to disk. Merging spilled run files would require another large multi-way merge. Instead, Taurus leaves runs unmerged and lets the graph reader reconstruct a sequential view on demand.
3.3
Graph Reader
Taurus uses a partially sequential graph reader that streams vertex embeddings once per pass (Fig. 4, orange). Input and output embeddings are stored as run files of sorted (vertex_id, embedding) records within contiguous vertex-ID ranges. Each run f stores IDs If , feature matrix Xf , and nf records, indexed by its ID range (If [0], If [nf −1]). Multiple run files may exist within a partition. The graph reader exposes a chunk -based iterator over contiguous vertex ranges. Each chunk contains CSR out-neighbors/offsets and embeddings in vertex order, and is enqueued into a reader queue for downstream processing. For chunk size C and feature size F , a chunk contains ⌊C/F ⌋ vertices. Chunks are partitioned by feature bytes rather than edge count. So high-degree vertices increase per-chunk edge work but not feature-read ordering. For a chunk spanning vertex IDs [s, e), the reader identifies overlapping runs. For each run f , two binary searches compute ℓf = min{i | If [i] ≥ s} and rf = min{i | If [i] ≥ e}, yielding rows [ℓf , rf ) within the chunk range. It issues one aligned pread per overlapping run using direct I/O (O_DIRECT), bypassing the OS page cache for single-pass embedding scans. Sorted rows from different runs are merged by vertex ID to reconstruct the chunk matrix. This merge-on-read avoids an external merge sort over all output embeddings at each layer. Run file descriptors are opened lazily, and a dedicated reader thread overlaps I/O with computation. 7
3.4
Orchestrator
The orchestrator preserves GNN semantics under broadcast execution (Fig. 4, green). It consumes reader chunks from the reader queue (§ 3.3), initializes the memory manager and graduation processor, and maintains compact O(|V |) arrays for each vertex’s pending-message count and execution state. For a GNN layer with mean aggregation, instead of gathering neighbors at destination (l) v, the orchestrator emits one normalized message per propagated edge (u, v): mu→v = (l) 1 |M(v)| hu , where M(v) is the full or sampled in-neighbor set propagated for v. It sends (l)
⟨v, mu→v , sv ⟩ to the memory manager, where sv is v’s current state. The orchestrator tracks partial aggregation using the state machine in Fig. 4. Vertices start in NOT_STARTED, move to IN_BUFFER when their first message creates a hot-store state, transition to EVICTED if spilled to the SSD cold store, or to ON_GPU if pinned in the GPU store (e.g., high in-degree vertices). After all expected messages arrive, they enter COMPLETED and become eligible for graduation. Valid transitions are NOT_STARTED → IN_BUFFER/ON_GPU, IN_BUFFER → EVICTED, EVICTED → IN_BUFFER, and IN_BUFFER/ ON_GPU → COMPLETED; ON_GPU vertices are pinned.
3.5
Memory Manager
The memory manager maintains partial aggregation states under a fixed (configurable) RAM budget (Fig. 4, red ) using a three-tier GPU–RAM–SSD hierarchy: GPU store, host-memory hot store and SSD-backed cold store. High-traffic states are pinned in GPU memory, active states occupy the hot store, and overflow states spill according to the eviction policy. 3.5.1
GPU Store
The highest tier is a GPU store for high in-degree vertices. Since power-law graphs route many messages to a small hub set, the memory manager pins their aggregation buffers (ON_GPU) in VRAM during initialization. Messages to these vertices bypass the hot store and accumulate directly in VRAM, reducing RAM pressure. 3.5.2
Hot Store
The hot store is a fixed-size host-memory slot array, with each slot holding one active vertex’s partial aggregation state (IN_BUFFER). Messages accumulate through a vertex-to-slot map; slots are allocated on entry to IN_BUFFER and released at COMPLETED. When full, selected states are evicted to the SSD-backed cold store. Since aggregation proceeds only for states in the hot or GPU store, evicted states must be reloaded before receiving further messages, and the memory manager reports all state transitions to the orchestrator. 3.5.3
Taurus Eviction Policy
When the hot store is full, the eviction policy selects spill victims. Random or recency-based choices may repeatedly evict vertices far from completion, causing eviction–reload cycles and extra SSD I/O. Taurus instead evicts vertices with the fewest pending messages, which are closest to graduation and least likely to be reloaded repeatedly from disk. Taurus’s eviction policy maintains an ordering of IN_BUFFER vertices by pending message count while supporting insertions, removals, score updates, and selection of the k lowestscoring vertices. Rather than a conventional Python heap, it exploits the bounded integer score range [1, max_in_degree], where a vertex’s score equals its pending message count (0
8
means COMPLETED). Vertices are organized into score-indexed buckets implemented as doubly linked lists, enabling O(1) insertion, removal, and score updates. Eviction proceeds by scanning the lowest non-empty buckets, yielding O(k) selection of k eviction candidates. 3.5.4
Cold Store
The cold store is a NumPy mmap file. Unlike single-pass feature and embedding scans, it uses buffered I/O because evicted states may be reloaded in later chunks, allowing the OS page cache to absorb eviction–reload traffic.
3.6
Graduation Processor
The graduation processor handles completed vertices (Fig. 4, blue). When a vertex’s pending message count reaches zero, the orchestrator instructs the memory manager to finalize its aggregation, append it to a configurable graduation buffer, and release the hot-store slot. Once full, a graduation buffer is enqueued for GPU offload. Double buffering lets one buffer collect newly completed vertices while the other is processed asynchronously. A dedicated GPU-offload thread applies the layer transformation using CUDA streams to overlap transfers and compute, preventing the orchestrator from blocking. Transformed embeddings are then enqueued to the write queue.
3.7
Embedding Writer
Transformed embeddings are materialized on disk for the next layer or final output. A dedicated writer consumes embeddings from the write queue and stages them in GPU-resident partition buffers. Since vertices arrive in graduation rather than vertex-ID order, the writer range-partitions them by vertex ID; when a partition buffer fills, it is sorted on the GPU and flushed sequentially as a sorted run file. Taurus writes run files through kvikio/cuFile, using GPUDirect Storage (GDS) when supported by the platform, and cuFile compatibility mode otherwise. On platforms where true GDS is unavailable, such as our consumer-GPU machine, we use cuFile-managed bounce buffers1 , but the interface remains unchanged. Even without true GPU-to-storage DMA, GPU-resident partition buffers avoid host-side sorting, reduce RAM pressure and leave more memory for the hot store.
3.8
Topology-aware Graph Reordering
Processing vertices in the original vertex-ID order can significantly increase hot-store residency. Vertices processed early may not complete until much later, leaving their partial states resident for extended periods. This increases memory pressure, eviction frequency, and ultimately execution time. While the ATLAS greedy reordering strategy maximizes completion rate, it does not explicitly minimize the span, i.e., the interval between a destination’s first and last incoming message. Next, we characterize vertex residency and I/O overhead in terms of vertex span and provide empirical evidence in Fig. 5. Let G = (V, E) be a directed graph, where π(v) denotes the processing rank of vertex v, and N − (v) its in-neighbors. A vertex becomes active when it receives a message from its first in-neighbor and completes after receiving the message from its last. Accordingly, its activation and completion ranks are a(v) = minu∈N − (v) π(u) and c(v) = maxu∈N − (v) π(u), respectively, yielding the span L(v) = c(v) − a(v). Let A(t) denote the number of active 1
https://docs.nvidia.com/gpudirect-storage/api-reference-guide/#cufile-compatibility-mode
9
Algorithm 1 Taurus topology-aware reordering Require: Graph G = (V, E), initial ordering π 1: Initialize r(v) ← π(v) 2: repeat P 3: µ(v) ← d−1(v) u∈N − (v) r(u), ∀v P − v∈N + (w) µ(v)/d (v) , ∀w 4: r(w) ← P − v∈N + (w) 1/d (v) 5: π ← ArgSort(r) 6: until π converges 7: return π
vertices at processing rank t. Then: A(t) =
X 1 a(v) ≤ t < c(v) , v∈V
and the cumulative active-state occupancy is X XX 1 a(v) ≤ t < c(v) M= A(t) = t
t
v∈V
XX X = 1 a(v) ≤ t < c(v) = L(v) = C(π) v∈V
t
v∈V
Thus, reducing C(π) lowers cumulative residency, eviction pressure and I/O. The objective is: X max π(u) − min π(u) . min C(π) = min π
π
v∈V
u∈N − (v)
u∈N − (v)
Bandwidth-minimization heuristics such as RCMK [31] optimize B = max(u,v)∈E |π(u) − π(v)|, a worst-case edge-length metric. Since all in-neighbors of v lie within [π(v) − B, π(v) + B], L(v) ≤ 2B and C(π) ≤ 2|V |B. This loose bound makes RCMK a useful bootstrap, but it does not directly minimize total span. Since optimizing C(π) over |V |! orderings is infeasible and bandwidth only bounds it loosely, Taurus instead minimizes in-neighborhood dispersion. Each vertex receives a realvalued position r(v), and we measure the spread of the in-neighbors of each destination around a center µ(v): J(r, µ) =
X v∈V
1 − d (v)
X
2 r(u) − µ(v) .
u∈N − (v)
The 1/d− (v) normalization factor gives each destination equal weight, preventing high-degree vertices from dominating the objective. In GNNs, each vertex aggregates a self message, so d− (v) ≥ 1 (normalization and µ(v) are well-defined); a vertex with only its self message has L(v) = 0. Since J is differentiable, we iteratively update the locality centers µ and the vertex positions r. For a fixed ordering r, setting ∂J/∂µ(v) = 0 gives µ∗ (v) =
1 d− (v)
10
X u∈N − (v)
r(u),
Reloads(×106)
Inf. Time(min)
120 100 80 60 40 20 0.5
R 2 = 0.95 R 2 = 0.83 1.0 1.5 2.0 C( )(×1015)
200
J (×1021)
(a) Inference time vs. C(π)
FS IL
5 0
FS IL
RC OG
R 2 = 0.93
R 2 = 0.96 1.0 1.5 2.0 C( )(×1015)
0 0.5
(b) Cold reloads vs. C(π)
RC OG
R 2 = 0.96
R 2 = 0.99 1.0 1.5 2.0 C( )(×1015)
0.5
(c) J vs. C(π)
Figure 5: Topology-based reordering on Friendster (FS) and IGB-large (IL) with RCMK (RC) and Original (OG) bootstraps, showing correlation between total span C(π) and (a) total inference time; (b) number of reloads from cold store, i.e., I/O cost; and (c) our proposed objective J. Each point corresponds to an intermediate ordering generated during Alg. 1 starting from OG and RC bootstraps. the average position of a vertex’s in-neighbors. Likewise, for fixed locality centers µ, setting ∂J/∂r(w) = 0 gives X r∗ (w) =
v∈N + (w)
X v∈N + (w)
µ(v) d− (v) 1 − d (v)
,
the weighted average of the locality centers of its out-neighbors (N + (w)). Sorting the updated positions yields the new ordering. These update rules naturally give rise to the iterative reordering procedure shown in Alg. 1. Fig. 5 validates this on FS and IL (Tab. 1) using OG and RCMK initializations; each point is an intermediate ordering from Alg. 1. Figs. 5a and 5b show that C(π) strongly predicts inference time and cold reloads (R2 =0.83–0.96), while Fig. 5c shows that the neighborhood dispersion objective J closely tracks C(π) (R2 =0.96–0.99). Thus, reducing J lowers span, I/O, and runtime.
3.9
Generalizability of Taurus
Taurus extends beyond GCN/GIN-style aggregation to other message-passing GNNs: SAGEConv and GATConv, and to sampled inference. SAGEConv SAGEConv can double hot-store footprint because its update concatenates a vertex’s aggregated neighborhood representation with its own embedding before applying a shared transformation: h(ℓ+1) = σ W [hagg,v ∥h(ℓ) v v ] , where hagg,v = Aggregateu∈N − (v) h(ℓ) u 11
Table 1: Graph datasets [11, 23, 32] used in experiments. Papers Friendster MAG-Cites IGB-Large IGB-Full Abbr. # Vertices # Edges Feat. Dim # Classes Top. Size (GiB) Feat. Size (GiB)
PA 111M 1.7B 128 172 14 54
FS 65M 3.6B 1024 64 28 251
MA 121M 1.4B 768 153 12 350
IL 100M 1.2B 1024 19 10 382
IF 269M 4B 1024 19 32 514 (FP16)
(ℓ)
Under broadcast execution, hv is immediately available when v is streamed, whereas hagg,v is ready only after all propagated messages arrive; retaining both vectors doubles per-vertex state (as seen in ATLAS [1]). Taurus eliminates this overhead by partitioning the neural′ ′ network matrix W ∈ Rd ×2d into W = [Wagg Wself ], where Wagg , Wself ∈ Rd ×d . Therefore, (ℓ) (ℓ) by linearity, W [hagg ∥hv ] = Wagg hagg + Wself hv , which is exactly equivalent to the original formulation. Thus, the hot store maintains only hagg ; once complete, Taurus writes Wagg hagg (ℓ) as a run file. A subsequent sequential pass computes Wself hv , sums the projected terms, applies bias/activation if present and emits the final embedding. This halves hot-store state and reduces eviction–reload cycles while preserving SAGEConv semantics, at the cost of one extra sequential pass. GATConv GATConv is more challenging because attention depends on both source and destination embeddings and is edge-specific: X h(ℓ+1) = σ αuv W hu , where v u∈N − (v)
αuv = softmax(euv ), and euv = LeakyReLU a⊤ [W hu ∥W hv ] Taurus supports GATConv using multiple sequential passes, avoiding O(|E|d′ ) edge-feature materialization. First, it materializes transformed features h′v = W hv for all vertices, reducing later passes from dimension d to d′ . This is followed by a topology-only pass over the CSR that computes edge attention scores and applies a neighborhood-wise softmax to obtain the attention coefficients αuv . Finally, Taurus streams the transformed embeddings again and aggregates using these attention weights. This preserves GAT semantics while trading random gathers and high-dimensional edge state for additional sequential passes. For multihead GAT, the same procedure applies per head, with proportional storage and pass costs. Sampling-based Inference Some applications require deterministic full-neighborhood inference, while others tolerate fanout-sampled approximation. Taurus supports both. For fanout k, Taurus constructs a layer-specific edge mask by sampling up to k incoming edges per vertex and suppressing the rest. Sampling reduces propagated messages, but not the sequential scan of vertex features/embeddings needed to compute updated representations. Thus, Taurus retains the same sequential read pattern while broadcasting only over sampled edges, reducing hot-store residency, eviction traffic and runtime.
12
I. GCN Inf. Time (s)
Extrapolated 16 0.1M 8 10K
10K
4
1K
256 1M
2560.1M
256 1M
256
64 0.1M
64 10K
64 0.1M
64
16 10K
16
16 10K
16 4
1K
4
1K
4 0.1K
4
1 10 16 1M
1 0.1K 256 1M
1 10 2560.1M
1 0.1K 256 1M
1 256
0.1M
8 0.1M
64 0.1M
64 10K
64 0.1M
64
10K
4
10K
16 10K
16
16 10K
16
1K
2
1K
4
1K
4 0.1K
4
4
0.1K 1M
1 0.1K 8 1M
1 0.1K 256 1M
1 10 2560.1M
1 0.1K 256 1M
1 256
0.1M
4 0.1M
64 0.1M
64 10K
64 0.1M
64
10K
2
10K
16 10K
16
16 10K
16
1K
4
4
1K
4 0.1K
4
0.5 0.1K 8 1M
1 0.1K 256 1M
1 10 2560.1M
1 0.1K 256 1M
1 256
0.1M
4 0.1M
64 0.1M
64 10K
64 0.1M
64
10K
2
10K
16 10K
16
16 10K
16
1K
1
1K
4
4 0.1K
4
4
0.1K
1K
L1 L2 L1 L2 0.5 0.1K L1 L2 L1 L2 1 0.1K L1 L2 L1 L2 1 GX DG DI TA GX DG DI TA GX DG DI TA A. PA B. FS C. MA
10
1K
Rel. Speedup
1
IV. GAT Inf. Time (s)
1K
0.1K 1M
1K
1K
Rel. Speedup
III. SAGE Inf. Time (s)
1K
1K
Rel. Speedup
2 0.1K
II. GIN Inf. Time (s)
1K
0.1K 1M
1K
1K
Rel. Speedup
1M 0.1M
L1 L2 L1 L2 1 0.1K L1 L2 L1 L2 1 GX DG DI TA GX DG DI TA D. IL E. IF
Figure 6: Sampled inference performance (fanout=10 per layer) of TAurus vs. GineX, DGL, and DGI for 2-layer GCN, SAGE, GIN and GAT on PA, FS, MA, IL and IF. Bars show total inference time (s, log scale); markers show speedup over Taurus. Hatched bars denote extrapolated runtimes.
4
Evaluation
We evaluate Taurus along four axes: end-to-end performance against OOC baselines, support for exact and sampled inference, impact of runtime components and resource usage.
4.1
Experimental Setup
We evaluate Taurus using four 2-layer vertex-classification GNNs: GraphConv (GCN) [2], SAGEConv (SAGE) [3], GINConv (GIN) [29], and single-head GATConv (GAT) [30], all with hidden dimension 128. Experiments use five open-source citation and social-network datasets (Tab. 1), with feature stores ranging from 54 GiB for Papers [23] to 514 GiB for IGB-Full [11]. While IGB-Full uses FP16, rest use FP32 precision. We report both exact fullneighborhood inference and fanout-sampled inference. For exact full-neighborhood inference, Taurus matches the output of an in-memory layer-wise DGL [27] implementation on PA using identical weights and precision, with mean per-vertex max absolute error 8 × 10−5 and mean relative error 2.8 × 10−6 . Unless otherwise stated, experiments run on a single workstation with a 12-core AMD Ryzen 9 9900X CPU (4.4 GHz), 128 GiB RAM, an NVIDIA RTX 5090 GPU with 32 GiB VRAM, a 2 TiB Samsung 990 PRO SSD, and Ubuntu 24.04.3 LTS. We clear the OS page cache before each run.
4.2
Taurus Implementation and Baselines
Taurus is implemented in Python with NumPy v2.0 and PyTorch v2.8; the graph reader and embedding writer are C++ PyTorch extensions compiled with ninja. Unless varied in 13
ablations, we use 8 MiB chunks, 256 MiB graduation buffers, queue size 20, 50 GiB hot store for PA/MA/IL, 100 GiB for FS/IF, and a 16 GiB GPU store after reserving VRAM for write buffers, intermediate tensors, and CUDA context. Framework overhead is 5–7 GiB. We compare Taurus (TA) with three OOC baselines: vertex-wise Ginex [18] (GX) and DGL GraphBolt OnDisk2 (DG), and layer-wise DGI [28] (DI). Ginex targets OOC GNN training with disk-resident neighbor and feature caches. We adapt it to inference by disabling backpropagation and using its default superbatch/batch sizes (2500–3300/1000), with 90 GiB feature cache and 10 GiB neighbor cache on our 128 GiB system. DGL GraphBolt OnDisk uses DGL’s OnDiskDataset abstraction for vertex-wise OOC execution. We use RCMK reordering [31] and batch size 32K to improve locality, retaining other defaults. DGI is a layer-wise inference framework using dynamic batching, RCMK ordering [31], and NumPy mmap files for features and CSC indices. We use the paper’s default settings. We report layer-wise results for DI and TA, and end-to-end results for GX and DG. Runs are capped at 4 h; incomplete runs are linearly extrapolated from completed vertex ranges/chunks, per layer for DI/TA and end-to-end for GX/DG. If a DI layer times out, later layers are measured with dummy inputs of matching dimensions and added to the extrapolated incomplete layer.
4.3
Comparison with Baselines
Fig. 6 compares sampled inference (fanout=10 ) across four 2-layer GNNs and five datasets. We use fanout 10 because smaller fanouts understate message-propagation costs, while larger fanouts make most baselines exceed the time budget. GX/DG bars are end-to-end times; DI/TA report per-layer times, and the total inference time is the sum across layers. Performance Improvements over Vertex-Wise Baselines TA outperforms GX by ≈ 40×, 62×, 57× and 140× on OOC FS, MA, IL and IF, respectively, because GX incurs repeated cache construction, sampled-batch materialization and feature I/O; Belady’s policy reduces but does not eliminate the latter. TA similarly outperforms DG by ≈ 60×, 57×, 90× and 96×, as DG gathers disk-resident features per batch. In contrast, TA streams embeddings sequentially: GCN/GIN require one read/write pass per layer, while SAGE/GAT add sequential passes but avoid repeated random gathers. Accordingly, DG and GX read up to ≈ 50× and ≈ 108× more data, respectively, in the evaluated workloads. Performance Improvements over Layer-Wise DGI Combined across both layers, TA achieves average speedups of 15×, 8.2×, and 7.2× over DI on FS, MA, and IL, respectively (Fig. 6, cols. B–D). On IF, DI’s first layer exceeded the 4 h cap; we estimate total DI time by extrapolating that layer and measuring later layers with matching dummy inputs. TA outperforms DI by ≈ 25× averaged across all models, while completing inference on the 514 GiB dataset in under 30 min. DI remains stronger than vertex-wise baselines because layerwise execution avoids redundant neighborhood expansion, but still reads up to 10× more data than TA. Performance When Features Fit in Memory On PA (Fig. 6, col. A), the graph topology and features fit entirely in memory, allowing the OS page cache to eliminate most disk I/O. Even so, TA outperforms the vertex-wise GX baseline by 3.7–8.1×, and DG by ≈ 1.9× for GCN and GIN, because it avoids repeated reads. Compared to the layer-wise DI baseline, TA remains up to ≈ 1.3× faster for GCN and GIN, but is ≈ 0.6× slower for SAGE 2
https://www.dgl.ai/dgl_docs/generated/dgl.graphbolt.OnDiskDataset.html
14
(a) FS/GCN2
C( ) (×1015)
Inf. Time (min)
120 RC OG C( ) 4 90 3 60 2 30 1 0 0 1 2 3 4 5 6 7 8 910 0 Iteration
C( ) (×1015)
Inf. Time (min)
90 RC OG C( ) 3 60 2 30 1 0 0 1 2 3 4 5 6 7 8 910 0 Iteration
(b) IL/GCN2
Figure 7: Impact of Taurus reordering bootstrap and iterations showing end-to-end inference time (bars, left Y axis) and span cost C(π) (markers, right Y axis). RC and OG mean RCMK and original ordering bootstraps, respectively. and GAT. Since I/O costs are mitigated by the OS page cache, the additional passes in TA to support SAGE and GAT dominate execution, while DI executes these operators directly in memory. Reduction in Layer Execution Time Across the FS, MA, IL, and IF datasets, Taurus’s execution time drops sharply from L1 to L2 by ≈ 4.8× across all models, because the first layer projects high-dimensional inputs (768/1024) to a 128-dimensional hidden space. Consequently, subsequent layers stream significantly less data from disk, reducing both I/O overhead and computational costs due to reduced dimensionality. PA exhibits the opposite trend because its input and hidden dimensions are 128, while the output dimension is 172. Hence, L2 is computationally costlier than L1 , resulting in an average ≈ 17% increase in execution time across all models. Single vs. Multi Pass GNNs SAGE and GAT require multiple sequential passes, increasing runtime relative to GCN/GIN. Compared to GCN, SAGE/GAT are 93%/112% slower on PA and 24–64%/0.7–17% slower on OOC datasets, consistent with reading 103%/51% more data. The increase is largest on PA because the graph and its features largely fit in memory; hence, the layer decomposition and repeated scans become costlier than out-of-the-box GNN operators. Lastly, the relatively small increase in GAT execution time compared to SAGE, despite the additional passes, is expected because GAT first projects the embeddings from d to d′ (d′ < d), allowing all subsequent passes to operate on the lower-dimensional representation. (§ 3.9).
4.4
Impact of Taurus Ordering
We compare OG, RD, RCMK (RC), ATLAS (AT) and Taurus (TA) orderings for exact 2layer GCN inference on IL and FS, using 50/80 GiB hot stores, a 10 GiB GPU store and TA eviction (Fig. 8). We also evaluate convergence under OG and RC bootstraps (Fig. 7). Convergence of Taurus reordering Fig. 7 shows the span objective C(π) (§ 3.8) decreases monotonically over 10 iterations across both OG and RC bootstraps for both datasets. This confirms that the iterative TA updates consistently improve the ordering. On FS (Fig. 7a), C(π) decreases by 36–52% (markers, right Y axis), translating to a 25–38% reduction in endto-end inference time (bars, left Y axis). The gains are larger on IL (Fig. 7b), where C(π) decreases by 65–75%, yielding a 51–64% reduction in runtime. RC is the stronger bootstrap, giving 6–60% lower final runtimes after convergence (≈ 50 vs. 53 min for FS; ≈ 22 vs. 55
15
(a) Inf. reloads.
time, I/O time, and #
NS GS HS EV Reuse %
Reuse %
100 100 80 90 60 80 40 70 20 60 0 OGRDRC AT TA OGRDRC AT TA 50 FS IL
Dst Mix (%)
Reloads (M, )
Inf. Time (min)
120 FS IL I/O Reloads320 90 240 60 160 30 80 0 OG RD RC AT TA 0
(b) Dst. distribution and hot-store reuse.
Figure 8: Impact of Taurus reordering on execution behavior. (a) E2E inference time (bars, left Y axis), I/O time (hatched bars, left Y axis), and # reloads (markers, right Y axis); (b) Avg. destination-state split (stacked bars, left Y axis) and hot-store reuse (markers, right Y axis) for Layer 1. min for IL), so Taurus uses RC by default. Lastly, over 84–90% of the total span reduction is achieved within the first 4–5 iterations, indicating rapid convergence. Based on this observation, we limit the number of iterations to 5 in Taurus. Impact on System Performance Fig. 8a compares the impact of different vertex orderings on system performance. TA uses RC as bootstrap with 5 iterations of refinement. By reducing span, TA consistently exhibits the lowest end-to-end inference time, reducing runtime by ≈ 31–48% over OG, RD, AT, and RC on FS to 56 min (green bars, left Y axis) and by ≈ 57–75% on IL to 21 min (blue bars, left Y axis). The improved span substantially lowers SSD traffic, reducing I/O time (both reads/writes) by 3.6–5.9× on FS (hatched green bars, left Y axis) and 73–158× on IL (hatched blue bars, left Y axis), while also decreasing eviction-reload cycles by 4.3× and 126×, respectively, on average (markers, right Y axis). On IL, TA reduces I/O overhead to near-negligible levels under the same GPU and hot-store budgets as the other orderings. Hot Store Utilization and Reuse Finally, we quantify this benefit at the chunk level in Fig. 8b, showing the average split of destination vertices for a chunk into HS (IN_BUFFER), GS (ON_GPU), NS (NOT_STARTED), and EV (EVICTED) states. TA increases the mean % of destination vertices in the hot store (HS) by ≈ 3.7% and 17% (green stack, left Y axis), with a corresponding decrease in the % in the cold store (EV) by ≈ 4% and 17% (orange stack, left Y axis) for FS and IL, respectively. Consequently, the hot-store reuse, i.e., the fraction of active destinations (vertices with ≥ 1 messages received) already resident in the HS hot store ( HS+EV ), increases from ≈ 90–93% to 98% on FS and from ≈ 67% to 99.8% on IL (markers, right Y axis). By keeping nearly all active states resident until graduation, TA reduces I/O thrashing.
4.5
Impact of Taurus Eviction Policy
Fig. 9 compares the proposed minimum pending-messages eviction policy (TA) against random (RD), first-in-first-out (FIFO), and least-recently-used (LRU) eviction for exact 2-layer GCN inference under TA ordering, using 80 GiB and 40 GiB hot stores for FS and IL, respectively, and a 10 GiB GPU store. Fig. 9a shows that prioritizing vertices for eviction with the fewest pending messages substantially improves execution time. Compared to RD and LRU, TA reduces inference time by ≈ 19–43% on FS and ≈ 20–26% on IL (solid bars, left Y axis). LRU performs poorly because it repeatedly evicts vertices far from completion, while FIFO improves upon LRU 16
40 20 30 15 RD 20 10 FIFO LRU 10 5 TA 01 2 4 8 163264 01 2 4 8 16 32 Reload Count Reload Count FS IL
# Verts x (M)
Inf. Time (min)
Reloads (M, )
FS I/O 120 IL Reloads 320 90 240 60 160 30 80 0 RD FIFOLRU TA 0
(a) Inf. time, I/O time, and # reloads.
(b) CDF of # reloads.
Figure 9: Impact of Taurus eviction on execution behavior. (a) E2E inference time (solid bars, left Y axis), I/O time (hatched bars, left Y axis), and # reloads (markers, right Y axis); (b) Cumulative reload distribution of the number of unique vertices reloaded from the cold store.
Inf. Time (min)
120 OG AT RC TA 90 60 30 0 30 40 50 60 70 Hot Store Size (GiB)
Inf. Time (min)
120 OG AT RC TA 90 60 30 0 60 70 80 90 100 Hot Store Size (GiB) (a) FS/GCN2
(b) IL/GCN2
Figure 10: Impact of hot-store capacity on E2E inference time under different vertex orderings for (a) FS and (b) IL. TA achieves identical performance with substantially lower hot-store budgets than OG, RC, and AT (dashed red line). by providing a minimum residency period after admission. TA outperforms both by evicting vertices closest to completion, allowing them to graduate soon after reload and minimizing eviction–reload cycles. Consequently, TA reduces I/O time by ≈ 31–80% to 9 min and 2.7 min (hatched bars, left Y axis), while reducing hot-store reloads by ≈ 1.4–6× to 40M and 21M on FS and IL, respectively (markers, right Y axis). This anti-thrashing behavior is directly validated in Fig. 9b showing the distribution of destination reload counts. Because TA selects near-complete victims, most vertices see only 2–4 reloads across FS and IL, versus 12–17 for RD. LRU has a long tail of 34–84 reloads because inactive but far-from-complete vertices are repeatedly swapped. By intelligently selecting victims that will rapidly graduate upon reload, TA minimizes repeated cold-store accesses, drastically reducing overall execution time.
4.6
Impact of Hot Store Size
Figs. 10 and 11 study the impact of varying hot-store sizes for exact 2-layer GCN inference on FS and IL, respectively, using a fixed 10 GiB GPU store. Larger hot stores only partially compensate for poor ordering. Figs. 10a and 10b show that increasing hot-store capacity reduces inference times by allowing more partially aggregated vertices to remain resident, thereby relieving eviction pressure, across all reorder strategies. However, the steep performance curves for OG, RC, and AT indicate that these orderings are fundamentally memory-starved. In contrast, TA exhibits a notably flatter curve, showing that it is consistently the least sensitive to hot-store capacity. While increasing the RAM budget improves TA by only ≈ 28–40%, the baseline orderings see improvements of ≈ 50–65% across FS and IL. Consequently, TA effectively decouples performance from memory scale, 17
(a) FS/TA/GCN2
Reloads (M)
I/O Time (min)
16 80 I/O Reloads 60 12 8 40 4 20 0 30 40 50 60 70 0 Hot Store Size (GiB)
Reloads (M)
I/O Time (min)
20 160 I/O Reloads 120 15 10 80 5 40 0 60 70 80 90 100 0 Hot Store Size (GiB)
(b) IL/TA/GCN2
Figure 11: Impact of hot-store capacity on TA. SSD I/O time (bars, left Y axis) and coldstore reloads (markers, right Y axis) under TA reorder for (a) FS and (b) IL.
Messages (%)
Reloads (M, )
Inf. Time (min)
(a) Inf. time and # reloads.
Verts (%, )
80 FS IL Verts (%) 8 60 6 40 4 20 2 0 0 2 4 6 8 10 12 14 16 0 GPU Store Size (GiB)
80 FS IL I/O Reloads 80 60 60 40 40 20 20 0 0 2 4 6 8 10 12 14 16 0 GPU Store Size (GiB)
(b) % of messages and vertices.
Figure 12: Impact of GPU-store capacity on Taurus. (a) E2E inference time (solid bars, left Y axis), I/O time (hatched bars, left Y axis) with corresponding cold-store reloads (markers, right Y axis); (b) % of messages aggregated (bars, left Y axis) and % of vertices (markers, right Y axis) on the GPU-store. requiring substantially less RAM to achieve identical throughput. For example, on FS, TA at just 60 GiB matches the peak performance of AT/OG at ≈ 90 GiB. Similarly, on IL, TA at 30 GiB matches RC at ≈ 60 GiB (dashed red line). By minimizing the average lifespan of active vertices (§ 3.8), TA ensures that partial states graduate quickly, drastically reducing the structural need for massive buffer capacity. TA approaches peak throughput with modest memory. Figs. 10 and 11 show that the 28– 40% reduction in end-to-end inference time achieved by TA is primarily due to lower SSD I/O overhead. As the hot-store size increases, reloads decrease from ≈ 120M to 10M on FS and from ≈ 75M to nearly zero on IL (circles, right Y axis), reducing I/O time from ≈ 17 min to 2 min on FS and from ≈ 12 min to negligible on IL (bars, left Y axis). For IL, reloads are almost completely eliminated beyond a 50 GiB hot store. At this point, the hot store can hold nearly all active vertex states, eliminating eviction–reload cycles. Consequently, SSD I/O becomes negligible, inference becomes compute-bound, and further increasing the hot-store size provides little additional speedup, demonstrating that Taurus can efficiently process a ≈ 400 GiB graph with only a 50 GiB hot-store budget (plus framework overheads).
4.7
Impact of GPU Store Size
Fig. 12 evaluates the impact of GPU-store capacity on 2-layer full-graph GCN inference for FS and IL using 80 GiB and 40 GiB hot stores, respectively. Pinning hub vertices reduces host-memory pressure. Without a GPU store, all intermediate aggregation states are managed by the hot store, causing high-degree vertices to occupy hot-store slots for extended durations and repeatedly evict lower-degree vertices. At 0 GiB, this results in ≈ 64M and 46M cold-store reloads (Figs. 12a, markers, right Y axis) on FS and IL, respectively, with end-to-end inference times of ≈ 66 min and 35 min and I/O times 18
Inf. Time (min)
Reloads (M, )
120 AT TA 160 90 120 60 80 30 40 0 100 100 100 0 +16 +0 HS (+GS) Size (GiB)
Reloads (M, ) Inf. Time (min)
I/O Default TA300 120 FS IL Reloads AT TA 90 225 60 150 30 75 0 50 60 70 80 90 100 50 50 0 +16 +0 HS (+GS) Size (GiB)
Figure 13: Taurus (TA) vs. ATLAS (AT). E2E time (bars, left Y axis), I/O time (hatched ), and # reloads (circles, right Y axis) for varying HS (+GS) budgets on IL (blue, left) and FS (green, right). of ≈ 12 and 9 min (Figs. 12a, bars, left Y axis). Increasing the GPU store permanently pins the highest-degree vertices, allowing them to directly absorb a disproportionate fraction of destination messages. At 16 GiB, the GPU store caches just 6.4% and 4.2% of vertices on FS and IL (Figs. 12b, markers, right Y axis), yet serves ≈ 51% and 39% of destination messages (Figs. 12b, bars, left Y axis), reducing reloads to 29.6M and 12.9M and inference time to 47 and 22 min. GPU-store gains track SSD I/O and saturate quickly. Fig. 12a shows that inference time closely tracks SSD I/O time as GPU-store capacity increases. At 16 GiB, inference time reduces by ≈ 27% (FS) and 36% (IL), while SSD I/O time decreases by ≈ 57% and 79%, respectively. The marginal benefit diminishes as capacity grows, since the highest in-degree vertices are pinned first. On IL, the first 2 GiB reduces inference time by ≈ 4.5 min (36% of the total 12.3 min gain), while the final 2 GiB saves only ≈ 0.5 min. FS exhibits a similar pattern (≈ 5 vs. 0.9 min). Consequently, the first 8 GiB captures 64–72% of the total improvement on both datasets.
4.8
Taurus Reordering Cost
Taurus’s reordering is a one-time preprocessing step and can be reused across all subsequent inference runs. The RCMK bootstrap takes ≈ 150–720s (largest on IF), after which the optimization loop itself is lightweight, with 5 iterations taking 56–177s (5–15%) of total time. The dominant cost is feature relabeling (125–1053s, 23–80%), an I/O-bound operation, while topology relabeling adds only 35–111s. Overall, preprocessing completes in 9–30 min across all datasets.
4.9
Comparison with ATLAS
Fig. 13 compares TA and AT for exact 2-layer GCN inference on IL and FS across hotstore (HS) and GPU-store (GS) budgets. On IL, the default TA setup (50+16 GiB) completes in ≈ 20 min, 1.65× faster than AT’s best at 100 GiB HS (33 min) while using two-thirds of total memory, and 5.7× faster than AT at an equal 50 GiB HS. Even without a GPU store (50+0), matching AT’s HS-only setup, Taurus completes in 24.4 min, 4.6× faster than AT at the same memory budget, isolating the benefit of Taurus’s topology-aware reordering. TA’s I/O is negligible (hatched, 1.3 min vs. 7–81 min for AT) as reloads fall from 265M to 11.5M (circles). On FS at an equal 100 GiB HS, TA is 1.9× faster (43.9 vs. 83.1 min), and 1.5× faster even with no GS. We omit additional AT memory configurations on FS due to prohibitively long runtimes. Overall, TA outperforms AT with lower memory and substantially lower I/O.
19
80 60 40 20 00
GPU Util. (%)
(a) FS/GCN2 CPU and Memory
(b) FS/GCN2 GPU and VRAM
0.8 0.6 0.4 0.2 0.0 600 1200 1800 2400 3000 Time (s)
Write IO (GB/s)
Read IO (GB/s)
2.0 1.5 1.0 0.5 0.0 0
32 24 16 8 0 600 1200 1800 2400 3000 Time (s)
GPU Mem. (GiB)
120 90 60 30 0 600 1200 1800 2400 3000 Time (s)
Memory (GiB)
CPU Util. (%)
2400 1800 1200 600 00
(c) FS/GCN2 Read/Write Bandwidth
Figure 14: Resource util. for a 2-layer GCN on the FS dataset. (a) CPU (blue, left Y axis) and memory (red, right Y axis) over time; (b) GPU util (sky blue, left Y axis) and memory (brick, right Y axis) usage over time; (c) SSD read (green, left Y axis) and write (orange, right Y axis) bandwidth over time.
4.10
Resource Utilization of Taurus
Fig. 14 shows CPU utilization, RAM usage, GPU utilization and memory, and SSD read/write throughput, sampled every second during exact 2-layer GCN inference on FS with a 100 GiB hot store and a 16 GiB GPU store. Taurus sustains average 1150% CPU utilization (Fig. 14a, blue lines, left Y axis) on our 12-core (2-way hyperthreaded) machine (§ 4.1) while average SSD reads remain low at 160 MiB/s (Fig. 14c, green lines, left Y axis), hence the streaming readers and writers do not stall compute, with an initial read burst (≈ 1.8 GiB/s) at layer starts reflecting filling of the read queue. CPU utilization drops (1200 → 1050%) while SSD writes increase between 900–1350 s (black arrows) due to vertex evictions. Resident memory (Fig. 14a, red lines, right Y axis) tracks the configured budget ≈ 105 GiB in Layer 1 (100 GiB hot store plus ≈ 5 GiB framework overhead) and ≈ 38 GiB in Layer 2 once 128-dim embeddings replace 1024-dim input features with the step at ≈ 2700 s marking the layer transition (red arrow ). GPU utilization (Fig. 14b, sky blue lines, left Y axis) is intermittent (peak 72%) because CPU-side aggregation and memory management dominate. GPU memory (Fig. 14b, brick lines, right Y axis) stays around ≈ 31 GiB reflecting the 16 GiB GPU store, 8 GiB write buffers, plus transient buffers. The periodic write spikes (peak 4.3 GiB/s) mark graduated vertices being flushed to the SSD.
5
Related Works
GNNs support large-scale recommendation and traffic applications [6,8,33]. As graphs evolve, systems periodically retrain models and refresh vertex representations to address drift in topology, features and labels [34, 35]. Taurus targets this refresh phase, accelerating billionscale inference on a single OOC machine.
20
5.1
Large-scale GNN Training and Inference
While many systems target scalable GNN training [8, 36, 37], comparatively little attention has been paid to large-scale inference. Systems such as AliGraph [8] and DistDGL [37] are designed for distributed mini-batch training, partitioning graphs across machines and optimizing sampling and gradient synchronization. Likewise, popular GNN libraries such as PyG [38] and DGL [27] provide extensive support for mini-batch training and neighborhood sampling, but are not optimized for single-machine full-graph inference. Consequently, billion-scale inference typically relies on distributed deployments (e.g., DistDGL-style clusters), incurring significant infrastructure overheads. Exact full-neighborhood inference creates larger working sets than sampled training, while sampled inference still suffers repeated feature movement under gather-based execution. Recent systems target GNN inference more directly. InferTurbo [24] uses GAS-style execution for scalable inference, but relies on large MapReduce clusters. DGI [28] converts training code to layer-wise inference with dynamic batching and graph reordering; however, its mmap-based OOC design remains gather-driven and incurs high read amplification beyond RAM. We therefore use DGI as the strongest layer-wise baseline.
5.2
Out-of-core GNN Training Systems
Several systems use SSD-backed single machines for large-scale OOC GNN training [14–18, 21, 39], primarily optimizing data movement, storage layout, and caching. DiskGNN [14] reduces read amplification by precomputing computation graphs, packing features contiguously on disk, and pipelining execution over a multi-level feature store, while MariusGNN [39] avoids precomputation by partitioning the graph and restricting aggregation to a vertex’s memory-resident neighbors within a partition. These strategies do not directly fit exact fullneighborhood inference: all vertices require outputs, precomputing every computation graph is prohibitive, and partition-restricted aggregation can drop cross-partition neighbors [14]. Capsule [16] uses partitioning, pruning, and optimized loading to fit training subgraphs in GPU memory, while Ginex [18] separates sampling from feature gathering to enable Beladyoptimal feature caching. These optimizations target sampled training; inference over all vertices, exact or sampled, still stresses feature movement and exposes random gather overheads. In contrast, Taurus makes sequential source broadcasts the primary abstraction for inference, then manages the resulting partial states through topology-aware ordering and tiered GPU–RAM–SSD storage.
6
Conclusion and Discussion
We introduced Taurus, a single-machine out-of-core system for billion-scale GNN inference, supporting both exact full-neighborhood and fanout-sampled execution. By replacing destination-centric gathers with pipelined source broadcasts, Taurus converts repeated random reads into sequential scans. Coupled with a tiered GPU-CPU-SSD hierarchy, topologyaware vertex reordering, and a pending-message eviction policy, it minimizes hot-store residency and I/O thrashing. Despite these gains, Taurus currently assumes static graph snapshots. We plan to address this by extending the broadcast execution model to support continuous, incremental updates on evolving topologies, and we also plan to explore multi-GPU scaling to accelerate compute-bound phases and support even larger embedding dimensions.
21
Acknowledgments The authors thank Roopkatha Banerjee and other members of the DREAM:Lab, Indian Institute of Science, for their assistance and insightful feedback. The authors used AI-assisted tools only for language refinement and clarity. All ideas, methods, experiments, results and conclusions are the authors’ own, and the final manuscript was reviewed and approved by the authors.
References [1] P. Naman and Y. Simmhan, “Atlas: Efficient out-of-core inference for billion-scale graph neural networks,” in ACM Symposium on High-Performance Parallel and Distributed Computing (HPDC), 2026. [2] T. N. Kipf and M. Welling, “Semi-supervised classification with graph convolutional networks,” in International Conference on Learning Representations (ICLR), 2017. [3] W. L. Hamilton, R. Ying, and J. Leskovec, “Inductive representation learning on large graphs,” in International Conference on Neural Information Processing Systems (NIPS), 2017. [4] Y. Dou, Z. Liu, L. Sun, Y. Deng, H. Peng, and P. S. Yu, “Enhancing graph neural network-based fraud detectors against camouflaged fraudsters,” in ACM International Conference on Information & Knowledge Management (CIKM), 2020. [5] B. Dasari, T. S. Dhiraj, G. Jambhrunkar, T. Kailasam, C. Vikram, S. Singla, P. Naman, and Y. Simmhan, “Billion-scale fintech analytics: Scalable data management and anomaly detection at npci,” in IEEE International Conference on Data Engineering (ICDE), 2026. [6] A. Derrow-Pinion, J. She, D. Wong, O. Lange, T. Hester, L. Perez, M. Nunkesser, S. Lee, X. Guo, B. Wiltshire et al., “Eta prediction with graph neural networks in google maps,” in ACM International Conference on Information & Knowledge Management (CIKM), 2021. [7] A. Sharma, P. Naman, R. Banerjee, P. Pansari, S. Gawali, M. Arya, S. Chandra, A. Josephraj, R. Ramesh, P. Rathore et al., “Scaling real-time traffic analytics on edgecloud fabrics for city-scale camera networks,” in TCSC SCALE Challenge, IEEE CCGRID Workshops, 2026. [8] H. Yang, “Aligraph: A comprehensive graph neural network platform,” in ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD), 2019. [9] D. Wu, Z. Li, and T. Mitra, “Inkstream: Instantaneous gnn inference on dynamic graphs via incremental update,” in IEEE International Parallel and Distributed Processing Symposium (IPDPS), 2025. [10] P. Naman and Y. Simmhan, “Ripple: Scalable incremental gnn inferencing on large streaming graphs,” in IEEE International Conference on Distributed Computing Systems (ICDCS), 2025. [11] A. Khatua, V. S. Mailthody, B. Taleka, T. Ma, X. Song, and W.-m. Hwu, “Igb: Addressing the gaps in labeling, features, heterogeneity, and size of public graph datasets for deep learning research,” in ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD), 2023. 22
[12] T. Kaler, A. Iliopoulos, P. Murzynowski, T. Schardl, C. E. Leiserson, and J. Chen, “Communication-efficient graph neural networks with probabilistic neighborhood expansion analysis and caching,” Proceedings of Machine Learning and Systems (MLSys), 2023. [13] S. Chen, X. Song, V. Theodore, and H. Liu, “Deal: distributed end-to-end gnn inference for all nodes,” arXiv preprint arXiv:2503.02960, 2025. [14] R. Liu, Y. Wang, X. Yan, H. Jiang, Z. Cai, M. Wang, B. Tang, and J. Li, “Diskgnn: Bridging i/o efficiency and model accuracy for out-of-core gnn training,” in Proceedings of the ACM on Management of Data (SIGMOD), 2025. [15] C. Su, H. Zhang, H. Zhao, W. Shen, B. Ai, Y. Li, K. Bian, and B. Cui, “Caliex: A disk-based large-scale gnn training system with joint design of caching and execution,” in International Conference on Data Engineering (ICDE), 2025. [16] Y. Xiang, Z. Ding, R. Guo, S. Wang, X. Xie, and S. K. Zhou, “Capsule: an out-of-core training mechanism for colossal gnns,” in Proceedings of the ACM on Management of Data (SIGMOD), 2025. [17] Z. Sheng, W. Zhang, Y. Tao, and B. Cui, “Outre: An out-of-core de-redundancy gnn training framework for massive graphs within a single machine,” in Proceedings of the VLDB Endowment, 2024. [18] Y. Park, S. Min, and J. W. Lee, “Ginex: Ssd-enabled billion-scale graph neural network training on a single machine via provably optimal in-memory caching,” in Proceedings of the VLDB Endowment, 2022. [19] K. Vora, “Lumos: Dependency-driven disk-based graph processing,” in USENIX Annual Technical Conference (USENIX ATC), 2019. [20] A. Roy, I. Mihailovic, and W. Zwaenepoel, “X-stream: Edge-centric graph processing using streaming partitions,” in ACM Symposium on Operating Systems Principles (SOSP), 2013. [21] J. Sun, M. Sun, Z. Zhang, Z. Shi, J. Xie, Z. Yang, J. Zhang, Z. Wang, and F. Wu, “Hyperion: Co-optimizing ssd access and gpu computation for cost-efficient gnn training,” in IEEE International Conference on Data Engineering (ICDE), 2025. [22] J. B. Park, V. S. Mailthody, Z. Qureshi, and W.-m. Hwu, “Accelerating sampling and aggregation operations in gnn frameworks with gpu initiated direct storage accesses,” in Proceedings of the VLDB Endowment, 2024. [23] W. Hu, M. Fey, M. Zitnik, Y. Dong, H. Ren, B. Liu, M. Catasta, and J. Leskovec, “Open graph benchmark: Datasets for machine learning on graphs,” in International Conference on Neural Information Processing Systems (NIPS), 2020. [24] D. Zhang, X. Song, Z. Hu, Y. Li, M. Tao, B. Hu, L. Wang, Z. Zhang, and J. Zhou, “Inferturbo: A scalable system for boosting full-graph inference of graph neural network over huge graphs,” in IEEE International Conference on Data Engineering (ICDE), 2023. [25] T. Kaler, N. Stathas, A. Ouyang, A.-S. Iliopoulos, T. Schardl, C. E. Leiserson, and J. Chen, “Accelerating training and inference of graph neural networks with fast sampling and pipelining,” Proceedings of Machine Learning and Systems (MLSys), 2022.
23
[26] P. Naman and Y. Simmhan, “A gpu is all you need: Rethinking distributed and out-ofcore gnn training,” in IEEE International Conference on High Performance Computing, Data and Analytics Workshop (HiPCW), 2025. [27] M. Y. Wang, “Deep graph library: Towards efficient and scalable deep learning on graphs,” in ICLR Workshop on Representation Learning on Graphs and Manifolds, 2019. [28] P. Yin, X. Yan, J. Zhou, Q. Fu, Z. Cai, J. Cheng, B. Tang, and M. Wang, “Dgi: An easy and efficient framework for gnn model evaluation,” in ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD), 2023. [29] K. Xu, W. Hu, J. Leskovec, and S. Jegelka, “How powerful are graph neural networks?” in International Conference on Learning Representations (ICLR), 2019. [30] P. Veličković, G. Cucurull, A. Casanova, A. Romero, P. Liò, and Y. Bengio, “Graph Attention Networks,” in International Conference on Learning Representations (ICLR), 2018. [31] W.-M. Chan and A. George, “A linear time implementation of the reverse cuthill-mckee algorithm,” BIT Numerical Mathematics, 1980. [32] J. Leskovec, “Stanford network analysis project,” https://snap.stanford.edu/, 2020. [33] R. Ying, R. He, K. Chen, P. Eksombatchai, W. L. Hamilton, and J. Leskovec, “Graph convolutional neural networks for web-scale recommender systems,” in ACM SIGKDD International Conference on Knowledge Discovery and Data Mining (KDD), 2018. [34] E. Rossi, B. Chamberlain, F. Frasca, D. Eynard, F. Monti, and M. Bronstein, “Temporal graph networks for deep learning on dynamic graphs,” arXiv preprint arXiv:2006.10637, 2020. [35] Y. Xia, Z. Zhang, H. Wang, D. Yang, X. Zhou, and D. Cheng, “Redundancy-free highperformance dynamic gnn training with hierarchical pipeline parallelism,” in International Symposium on High-Performance Parallel and Distributed Computing (HPDC), 2023. [36] S. Gandhi and A. P. Iyer, “P3: Distributed deep graph learning at scale,” in USENIX Symposium on Operating Systems Design and Implementation (OSDI), 2021. [37] D. Zheng, C. Ma, M. Wang, J. Zhou, Q. Su, X. Song, Q. Gan, Z. Zhang, and G. Karypis, “Distdgl: Distributed graph neural network training for billion-scale graphs,” in IEEE/ACM Workshop on Irregular Applications: Architectures and Algorithms (IA3), 2020. [38] M. Fey and J. E. Lenssen, “Fast graph representation learning with pytorch geometric,” arXiv preprint arXiv:1903.02428, 2019. [39] R. Waleffe, J. Mohoney, T. Rekatsinas, and S. Venkataraman, “Mariusgnn: Resourceefficient out-of-core training of graph neural networks,” in European Conference on Computer Systems (EuroSys), 2023.
24