Exact sampling of graphs with prescribed degree correlations

Many real-world networks exhibit correlations between the node degrees. For instance, in social networks nodes tend to connect to nodes of similar degree. Conversely, in biological and technological networks, high-degree nodes tend to be linked with low-degree nodes. Degree correlations also affect the dynamics of processes supported by a network structure, such as the spread of opinions or epidemics. The proper modelling of these systems, i.e., without uncontrolled biases, requires the sampling of networks with a specified set of constraints. We present a solution to the sampling problem when the constraints imposed are the degree correlations. In particular, we develop an efficient and exact method to construct and sample graphs with a specified joint-degree matrix, which is a matrix providing the number of edges between all the sets of nodes of a given degree, for all degrees, thus completely specifying all pairwise degree correlations, and additionally, the degree sequence itself. Our algorithm always produces independent samples without backtracking. The complexity of the graph construction algorithm is O(NM) where N is the number of nodes and M is the number of edges.


Introduction
Complex systems often consist of a discrete set of elements with heterogeneous pairwise interactions. Networks, or graphs have proven to be a useful representational paradigm for the study of these systems [1][2][3][4]. The nodes, or vertices, of the graphs represent the discrete elements, and the edges, or links, represent their interaction. In empirical studies of real-world systems, however, for reasons of methodology, privacy, or simply lack of data, frequently there is only limited information available about the connectivity structure of a network. When this is the case, one has to take a statistical approach and study ensembles of graphs that conform to some structural constraints. This statistical approach enables the computation of ensemble averages of network observables as determined solely by the constraints, i.e., by the specified structural properties of the graphs. Ensemble modeling of this type is necessary to determine the relationship between the given structural constraints and the behavior of the complex system as a whole. Calculating ensemble averages, though, requires the ability to construct all the graphs that are consistent with the required structural constraints, a highly non-trivial problem.
Perhaps one of the simplest examples of structural constraints that occur in data-driven studies of real-world systems is to fix the degree of each node, which is the number of edges that are connected to, or are incident on the node. For an undirected graph with N nodes this information is specified by a degree sequence  = ⋯ Graph construction and sampling becomes even more difficult when there are structural constraints of higher order, such as correlations amongst the node degrees. Degree correlations can be expressed in several ways, for example with the help of the conditional probability P d d ( ) ′| that a node of degree d will have a neighbor of degree d′, or more simply, by the average degree of the neighbors of a node with degree d, [36]. The properties of d d ( ) ′ characterize the so-called assortativity of a graph, which is a measure of the tendency of a node to connect to nodes of similar degree. If d d ( ) ′ is increasing in d, the graph is degree-assortative, if it is decreasing the graph is degree-disassortative, and if it is constant, the graph is degreeuncorrelated. Even more coarse-grained measures of degree correlations are possible, including the Pearson coefficient [37], the Spearman coefficient [38] and the Kendall coefficient [39]. These coefficients assume values ranging from 1 − , for highly disassortative graphs, to 1, for highly assortative ones. A more precise way to express degree correlations is via the use of a JDM. The JDM of a given undirected simple graph is a symmetric matrix whose ( , ) α β element is the number of edges between nodes of degree α and nodes of degree β. The dimensions of the JDM are Δ × Δ, where Δ is the largest degree of a node in the graph. The degree correlation measures discussed above specify the correlations only statistically, but they do not fix the number of edges between nodes of given degrees, whereas the JDMs do. In this sense, the relationship between JDMs and the statistical degree correlation measures is similar to the relationship between degree sequences and degree distributions. Degree correlations have generated considerable interest, as they are known to affect many structural and dynamical properties of graphs and the processes they support [40][41][42][43][44][45][46][47]. Nevertheless, even though their importance is well established, it has heretofore not been possible to perform ensemble modeling of graphs with prescribed JDMs. In this Article, we solve this problem by developing an algorithm based on the stub-matching method to construct and sample ensembles of graphs with a specified JDM.

Mathematical foundations 2.1. Graphicality of JDMs
The problem of graphicality for JDMs asks whether a specified symmetric matrix can be the JDM of a simple graph. Our starting point is an Erdős-Gallai-like theorem that gives the requirements for a JDM to be graphical [48][49][50].
Before stating the theorem, though, note that a JDM specifies uniquely the degree sequence of the graphs that realize it [48]. Given a JDM J, the number of nodes with degree α is where V α is the set of nodes, or degree class, with degree α. As a general rule of notation we will use lowercase Greek letters to indicate degree values and lowercase Latin letters for node indices. In the equation above the sum of each row α of J is the number of connections involving nodes of degree α (i.e., all nodes in class V α ). As each node of degree α has exactly α stubs the total number of nodes of degree α is given by the notal number of stubs from all nodes in class V α divided by α. Moreover, each edge between nodes of the same degree involves 2 stubs. Thus, the diagonal elements must be double-counted. Note that multiple JDMs can specify the same degree sequence and thus prescribing a JDM is more constraining than only prescribing a degree sequence. With the definitions above, the necessary and sufficient conditions for a JDM to be graphical can be stated as follows [48][49][50]: Theorem 1. JDM graphicality. A symmetric Δ × Δ matrix J with non-negative integer elements is a graphical JDM if and only if: It is important to observe that any graphical realization of a JDM can be decomposed into the disjoint union of a set of subgraphs G αβ that are bipartite (α β ≠ ) with node sets V α and V β and J αβ edges between them or unipartite (α β = ) with node set V α and J αα edges within that set. We are going to call such representation of a graphical realization a degree class representation. A simple example of a graphical JDM with N = 10 and 4 Δ = is given by the matrix: ( 1 ) Panels (a) and (b) of figure 1 show a graphical realization of J in degree-class representation and regular representation, respectively. Panels (d) and (e) of the same figure show another realization of J in the two representations. The color of the edges indicate the subgraph they belong to. For example, G 24 is a bipartite graph between nodes of degree 2 (V 2 ) and 4 (V 4 ), respectively, having J 4 2,4 = edges drawn in green color, whereas G 33 is unipartite with a single J 1 33 = edge drawn in blue. Note that while both graphical realizations have the same JDM, they are very different graphs. To see this, consider the counts n ℓ of cycles C ℓ of length ℓ (a cycle is a closed path without repeated nodes). The graph in figure 1(b) has n 1 3 = , n 2 4 = , n 1 5 = , n 2 6 = , n 3 7 = and n 3 8 = , whereas the one in figure 1(e) has n 1 3 = , n 1 4 = , n 2 5 = , n 3 6 = , n 4 7 = and n 1 8 = .
Theorem 1 is an existence theorem, just like the Erdős-Gallai theorem for the case of degree sequences, and as such it does not provide an algorithm that can generate simple graphs with a given JDM. More importantly, we also need an algorithm that does not exclude classes of graphical realizations of a given JDM, but that can construct in principle any such realization. The situation is similar to that of degree sequences. In that case the Havel-Hakimi method [6,7] is always able to create a graphical realization of a graphical degree sequence, but cannot construct them all, i.e., there will be some realizations that can never be built by this algorithm. This was the reason for the introduction of the notion of star-constrained graphicality in [32,33] and the subsequent construction algorithms in [9,11]. Here as well, we want to have a direct construction algorithm and ultimately an exact sampler that does not exclude any realization of a JDM. Due to the different nature of the constraints from the degree-sequence-based case, we need to develop a novel approach.
The idea of the approach is based on the degree-class representation above. Since the edges of the subgraphs G αβ are disjoint, we could build a graphical realization G of the JDM J by building all these subgraphs, while respecting the constraints. For a G αβ subgraph we know its node set(s) and its total number of edges J αβ . Consider then a node v V ∈ α . We are not given its degree in G αβ for any β, but we know that the sum of its degrees within every one of these subgraphs must add up to α. For example, the sum of the number of the purple, green and red edges coming out of node 2 in figure 1(b) must add to 4. In addition, we also have the constraints that the sum of the degrees of one color of all nodes within V α must equal to the corresponding given JDM entry. Indeed, for example, the sum of all green edges in figure 1(a) or figure 1(b) is J 4 2,4 = , for orange is 4, red is 3, etc. Thus, the idea of the algorithm is to first determine the degree of a given color respecting the constraints for all nodes and all colors, then use our methods introduced earlier [9, 11] (see A) to build the G αβ subgraphs based on the corresponding degree sequences of their nodes. Different graphical realizations will be obtained from different assignments of color degrees and, of course, from the different graphical realizations of the same set of degrees. Note that for the bipartite subgraphs G αβ we are specifying degree sequences for nodes in both partitions V α and V β and thus we can use our graph construction method for directed graphs [11], because a bipartite graph can be represented as a directed graph if nodes in one partition have only outgoing edges and in the other only incoming edges. In the following it will be useful to introduce the notion of degree spectra, representing the degrees of different colors of a node, as described above.

Degree spectra
Consider a single row α of a graphical JDM J. The information contained in the row determines the precise number of edges needed between nodes of degree α and nodes of every degree. In other words, of all the stubs coming from V α , J ,1 α of them must end in a node of degree 1, J ,2 α of them must end in a node of degree 2, and so on. However, these matrix elements do not specify how to distribute these edges within and between the degree classes. To better specify these connections one introduces the notion of the degree spectra, which can be conveniently represented as a matrix. The degree spectrum of a node is the sequence of its degrees towards all the degree classes, including its own degree class. A degree-spectra matrix S is a N Δ × matrix whose i ( , ) α element S i α is the number of edges between node i and degree class α (the set of nodes of degree α). The ith column of S defines the degree spectrum of node i. Panels (c) and (f) of figure 1show two representations of the same JDM given in equation (1). In general, there are many degree-spectra matrices that correspond to the same JDM. As described in the previous section, we employ a two-step process in order to randomly sample graphs that realize a given JDM. First, we generate a random degree-spectra matrix from the JDM. Second, we construct a random graph that realizes the JDM and that obeys or is consistent with the chosen spectra matrix. This approach creates the need for a method to guarantee that the spectra generated from a JDM are graphical.
The generation of a graphical degree-spectra matrix proceeds systematically, node by node. Therefore, at each step, some nodes will have an already fixed number of links within some of the subgraphs (links of a given color), while for the rest these numbers will not have been determined yet. Thus, at any time during this process we have a partial degree sequence of a bipartite graph. As the subgraphs must be simple graphs (realizable), one must be able to decide whether a partial bipartite degree sequence is graphical. The sufficient and necessary criterion for the graphicality of a partial bipartite degree sequence will be given in theorem 2 below. However, that will not necessarily mean that the whole JDM J is still realizable, in other words, how do we know that by guaranteeing the graphicality for a subgraph G αβ we have not precluded graphicality of some other subgraph G γδ , and ultimately of J? The answer to this question will be given by theorem 3, later on. Together, these theorems form the basis for our algorithm to generate graphical degree spectra.
Before proving a theorem that provides a graphicality test for partial bipartite degree sequences, we need to set some notations. Let A, B, H and K be four sets of nodes: H K figure 2). The sets can be of different size, but neither U nor V can be empty. Now, let be two given sequences of integers. They will represent the partial bipartite degree sequences that have already been fixed by the algorithm up to that point. The degrees of the other nodes, specifically those in the sets B and K, are not yet specified. What is specified is the total number of edges ε in the bipartition, i.e., the total number edges running between the sets U and V. Then, the partial bipartite degree sequence triplet ( , , )   ε , hereafter simply called a triplet, is graphical if there exists a bipartite graph on U and V with ε edges and degree sequences In other words, the bipartite graph must be such that the nodes in A have degree sequence  and those in H have degree sequence . The partial degree sequence problem is to decide whether one can choose the degrees of the nodes in the sets B and K such that the above constraints are satisfied and the bipartite degree sequence  is graphical.
Since the graph realizing a triplet is bipartite, the number of edges ε equals the number of stubs in either set of nodes: Figure 2. Schematic for the partial degree sequence problem.
The imposed partial sequences  and  prescribe a certain number of stubs in the first A | |nodes of U and in the first H | |nodes of V. Let these be P p , respectively. Then, the set B must contain exactly P ε − stubs; similarly, the set K must contain exactly Q ε − stubs. With these considerations, we first define the concept of a balanced realization of a triplet. Let A realization of a triplet is defined to be balanced if and only if the degree of any node in B is either μ ⌊ ⌋or μ ⌈ ⌉, and the degree of any node in K is either ν ⌊ ⌋or ν ⌈ ⌉. Notice that this means that if μ or ν are integers, then all the nodes in B or K must have exactly degree μ or ν, respectively. Conversely, if they are not integers, then the degrees of any two nodes in B or in K, respectively, can differ at most by 1. That is, a realization is balanced if and only if all the degrees of the nodes that one is free to choose (those in B and K) are as close as possible to their averages μ and ν. The definition can be equivalently formalized by introducing a functional f acting on B and K: Then, a realization of a triplet is balanced if and only if both f B ( )and f K ( )vanish. An important theorem about the graphicality of triplets can now be proven.
  ε is graphical if and only if it admits a balanced realization.
Proof. Sufficiency is obvious. If the triplet admits any realization, balanced or not, it is graphical by definition.
To prove necessity, suppose the triplet is graphical. Then, it admits a realization G. If G is balanced, then there is nothing to do. Conversely, if G is not balanced, then f B ( ), f K ( ), or both, are greater than 0. Without loss of generality, assume that f B Again without loss of generality, assume that d b i μ < ⌊ ⌋(the other cases are treated analogously). Then, since the number of stubs within B is fixed, there must exist a node b This yields a different realization with the same degrees for the nodes in V, and in which f B ( )is decreased by at least 1, as the degrees of B moved towards the balanced condition. The procedure can be repeated until f B ( ) 0 = , resulting in a balanced realization. □ A key consequence of this theorem is the following.
Corollary 1. Let ( , , )   ε be a graphical triplet, and let x be a node in B or in K. If there is a realization of the triplet in which d x α = and another in which d x β = , with α β < , then for all γ with α γ β ⩽ ⩽ there exists a realization in which d x γ = .
Proof. Without loss of generality, assume x B ∈ . Then, there are several cases, each determined by the relative values of α, β and μ ⌊ ⌋. The most general case is α μ β < ⌊ ⌋ < , so consider only this situation. Start from the realization with d x β = . Repeated applications of the method in the proof of theorem 2 will eventually yield a realization in which d x μ = ⌊ ⌋. For each step, the degree of x will have decreased by 1. Therefore, one realization of the triplet will have been found with d x γ = for all μ γ β ⌊ ⌋ ⩽ ⩽ . Now, start from the realization with d x α = . Applying the same step from the proof of theorem 2 repeatedly will eventually yield a realization in which d x μ = ⌊ ⌋. For each of these steps, the degree of x will have increased by 1. Therefore, one realization of the triplet will have been found with d Notice that, given a graphical triplet, corollary 1 also implies the existence of minimum and maximum allowed degrees for each node whose degree has not yet been fixed in that triplet (namely, in B and K). That is, a realization of the triplet exists with a node having either its minimum or maximum degree, or any degree between these two values. Of course, the value of the minimum and maximum degree will depend on which degrees have been fixed up to that point, so these need to be computed on the fly. How to calculate these degree bounds will be explained in section 3.1.

Building a degree-spectra matrix
Corollary 1 suggests the possibility of a direct, sequential way to build a degree-spectra matrix from a JDM. However, building the degree-spectra matrix node by node is a local process, which guarantees via theorem 2 only that the bipartite graph in which the node whose degree spectrum is being set resides is graphical. There is a global constraint, however, on every node, namely that the sum of their degree-spectrum elements must add up to the degree of the class they belong to. We have to make sure that the local construction process also respects the global constraints, i.e., it is feasible with it. The theorem below will show that this sequential construction process is feasible, and just as importantly, all graphical realizations of a JDM J can be constructed in this way, i.e., all graphical degree-spectra matrices can be obtained by this sequential construction process.
Theorem 3. Let  be the subset of all the nodes with fixed spectra; then, there exists a realization of a JDM J consistent with the fixed spectra if and only if for every ( , ) α β pair with , {1, , } α β ∈ … Δ there exists a graph G αβ with J , α β edges also satisfying the fixed spectra of .
Proof. Necessity is obvious. If there exists a realization of J satisfying the spectra, then each subgraph between any pair of degree classes both satisfies the spectra and has the right number of edges.
To prove sufficiency, assume that we have a fixed degree spectrum for all the nodes in  and we have guaranteed the graphicality of all the subgraphs G αβ . They have the right number of edges J , α β and their nodes satisfy the fixed spectra specified in the subset . Since we have guaranteed graphicality for all the G αβ subgraphs with these constraints, let us consider some graphical realization for each such subgraph and consider their union graph G. If the 'free' nodes, i.e., those without a fixed spectrum, have all the correct degree in G (i.e., every node v V ∈ α has d v α = for all α), then there is nothing to do. Now, assume they don't. Since the total number of edges in each G αβ is correct by hypothesis, there must exist a degree α and two free nodes v and w belonging to V α such that d v α < and d w α > . Thus, there must exist a node u connected to w but not to v. Then, erase the edge (u, w), and replace it with (u, v). This leaves the numbers of edges in all G αβ unchanged, and does not change the degree spectrum of u, because v and w belong to the same degree class. Repeating this procedure results eventually in all the nodes having the correct degree. □ Theorem 3 is fundamentally important as it justifies a systematic, node-by-node approach in building a graphical degree-spectra matrix. In fact, so long as one guarantees the possibility of subgraphs with the correct number of edges, a partial degree-spectra matrix maintains the graphicality of the JDM.
The only detail left is specifying how to choose the numbers that form the degree spectra. Fortunately, this is straightforward. As mentioned in the previous section, an implication of corollary 1 is the existence of minimum and maximum allowed degrees for nodes in partial degree sequences. Let them be m (minimum) and M (maximum). But a partial degree sequence is nothing else than a partially built degree spectrum, if one recognizes the node sets U and V as two degree classes. Then, a condition that must be satisfied in building a degree-spectra matrix is that any new number chosen to augment a partially built degree spectrum has to be within these bounds. However, one must also consider that if a node belongs to a certain degree class, it must have the correct total degree.
To state both conditions, assume the degree spectrum of node v V ∈ α is being built. Let Γ be the set of degree classes for which a spectrum element has already been chosen, and let S v β be the element to determine next. Then, a valid value k for S v β must satisfy the two conditions m k M , Below, in section 3.1 we describe how to compute the min and max values for degree spectra elements.

Description
We are now ready to describe our JDM sampling algorithm. The algorithm is composed of two parts. The first is a spectra sampler that randomly generates degree-spectra matrices from a graphical JDM J: (i) Initialize i = 1.
(iii) Let l be the number of the residual, unallocated stubs of node i. If l 0 ≠ : (4) Extract an integer S i , α uniformly at random between r and R.
(b) Increase α by 1, and go to step (iii).
(iv) Increase i by 1. If i N ⩽ , go to step (ii).
To find the values of m and M in step (iii).a.1 above, consider the degrees of the nodes belonging to V α and V β in G αβ . In the formalism of section 2.2, the already fixed spectra elements are equivalent to the sequences  and . Then, to test the viability of a given value as a degree-spectrum element, assign it to the element being determined, complete the degree sequence making it balanced, and test it for graphicality, see figure 3. If the sequence is graphical, then the triplet has a balanced realization, which by theorem 2 is a necessary and sufficient condition for the existence of a subgraph corresponding to the spectrum element being determined. If G αβ is unipartite, the graphicality test can be done using the fast method described in [9]. The situation is marginally different if G αβ is bipartite. In this case, as previously mentioned, the degree sequence can be built as a BDS in which nodes of degree α only have incoming edges, and nodes of degree β only have outgoing ones. This sequence can then be tested with the fast directed graphicality test described in [11].
Thus, to find the minimum value m one can simply run a sequential test, checking for valid spectrum values from 0 onwards. The first successful value is m. Then, to find M, use bisection to test all the values from m 1 + to the theoretical maximum, looking for the largest number allowed. Clearly, the theoretical maximum at that stage is the degree of the class the node belongs to minus the sum of the already fixed spectra values for that degree.
These considerations also clarify the nature of the second part of the algorithm, which samples realizations of the JDM from an extracted degree spectra matrix. Summarizing, • JDM realizations can be decomposed into a set of independent unipartite and bipartite graphs.
• The degree spectra define the degree sequences of the component subgraphs.
Then, to accomplish the actual sampling, extract the degree sequences from the degree spectra and use them in the graph sampling algorithms for undirected and directed graphs presented in [9,11] and in here in A. Every time a sample is generated, it constitutes a subgraph of a JDM realization. All that is needed in the end is simply to list the edges correctly, since the graph realizing the JDM is the union of all the unipartite and bipartite subgraphs into which it has been decomposed.

Sampling weights
Our algorithm does not extract all degree-spectra matrices from a JDM with the same probability. However, the relative probability for the extraction of each spectra matrix is easily computed, and it can be used to reweight the sample and obtain unbiased sampling. If every new element of a degree-spectra matrix is extracted uniformly at random between r and R, its probability of being chosen is simply  In the expression above, Q i is the value that Q assumes on the ith sampled matrix. Indicating by r j and R j the values that r and R assume for the jth matrix element extracted, the weights are Of course, besides the spectra matrix, every subgraph has its own sampling weight. Thus, the total weight of a single JDM sample is the product of the corresponding spectrum weight and all the subgraph weights. To describe the distribution of the sample weights, first recall that the individual subgraph weights are log-normally distributed [9,11]. Thus, as the sample weights are their product, we expect them to be log-normally distributed too. Also, for large JDMs, where 1 2 Δ ≫ , the m factors in equation (5) are effectively random. Thus, our expectation is that the spectra weights are log-normally distributed as well. To verify this, we extracted the JDM of a random scale-free network with 1000 nodes and power-law exponent of 2.5, and used it to generate an ensemble of 10 5 degree-spectra matrices and one of 10 8 JDM samples of a single spectra matrix. Figure 4 shows that the histograms of the logarithms of spectra matrix weights and sample weights are well approximated by a Gaussian fit, supporting our assumptions.
A simple and small example is provided in B. There, we analytically compute the JDM ensemble averages of the local clustering coefficients of nodes of all degrees, based on unweighted sampling and also based on weighted sampling, with the weights provided by the algorithm. In table B.1 , we show the results of simulations using our algorithm, taking into account the sample weights (as described above), and simply computing the averages of the clustering coefficients over the samples generated. The results between theoretical and simulated measures agree very well. The differences between weighted and unweighted versions can be also appreciated, and while they are small in this example, they are measurable and need to be taken into account in general.

Computational complexity
To determine the computational complexity of the algorithm, first note that the main cost in creating a spectra matrix comes from the repeated graphicality tests. Let A be the number of non-empty degree classes in the JDM
Then, for each of the NA non-trivial elements in the degree-spectra matrix, A tests are needed, each with a computational complexity of the order of the number of nodes in the corresponding degree class. Thus, the total computational complexity for the spectra construction part of the algorithm is Notice that in our treatment one is free to choose the order of the degree classes. Thus, to minimize the complexity, one can simply determine the degree-spectra elements in descending order of degree-class size. Then, the worst case corresponds to the equipartition of the nodes amongst degree classes, V if the number of degree classes is of the same order as the number of nodes.
A more precise estimate for a given JDM can be obtained by rewriting equation (6) as where the degree distribution P d V N ( ) d = | | is the probability that a randomly chosen node has degree d. It is easy to see, then, that the worst case is unlikely to occur. Consider for instance systems of widespread insterest, such as scale-free networks, for which P d d ( ) ∼ γ − with 2 γ > . Then, in the limit of large networks, the equation above becomes ( ) Thus, in this case, the complexity leading order for spectra matrix extraction is only quadratic. Given a degree-spectra matrix, to construct a JDM realization one then needs to build A ( ) . Therefore, the total complexity of the graph construction part of our method is N ( ) 2  for sparse networks, and N ( ) 3  for dense ones. Once more, we do not expect the worst case complexity to occur often. For example, in the already mentioned case of scale-free networks, which are always sparse [51], the total complexity of our algorithm would only be quadratic. A less efficient sampling method has been developed recently [52], but it is based on backtracking, producing results containing biases that are uncontrolled and that cannot be estimated. An alternative algorithm, that does not involve backtracking, has been proposed [53]. However, despite having a computational complexity that is comparable to ours, being (N 2 ) on highly heterogeneous networks, it is not an exact sampling algorithm.

Conclusions
In summary, we have solved the problem of constrained graphicality when degree correlations are specified, developing an exact algorithm to construct and sample graphs with a specified JDM. A JDM specifies the number of edges that occur between degree classes of nodes (nodes of given degrees), and thus completely determines all pairwise degree correlations in its realizations. Our algorithm is guaranteed to successfully build a random JDM sample in polynomial time, systematically, and without backtracking. It is also guaranteed to be able to build any of the graphical realizations of a JDM. Each graph is constructed independently and thus there are no correlations between samples. Although the algorithm does introduce a sample bias, the relative probability for the construction of each sample is computable, which allows the use of weighted averages to obtain unbiased sampling (importance sampling). However, importance sampling is only exact in the limit of an infinite number of samples. This raises the issue of convergence. The log-normal distribution of weights makes convergence slow, but for small-to medium-sized networks good accuracy can be achieved, and quantities computed as if from uniform sampling. Improving the speed of convergence is a challenging problem, partly because it depends on the constraining JDM, and will be addressed in future publications. Degree correlations in real-world systems have been widely observed. Social networks are known to be positively correlated, and the concept of assortativity was known to the sociological literature before it was employed in applied mathematics. Technological networks are also characterized by particular correlation profiles. Moreover, correlations significantly affect the dynamics of spatial processes, such as the spread of epidemics [3]. Thus, with our algorithm, one can model correctly complex systems of general interest with desired degree assortativity. For the first time, this enables the study of networks in which the correlations are not determined solely by the nodes' degrees. For instance, there exist many studies about social networks, consisting of a comparison between some specific real-world network and a randomized ensemble of networks with the same degree sequence or degree distribution. As social networks are scale-free, these studies often just sample the same sequence or the same type of power-law sequences to produce null-model results. However, social networks are assortative, while random scale-free networks are on average disassortative. Thus, the average correlations of scale-free networks make degree-sequence and degree-distribution sampling problematic if one is trying to consider a random model of a social network. Our method allows one to avoid this problem by directly imposing the correlations, rather obtaining only those imposed by the degree sequence.
Upper bounds on the computational complexity of our algorithm show that in the worst case it is cubic in the number of nodes. However, we provide a way to compute the expected worst-case complexity if the degree distribution of the networks considered is known. This shows that, for commonly studied cases such as scalefree networks, the maximum complexity is only of the order of N 2 , making the algorithm even more efficient.  Given a non-increasing graphical degree sequence , a random undirected graph that realizes  can be constructed by: (i) To each node, assign a number of stubs equal to its degree.
(ii) Choose a hub node i. Any node can in principle be chosen, for example, the node with the largest degree.
(iii) Create a set of forbidden nodes X, which initially contains only i.
(iv) Find the set of allowed nodes A to which i can be linked preserving the graphicality of the remaining construction process. To find A, first determine the fail degree κ using the method described below. Then A will consist of all nodes j X ∉ that have remaining degree greater than κ.
(v) Choose a random node m A ∈ and connect i to it.
(vi) Reduce the value of d i and d m in  by 1, and reorder it.
(vii) If m still has unconnected stubs, add it to the set of forbidden nodes X.
(viii) If i still has unconnected stubs, return to step (iv).
(ix) If nodes still have unconnected stubs, return to step (ii).
To determine the fail degree in a degree sequence  being sampled, build the residual degree sequence ′, by connecting the hub node i with remaining degree d i to the d 1 i − nodes with the largest degrees that are not in the forbidden set X and reducing the elements of  accordingly. Then, compute the graphicality test inequalitites. Each inequality potentially yields a fail-degree candidate, depending on the values of L k and R k . For each value of k there are only 3 possibilities: In case (a), the degree of the first non-forbidden node whose index is greater than k is the fail-degree candidate. In case (b), the degree of the first non-forbidden node whose index is greater than k and whose degree is less than k 1 + is the fail-degree candidate. In case (c), there is no fail-degree candidate. The sequence of candidate nodes is non-decreasing until the fail-degree is found. Thus, one can stop the calculation when either the current fail-degree candidate is less than the previous one, or when a case (a) happens.
This algorithm generates graph samples biasedly. However, the relative probability of generating a particular sample μ is where d i is the residual degree of node i when it is chosen as a hub, m is the total number of hubs used, and A i j is the allowed set for the jth link of hub i. Thus, an unbiased estimator for a network observable Q for any target distribution P is the weighted average and G k and Ḡ k are defined as follows. Let where δ is the Kronecker delta, and Ḡ is given by the recurrence relation To efficiently test the graphicality of a BDS :  (iv) Test if L R 1 1 ⩽ . If false, then stop;  is not graphical. If true, set k = 2 and continue.
(v) Test if L R k k ⩽ . If false, then stop;  is not graphical. If true, increase k by one and repeat. Continue until k N 1 = − , then stop;  is graphical.
Given a graphical BDS of integer pairs  in lexicographic order, a random directed graph that realizes  can be constructed by: (i) Assign in-stubs and out-stubs to each node according to its degrees.
(ii) Define as current hub the lowest-index node i with non-zero out-degree.
(iii) Create a set of forbidden nodes X, which initially contains i and all nodes with zero in-degree.
(iv) Find the set of allowed nodes A to which an out-stub of i can be connected without breaking graphicality.
To find A, first determine the fail in-degree κ using the method described below. Then A will consist of all nodes j X ∉ that have remaining in-degree greater than κ.
(v) Choose a random node m A ∈ and connect an out-stub of i to one of its in-stubs.
(vi) Reduce the value of d i + and d m − in  by 1, and reorder it accordingly.
(vii) Add m to the set of forbidden nodes X.
(viii) If i still has unconnected out-stubs remaining, return to step (iv).
(ix) If nodes still have unconnected out-stubs, return to step (ii).
The following simple procedure can be used to efficiently find the fail-in-degree in step (iv) of the sampling algorithm. (v) Find the first non-forbidden node in ′ whose index is greater than k.
(vi) Identify this node in the original BDS . Its in-degree is the fail-in-degree. Stop.
As in the case of the sampling algorithm for undirected graphs, this algorithm generates directed graph samples biasedly. However, an unbiased estimator for a network observable Q for any target distribution P is the weighted average given by equation (A.6). In this case the weights are where ν is the total number of hubs used, A i j | | is the size of the allowed set immediately before placing the jth connection coming from the ith hub, and d i + is the out-degree of the ith node chosen as a hub. Note that, unlike the case for undirected networks, there is no factorial combinatorial factor in the weights. This is because while the particular sequence of hub nodes chosen depends on the links placed, every node with non-zero out-degree will be selected, sooner or later, as the hub. Therefore, all the samples produced would have an extra, identical, multiplicative factor of As only the relative probabilities are needed for estimating an observable, and this factor is the same for every possible sample, it is eliminated from the formula for the weights.

Appendix B. An explicit example
To illustrate the sampling mechanism and the difference between weighted and unweighted estimation, we consider the realizations of the JDM B.1. Unweighted estimate To calculate the theoretical results for the unweighted case, we need to consider the probability with which our algorithm generates each degree-spectra matrix from J. To this purpose, first note that there are several degreespectra matrices whose realizations are all pentagon graphs. Also, all the hexagon and bow tie graphs have the same degree-spectra matrix This allows us to compute just the probability of generating S, as all the other matrices will yield the same contribution to c P 2 〈 〉 and c P 3 〈 〉 . Our method chooses the elements of the degree-spectra matrix S being created in a systematic way, node by node. As there are no nodes of degree 1, all the element in the first row of the matrix are fixed to 0. Then, the first element to choose is S 2,1 , that is, the number of edges between node 1 and nodes of degree 2. The possible choices for this element are 0, 1, and 2. Choosing 0 or 2 will result necessarily in a degree-spectra matrix whose realizations are all pentagon graphs. In fact, from figure B1 one can see that, amongst the realizations of J, the pentagon graphs are the only ones in which a node of degree 2, such as node 1, has either no edges or 2 edges with nodes of degree 2. Thus, choosing the value of S 2,1 with uniform probability, at this stage one generates pentagon graphs with probability 2 3.
The remaining choice, S 1 2,1 = , happens with probability 1 3. In this case, S 3,1 is forced to be 1, since the elements in the first column of S must sum up to the degree of the node 1, which is 2. The next element to determine is then S 2,2 . Similarly to the previous case, the possible values are 0, 1, and 2. Choosing 0 or 2 will always result in pentagon graphs, whose probability of being generated increases by 1 3 · 2 3 2 9 = . Choosing S 1 2,2 = , which occurs with total probability 1 3 · 1 3 1 9 = , forces S 1 3,2 = . The next value to determine is that of S 2,3 . As before choosing 0 or 2 yields pentagon graphs, whose total probability of being generated increases by 1 3 · 1 3 · 2 3 2 27 = . The choice of S 1 2,3 = , which has a total probability 1 3 · 1 3 · 1 3 1 27 = of happening, implies that S 1 2,3 = . Then, the degree-spectra matrix being built can only be S HB . In fact, as it is evident from figure B1, the only graphs realizing J in which at least 3 nodes of degree 2 are linked exactly to one other node of degree 2 and one of degree 3, are hexagon and bow tie graphs.

= + =
Then, knowing that S HB is sampled with probability 1 27, and the remaining degree-spectra matrices always yield pentagon graphs, it is

B.2. Weighted estimate
In order to obtain an analytical result for the weighted estimate, rather than computing the probability of occurrence of each degree-spectra matrix as sampled by our algorithm, we need to compute their actual number. From the previous section, we already know that all the hexagon and bow tie graphs come from the same degreespectra matrix S HB , which is unique. Then, we only need to compute the number of degree-spectra matrices corresponding to pentagon graphs. To do so, remember that the first choice in the construction of a degree-spectra matrix from J is the value of the element S 2,1 . If S 0 2,1 = or S 2 2,1 = , then we are guaranteed to get a pentagon graph. However, while each of these two choices fixes the value of S 3,1 , we are still free to select a value for the next 'free' element, S 2,2 .
If S 0 2,1 = , the allowed values for S 2,2 are 1 and 2. Choosing 2 fixes all the other elements of the degreespectra matrix. Conversely, choosing 1 results in S 2,3 still to be determined. Its possible values are 1 and 2. Thus, there are 3 different degree-spectra matrices with S 0 2,1 = . If, instead, S 2 2,1 = , the situation is very similar to the first case. The possible choices for S 2,2 are 0 and 1. Choosing 0 fixes the entire matrix; choosing 1 requires to select a value for S 2,3 , which can be either 0 or 1. Thus, there are 3 matrices with S 2 2,1 = . The third possibility of S 1 2,1 = still allows degree-spectra matrices corresponding to a pentagon graph. Similarly to the previous case, the simplest way to construct one is to impose S 0 2,2 = or S 2 2,2 = . In both cases, one must then choose a value for S 2,3 . The possibilities are 1 and 2, if S 0 2,2 = , or 0 and 1, if S 2 2,2 = . Any choice for S 2,3 fixes all the remaining elements of the matrix. Finally, it is still possible to choose S 1 2,1 = and S 1 2,2 = , and still construct matrices corresponding to a pentagon graph. The choice is again on S 2,3 . Choosing S 0 2,3 = or S 2 2,3 = fixes all the other elements of the matrix, whose realizations will be pentagon graphs. Imposing S 1 2,3 = , instead results in the matrix S HB , exhausting all possibilities. This shows that there are 6 different matrices with S 1 2,1 = that generate pentagon graphs.
The decisional tree we just described is shown in figure B3 as a visual aid. In summary, there are 12 different degree-spectra matrices that realize J and whose realizations are always pentagon graphs. Knowing c P

B.3. Numerical verification
To validate our algorithm against the analytical results presented in the two sections above, we performed extensive numerical simulations, generating 10 4 degree-spectra matrices, and 10 4 samples per matrix, for a total of 10 8 graphs. For each graph generated, we saved the average local clustering coefficients for nodes of both degrees. Then, we obtained both weighted and unweighted results by averaging the data first naively, and then with a proper use of the weights according to equation (4). The results, shown in table B.1, show that the weighted averages obtained using our algorithm converge to the correct result. Also, the difference between weighted and unweighted results can be appreciated even when it is quite small, as in our example. This illustrate the sensitivity of our method, as well as the necessity of using proper sampling when performing this kind of studies.