# How to implement the directed Barabasi-Albert algorithm in Mathematica?

Does anyone know how to implement an efficient version of the directed Barabasi-Albert algorithm in Mathematica that scales until hundred of thousands nodes? The directed network mechanism is explained in the following pointer; it suggests that the scaling exponent for the directed case is alpha = 2. I tried to implement the following version:

(* BARABASI-ALBERT IMPLEMENTATION *)
(* Script to generate a \
directed-syntetic graph based on the Barabasi-Albert model with \
parameter m *)

(* Represent the directed graph by using adjacency lists *)

directedInLinks = {}; (* List of lists of incoming neighbors *)
\
directedOutLinks = {}; (* List of lists of outgoing neighbors *)

inDegree = {}; (*List of lists of in-degrees*)

outDegree = {}; (*List of lists of out-degrees*)

m = 2;  (* Number of nodes that every incoming node try to attach *)

T = 1 10^4; (*Total number of nodes to add*)

AbsoluteTiming[

Do[(* For every node *)

AppendTo[inDegree, 1]; (*Assume that every incoming node has in-
degree equals to 1. This is for implememting the preferential \
attachment step*)
AppendTo[outDegree, 0];

parentsPossibilities =
Range[index - 1]; (*Possible existing nodes to attach*)

If[Length[parentsPossibilities] >= m, (*
If -- Verify if there is enough nodes *)

(*Preferential atachment step,
not sure how to made it more efficient*)

parents =
RandomSample[
inDegree[[1 ;; Length[parentsPossibilities]]] ->
parentsPossibilities, m];

Do[(*For every parent*)

inDegree[[mr]]++;
outDegree[[index]]++;

, {mr, parents}];
] (* If -- Verify if there is enough nodes *)

If[Mod[index, T/10] == 0,
Print[index]]; (*Track the growing process*)

, {index, 1, T}];

](*AbsoluteTiming*)

(* To have the exact in-dergee in evey node*)
inDegree = inDegree - 1;

(*PLOT*)
(*Get the empirical distirbution*)

emdist = EmpiricalDistribution[inDegree];

(*Analytical CDF*)
analyticalDegree[inDegree_] := (inDegree + 0.1)^-1;

(*Plot in-Degree*)
Show[
ListLogLogPlot[{Table[{x, 1 - CDF[emdist, x]}, {x,
DistributionDomain[emdist]}]}, PlotRange -> All,
PlotMarkers -> {Graphics[{Red, Circle[{0, 0}, 1]}], 0.03}],
ListLogLogPlot[Table[{i, analyticalDegree[i]}, {i, Union[inDegree]}],
Joined -> True, PlotStyle -> Black]]


Although, it seems to work, I was wondering a way to do it more efficiently (I need to run orders of magnitude close to half million nodes) and it is still too slow, I guess that the bottle neck is in the preferential attachment condition. Any suggestion? Thank you.

• Please at least provide a good link to the algorithm description Jan 14 '15 at 19:33
• What do you mean by "Barabási-Albert algorithm"? Do you mean preferential attachment? Are you aware of BarabasiAlbertGraphDistribution? RandomGraph@BarabasiAlbertGraphDistribution[100000, 1] takes no perceptible time. Jan 14 '15 at 19:38
• Voting to close, as it's already built-in. Just use RandomGraph[BarabasiAlbertGraphDistribution[100000, 5]]; or similar. Generating a Barabasi-Albert graph with one hundred thousand notes executes in a couple milliseconds on my computer, and easily scales to tens of millions (or billions, if you have enough RAM) nodes. Jan 14 '15 at 19:41
• I mean for a directed network
– Paul
Jan 14 '15 at 19:49
• Hmm, it looks like DirectedEdges -> True is not supported in RandomGraph for BarabasiAlbertGraphDistribution. In that case I'll remove the close vote. In the mean time, could you edit your question to explain what you mean by a directed Barabasi-Albert distribution? All the times I've seen Barabasi-Albert distributions, they've been undirected, and the Wikipedia page only covers the undirected case. Jan 14 '15 at 20:13

Here is an attempt with roughly $O(n^{1.6})$ runtime complexity. First, let n be the number of nodes, m be the number of out-edges from each node, p the randomness parameter in the paper you linked, and iD be the in-degrees of each node. For initialization, I'll let the first m nodes have in-degree 1:

SetSystemOptions["SparseArrayOptions" -> {"TreatRepeatedEntries" -> 1}];
n = 10000;
m = 2;
p = 0.5;
iD = ConstantArray[0, n];
c = ConstantArray[1, n];
iD[[;; m]] = 1;


Now we can construct the sparse adjacency matrix:

A = SparseArray[Flatten[Reap[Do[
If[RandomReal[] < p, w := c, w := iD];
t = RandomChoice[w[[;; i - 1]] -> Range[i - 1], m];
iD[[#]]++ & /@ t;
, {i, m + 1, n}]][[-1]]], {n, n}];


On my machine, a 10,000-node system constructs in 0.36 seconds, which is about 45 times faster than the initial implementation. The main reason it's faster than your version is because it avoids the use of AppendTo, which is slow when used to iteratively grow a list (because it constructs a new list at each step, incurring an $O(n^2)$ cost), and instead uses the Reap and Sow functions (which have $O(n)$ cost).

Note that by construction, the adjacency matrix is lower triangular:

LowerTriangularize[A] == A


which produces

True


### Performance Concerns

While this is much faster than the initial attempt, it is still orders of magnitude slower than the built-in undirected Barabasi-Albert algorithm. For example,

RandomGraph[BarabasiAlbertGraphDistribution[500000, 2]];


executes in 0.18 seconds, whereas the directed Barabasi-Albert $n=500000$ case took 20 minutes on my computer. The bottleneck is the following step:

t = RandomChoice[w[[;; i - 1]] -> Range[i - 1], m];


which makes a weighted choice of nodes that node i will connect to. All other steps in the Do loop proceed in $O(1)$ time. Does anyone have suggestions as to how to improve this? Unfortunately, w changes after each step, which makes things difficult.

I'm not sure if I understood the question properly, but I was trying to generate a random graph from the Albert-Barabasi distribution, with directed edges that would differ in the IN and OUT degrees for each node. It seems that the option for DirectedEdges -> True is not yet implemented for a random graph with Albert-Barabasi distribution. So, at the end, I tried the following, which it might not be what you were looking for, but maybe some people will find it useful:

g0 = RandomGraph[BarabasiAlbertGraphDistribution[50000, 2]];
g = DirectedGraph[g0, "Random"];
VertexInDegree[g]
VertexOutDegree[g]


So, in my case, this works because I just wanted to randomly choose which edges would be in/out for each node, so I hope this is of some use.

Cheers, Pedro

• IGraph/M includes a directed version through the IGBarabasiAlbertGame function. It supports the DirectedEdges options. Nov 3 '17 at 7:31
• Cool! I'll check it out :) Nov 4 '17 at 16:25