Large deviation and anomalous fluctuations scaling in degree assortativity on configuration networks

TL;DR

This paper investigates the probability distribution of degree assortativity in networks, revealing large deviation principles and anomalous fluctuation scaling, especially in heterogeneous scale-free networks.

Contribution

It introduces a multicanonical Monte Carlo method to analyze the full distribution of degree assortativity and establishes a large deviation principle with a novel scaling exponent.

Findings

The distribution obeys a large deviation principle with a convex rate function.

The scaling exponent $\xi$ equals 1 for Poisson graphs and varies for scale-free networks.

Fluctuations exhibit anomalous scaling in highly heterogeneous networks.

Abstract

By constructing a multicanonical Monte Carlo simulation, we obtain the full probability distribution $\rho_N(r)$ of the degree assortativity coefficient $r$ on configuration networks of size $N$ by using the multiple histogram reweighting method. We suggest that $\rho_N(r)$ obeys a large deviation principle, $\rho_N \left(r-r_N^* \right) \asymp {e^{ - {N^\xi }I\left( {r- r_N^* } \right)}}$, where the rate function $I$ is convex and possesses its unique minimum at $r=r_N^*$, and $\xi$ is an exponent that scales $\rho_N$'s with $N$. We show that $\xi=1$ for Poisson random graphs, and $\xi\geq1$ for scale-free networks in which $\xi$ is a decreasing function of the degree distribution exponent $\gamma$. Our results reveal that the fluctuations of $r$ exhibits an anomalous scaling with $N$ in highly heterogeneous networks.

Equations

$$r=\frac{M^{-1}{\sum }_i{j}_i{k}_i-{\left[M^{-1}{\sum }_i\frac{1}{2}\left(j_i+k_i\right)\right]}^2}{M^{-1}{\sum }_i\frac{1}{2}\left(j_i^2+k_i^2\right)-{\left[M^{-1}{\sum }_i\frac{1}{2}\left(j_i+k_i\right)\right]}^2},$$

$$F_{ijk\ell ;[1]}:(i,j)$$

$$F_{ijk\ell ;[2]}:(i,j)$$

$$F_{ijk\ell ;[3]}:(i,k)$$

$$P_{acc}=\min \left\{1,e^{-\beta \Delta r}\frac{|{\Phi }_{\mathbf{A}}|}{|{\Phi }_{F\mathbf{A}}|}\right\},$$

$$|{\Phi }_{\mathbf{A}}|$$

$$|{\Phi }_{F\mathbf{A}}|=|{\Phi }_{\mathbf{A}}|+\Delta |\Phi |,$$

$$p_i\left(r\right)={\rho }_N\left(r\right)\frac{e^{-{\beta }_{i}r}}{Z_i},$$

$$p_i\left(r\right)\mathrm{d}r=\frac{N_i(r)}{n_i}.$$

$${\rho }_N\left(r\right)\mathrm{d}r=\frac{N_i(r)Z_i}{n_i{e}^{-{\beta }_{i}r}}.$$

$${\rho }_N\left(r\right)\mathrm{d}r=\frac{{\sum }_{i=1}^R{N}_i\left(r\right)}{{\sum }_{j=1}^R{n}_j{Z}_j^{-1}{e}^{-{\beta }_{j}r}},$$

$$Z_k={\int }_r{\rho }_N\left(r\right)e^{-{\beta }_{k}r}\mathrm{d}r={\int }_r\frac{{\sum }_{i=1}^R{N}_i\left(r\right)}{{\sum }_{j=1}^R{n}_j{Z}_j^{-1}{e}^{\left({\beta }_k-{\beta }_j\right)r}}\mathrm{d}r.$$

$$⟨r^n⟩=\frac{\int r^n\rho \left(r\right)\mathrm{d}r}{\int \rho \left(r\right)\mathrm{d}r}.$$

$$P\left(G\right)=\frac{1}{Z}{e}^{-H\left(G\right)},$$

$$H\left(G\right)=\sum _i{\theta }_i{k}_i\left(G\right)=\sum _{i<j}\left({\theta }_i+{\theta }_j\right)A_{ij},$$

$$Z=\sum _G{e}^{-H\left(G\right)}=\prod _{i<j}\left(1+e^{-{\theta }_i-{\theta }_j}\right).$$

$$P\left(G\right)=\prod _{i<j}{p}_{ij}^{A_{ij}}{\left(1-p_{ij}\right)}^{1-A_{ij}},$$

$$p_{ij}=\frac{x_i{x}_j}{1+x_i{x}_j},$$

$$⟨k_i⟩=\sum _{j\ne i}{p}_{ij}=\sum _{j\ne i}\frac{x_i{x}_j}{1+x_i{x}_j}.$$

Full text

I Introduction

Over the past two decades, we have witnessed the success of complex networks in describing the pattern discovered ubiquitously in real world [1], such as community structure and scale-free structure, and modelling many dynamical processes in nature [2], such as synchronization [3, 4], epidemic spreading [5], opinion formation [6], etc [7]. In particular, how to characterize the structural features of complex networks is essential not only for uncovering the organizational principles of real systems, but also for understanding and controlling the dynamical processes on them [8, 9, 10].

An important feature in complex networks is so-called degree assortativity, which quantifies the tendency of nodes to be connected to other nodes of similar degree. A networks is called assortative if nodes with high degree preferably connect to other nodes with high degree, and dissortative if nodes with high degree are linked to nodes with low degree. Technical and biological networks have been found to be dissortatively mixed, while social networks show assortative correlations [11, 12, 13]. It was shown that, on the one hand, degree correlations are key to many structural properties of networks, such as percolation [12, 14], mean distance [14], and robustness [13, 15]. On the other hand, degree correlations affect the properties of dynamical processes taking place on networks, such as epidemic spreading [16, 17, 18], stability against stimuli and perturbation [19, 20], and synchronization of oscillators [21].

In his seminal papers [12, 13], Newman introduced the assortativity coefficient $r$ to measure the degree correlation, which is defined as

[TABLE]

where $M$ is the number of edges, and $j_{i}$, $k_{i}$ are the degrees of the nodes at the ends of the $i$th edge, with $i=1,\cdots,M$. The assortativity coefficient $r$ is actually the Pearson’s correlation coefficient between the degrees of neighboring nodes, which is supposed to have natural bounds $r\in[-1,1]$. A network is assortative when $r>0$ and disassortative when $r<0$.

Most of previous works on this subject were performed on scale-free networks with power-law degree distributions $P(k)\sim k^{-\gamma}$ [22, 23, 24, 25, 26, 27, 28, 29]. It has been shown that, on the one hand, for degree distribution exponent $2<\gamma<4$ the assortativity coefficient $r$ is usually negative in finite-size networks. On the other hand, $r$ always decreases in magnitude as network size increases, and $r$ equals to zero in the infinite networks. Maslov et al. [22] have shown by using computer simulations that the degree dissortativity results from the restriction of at most one edge between any pair of nodes. Furthermore, Park and Newman [23] verified this result in theory. They proposed a grand canonical ensemble of graphs such that analytical calculation of degree correlations becomes feasible. Johnson et al. [24] proposed an alternative explanation for the phenomenon by information entropy, and they showed that the Shannon entropy is maximized at some negative value of assortativity coefficient $r$ for highly heterogeneous scale-free networks. Menche et al. [25] analyzed the maximally disassortative scale-free networks and found that the lower bound of $r$ approaches to zero as network size increases in a power-law way. Dorogovtsev et al. [26] also found the results in a specific class of recursive trees with power-law degree distribution. Yang et al. [27] derived analytically the lower bound of assortativity coefficient in scale-free networks. Similar phenomenon was also discussed in some related works [28, 29], although the authors therein argued the availability of the Pearson’s coefficient for measuring degree correlations in large-size heavy-tailed networks, and alternatively they proposed other measurements such as Kendall-Gibbons’ $\tau_{b}$ [28] and Spearman’s $\rho$ [29].

Previous works mainly focused on either the typical behavior of $r$, such as how the expected value of $r$ changes with network size and degree heterogeneity [22, 23, 24, 26], or how to obtain a class of specific networks with some atypical value of $r$ [14, 25, 27]. For an ensemble of random networks with a given degree sequence (i.e. configuration model), it is known that the assortativity coefficient $r$ varies from one network realization to another. An interesting question arises: what is the probability of generating a configuration network whose assortativity coefficient $r$ falls in an interval $[r,r+\mathrm{d}r)$? The question is equivalent to finding the probability distribution function $\rho_{N}$ of $r$ with network size $N$. For the purpose, we shall employ a statistical-mechanics inspired Monte Carlo (MC) method, multiple histogram reweighting (MHR) [30, 31], to fully sample $\rho_{N}$ over a wide range of $r$. The method is computationally efficient and enable us to cover rare-event tails with very low probabilities of $r$. Recently, the MHR method was applied to investigate the large deviation properties of the largest connected [32] or biconnected component [33], the diameters [34] for random graphs, and resilience of transportation networks [35] as well as power grids [36]. Related algorithms [37], for example, Wang-Landau algorithm [38], has been used to efficiently sample large spectral gap [39] and prescribed motif densities in networks [40], and rare trajectories in chaotic systems [41, 42].

To that end, we first build a canonical ensemble MC sampling by a random edge-swapping scheme [43] and then collect a series of histograms of $r$ at different inverse temperatures. Finally, $\rho_{N}(r)$ is obtained by using the MHR method. By implementing the method on the configuration models with Poisson degree distributions and power-law degree distributions, we find that for all the cases under consideration $\rho_{N}(r)$ is unimodal and its width becomes narrower as $N$ increases. The expected value of $r$ is negative and decays in magnitude as $N$ increases in a power-law way, as reported in previous literatures. The variance $\sigma_{r}^{2}$ of $r$ decreases in power-law form, $\sigma_{r}^{2}\propto 1/N^{\xi}$, with the increase of $N$ as well. For homogeneous networks such as Poisson random graphs, $\xi=1$ such that the fluctuation in $r$ is standard. Strikingly, for highly heterogeneous networks such as scale-free networks with $\gamma<3$, we have $\xi>1$ and thus the fluctuation scaling of $r$ with $N$ is anomalous. Moreover, we suggest that $\rho_{N}(r)$ obeys a large deviation principle [44], $\rho_{N}\left(r-r_{N}^{*}\right)\asymp{e^{-{N^{\xi}}I\left({r-r_{N}^{*}}\right)}}$, where ${I\left({r}\right)}$ is so-called large deviation rate function which plays a role of microcanonical entropy of the network configuration model [45, 46, 47, 48]. $r_{N}^{*}$ is the most probable value of $r$, and $\xi$ is just mentioned that is the exponent scaling the $\sigma_{r}^{2}$’s with $N$.

II Multi-Canonical ensemble Monte Carlo sampling

The configuration model is an ensemble of random graphs with a given degree sequence $\{k_{1},\cdots,k_{N}\}$, where $k_{i}$ is the degree of node $i$ and $N$ is the number of nodes. The model was formulated by Bollobás [49], inspired by Ref.[50]. It was popularized by Newman, Strogatz, and Watts [51], who realized that it is a useful and simple model for real-world networks. The configurations networks are generated as follows. Firstly, each node $i$ is assigned a given number of half-edges equal to its observed degree $k_{i}$, with $\sum\nolimits_{i=1}^{N}{k_{i}}$ assumed to be even. Each half-edge is then connected to a randomly chosen other half-edge to form an edge in the graph. Finally, all the self-loops and all the parallel edges between two different nodes are removed by an algorithm to reshuffle edges that ensures the degree distribution unchanged. It was pointed out that the algorithm produces a bias in resulting network configurations [52]. Such a bias can be eliminated by a refusal algorithm [53], but the latter is more computationally time-consuming. However, it does not produce any effect in our model whether algorithm is applied. This is because that the first generated network is only used as the starting point of Monte Carlo sampling introduced below. In the long time, the results do not sensitive to the initial configuration.

We consider a Markov Chain Monte Carlo (MCMC) algorithm in which we weight each network configuration ${\bf A}$ with a Boltzmann weight $r({\bf A})$, where ${\bf A}$ is the adjacency matrix of the underlying network whose entries are defined as $A_{ij}=A_{ji}=1$ if nodes $i$ and $j$ are connected and $A_{ij}=A_{ji}=0$ otherwise, and $r({\bf A})$ is the assortativity coefficient of the network ${\bf A}$. To perform the MCMC, we consider the elementary edge swap moves that preserve the degree distribution of the network. We consider four different nodes, $i,j,k,\ell$, and one of the three invertible moves ${\bf A}\to F{\bf A}$ in which the following edge swaps are performed [54, 55]

[TABLE]

In order to perform any of these three moves the initial two links between the four nodes must be present in the network while the final two links must be absent or vice versa as multiple edges between two different vertices are forbidden. Since not all moves are accepted by the algorithm the MCMC algorithm should take into account the fact that some network configurations might allow more moves than others. In Fig. 1(a), we show a simple graph of four nodes with two edges. There exist two possible configurations to move by edge swaps. However, if an additional edge is introduced between nodes 1 and 3 (see Fig. 1(b)), one of the resulting move configurations is forbidden since the parallel edges are present.

We indicate with $|\Phi_{\bf A}|$ the number of edge swaps allowed if starting from adjacency matrix ${\bf A}$. Each single allowable edge-swap move ${\bf A}\to F{\bf A}$ is accepted by the Metropolis probability which ensure unbiased sampling of the network configurations

[TABLE]

where $\beta$ is the inverse temperature played a role of a conjugated field acting on the assortativity coefficient $r$. Generally, for a larger $\beta$, the sampled networks prefers to smaller values of $r$, and thus $\beta$ can be used to adjust the bias on sampling assortativity coefficient. $\Delta r$ is the change in the assortativity coefficient $r$ due to the edge-swapping trial. Therefore at each step the algorithm selects a value of $\alpha\in\{1,2,3\}$ with uniform probability and draws four nodes $i<j<k<\ell$ until the move $F_{ijk\ell;[\alpha]}$ is allowed. It then accepts the allowed move with probability $P_{acc}$.

We note that $|\Phi_{\bf A}|$ admits the following expression

[TABLE]

This expression can be used to calculate $|\Phi_{\bf A}|$ at the beginning of the MCMC algorithm. In order to calculate how $|\Phi_{\bf A}|$ changes at each step of the MonteCarlo step it is more convenient to consider the expression

[TABLE]

Indeed using this expression one can just write

[TABLE]

where $\Delta|\Phi|$ can be calculated by considering only the terms that change in Eq. (5)

In fact, the term ${|\Phi_{\bf A}|}/{|\Phi_{F{\bf A}}|}$ in Eq.(3) is very close to one since ${|\Phi_{\bf A}|}$ is the order of square of the number of edges, $M^{2}$, and thus the deviation of ${|\Phi_{\bf A}|}/{|\Phi_{F{\bf A}}|}$ from one is the order of $1/M^{2}$. Therefore, dropping such a term is expected to generate not much effect to the results, but it is bound to improve computing efficiency significantly. We have tested several networks and found that the results are consistent whether the term exists or not.

Similar procedure was also used to study the relation between degree correlations and other topological features such as clustering coefficient [56] and percolation property [57]. For a given inverse temperature $\beta_{i}$, the probability density $p_{i}(r)$ of generating a network with the assortativity coefficient $r$ follows the Boltzmann distribution [58, 59, 60],

[TABLE]

where $\rho_{N}(r)$ is probability density function of $r$ we want to obtain, and ${Z_{i}}=\int{\rho_{N}\left(r\right){e^{-{\beta_{i}}r}}\mathrm{d}r}$ is the partition function (normalized factor) at the inverse temperature $\beta_{i}$. In practice, $p_{i}$ can be obtained by performing MC simulations at $\beta_{i}$. To that end, we build a histogram $N_{i}(r)$ of the number of times out of $n_{i}$ that an interval $\left[r,r+\mathrm{d}r\right)$ is observed, and thus we have

[TABLE]

In simulations, we have performed $N\times 10^{5}$ (with $N$ being the size of the underlying network) trials for edge swaps and the last $5N\times 10^{4}$ trials are used to count bins of histogram of $r$. Using Eq. (8), Eq. (7) can be rewritten as

[TABLE]

The MHR method takes advantage of collecting a series of histograms at nearby temperature overlap. We perform a series of $R$ MC simulations in the canonical ensemble corresponding to $R$ different inverse temperature $\beta_{i}$ with $i=1,\cdots,R$, where $\beta_{i}$ is chosen uniformly from the interval $\left[\beta_{min},\beta_{max}\right]$. The improved estimate for $\rho_{N}(r)$ is given by [61]

[TABLE]

where the partition function $Z_{j}$ can be found self-consistently by iterating the following equations,

[TABLE]

During the iterations for Eq.(II), we have used a rescaling of $Z$-values (divided all by the smallest) after each step to avoid an overall growth.

Once the $\rho_{N}(r)$ is obtained, we can compute the $n$th moment of the assortativity coefficient $r$,

[TABLE]

In particular, $\left\langle r\right\rangle$ is the expected value of $r$, and $\sigma_{r}^{2}=\left\langle{{r^{2}}}\right\rangle-{\left\langle r\right\rangle^{2}}$ is the variance of $r$.

III Poisson random graphs

We first consider the Poisson random graphs whose degree distribution follows $P\left(k\right)={{{e^{-\left\langle k\right\rangle}}{{\left\langle k\right\rangle}^{k}}}\mathord{\left/{\vphantom{{{e^{-\left\langle k\right\rangle}}{{\left\langle k\right\rangle}^{k}}}{k!}}}\right.\kern-1.2pt}{k!}}$ with average degree ${\left\langle k\right\rangle}=6$. In Fig. 2, we show the logarithm values of $\rho_{N}(r)$ for several different $N$. Using the MHR method, the probabilities as small as $e^{-200}\simeq 10^{-87}$ are easily accessible. As $N$ increases, the width of the distribution of $\rho_{N}(r)$ becomes narrower. The typical value $r_{N}^{*}$ of $r$, i.e. the most probable value of $r$ corresponding to the maximum in $\rho_{N}(r)$, is very close to zero. To investigate the size effect of $\rho_{N}(r)$ in more detail, we have computed the expected value $\left\langle r\right\rangle$ and the variance of $\sigma_{r}^{2}$ of $r$ as a function of $N$. We find that $\left\langle r\right\rangle$ is always negative for all the $N$’s and decays in magnitude with $N$. As shown in Fig. 3(a), the minus $\left\langle r\right\rangle$ can be well fitted linearly with $N$ in the log-log plot, $-\left\langle r\right\rangle\sim{N^{-\nu}}$, with the exponent $\nu=1.08$. In Fig. 3(b), we show that $\sigma_{r}^{2}$ decreases with $N$ in a power-law way as well, $\sigma_{r}^{2}\sim{N^{-\xi}}$, with the exponent $\xi=0.99$ that is very close to one. This implies that the fluctuation of $r$ on Poisson random graphs is inversely proportional to the system size $N$, in accordance with the central limit theorem.

Next, we want to check whether the $\rho_{N}(r)$ obeys a large deviation principle. To that end, we first make a shift $r_{N}^{*}$ in $r$ such that the locations of the maximum in $\rho_{N}(r-r_{N}^{*})$ coincide for all the $N$’s. We then scale the logarithm of $\rho_{N}(r-r_{N}^{*})$’s with $N^{-\xi}$ providing that the $\rho_{N}(r)$ obeys a Gaussian form around $r=r_{N}^{*}$. Thus, we suggest a form of $\rho\left(r-r_{N}^{*}\right)\asymp{e^{-{N^{\xi}}I\left({r-r_{N}^{*}}\right)}}$, where $I$ is the large deviation rate function that is convex and possesses its unique minimum at $r=r_{N}^{*}$. Finally, we make a shift on $I$ so that $I_{\min}=0$ at $r=r_{N}^{*}$, which is often done because only $I_{\min}=0$ makes sense for $N\to\infty$. This suggestion is verified in Fig. 4, in which one can see that all the curves for each $N$ coincide not only near $r_{N}^{*}$, but also far from $r_{N}^{*}$.

IV Scale-free networks

We now consider the case of scale-free networks whose degree distribution follows a power-law function, $P(k)=(\gamma-1)k_{0}^{\gamma-1}k^{-\gamma}$, where $k_{0}$ is the minimal degree, and $\gamma$ is degree distribution exponent. Here we focus on the range $2<\gamma<4$. The maximal degree $k_{\max}$ is chosen by a natural cutoff, ${k_{\max}}=\min\left({{k_{0}}{N^{1/(\gamma-1)}},N-1}\right)$ such that $\int_{k_{0}}^{\infty}P(k)\mathrm{d}k=1/N$. In Fig. 5, we shows the logarithm of $\rho_{N}(r)$ for three different $\gamma=2.3$ (a), 2.5 (b), 3.0 (c) and for five different $N$’s. It can easily seen that for all cases $\rho_{N}(r)$ are always unimodal. All the expected value of $r$ are negative, $\left\langle r\right\rangle<0$. This is especially obvious for smaller $\gamma$. With the increment of $N$, $\left\langle r\right\rangle$ moves to zero gradually. In Fig. 6(a), we show that $\left\langle r\right\rangle$ can be well fitted by the form of $-\left\langle r\right\rangle\sim{N^{-\nu}}$. The exponent $\nu$ is dependent on $\gamma$, which is $\nu=0.167$, 0.214, and 0.443 for $\gamma=2.3$, 2.5, and 3.0, respectively. The fluctuations of $r$, $\sigma_{r}^{2}$, obey the scaling law as well, $\sigma_{r}^{2}\sim{N^{-\xi}}$, as shown in Fig. 6(b). The exponent $\xi$ decreases as $\gamma$ increases, which is $\xi=1.59$, 1.28, and 0.99 for $\gamma=2.3$, 2.5, and 3.0, respectively. That is to say, for highly heterogeneous networks, they exhibit anomalously small fluctuations in $r$, since $\xi>1$ implies that the fluctuations decay with $N$ faster than the standard $1/N$ scaling.

In Fig. 7, we show the large deviation functions for scale-free networks. As mentioned before, the large deviation functions are obtained by $I\asymp-{{N^{-\xi}}}\ln\rho\left(r-r_{N}^{*}\right)$. As expected, all the data coincide for different $N$.

V Configuration network model with soft constraints

Finally, we shall compare the scaling behavior of the assortativity coefficient $r$ between two different ensembles of configuration model. The first one, as we studied before, is microcanonical, in which degree sequence $\{k_{1},\cdots,k_{N}\}$ are fixed. The second one is canonical ensemble that is easier to handle mathematically, and it is called the exponential random graph model in network science [58, 59, 60]. In the canonical ensemble, the hard constraints in microcanonical ensemble are softened by enforcing only as expected values, i.e. $\left\langle{{k_{i}}}\right\rangle={\bar{k}_{i}}$ for $i=1,\cdots,N$. The canonical probability of a graph ${G}$ is written as [58, 59, 60, 62, 63, 64]

[TABLE]

where $H$ is the graph Hamiltonian defined as

[TABLE]

and the normalizing quantity $Z$ is partition function that can be calculated exactly,

[TABLE]

Substituting Eq. (14) and Eq. (15) into Eq. (13), $P(G)$ can be written as the mass probability function of a Bernoulli-distributed binary random variable $A_{ij}$ (adjacency maxtrix),

[TABLE]

with success probability

[TABLE]

where $x_{i}=e^{-\theta_{i}}$ is called fugacity that can be obtained numerically by solving constraint equations,

[TABLE]

In Fig. 8, we compare the results of canonical scale-free model with those of the microcanonical scale-free model for three different values of $\gamma$. For each $\gamma$ and each $N$, we generate at least 5000 realizations of canonical configuration networks according to Eq. (17) to obtain mean value and variance of $r$, in which the expected values of node degrees are the same as the degree sequence in microcanonical configuration networks. In canonical model, one can see that both $-\left\langle r\right\rangle$ and $\sigma_{r}^{2}$ decay with power-law as $N$ increases. On the one hand, the values of $\left\langle r\right\rangle$ are almost independent of specific ensemble and share the same scaling exponent $\nu$. On the other hand, the values of $\sigma_{r}^{2}$ in canonical model are always larger than those in microcanonical model. This is especially obvious for smaller values of $\gamma$. The result is as expected because in canonical model the degree of each node is fluctuating from one network realization to another. For $\gamma=2.5$ and $\gamma=3$, the scaling exponents $\xi$ are almost the same in the two ensembles. However, for $\gamma=2.3$, $\xi\simeq 1.2$ in canonical model is less than 1.59 in the microcanonical model.

In Fig. 9(a) and Fig. 9(b), we show the scaling exponents $\nu$ and $\xi$ as a function of $\gamma$, respectively. In the two ensembles $\nu$ increases monotonically as $\gamma$ increases. When $\gamma<3$, $\nu$ are almost the same, and when $\gamma>3$, $\nu$ in canonical ensemble is slightly larger. However, $\xi$ changes with $\gamma$ in two different trends. When $\gamma>2.5$, $\xi$ in the two ensembles are almost the same, and remains constant around one when $\gamma>3$. For $\gamma<2.5$, $\xi$ in microcanonical model are obviously larger than those in canonical model. For example, for $\gamma=2.1$ we have $\xi=2.38$ in microcanonical model and $\xi=1.22$ in canonical model. From Fig. 9(b), one can see that when $\gamma>3$ the scale-free networks start to share the same scaling exponents as the Poisson-distributed random graphs. Intuitively, it seems to be relevant to the divergence of the second moment of the degree distribution on scale-free networks with $\gamma>3$. It may be hopeful to establish this possible connection in the exponential random graph models as it is easier to handle mathematically in the canonical ensemble. We have realized that in a recent paper [64], the authors used the two-star model [65] to study degree correlations between the nearest and next nearest neighboring nodes. They analytically calculated the degree assortativities and showed that they are nonmonotonic functions of the model parameters, with a discontinuous behavior at a first-order transition. However, in the work the authors did not observe a broad degree distribution such as power law form that are properties of many empirical networks. Therefore, it is still a challenging problem.

VI Conclusions

In summary, we have used the MHR method to obtain the probability distribution $\rho_{N}$ of the assortativity coefficient $r$ on configuration networks. This method enable us to obtain the rare-probability tails of $\rho_{N}(r)$ within the allowable computational time. We show that $\rho_{N}(r)$ satisfies a large deviation principle after a shift $r_{N}^{*}$ in $r$, $\rho_{N}\left(r-r_{N}^{*}\right)\asymp{e^{-{N^{\xi}}I\left({r-r_{N}^{*}}\right)}}$, in which $I(r)$ is the large deviation rate function that is convex and possesses its unique minimum at $r=r_{N}^{*}$. We find that $\xi=1$ in Poisson random graphs and scale-free networks with $\gamma>3$, indicating a normal fluctuations scaling of $r$ with $N$ in such networks, $\sigma_{r}^{2}\propto 1/N$. Interestingly, $\xi>1$ for $\gamma<3$, showing an anomalously fast decay in the fluctuation of $r$ as $N$ increases. Such an anomalous phenomenon in time-consuming observables have also been found in some other systems [66, 67, 68, 69, 70]. Furthermore, we show that in the canonical ensemble $\xi$ is slightly greater than one for $\gamma<2.5$ but is obviously less than that in the microcanonical model. This suggests that the anomaly in fluctuations of $r$ is not very significant in the canonical ensemble.

In the future, it is worthy investigating the joint distribution of assortativity coefficient $r$ and other topological observables, such as the average shortest path length, the largest eigenvalue of adjacency matrix or the second smallest eigenvalue of the Laplacian matrix, using the MHR method. This will surely deepen the understanding of the role of degree assortativity on dynamical precesses on configuration networks [16, 17, 18, 19, 20, 21].

Recently, we have noticed that large deviation theory has been used to uncover atypical structural and dynamical characteristics of complex networks, such as a first-order percolation transition subject to a rare initial damage [71, 72], a first-order phase transition in the condensation of node degrees [73], localization transitions [74, 75, 76] and optimal paths [77] of dynamical observables in random walk model , and epidemic extinction [78, 79, 80] and spin model [80, 81]. In the future, we believe that large deviation theory and related rare-event simulation methods may inspire more research works in network science.

Acknowledgements.

We acknowledge supports from the National Natural Science Foundation of China (Grant Nos. 11875069, 11975025, 12011530158, 61973001) and the Key Scientific Research Fund of Anhui Provincial Education Department (Grant No. KJ2019A0781))