# Optimising a routine which contains calls to a very large list

I have a large graph called ndgraph taken from some network data.

The following code builds a model of the data (the nd.edu graph) by 'growing' a graph in $n$ steps:

nd = Import["nd.txt", "Table"];
For[i = 1, i <= Length[nd], i++,
nd[[i]][[1]] = nd[[i]][[1]] + 1;
nd[[i]][[2]] = nd[[i]][[2]] + 1;
];(*Shifting by one, as a node 0 exists*)
ndgraph = Graph[nd];
edgelist = ConstantArray[0, 1];
Ag = Graph[edgelist];
n = Length[VertexList[ndgraph]];
lengthvertexlist = Length[VertexList[ndgraph]];
vd = VertexDegree[ndgraph];
sumofvertexdegrees = Sum[vd[[i]], {i, 1, lengthvertexlist }];
theta = RandomSample[Table[vd[[i]], {i, 1, lengthvertexlist }]];
norm = Abs[theta[[1]]];
k = RandomVariate[PoissonDistribution[#/2]] & /@ theta;

For[i = 2, i <= n, i++,
j = RandomChoice[
1/norm theta[[Range[1, i - 1]]] -> Range[1, i - 1], {k[[i]]}];
norm = norm + theta[[i]];
];


For some reason, asking for certain values of theta, which has over a million elements from the beginning, is very slow.

Is there some way I can speed this up by removing calls to this list? Perhaps using part isn't a good idea?

• Instead of theta[[Range[1, i - 1]]], try theta[[1;;i-1]]. Apart from that, I don't think that the data requests to theta are the bottleneck. It is probably VertexAdd and EdgeAdd. Actually, you don't need to do that; you can also collect new vertices and edges in a List or in a InternalBag and add them in the end at once. Moreover, I would suggest consider to replace For by Do. – Henrik Schumacher Feb 12 '18 at 9:39
• Btw.: a example graph that highlights the performance issues is missing. – Henrik Schumacher Feb 12 '18 at 9:39
• I’ll add it, hang on – Alexander Kartun-Giles Feb 12 '18 at 9:40
• Minor improvement: norm += theta[[i]]; – Henrik Schumacher Feb 12 '18 at 9:41
• No problem. I’ll try using a replacement to edge add – Alexander Kartun-Giles Feb 12 '18 at 9:42

In the meantime, I was able to create a CompiledFunction for the body of your long For loop (surprisingly, RandomChoice is compilable!):

cf = Compile[{{i, _Integer}, {norm, _Real}, {k, _Integer}, {theta, _Real, 1}},
Transpose[{
Table[i, {l, 1, k}],
RandomChoice[1./norm theta[[1 ;; i]] -> Range[1, i], {k}]
}],
CompilationTarget -> "C",
RuntimeAttributes -> {Listable},
Parallelization -> True
]


The following shows how you can retrieve, create, and process data a bit more efficiently. DeveloperToPackedArray takes rectangular arrays of integer, real, or complex numbers and tries to store them in a packed format that can be accessed quicker.

ndgraph = Graph[DeveloperToPackedArray[Import["nd.txt", "Table"]] + 1]; // AbsoluteTiming
theta = RandomSample[VertexDegree[ndgraph]]; // AbsoluteTiming
norms = Accumulate[theta]; // AbsoluteTiming
k = RandomVariate[PoissonDistribution[# 0.5]] & /@ theta; // AbsoluteTiming
results = cf[Range[VertexCount[ndgraph]], norms, theta, k + 1]; // AbsoluteTiming
edgelist = UndirectedEdge @@@ (Join @@ results); // AbsoluteTiming


{8.70484, Null}

{0.082112, Null}

{0.002379, Null}

{9.34685, Null}

{138.481, Null}

{1.04823, Null}

I haven't run the original code, yet, so I cannot tell how these timings compare to it.

Some remarks:

• Adding 1 to a List list can quicker be done with list + 1 than with a loop.

• I have to add 1 to the list k since compiled functions do not like empty lists {}. I do not know what the best way to get the original functionality would be.

• cf is parallelized. However, it manages to leverage only half the computational power of my quad core; proabably because the code in cf gets memory bound for large values of i.

• Joining of lists is done once in the end as well as the conversion to UndirectedEdges: Each Join means a copy; UndirectedEdge unpacks arrays and leads to data structures that cannot be digested by CompiledFunctions.

• Thank you, I will try and compile things and see how much faster it runs. – Alexander Kartun-Giles Feb 12 '18 at 13:29
• You're always welcome! – Henrik Schumacher Feb 12 '18 at 15:16
• Yes, it now runs in a matter of seconds what previously took 6 hours, have edited my own answer to show how to achieve the original functionality. – Alexander Kartun-Giles Feb 23 '18 at 10:52

By creating a list of edges, which is the built with Graph[edgelist], its an order of magnitude faster.

(*nd= Import["nd.txt","Table"];
For[i=1,i\[LessEqual]Length[nd],i++,
nd[[i]][[1]]=nd[[i]][[1]]+1;
nd[[i]][[2]]=nd[[i]][[2]]+1;
];
ndgraph=Graph[nd];*)
edgelist = {};
n = Length[VertexList[ndgraph]];
lengthvertexlist = Length[VertexList[ndgraph]];
vd = VertexDegree[ndgraph];
sumofvertexdegrees = Sum[vd[[i]], {i, 1, lengthvertexlist }];
theta = RandomSample[Table[vd[[i]], {i, 1, lengthvertexlist }]];
norm = Abs[theta[[1]]];
k = RandomVariate[PoissonDistribution[#/2]] & /@ theta;
For[i = 2, i <= n, i++,
If[Mod[i, 1000] == 0, Print["Time = ", i];];
j = RandomChoice[
1/norm theta[[1 ;; i - 1]] -> Range[1, i - 1], {k[[i]]}];
edgelist = Join[edgelist, Thread[i <-> j]];
norm += theta[[i]];
];


EDIT: To modify using compile, as in Henrik's answer, we have

cf = Compile[{{i, _Integer}, {j, _Integer}, {norm, _Real}, {k, \
_Integer}, {theta, _Real, 1}},
Transpose[{Table[i, {l, 1, k}],
RandomChoice[1./norm theta[[1 ;; j]] -> Range[1, j], {k}]}],
CompilationTarget -> "C", RuntimeAttributes -> {Listable},
Parallelization -> True]

SampleCount = 10000;
ndgraph =
Graph[DeveloperToPackedArray[Import["nd.txt", "Table"]] +
1]; // AbsoluteTiming
theta = RandomSample[VertexDegree[ndgraph],
SampleCount]; // AbsoluteTiming
norms = Accumulate[theta]; // AbsoluteTiming
k = RandomVariate[PoissonDistribution[# 0.5]] & /@
theta; // AbsoluteTiming
results =
cf[Range[Length[theta]] + 1, Range[Length[theta]], norms, k + 1,
theta]; // AbsoluteTiming
results =
Delete[results[[#]], 1] & /@ Range[Length[theta]]; // AbsoluteTiming
edgelist = UndirectedEdge @@@ (Join @@ results); // AbsoluteTiming
nd100 = Graph[edgelist];


which outputs:

where the original functionality is achieved by removing the 1st added edge. The speed, as Henrik points out, is a matter of seconds, compared to the original 6 hours.

• That's good to hear! – Henrik Schumacher Feb 23 '18 at 11:56