Scale-Free Networks, in depth: preferential attachment in linear time, the exact Barabási-Albert degree law, fitting and testing power laws properly, the rarity debate, and robustness to random versus targeted failure

Preferential attachment and a linear-time generator

The Barabási-Albert (BA) model has one parameter m. Start with a small seed graph. Add nodes one at a time; each new node creates m edges to existing nodes, choosing node i with probability proportional to its current degree k_i. Rich nodes get richer, and the oldest nodes, which have been collecting links the longest, become the hubs.

Naively, each attachment means computing a degree-weighted random choice over all n nodes, so the generator costs O(n²). The standard trick removes that: keep a list containing every edge endpoint, once per edge end. A node with degree k appears exactly k times, so choosing a uniformly random element of the list is choosing a node with probability proportional to degree. Each step appends 2m entries, and the whole generator is O(nm).

import random
from collections import Counter

def ba_graph(n, m, seed=0):
    rng = random.Random(seed)
    edges = []
    targets = list(range(m))      # first new node links to m seed nodes
    repeated = []                 # every edge endpoint, once per edge end
    for v in range(m, n):
        for t in targets:
            edges.append((v, t))
            repeated += [v, t]
        chosen = set()
        while len(chosen) < m:    # m distinct, degree-proportional picks
            chosen.add(rng.choice(repeated))
        targets = list(chosen)
    return edges

def degrees(n, edges):
    d = [0] * n
    for a, b in edges:
        d[a] += 1
        d[b] += 1
    return d

networkx's barabasi_albert_graph(n, m, seed) uses the same repeated-endpoint trick. Note the set: sampling without duplicates means the probabilities are only approximately proportional to degree within a single step, which is the usual, accepted approximation. Running it with n = 100,000 and m = 2 gives 199,996 edges, a mean degree of 4.0 (always 2m), and a maximum degree of 1,029, roughly 250 times the mean. In an Erdős-Rényi graph of the same size and mean degree the largest degree would be around 15.

The exact degree distribution

For the BA model the stationary degree distribution is known exactly (Dorogovtsev, Mendes and Samukhin; Bollobás, Riordan, Spencer and Tusnády for the rigorous version):

P(k) = 2m(m+1) / (k(k+1)(k+2)) for k ≥ m

For large k this behaves like 2m²k^-3, so γ = 3. The formula also telescopes: summing over all k ≥ m gives exactly 1, a good sanity check if you implement it. With m = 2 it is 12 / (k(k+1)(k+2)): P(2) = 0.5, P(3) = 0.2, P(4) = 0.1, P(10) = 0.0091. The simulated graph above measured 0.4995, 0.2001, 0.0997 and 0.0088, which is the kind of agreement you should demand of your own generator before trusting anything built on it.

Two consequences follow from γ = 3. The mean degree is finite (it is 2m), but the second moment ⟨k²⟩ grows with network size, logarithmically for γ = 3 and as a power of n for 2 < γ < 3. Many processes depend on ⟨k²⟩, which is why hubs dominate spreading and percolation. And the largest degree grows roughly like m√n, so hubs keep getting bigger as the network grows rather than settling at a typical size.

Fitting a power law without fooling yourself

The most common mistake in this area is to plot a degree histogram on log-log axes, fit a straight line by least squares and report the slope as γ. That procedure is biased, gives no measure of whether a power law is plausible at all, and is sensitive to binning. Lognormal and stretched-exponential distributions also look nearly straight over one or two decades. Clauset, Shalizi and Newman (2009) set out the method now considered standard:

1. Estimate γ by maximum likelihood for each candidate lower cutoff k_min, since power-law behaviour usually holds only in the tail.
2. Choose k_min as the value that minimises the Kolmogorov-Smirnov distance between the data above it and the fitted model.
3. Test goodness of fit with a bootstrap: generate synthetic data from the fitted model, refit, and see how often the synthetic KS distance exceeds the observed one. A small p-value (they suggest below 0.1) rules out the power law.
4. Compare with alternatives (lognormal, exponential, power law with cutoff) using a likelihood-ratio test, which tells you which is favoured and whether the difference is significant.

import math

def mle_gamma(ks, kmin):
    """Discrete power-law MLE, using the continuous approximation of
    Clauset et al.: good when kmin is not too small (roughly >= 6)."""
    tail = [k for k in ks if k >= kmin]
    return 1 + len(tail) / sum(math.log(k / (kmin - 0.5)) for k in tail)

d = degrees(100_000, ba_graph(100_000, 2))
for kmin in (4, 6, 10, 20):
    print(kmin, round(mle_gamma(d, kmin), 2))
# 4 2.68   6 2.79   10 2.88   20 2.96

On the BA graph whose true asymptotic exponent is 3, the estimate climbs from 2.68 at k_min = 4 to 2.96 at k_min = 20. That is not noise: the exact law 12/(k(k+1)(k+2)) only approaches k^-3 for large k, and the approximation in the estimator is poorest at small k_min. It is a concrete demonstration of why k_min must be chosen by a principled rule and reported alongside γ. In practice, use the Python powerlaw package (Alstott, Bullmore and Plenz), which implements this procedure: powerlaw.Fit(data, discrete=True) chooses k_min and estimates the exponent, and fit.distribution_compare('power_law', 'lognormal') returns the log-likelihood ratio and its p-value.

Variants: what changes the exponent

The BA model fixes γ = 3, but real heavy-tailed networks report a range of exponents, and small changes to the attachment rule explain much of that range. Three variants are worth knowing because they show which ingredient produces which effect.

• Initial attractiveness. Attach with probability proportional to k_i + A instead of k_i. Dorogovtsev, Mendes and Samukhin showed the exponent becomes γ = 3 + A/m, so the tail can be made steeper by giving every node a baseline appeal; with negative A (down to -m) it becomes shallower, towards 2.
• Non-linear attachment. With probability proportional to k_i^α, Krapivsky, Redner and Leyvraz found that only α = 1 gives a pure power law. Sub-linear α < 1 gives a stretched-exponential tail; super-linear α > 1 gives a winner-takes-all network in which one node gathers a finite fraction of all links.
• Fitness. In the Bianconi-Barabási model each node has a fitness η_i and attracts links in proportion to η_ik_i, so a late but fit node can overtake old hubs, which pure BA never allows.

The lesson for modelling is that a power law is a fragile outcome: it needs growth and attachment that is linear in degree. If either is missing, expect a heavy tail of some other shape, which is one reason real data so often fits a lognormal at least as well.

Are scale-free networks rare?

When Broido and Clauset (2019) applied these tests to nearly a thousand real-world networks from many domains, they found strong evidence for a power-law degree distribution in only a small minority; lognormals fitted as well or better in many cases, and social networks were at best weakly scale-free. Their paper was titled Scale-free networks are rare and it prompted a lively response from Barabási and others, who argued that the criteria were too strict for finite, noisy data and that heavy tails, not exact power laws, are what matter.

The practical lesson does not depend on who wins that argument. Heavy-tailed degree distributions are very common, and many of the important consequences (hubs, large ⟨k²⟩, fragility under targeted attack) come from the heavy tail itself. Whether the tail is exactly a power law with a specific γ is a much stronger claim, and you should only make it if the tests above support it. Say "heavy-tailed, consistent with a power law above k_min = 12 with γ = 2.4 ± 0.1, lognormal not excluded" rather than "scale-free".

Hubs, robustness and spreading

Albert, Jeong and Barabási (2000) pointed out a striking asymmetry: scale-free networks tolerate random failures well, because a random node is almost always a low-degree one, but break quickly when the hubs are removed deliberately. The generator above makes this easy to check with a union-find over surviving edges, measuring the fraction of all nodes in the largest connected component:

Removing a fifth of the nodes at random leaves a giant component holding 78% of all nodes, nearly all of the 80% that survive (about 97% of them). Removing the top fifth by degree shatters it entirely. The top 1% of nodes in this graph hold 11.8% of all edge endpoints, which is the leverage an attacker, or a well-targeted intervention, exploits. The same arithmetic drives spreading: with γ ≤ 3 the epidemic threshold of the SIS model vanishes as the network grows (Pastor-Satorras and Vespignani, 2001), covered in epidemic models on networks. For the percolation theory behind giant components, see percolation theory, and for engineering-style reliability analysis, network reliability.

Operational consequences

• Graph processing. Hubs make work skewed: one vertex's adjacency list can be millions of entries, so vertex-partitioned systems get stragglers. Split high-degree vertices (edge partitioning, vertex cuts) and process hub neighbourhoods in chunks.
• Sampling. Random-walk and BFS samples over-represent hubs; random-node samples miss them. Correct for degree bias or report it.
• Ranking. Degree-based and eigenvector-style rankings are dominated by hubs; PageRank is the classic way to weight links by the importance of their source rather than raw count.
• Modelling. If you need a null model with the same degree sequence as your data, use the configuration model or degree-preserving edge swaps rather than BA, which fixes γ = 3 and has almost no clustering. Real networks also tend to be small worlds with high clustering; small-world networks covers that property.
• Testing. Validate any generator against the exact P(k) above and check 2×edges = sum of degrees before trusting downstream results.

What to do next

1. Run the generator at n = 100,000, m = 2 and reproduce P(2), P(3), P(4) and the maximum degree.
2. Fit your own network's degree sequence with powerlaw.Fit(..., discrete=True), report k_min, γ and the lognormal comparison, and stop calling it scale-free if the tests do not support it.
3. Reproduce the robustness table, then try removing nodes by betweenness instead of degree.
4. Profile one graph job for degree skew and add hub splitting if a few partitions dominate.
5. Replace any BA null model in your analysis with a degree-preserving randomisation.