Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

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.

share|improve this question
    
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 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. –  popffabrik Mar 8 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 at 0:46

2 Answers 2

up vote 4 down vote accepted

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.

share|improve this answer
    
Is the List@@ that does nothing there for a reason? ;-} Outer to the rescue again! Nice. +1 –  rasher Mar 8 at 11:43
1  
@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 at 11:55
    
ah, thanks for link! –  rasher Mar 8 at 11:58
    
Thanks! I ticked this answer because it gives me the best performance (and the code is so short). –  popffabrik Mar 10 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....

share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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