Fast code to generate large random graph with fixed degree sequence

Question

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[[1]] -> 0, 
 pick[[2]] -> 0, {pick[[1, 1]], pick[[2, 2]]} -> 
 partM[old, pick[[1]]], {pick[[2, 1]], pick[[1, 2]]} -> 
 partM[old, pick[[2]]]}]]

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

Answer 1

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[25]= {{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).

Answer 2

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:

IGRewire[AdjacencyGraph[mat], 10000]

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.

Answer 3

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[#[[2]], n - 1], 
         Join[ConstantArray[1, #[[1]]], ConstantArray[0, n - 1 - #[[1]]]]]] & /@ obj, 1]
   ];
m = {4, 2, 2, 1, 3, 2, 2};
AdjacencyGraph[randomDegreeGraph@m, VertexLabels -> "Name", DirectedEdges -> True]
Tr /@ AdjacencyMatrix@g == m