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

Question

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.

Answer 1

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.

Answer 2

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....