Betweenness centrality quantifies the importance of a node as an intermediary on shortest paths. For a node v, the standard definition is:
where σ(s, t) is the total number of shortest paths from s to t and σ(s, t | v) is the number of those paths that pass through v. Naïvely computing BC requires enumerating all-pairs shortest paths in O(V²(V + E)) time. Brandes' algorithm [1] reduces this to O(V(V + E)) by computing, for each source s, a recursive pair-dependency score via backward accumulation over the BFS tree:
The final BC score is then:
where the factor of ½ corrects for double-counting in undirected graphs. This two-phase BFS approach is the standard algorithm for exact BC computation on unweighted graphs.
For each source vertex s, Brandes executes two phases:
Per-source working arrays used by Brandes' algorithm. Each must be allocated, initialized, and reset for every one of the n source iterations.
| Structure | Type | Purpose |
|---|---|---|
| adj[v] | vector<vector<int>> | Adjacency list (0-based node IDs) |
| sigma[v] | double[n] | Shortest-path count from s to v |
| dist[v] | int[n] | BFS distance from s |
| pred[w] | vector<int>[n] | Predecessors of w on shortest paths |
| delta[v] | double[n] | Accumulated dependency for v |
| stk[] | int[n] | BFS-visited nodes in discovery order |
| bc[v] | double[n] | Final betweenness centrality score |
Parallelizing Brandes efficiently on this graph is non-trivial for several compounding reasons.
Memory pressure scales with thread count. Each thread requires private copies of five working arrays of size n, totaling ≈750 KB per thread. At 8 threads this footprint reaches ≈6 MB — at or above the L3 cache capacity on most consumer CPUs — meaning threads actively evict each other's working data. Simply adding threads increases total work rather than reducing it, as we confirm experimentally.
Severe and irregular load imbalance. Per-source execution times span a 255× range (21 µs to 5.4 ms), driven by the PPI graph's extreme degree skew. A small number of hub proteins (BRCA1, TSG101, CASP8) are connected to thousands of other nodes and participate in millions of shortest paths; sources adjacent to these hubs require orders of magnitude more work than isolated proteins. Static work distribution strategies risk leaving threads idle while others process the heaviest sources.
Shared state and accumulation races. The final bc[] array is updated by every source iteration for every node on its shortest paths. Concurrent writes require either fine-grained synchronization — which serializes updates on the high-degree hub nodes that are touched by all threads simultaneously — or private per-thread accumulators with a reduction phase, which introduces its own memory and bandwidth overhead.
Fine-grained parallelism is bottlenecked by the graph's small diameter. An alternative to source-level parallelism is to parallelize within each source's BFS. However, the PPI graph has a diameter of only ≈6 levels, yielding very small per-level frontiers. Any fine-grained parallelism strategy must contend with this: the work per parallel section is smaller than the synchronization overhead required to manage it.
We use the BioGRID v5.0.256 human PPI dataset (Homo sapiens, organism ID 9606) [2]. The raw TAB3 file contains interactions across multiple experimental system types and species. We filter to physical interactions only, discarding genetic interactions, non-human partners, self-loops, and duplicate undirected edges. After deduplication the graph has 18,727 nodes; 2,113 have zero betweenness (isolated proteins or nodes in small components). Entrez Gene IDs are remapped to contiguous integers 0…n−1 at load time.
The serial implementation pre-allocates all working arrays once and reuses them across all 18,727 source iterations. The BFS uses a hand-rolled queue (vector<int> with head/tail pointers) to avoid repeated allocation overhead. Before each source, all n entries of sigma, dist, delta, and pred are reset via a full O(n) sweep. This reset is necessary for correctness but constitutes a significant fraction of total runtime, as each sweep touches 750 KB of memory cold. Serial total time: 52.61 s.
The most direct approach to parallelism is to distribute the source loop across threads. The key design decision is how to handle the shared bc[] accumulator: concurrent writes from multiple threads require either synchronization or replication. We chose replication: each thread maintains a private bc_local[tid][·] array and accumulates into it exclusively. After the parallel region, the per-thread arrays are reduced serially into the final bc[]. This eliminates all synchronization during the parallel region at the cost of O(nT) memory and a serial reduction phase.
We use schedule(static), assigning each thread a contiguous block of ≈n/T sources. Each thread independently allocates its own private working arrays.
#pragma omp parallel num_threads(T) { int tid = omp_get_thread_num(); // private working arrays, size n each (~750 KB per thread) ... #pragma omp for schedule(static) for (int s = 0; s < n; s++) { // BFS forward + backward on private arrays // write results to bc_local[tid] } } // serial reduction: bc[v] = sum_t bc_local[t][v]
To address load imbalance from the 255× spread in per-source times, we replace static scheduling with schedule(dynamic, 1), handing sources out one at a time from a shared work queue. This ensures that no thread sits idle while another processes a heavy hub-adjacent source. The tradeoff is synchronization overhead: every source assignment requires an atomic increment on a shared counter, adding 18,727 atomic operations to the critical path.
Rather than replicating bc[] per thread, this variant uses a single shared array with #pragma omp atomic for every backward-accumulation write. This eliminates the O(nT) memory overhead of private accumulators and the serial reduction phase, but introduces fine-grained synchronization on bc[]. The critical problem is that hub proteins appear on millions of shortest paths, so all threads simultaneously write to the same small set of bc[v] entries during the backward pass. This serializes updates on the hottest entries and generates false sharing across nearby cache lines.
An alternative parallelization strategy is to keep the source loop serial and instead parallelize within each source's accumulation phase. We perform BFS serially (building a levels[] array grouping nodes by BFS depth), then process the backward pass level-by-level in parallel:
for (int lvl = levels.size()-1; lvl >= 1; lvl--) { #pragma omp parallel for num_threads(T) schedule(static) for (int idx = 0; idx < levels[lvl].size(); idx++) { int w = levels[lvl][idx]; for (int v : pred[w]) { double contrib = (sigma[v]/sigma[w]) * (1.0 + delta[w]); #pragma omp atomic delta[v] += contrib; } } for (int w : levels[lvl]) if (w != s) bc[w] += delta[w]; }
Nodes at the same BFS level write to predecessors one level shallower, so a barrier between levels ensures correctness. Atomic adds to delta[] are necessary because multiple nodes at depth ℓ+1 can share a predecessor at depth ℓ, creating write conflicts. The PPI graph's small diameter limits the available parallelism at each level, as discussed in Section 3.
Profiling the serial baseline revealed that the O(n) reset before each source represents a disproportionate share of total runtime for sparse sources. Critically, this reset precedes BFS, meaning it touches memory cold — dirtying cache lines before any useful work has been done. For isolated nodes and small components, the reset processes all 18,727 entries despite BFS visiting only a handful.
The lazy-reset optimization defers the reset to after the accumulation phase, resetting only the nodes actually visited during BFS:
// Track visited nodes during BFS visited.push_back(w); // whenever a new node is discovered // After accumulation: reset only visited nodes (cache-warm) for (int v : visited) { pred[v].clear(); dist[v] = -1; sigma[v] = 0.0; delta[v] = 0.0; }
This reduces the reset from O(n) to O(|visited|) and ensures the reset memory is warm in cache (having just been written during BFS and accumulation), reducing effective memory bandwidth.
The final variant combines static inter-source parallelism (Variant 1) with lazy reset (Variant 5). Each thread maintains its own private visited[] vector alongside its other working arrays. This directly addresses the memory bandwidth saturation diagnosed in Variant 1: by reducing each thread's active working-set size, the combined footprint of 8 simultaneous threads is more likely to fit within the shared L3 cache, reducing the mutual eviction that inflates per-source time.
All experiments ran on the GHC machines, specifically machine 78. Wall time for parallel variants is approximated as (sum of per-source times) / T; actual wall time equals or slightly exceeds this estimate depending on load balance.
intra-bfs is 26% slower than serial; omp-atomic performs comparably to omp-static despite its synchronization overhead.Figure 1 shows static and dynamic scheduling performing nearly identically (3.57× vs. 3.48× at 8 threads). This is a meaningful result: despite the 255× spread in per-source times, dynamic scheduling offers no benefit. Figure 5 explains why. The apparent imbalance is almost entirely due to 6 isolated nodes; 99.97% of sources fall within the narrow 2–4 ms band of the large connected component. Under static scheduling these sources are distributed evenly across threads, producing minimal thread imbalance. Dynamic scheduling pays 18,727 atomic work-counter increments to solve an imbalance problem that does not materially exist, resulting in slightly worse performance than static.
Figure 3 shows the intra-BFS variant at 0.79× serial speed — a 26% slowdown. This is the most instructive failure of the project. The core difficulty is that the PPI graph's ≈6-level diameter means each parallel section contains only ~800 nodes per thread. Thread-launch and barrier overhead (5–20 µs per level) dominates the actual computation time, so threads spend more time synchronizing than computing.
Compounding this, the atomic adds to delta[] required for correctness — multiple nodes at depth ℓ+1 can share a predecessor at depth ℓ, creating write conflicts — serialize the hottest updates within each parallel section. The result is that intra-BFS parallelism on this graph pays synchronization costs twice: once at the level boundary via implicit barriers, and once within each level via atomic contention. Graphs with large diameter and wide frontiers would not suffer this problem, but the PPI graph's biological small-world structure makes this approach fundamentally unsuitable here.
The omp-atomic variant (3.61×) performs comparably to private-array static scheduling (3.57×). This near-tie occurs because the two approaches exchange one bottleneck for another of roughly equal cost: the private-array variant pays an O(nT) serial reduction phase after the parallel region, while the atomic variant pays contention on hub nodes during the backward pass. At n = 18,727 and T = 8, these costs approximately cancel. At larger n or denser graphs, hub contention would scale worse than the linear reduction, making the private-array approach strictly preferable.
The central finding of this project is that memory-bandwidth saturation, not load imbalance or synchronization, is the binding constraint on scaling. Figure 4 provides direct evidence: if threads operated without mutual interference, the sum of per-source times would remain flat at 52.61 s regardless of thread count. Instead, at 8 threads it reaches 118 s — a 2.24× inflation over serial. Each thread's work takes 2.24× longer simply because 8 threads are running simultaneously.
The culprit is the O(n) array reset that begins every source iteration. Each reset sweeps ≈750 KB of working arrays cold, before any BFS work has touched those cache lines. At 8 threads, eight simultaneous cold sweeps continuously evict each other's data from shared L3 cache, saturating the memory bus. The consequence is visible in Figure 2: parallel efficiency drops from 91% at 2 threads to 45% at 8 threads. Each additional thread adds compute capacity but also adds to the collective memory pressure, and beyond 4 threads the pressure grows faster than the compute benefit.
The lazy-reset optimization (Variants 5 and 6) directly targets this bottleneck. By deferring the reset to after accumulation and limiting it to visited nodes, threads reset only warm cache lines rather than sweeping cold memory. For the 2,113 isolated-node sources, the reset shrinks from O(18,727) to O(1); for connected sources the cache-warm effect reduces effective bandwidth regardless of visited count.
The observed speedup plateaus near 3.5× at 8 threads, well short of the ideal 8×. The primary cause is memory-bandwidth saturation — each thread must reset ≈750 KB of working arrays before every source iteration, and at 8 threads these concurrent cold sweeps saturate the shared memory bus and continuously evict each other's cache lines from L3. A secondary contributor is load imbalance: even though 99.97% of sources complete within the narrow 2–4 ms band, the few hub-adjacent sources that take significantly longer cause the slowest thread to finish after the others have idled, wasting available compute. Together, these two effects — memory contention growing super-linearly with thread count, and a long tail of heavy sources — cap practical speedup at roughly T/2 for this graph and problem size.
| Rank | Node ID | Entrez ID | Gene | BC Score |
|---|---|---|---|---|
| 1 | 238 | 672 | BRCA1 | 991,816 |
| 2 | 141 | 7182 | TSG101 | 986,648 |
| 3 | 236 | 330 | BIRC2 | 927,766 |
| 4 | 872 | 821 | CASP8 | 884,295 |
| 5 | 336 | 9531 | CFLAR | 836,256 |
| 6 | 579 | 3329 | HSPD1 | 833,134 |
| 7 | 6 | 6472 | SHMT1 | 828,714 |
| 8 | 43 | 2885 | GRB2 | 807,849 |
BRCA1 and TSG101 ranking first is biologically meaningful: both are hub proteins with roles in DNA repair, cell signaling, and cancer pathways, connected to a large and diverse set of other proteins. These nodes serve as critical bridges between otherwise-distant parts of the network — precisely the condition that maximizes betweenness centrality — and their high degree is also what makes them the primary bottleneck for parallel accumulation, as all threads contend on their bc[] entries simultaneously.
A multicore CPU is the appropriate target for this problem. BFS on irregular adjacency lists causes severe warp divergence on GPUs due to non-uniform neighbor counts and unpredictable memory access patterns. The strict sequential dependencies within each source's backward pass are difficult to map to SIMD execution, and the high-precision double arithmetic required by BC scores exceeding 900,000 is better served by CPU floating-point units. A GPU would become competitive for graphs with millions of nodes, where wide BFS frontiers could fill warps efficiently and amortize launch overhead; for an 18,727-node graph with 6-level BFS the granularity is too small to benefit from GPU-scale parallelism.
std::chrono timing and per-source profiling instrumentation for the serial baseline