Understanding performance: graph connected components

I need to analyze directed graphs with 10M edges, 1M vertices, and 300K strongly connected components, so that the largest one contains 400K vertices.

I read some explanations of Leonid Shifrin here, here and here. Although I don't quite understand why his code merge is so blazingly fast, I learned two things: no recursions only iterations, linked lists are cool. Trying to mimic his approach, I prepared two function for finding strongly connected components. The first one is Kosaraju's algorithm:

ClearAll[readEdge, SowDFS, DFS, kosaraju]

(* Sow as a postvisit vertex function *)
SowDFS[q[v_?newQ, rest_]] := SowDFS[newQ[v] = False; adjIn[v] /. q[] -> q[Sow[v], rest]];
SowDFS[q[v_, rest_]] := SowDFS[rest];

(* deep first scan *)
DFS[q[v_?newQ, rest_]] := DFS[Sow[v, tag]; newQ[v] = False; adjOut[v] /. q[] -> rest];
DFS[q[v_, rest_]] := DFS[rest];

(* finding connected components *)
kosaraju[graphList_] := Block[{adjOut, adjIn, order, newQ, q, tag = 1, $IterationLimit = Infinity}, (* construction of adjacency lists *) adjOut[v_] = q[]; adjIn[v_] = q[]; Scan[readEdge, graphList]; SetAttributes[q, HoldAllComplete]; (* the first scan: topological sort of the reverse graph *) newQ[v_] = True; order = Reverse[Reap[Scan[SowDFS[q[#, q[]]] &, DeleteDuplicates@Flatten[graphList]]][[2, 1]]]; (* the second scan: finding components *) Clear[adjIn, newQ]; newQ[v_] = True; Last@Reap[Scan[(tag++; DFS[q[#, q[]]]) &, order], _, #2 &] ] The second one is based on Tarjan's algorithm, as it is described in Wikipedia. I only adapted it in non-recursive way: (* reading edge *) readOut[{a_, b_}] := (adj[a] = q[b, adj[a]]); (* test to start component *) strongQ[v_] := index[v] == link[v]; (* pop from stack *) pop[v_, q[v_, rest_]] := (inQ[Sow[v, v]] = False; rest); pop[v_, q[a_, rest_]] := pop[v, inQ[Sow[a, v]] = False; rest]; (* deep first scan with two stacks *) biDFS[q[v_?newQ, rest_], stack_] := biDFS[ (* previsit function *) newQ[v] = False; index[v] = link[v] = ++idx; {p[v], root} = {root, v}; adj[v] /. q[] -> q[h[v], rest], inQ[v] = True; q[v, stack]]; biDFS[q[h[v_?strongQ], rest_], stack_] := biDFS[rest, pop[v, stack]]; biDFS[q[v_?inQ, rest_], stack_] := biDFS[link[root] = Min[link[root], index[v]]; rest, stack]; biDFS[q[v_, rest_], stack_] := biDFS[rest, stack] (* postvisit vertex function *) h[a_] := (root = p[a]; link[root] = Min[link[root], link[a]]; a); (* start scan *) start[v_?newQ] := Block[{p, root = v, inQ, link, index, idx = 0}, (* p is for parent *)p[a_] := a; inQ[a_] = False; biDFS[q[v, q[]], q[]]] (* finding connected components *) tarjan[graphList_] := Block[{adj, newQ, q,$IterationLimit = Infinity},
(* construction of adjacency out-lists *)
adj[v_] = q[]; Scan[readOut, graphList]; SetAttributes[q, HoldAllComplete];
(* finding components *)
newQ[v_] = True;
Last@Reap[Scan[start, DeleteDuplicates@Flatten[graphList]], _, #2 &]
]

I try not to use WM build in functions for graph analysis in order to have a clear comparison.

To test the performance I chose Google+ graphs from Stanford Large Network Dataset Collection. The renamed version of one of them is available here.

SetDirectory[NotebookDirectory[]];
graphList = Import["edges.dat"];

This graph contains 5126172 edges but only 16405 vertices. The structure of SCC is not completely trivial: a single connected component with 11064 vertices, another one — 11 vertices, one more — 5 vertices, three ones — 3 vertices each, ten components with 2 vertices, and 5296 vertices which are not strongly connected.

On the one hand, WM somehow requires a lot of time to form the graph:

result2 = ConnectedComponents[
Graph[DirectedEdge @@@ graphList]]; // AbsoluteTiming
(* {139.523577, Null} *)

and only half of a second to find the strongly connected components.

On the other hand, reading the graph takes 29 seconds in my function tarjan, and 12 seconds to find the components, thus 42 second in total:

result1 = tarjan[graphList]; // AbsoluteTiming
(* {42.207118, Null} *)

Hence, my question is it possible to speed up the tarjan function so that it will be at least 10 times slower than the built-in function ConnectedComponents?

My nb-file is available here.

• The merge function is fast because it is type-specialized on numerical lists, on which it gets compiled. No time to look at your code in detail right now, but any top-level code with some sort of iteration will likely be much slower than the built-in function. You may want to check out my answer here, for a semi-compiled implementation of connected components which might rival a built-in one, and actually the entire thread to which it belongs, for some extra context / ideas. – Leonid Shifrin Jun 10 '14 at 19:48
• If you don't want to restrict yourself to Mathematica you can take a look here and here – Sektor Jun 10 '14 at 20:04
• @LeonidShifrin Thanks for the quick response. Sure, I will read the thread, but at a glance it is about connected components of an undirected graph, that is a bit easier. – Grisha Kirilin Jun 10 '14 at 20:22
• @Sektor looks interesting, but I have a problem for one time only. I use it to learn something new about Mathematica. – Grisha Kirilin Jun 10 '14 at 20:26
• There might be methods of interest at this link. My guess is mostly they will be too slow though (except maybe the built-in ConnectedComponents). – Daniel Lichtblau Jun 11 '14 at 13:47