5

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.

matimo2
  • 493
  • 2
  • 9
  • Have you tried the solutions from here? http://mathematica.stackexchange.com/q/608/12 You'll notice that Mma has builtin algorithms for this. – Szabolcs Mar 08 '14 at 00: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 08 '14 at 00: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 08 '14 at 00:46

2 Answers2

7

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

ciao
  • 25,774
  • 2
  • 58
  • 139
5

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.

Simon Woods
  • 84,945
  • 8
  • 175
  • 324