Fast code to generate large random graph with fixed degree sequence

I'm trying to create a function so to generate large random graph with a fixed degree sequence $\{v_1,v_2,\dots,v_N\}$, where $N$ is the number of vertices in the graph. The Mathematica function that should do the work is:

RandomGraph[DegreeGraphDistribution[{v1,v2,...vN}]]

but unfortunatly, it is useless for large graph, as it is extremely slow as the number of nodes increases.

I have tried to build my own function, where instead of generating a new graph with a given degree sequence, I start from a non random matrix m, that has the desidered degree sequnce $\{v_1,v_2,\dots,v_N\}$ and then I rewire the links so to keep the degree fixed. Here is my code:

partM[m_, {a_, b_}] := m[[a, b]]

confmodel0[m_] :=
Module[{x, y, old, pos, pick}, {x, y} = Dimensions[m];
pos = Position[m, x : _ /; x > 0.]; pick = RandomChoice[pos, 2];
old = m;
While[m[[pick[[1, 1]], pick[[2, 2]]]] != 0 ||
m[[pick[[2, 1]], pick[[1, 2]]]] != 0 ||
pick[[1, 2]] - pick[[2, 2]] == 0 ||
pick[[1, 1]] - pick[[2, 1]] == 0, pick = RandomChoice[pos, 2]];
ReplacePart[
m, {pick[] -> 0,
pick[] -> 0, {pick[[1, 1]], pick[[2, 2]]} ->
partM[old, pick[]], {pick[[2, 1]], pick[[1, 2]]} ->
partM[old, pick[]]}]]

confmodel[m_] := Nest[confmodel0, m, 10*Total[Boole[Positive[m]], 2]];

However, also this function is too slow for large graph and I cannot generate real randomized graphs.

In network theory, this model is called configuation model. It would be nice to have a function for both directed and undirected graph.

Here you find a test matrix. I would like to generate a random matrix with the same degree sequence of this matrix: https://dl.dropboxusercontent.com/u/62056077/TestMatrix.m

Thanks for any suggestion

• Are you aware of the work by Fabien Viger (Google dude) in this area? I've used his algo to generate graphs, very fast. Perhaps a peek at the source might give you some ideas. – ciao Apr 1 '14 at 21:35
• How large do the graphs you're generating need to be? And do you have a test data set we could use for benchmarking? – Pillsy Apr 1 '14 at 22:10
• @Pillsy. Good idea. I added a link so you can download test data. Thanks – sam84 Apr 2 '14 at 9:13
• DegreeGraphDistribution implements the configuration model as well. Do you have any reason to believe that its implementation is not efficient enough and it is at all possible to do better in pure Mathematica, not in C. – Szabolcs Apr 6 '14 at 21:18

This comes with several caveats:

(1) It is also often slow. I have reason to believe it gets careless about certain "painted into a corner" situations, and thus might simply fail.

(2) It gives results that are in no sense uniformly random, across the range of possible graphs that meet the requirements.

(3) I wrote it some time ago and no longer understand any of it. So it comes "as is", with no guarantees of anything beyond "you get what you get".

(4) The form of result might not be the most desirable.

insufficient[degs_, edgelist_, donelist_] := Catch[Module[
{j, k, n = Length[degs], diff},
For[j = 1; j <= n, j++,
If[degs[[j]] == 0, Continue[]];
diff = degs[[j]] - Length[Union[donelist, edgelist[[j]]]];
If[diff <= 0, Throw[True]];
];
False
]]

decode[val_, degs_] :=
Module[{j = 0, k = 0}, While[val > k && j < Length[degs], j++;
k += degs[[j]];];
j]

decode[val_, degs_, skip_] :=
Module[{j = 0, k = 0}, While[val > k && j < Length[degs], j++;
If[! MemberQ[skip, j], k += degs[[j]]];];
j]

randomDegreeGraph[degrees : {_Integer ..}, tries_: Infinity] :=
Catch[Module[
{edgelist, degs, n = Length[degrees], nedges, npairs, j, v1, v2,
subtot, try = 0, donelist = {}, prob},
If[(tries =!= Infinity && (! IntegerQ[tries] || tries < 0)),
Throw[$Failed]]; If[! EvenQ[Total[degrees]], Throw[$Failed]];
While[try < tries, try++;
npairs = n*(n - 1)/2;
nedges = Total[degrees]/2;
edgelist = Table[{j}, {j, n}];
degs = degrees;
While[npairs >= nedges > 0,
j = RandomInteger[{1, 2*nedges}];
v1 = decode[j, degs];
subtot = 2*nedges - Total[degs[[edgelist[[v1]]]]];
If[subtot <= 0, Break[]];
j = RandomInteger[{1, subtot}];
v2 = decode[j, degs, edgelist[[v1]]];
nedges--;
npairs--;
edgelist[[v1]] = Append[edgelist[[v1]], v2];
edgelist[[v2]] = Append[edgelist[[v2]], v1];
degs[[v1]] = degs[[v1]] - 1;
degs[[v2]] = degs[[v2]] - 1;
If[(degs[[v1]] == 0 || degs[[v2]] == 0) &&
insufficient[degs, edgelist, donelist], Break[]];
If[degs[[v1]] == 0, donelist = Append[donelist, v1];
npairs = npairs - (n - Length[Union[donelist, edgelist[[v1]]]])];
If[degs[[v2]] == 0, donelist = Append[donelist, v2];
npairs =
npairs - (n - Length[Union[donelist, edgelist[[v2]]]])];];
If[nedges == 0, Throw[Map[Sort, edgelist[[All, 2 ;; -1]]]]];];
Throw[\$Failed["too many iterations"]];]]

Example:

randomDegreeGraph[{4, 2, 2, 1, 3, 2, 2}]

(* Out= {{2, 4, 5, 7}, {1, 5}, {6, 7}, {1}, {1, 2, 6}, {3, 5}, {1,
3}} *)

The resulting graph has vertices adjacent to vertex j given as a list in the jth position. So vertex 1 is connected to vertices 2, 4, 5, and 7, and vertex 2 is connected to vertices 1 and 5, etc.

Feel free to post any improvements, especially on item (1).

• unfortunatly the code seems now working for large matrix (say 100 nodes). I will open a bounty, and also post some updates if I can find a better solution. Thanks – sam84 Apr 6 '14 at 21:06

Not completely sure I'm understanding. Anyway, for Directed Graphs without self loops:

randomDegreeGraph[m_List] := Module[{n, t, t1, obj, s, x},
n = Length@m;
t = Table[If[i != j, x[i, j], Sequence @@ {}], {i, n}, {j, n}];
t1 = Table[If[i != j, x[i, j], 0], {i, n}, {j, n}];
obj = Transpose[{List /@ m, t}];
t1 /. Flatten[Thread[Rule[RandomSample[#[], n - 1],
Join[ConstantArray[1, #[]], ConstantArray[0, n - 1 - #[]]]]] & /@ obj, 1]
];
m = {4, 2, 2, 1, 3, 2, 2};
AdjacencyGraph[randomDegreeGraph@m, VertexLabels -> "Name", DirectedEdges -> True]
Tr /@ AdjacencyMatrix@g == m • The function you propose works very well. It would be nice to modify the function so to conserve both out and in-degrees (now the control is only on the out degree). – sam84 Apr 8 '14 at 7:22

Update: Now I recommend using IGraph/M instead of IGraphR when the desired functionality is already implemented.

IGDegreeSequenceGame[Total[mat], Method -> "VigerLatapy"]

The Viger–Latapy method always generated connected graphs, but doesn't support directed ones.

We can also take the directed graph corresponding to mat and randomly rewire its edges while keeping the degree sequence:

This does 10000 rewiring trials, not all of which may be successful.

Original post:

You can use igraph through my RLink-based package IGraphR.

mat = Import["https://dl.dropboxusercontent.com/u/62056077/TestMatrix.m", "Package"];

mat == Transpose[mat]

(* ==> False *)

Your adjacency matrix is not symmetric. Are you looking for directed or undirected graphs?

degs = Total[mat];

This degree sequence happens to be graphical, so let's do undirected. Please set up and load IGraphR` as described in the linked post first.

IGraph["degree.sequence.game"][degs, Null, "vl"]

This takes less than a second on my machine.

This igraph function is documented here.

The "vl" method generates undirected connected graphs. It has methods that work for directed graphs too but those do not sample with uniformly from the set of allowed graphs.