Optimization of a graph generator algorithm

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.

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]}]


• I think by preallocating to the max size you are taking a big hit on memory, and maybe also on speed. If you preallocate to say 2^Ceiling[Log[2*m]], and then double all sizes whenever any node hits that many vertices, you will both save on the memory/speed and still be able to use Compile. Hugely faster than my Append-based neighbor lists. – Daniel Lichtblau Oct 31 '17 at 14:23
• @DanielLichtblau tried that. Worse top-level time. I think the issue might actually be array unpacking. Need to look into it. – b3m2a1 Oct 31 '17 at 22:52
• @DanielLichtblau The main issue seems to be some form of pre-factor, based on a basic FindFit comparison of model2, model3, and model3c. So I guess the Appends are being handled pretty well internally. Either that or there simply aren't enough of them occurring at the n I generated data for. – b3m2a1 Oct 31 '17 at 23:17

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.

• You could escape using Append by preallocating an array for each element of g and then using vdeg as an index tracker for insertion. – b3m2a1 Oct 31 '17 at 6:26
• @b3m2a1 I also had that idea, somewhat belatedly. Every time we hit the size limit, double the preallocation, pad with zeros, and keep an array of the max lengths (or just use powers of 2). It would give a considerable speed improvement, at modest cost in memory. Also might lend itself to use of Compile, if we are willing to make all the neighbor lists as large as the largest one. Instead f attempting any of this (which would be real work) I think I'll just take the expedient route, and upvote yours. – Daniel Lichtblau Oct 31 '17 at 14:19
• Thanks for your answer, it really did help me and saved me a lot of time. I just have one doubt, because I'm trying to add some modifications in the code: is there a way to calculate the LocalClusteringCoefficient at every step inside the Do loop without having to generate the graph and apply the built-in function to return the results? – Augusto Bellinaso Nov 7 '17 at 20:34
• Yes it can be done, but also it can be done with the method by @b3m2a1 and that approach is much faster. Either way will require a bit of coding though. At every step you have to determine the clustering for the new vertex (easy in case 1, harder in case 2), and also update for all neighbors of the new vertex. – Daniel Lichtblau Nov 8 '17 at 0:02
• The first lines of your code (before the 'Do' loop), generate a 'CompleteGraph'. I'm not sure if it is possible, but could I use a similar structure to generate a 'CycleGraph' instead? – Augusto Bellinaso Nov 8 '17 at 16:01