Optimization of a graph generator algorithm

Question

I'm working with complex networks and I have this algorithm that I'm using to generate my networks, which I use to compare with empirical results. The idea of the algorithm is kind simple: I started from the Barabasi-Albert algorithm, which I used as a base for my model. The model starts with a complete graph with (m+1) nodes, where m is the same parameter used in the BA Model (number of links that the new node will have), and then I create a new node. This node will be add to the network following one of the two processes:

• duplication: where I choose a node from the network to duplicate, and the new node will have the same neighbours as the parent node; the probability used to choose this node is similar to the preferential attachment probability defined by Barabasi.
• preferential attachment: new node will follow the rule defined by Barabasi, choosing one of the most connected nodes in the network and connecting to m different nodes.

model[n_, m_, r_] /; n >= 3 :=     
  Module[{g = CompleteGraph[m + 1], vc, vl, el, el2, node, newnode},
    Do[
      vc = VertexCount[g];
      vl = VertexList[g];
      el = EdgeList[g];
      If[RandomReal[] <= p, 
        node = RandomChoice[(VertexDegree[g]^(r)) -> vl];
        newnode = AdjacencyList[g, node];
        el = ParallelMap[UndirectedEdge[vc + 1, #] &, newnode];
        ,
        el = ParallelMap[
               UndirectedEdge[vc + 1, #] &, 
               RandomSample[VertexDegree[g] -> vl, m]
             ]
      ];
      g = Graph[Join[EdgeList[g], el]];
      , 
      {n - m - 1}
    ];
    g
  ]

In the code, n is the size of the network, m is the BA model parameter, r is the expoent I'm using in the probability to choose the node to duplicate and p is tha probability that the new node will either duplicate an existing one or follow the preferential attachment rule.

Since I'm working with large networks (around 10000 nodes) and with high values of p (from 0.5 to 0.9), I'm trying to find a way, if there is, to write this code more efficient, because it takes almost a week to genereta a network with n=11711, m=5, r=1 and p =0.5.

Answer 1

Here's an adaptation of Daniel Lichtblau's version:

modelGraph[g_] :=
  Graph[
   Union[
    Map[Sort,
     Flatten[MapIndexed[Thread[{#2[[1]], #1}] &, g], 1]
     ]
    ]
   ];

model3Core[n_, m_, r_, p_] :=

  Module[{g, vdeg, el, el2, node, nodelist, newnode,
    sampPadding = Table[0, {i, n}]
    },
   g =
    With[{base = Table[0, {i, m + n}]},
     Table[base, n]
     ];
   Do[
    With[{l = Delete[Range[m + 1], j]},
     g[[j, 1]] = Length@l;
     g[[j, 2 ;; Length@l + 1]] = l;
     ],
    {j, m + 1}
    ];
   vdeg = Table[0, {i, n + m}];
   Do[vdeg[[i]] = m, {i, m + 1}];
   Do[
    If[RandomReal[] <= p,
      node = RandomChoice[vdeg[[1 ;; j - 1]]^r -> Range[j - 1]];
      g[[j]] = g[[node]];
      vdeg[[j]] = vdeg[[node]];
      Do[
       With[{k = g[[node, k]]},
        vdeg[[k]] += 1;
        g[[k, vdeg[[k]] ]] = j;
        ],
       {k, vdeg[[node]]}
       ];,
      nodelist = RandomSample[vdeg[[1 ;; j - 1]] -> Range[j - 1], m];
      g[[j]] = Join[nodelist, sampPadding];
      vdeg[[j]] = m;
      Do[
       vdeg[[k]] += 1;
       g[[k, vdeg[[k]] ]] = j;,
       {k, nodelist}
       ]
      ];,
    {j, m + 2, n}
    ];
   Append[g, vdeg]
   ];
model3[n_, m_, r_, p_] /; n >= 3 :=

 MapThread[Take, {Most@#, Drop[Last@#, -m]}] &@
  model3Core[n, m, r, p]

All I did really was remove the Appends

If we redefine model2 to remove the Graph call we can directly compare them:

Map[
 First@AbsoluteTiming@model2[#, 5, 1, .5] &,
 {100, 1000, 5000}
 ]

{0.009174, 0.405928, 12.1925}

Map[
 First@AbsoluteTiming@model3[#, 5, 1, .5] &,
 {100, 1000, 5000}
 ]

{0.010279, 0.627182, 20.7768}

And... somehow I made it worse.

But my version can be Compiled:

model3Comp =
  With[{t = Extract[DownValues[model3Core], {1, 2}, Unevaluated]},
   Compile @@ Hold[{
      {n, _Integer},
      {m, _Integer},
      {r, _Real},
      {p, _Real}
      },
     t
     ]
   ];
model3c[n_, m_, r_, p_] /; n >= 3 :=

 MapThread[Take, {Most@#, Drop[Last@#, -m]}] &@
  model3Comp[n, m, r, p]

Map[
 First@AbsoluteTiming@model3c[#, 5, 1, .5] &,
 {100, 1000, 5000}
 ]

{0.0009, 0.048616, 1.32722}

And that brings us into a good domain to work with.

Now we can do the real call in reasonable time:

model3c[11711, 5, 1, .5] // AbsoluteTiming // First

10.096

And just to check that it lines up on a smaller system:

GraphicsRow@
 Map[Rasterize@*modelGraph, {model2[50, 5, 1, .5], 
   model3c[50, 5, 1, .5]}]

Answer 2

One issue is the creation and part extraction of new graphs at every step. Another is that parallelizing might not be useful in this situation since the operations are fast, there are not many to be done in each step, and the data transfer cost might thus dwarf the actual iteration.

The code below seems to produce identical results. At least they look similar. I have not made IsomorphicGraphQ behave well enough to be sure though.

model2[n_, m_, r_, p_] /; n >= 3 := Module[
  {g, vdeg, el, el2, node, newnode},
  g = ConstantArray[{}, n];
  g[[1 ;; m + 1]] = Table[Delete[Range[m + 1], j], {j, m + 1}];
  vdeg = ConstantArray[0, n];
  vdeg[[1 ;; m + 1]] = m;
  Do[
   If[RandomReal[] <= p,
     node = RandomChoice[vdeg[[1 ;; j - 1]]^r -> Range[j - 1]];
     g[[j]] = g[[node]];
     vdeg[[j]] = vdeg[[node]];
     Do[vdeg[[k]] += 1; g[[k]] = Append[g[[k]], j], {k, g[[node]]}];
     ,
     node = RandomSample[vdeg[[1 ;; j - 1]] -> Range[j - 1], m];
     g[[j]] = node;
     vdeg[[j]] = m;
     Do[vdeg[[k]] += 1; g[[k]] = Append[g[[k]], j], {k, node}]
     ];
   , {j, m + 2, n}];
  Graph[Union[
    Map[Sort, Flatten[MapIndexed[Thread[{#2[[1]], #1}] &, g], 1]]]]
  ]

At n=1000 the original takes 87 seconds on my laptop. The variant above takes 2.5 seconds. I still am not thrilled with the quadratic+ complexity caused (mostly) by the use of Append. Might be room for more improvement there.