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.

  • $\begingroup$ Have you tried the solutions from here? mathematica.stackexchange.com/q/608/12 You'll notice that Mma has builtin algorithms for this. $\endgroup$
    – Szabolcs
    Mar 8, 2014 at 0:17
  • $\begingroup$ 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. $\endgroup$
    – matimo2
    Mar 8, 2014 at 0:30
  • $\begingroup$ I am sorry, you are correct. I realized this, but then I had to leave suddenly so I couldn't remove my comment. $\endgroup$
    – Szabolcs
    Mar 8, 2014 at 0:46

2 Answers 2


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.

  • $\begingroup$ Is the List@@ that does nothing there for a reason? ;-} Outer to the rescue again! Nice. +1 $\endgroup$
    – ciao
    Mar 8, 2014 at 11:43
  • 1
    $\begingroup$ @rasher, the List@@ is there to unpack at level 1, so that a and b will remain packed. See this for discussion. $\endgroup$ Mar 8, 2014 at 11:55
  • $\begingroup$ ah, thanks for link! $\endgroup$
    – ciao
    Mar 8, 2014 at 11:58
  • $\begingroup$ Thanks! I ticked this answer because it gives me the best performance (and the code is so short). $\endgroup$
    – matimo2
    Mar 10, 2014 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....


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.