Kosaraju’s SCC Algorithm, in depth: CSR arrays, a counting-sort transpose, memory budgets and the parallel forward-backward variant

The algorithm in two passes

Two vertices are strongly connected if each can reach the other; that relation partitions the vertices into SCCs, and contracting each SCC gives the condensation, a DAG. Pass 1 runs DFS over the whole graph and appends each vertex to a list when it finishes. The vertex that finishes last lies in a source component of the condensation. Pass 2 walks the transposed graph, starting from vertices in reverse finish order and skipping assigned ones. In the transpose, a source component becomes a sink, so a search started inside it cannot leak out. Everything it reaches is that one SCC. Remove it and repeat.

One detail matters for performance: only pass 1 needs depth-first order. Pass 2 only needs the set of unassigned vertices reachable in the transpose, so any traversal works, including breadth-first search, which is easier to vectorise and parallelise.

CSR: the layout that makes it fast

Adjacency lists of Python objects or of std::vector cost tens of bytes per edge and scatter memory. CSR stores the graph in two arrays: dst holds all edge targets grouped by source, and off holds V+1 offsets so that the out-neighbours of u are dst[off[u]:off[u+1]]. Building it is a counting sort on the source column; building the transpose is the same counting sort on the target column. Both are O(V + E) and need no comparison sort.

import numpy as np

def build_csr(n, src, dst, dtype=np.int64):
    """Counting sort of edges by src. Returns (off, adj) with adj grouped by source."""
    counts = np.bincount(src, minlength=n)
    off = np.zeros(n + 1, dtype=dtype)
    np.cumsum(counts, out=off[1:])
    order = np.argsort(src, kind="stable")   # stable sort on small ints; a manual
    adj = dst[order].astype(np.int32)        # counting-sort scatter is O(E) in C
    return off, adj

def build_graph(n, src, dst):
    fwd = build_csr(n, src, dst)
    rev = build_csr(n, dst, src)              # transpose: same routine, columns swapped
    return fwd, rev

Self-loops and parallel edges need no special handling: neither changes reachability. Vertex ids must be dense integers 0..V-1, so map external keys such as URLs or symbol names to ids first and keep the mapping for the output.

Both passes as array loops

Pass 1 is an iterative DFS. A recursive one overflows the call stack on a long path, which real graphs contain. Each stack frame is a vertex plus the index of its next unexplored edge, held in two parallel integer arrays rather than tuples.

def pass1_postorder(n, off, adj):
    seen = np.zeros(n, dtype=np.bool_)
    order = np.empty(n, dtype=np.int32)       # vertices in finish order
    stack_v = np.empty(n, dtype=np.int32)
    stack_i = np.empty(n, dtype=np.int64)     # next edge index per frame
    k = 0
    for s in range(n):
        if seen[s]:
            continue
        seen[s] = True
        top = 0; stack_v[0] = s; stack_i[0] = off[s]
        while top >= 0:
            v = stack_v[top]; i = stack_i[top]
            if i < off[v + 1]:
                stack_i[top] = i + 1
                w = adj[i]
                if not seen[w]:
                    seen[w] = True
                    top += 1; stack_v[top] = w; stack_i[top] = off[w]
            else:
                order[k] = v; k += 1; top -= 1   # v finishes
    return order

def pass2_label(n, roff, radj, order):
    comp = np.full(n, -1, dtype=np.int32)
    queue = np.empty(n, dtype=np.int32)
    c_id = 0
    for s in order[::-1]:                       # reverse finish order
        if comp[s] != -1:
            continue
        comp[s] = c_id; head = 0; tail = 1; queue[0] = s
        while head < tail:                      # BFS is enough here
            v = queue[head]; head += 1
            for w in radj[roff[v]:roff[v + 1]]:
                if comp[w] == -1:
                    comp[w] = c_id; queue[tail] = w; tail += 1
        c_id += 1
    return comp, c_id

Component ids come out in topological order of the condensation, sources first, because pass 2 starts in a source component each time. That is free information: a build system can schedule components in id order without a separate topological sort. In pure Python with NumPy arrays, both passes over a random million-edge graph take under two seconds; compiled with Numba or written in C, they are limited by memory latency, which the next section budgets.

Memory budget for a billion edges

Take a graph with V = 100 million vertices and E = 1 billion edges, such as a web crawl or a whole-organisation call graph. With 32-bit vertex ids and 64-bit offsets (E exceeds 2³¹):

The transpose is the price of Kosaraju: it doubles edge storage. Tarjan's and Gabow's algorithms need only the forward graph and one DFS, which is why they are the default for large single-machine runs. Kosaraju earns its place when the transpose already exists (many graph stores keep both directions for in-edge queries), when pass 2's freedom to use BFS matters, or when simplicity and auditability matter more than a 2x memory difference. If the edge arrays do not fit in RAM, memory-map them: pass 2 in BFS form streams far better from disk than a DFS does.

Both passes are dominated by random reads of seen or comp for each neighbour. Renumbering vertices so that neighbours have nearby ids, for example in BFS order or by a graph-reordering heuristic, often speeds both passes by a factor of two or more on large graphs, at the cost of a permutation pass up front.

Trimming trivial components

Real graphs contain huge numbers of trivial SCCs: a vertex with no incoming or no outgoing edges cannot be on any cycle, so it is an SCC on its own. Removing it can expose more such vertices, so the rule is applied repeatedly. Trimming is embarrassingly parallel and in web and social graphs often removes a large share of vertices before any search runs.

def trim(n, off, adj, roff, radj, alive):
    """Peel vertices with zero live in-degree or out-degree. Returns peeled ids."""
    indeg = np.array([alive[radj[roff[v]:roff[v+1]]].sum() for v in range(n)])
    outdeg = np.array([alive[adj[off[v]:off[v+1]]].sum() for v in range(n)])
    frontier = [v for v in range(n) if alive[v] and (indeg[v] == 0 or outdeg[v] == 0)]
    peeled = []
    while frontier:
        nxt = []
        for v in frontier:
            if not alive[v]:
                continue
            alive[v] = False; peeled.append(v)
            for w in adj[off[v]:off[v+1]]:          # v no longer feeds w
                if alive[w]:
                    indeg[w] -= 1
                    if indeg[w] == 0: nxt.append(w)
            for u in radj[roff[v]:roff[v+1]]:       # u no longer reaches v
                if alive[u]:
                    outdeg[u] -= 1
                    if outdeg[u] == 0: nxt.append(u)
        frontier = nxt
    return peeled

Forward-backward: SCC without DFS

DFS is inherently sequential, so parallel SCC codes replace it with reachability, which parallelises well as level-synchronous BFS. The forward-backward (FW-BW) algorithm of Fleischer, Hendrickson and Pinar (2000) picks a pivot p, computes F, the set reachable from p, and B, the set that reaches p (reachability in the transpose, the same array Kosaraju builds). Then F ∩ B is exactly p's SCC. Every other SCC lies entirely inside F \ B, B \ F, or the remainder, since two vertices on a common cycle are either both reachable from p or both not, and likewise for reaching p. The three sets are processed independently, in parallel.

def fwbw(vertices, pivot_choice, reach_fwd, reach_bwd, emit):
    tasks = [vertices]
    while tasks:                     # in a real code: a parallel task pool
        S = tasks.pop()
        if not S:
            continue
        p = pivot_choice(S)          # e.g. highest in-degree x out-degree
        F = reach_fwd(p, S)          # BFS restricted to S
        B = reach_bwd(p, S)          # BFS on the transpose restricted to S
        emit(F & B)                  # one SCC
        tasks += [F - B, B - F, S - F - B]

The worst case is poor: on a long chain of singleton SCCs each step removes one vertex, giving O(V(V + E)) work. Production codes combine three steps, as in the Multistep method of Slota, Rajamanickam and Madduri (2014). First, trim. Second, run one FW-BW from a high-degree pivot, which in real graphs usually lands in the single giant SCC and removes it with two parallel BFS runs. Third, handle the many small leftover SCCs with colouring: every vertex repeatedly takes the largest id among itself and its live in-neighbours until nothing changes, so each vertex ends up with the largest id that can reach it. A vertex whose colour equals its own id is a root, and its SCC is the set of same-coloured vertices that can reach it, found by a backward BFS restricted to that colour. Remove those SCCs and repeat.

Worked example: eight vertices, three methods

Take eight vertices and the edges 0→1, 1→2, 2→0, 2→3, 3→4, 4→5, 5→3, 5→6, 7→0. Trim: vertex 7 has in-degree 0 and vertex 6 has out-degree 0, so both peel as singleton SCCs; nothing else drops. Sequential Kosaraju on the remaining six: pass 1 from 0 visits 0, 1, 2, 3, 4, 5 and finishes in the order 5, 4, 3, 2, 1, 0. Pass 2 starts at 0 in the transpose, where 0's in-edges come from 2, and 2's from 1. It collects {0, 1, 2} and stops: the original edge 2→3 is reversed to 3→2, so nothing leads out toward 3. Next unassigned in reverse order is 3, giving {3, 4, 5}. Ids 0 and 1 are in topological order: {0,1,2} feeds {3,4,5}.

FW-BW with pivot 3: F = {3, 4, 5} (6 was trimmed), B = {3, 4, 5, 2, 1, 0}. F ∩ B = {3, 4, 5}. F \ B is empty, B \ F = {0, 1, 2}, and the remainder is empty. Recursing on {0, 1, 2} with pivot 0 gives F = B = {0, 1, 2}. Same answer, and the two recursive calls could have run on different cores.

Choosing an SCC algorithm

Failure modes

• Recursion. A recursive pass 1 crashes on a path of a few hundred thousand vertices.
• Int32 offsets. Offsets overflow silently past 2³¹ edges; use 64-bit offsets.
• Pass 2 on the wrong graph. Searching the forward graph in reverse finish order merges components. A random-graph test against a reference catches it at once.
• Unrestricted reachability in FW-BW. Each BFS must stay inside its current set S, or SCCs leak across tasks and are emitted twice.
• Pivot on a chain. FW-BW without trimming and colouring can degrade to quadratic work on path-like graphs.
• Id mapping lost. Dense ids are internal; keep the map back to external keys and test it.

Testing it

Test the scaled code against a tiny reference. On thousands of random graphs (sparse, dense, with self-loops, with duplicate edges, a single long cycle, a long path), compare the partition from Kosaraju-CSR, from FW-BW and from a simple Tarjan, as sets of frozensets so that id order does not matter. Additionally check that component ids are a valid topological order of the condensation: for every edge u→v, comp[u] <= comp[v].

What to do next

1. Implement build_csr and both passes, then test them against a reference SCC on random graphs.
2. Add the topological-order assertion on component ids and run it in CI.
3. Measure edges per second on a million-edge graph in Python, then compile the same loops with Numba and compare.
4. Add trimming, and report what fraction of vertices it peels on a real graph you own, such as an import or call graph.
5. Implement FW-BW with BFS restricted to a set, and confirm on a long path that it degrades without trimming.
6. Keep learning: Tarjan's one-pass SCC, Gabow's path-based SCC, iterative DFS with explicit stack frames and topological sort for scheduling the condensation.