# Improving performance in the generation of random (directed acyclic) graphs

I need to generate large (directed, acyclic, transitive) graphs. My code in Mathematica gets very slow for large numbers of vertices. I've tried a number of approaches and I can't find any significant improvements on speed.

The graphs are generated as follows. You are given a list aList of size nS of randomly generated points (tuples of reals) and a boolean "edge function" f acting on pairs of such points. You then generate a graph with nS vertices labelled by integers i = 1,...,nS whose edges are found as follows: f(aList[[i]],aList[[j]]) == True $\iff$ i->j. The properties of f are always such that the resulting graph is a transitive DAG.

Here an explicit example. My list of points (2-tuples in this case) is

nS = 1000;
aList = Sort[Table[{RandomReal[], RandomReal[]}, {nS}]];


and my edge function is

f[i_, j_] := aList[[j, 1]] - aList[[i, 1]] > Abs[aList[[j, 2]] - aList[[i, 2]]]


Now I want to generate my graph. I have thought of either finding all the Rules (edges) using Reap/Sow or computing an adjacency matrix, but all the approaches I have tried are very slow for large nS. Below my different approaches:

r1 =       Reap@Do[If[f[i, j], Sow[Rule[j, i]]],    {i, nS}, {j, i}]; // Timing
r2 = ParallelTable[If[f[i, j], Rule[j, i], ## &[]], {i, nS}, {j, i}]; // Timing
r3 = ParallelTable[If[f[i, j], 1., 0.],             {i, nS}, {j, i}]; // Timing

(*
{4.626997, Null}
{2.023294, Null}
{1.397499, Null}
*)


Each of these approaches is very slow for nS = 1000 already. To generate a boolean adjacency matrix this way in C++ takes me ~0.01 seconds instead of ~1.0 second. Is there some way to improve performance in Mathematica? Perhaps using Compile? I have tried tweaking different things without much success.

I haven't given any background/motivation but if you'd like me to I'd be happy to explain.

• Have you tried the solutions from here? mathematica.stackexchange.com/q/608/12 You'll notice that Mma has builtin algorithms for this. – Szabolcs Mar 8 '14 at 0:17
• Thanks, yes I have seen the solution in the link but I'm afraid I don't think it helps in my case. The solution addresses one particular random process for generating DAGs (from uniform random adjacency matrices). I have also seen the built-in algorithms for some particular types of random DAGs, but no algorithms that generate DAGs from underlying sets and an "edge function" as described above. The types of random DAGs that I am generating are qualitatively different. – matimo2 Mar 8 '14 at 0:30
• I am sorry, you are correct. I realized this, but then I had to leave suddenly so I couldn't remove my comment. – Szabolcs Mar 8 '14 at 0:46

For the specific example in the question, this generates the adjacency matrix in about 20ms for nS = 1000:

{a, b} = List @@ Transpose[aList];
adjMat2 = UnitStep[Abs@Outer[Plus, -b, b] ~Subtract~ Outer[Plus, -a, a]] ~BitXor~ 1


If you have lots of different edge functions to play with and don't want to work out a different optimised algorithm for each one, I suspect using Compile will be the easiest route.

• Is the List@@ that does nothing there for a reason? ;-} Outer to the rescue again! Nice. +1 – ciao Mar 8 '14 at 11:43
• @rasher, the List@@ is there to unpack at level 1, so that a and b will remain packed. See this for discussion. – Simon Woods Mar 8 '14 at 11:55
• ah, thanks for link! – ciao Mar 8 '14 at 11:58
• Thanks! I ticked this answer because it gives me the best performance (and the code is so short). – matimo2 Mar 10 '14 at 15:44

Here's a start, about fifteen times faster in my limited tests on a case of 2000... (aList, etc. setup same as yours, you can throw the pieces to parallel kernels, of course, particularly the mapping from tuples to rules which is >80% of the time used in below).

i1 = Join @@ Table[ConstantArray[i, i], {i, 1, nS}];
i2 = Join @@ Table[Range[1, j], {j, 1, nS}];
tt = Subtract[Subtract[aList[[i1, 1]], aList[[i2, 1]]],
Abs[Subtract[aList[[i1, 2]], aList[[i2, 2]]]]];
pp = Pick[Transpose[{i2, i1}], Sign[tt], 1];
result = Rule @@@ pp;


If you don't need the tuples as rules, i.e., you just care about getting the adjacency matrix, simply change the last piece to

adjMat= SparseArray[pp -> 1, {nS, nS}];


Should cut time in half again. Break tuples list, spread array creation, BitOr results, who knows how fast....