Barabási-Albert Model#

The Barabási-Albert (BA) model generates scale-free networks through preferential attachment: newly added nodes connect preferentially to nodes that are already highly connected. This rich-get-richer mechanism produces degree distributions following a power law, observed empirically in the Web, citation networks, and biological interaction networks.

Formal Definition#

Definition (BA Model). Start with a seed graph \(G_0\) on \(m_0 \geq m\) nodes. Construct \(G_1, G_2, \ldots, G_n\) iteratively:

At each step \(t = 1, 2, \ldots\):

  1. Add a new node \(v_t\).

  2. Connect \(v_t\) to exactly \(m\) distinct nodes in \(G_{t-1}\). Node \(i\) is selected with probability $\(\pi_i = \frac{k_i}{\sum_{j \in G_{t-1}} k_j}\)\( where \)k_i\( is the current degree of node \)i$.

At time \(t\), the total degree is \(\sum_j k_j = 2(|E(G_0)| + mt) \approx 2mt\) for large \(t\), so: $\(\pi_i = \frac{k_i}{2mt}\)$

This rule is called preferential attachment: nodes with higher degree attract disproportionately more new links.

Degree Distribution via Mean-Field Theory#

Treating \(k_i(t)\) as a continuous quantity, the expected increment per step is: $\(\frac{\partial k_i}{\partial t} = m \cdot \pi_i = \frac{k_i(t)}{2t}\)$

Node \(i\) enters at time \(t_i\) with degree \(k_i(t_i) = m\). Solving the ODE: $\(k_i(t) = m \left(\frac{t}{t_i}\right)^{1/2}\)$

Deriving \(P(k)\). Assuming arrival times \(t_i\) are uniform on \([m_0, t]\): $\(P(k_i(t) < k) = P\!\left(t_i > \frac{m^2 t}{k^2}\right) = 1 - \frac{m^2}{k^2}\)$

Differentiating yields the stationary degree distribution: $\(\boxed{P(k) = 2m^2\, k^{-3}, \qquad k \geq m}\)$

The exponent \(\gamma = 3\) is universal for the BA model. Networks with \(P(k) \propto k^{-\gamma}\) are called scale-free; real-world networks typically have \(\gamma \in [2, 3]\).

import numpy as np
import networkx as nx
import matplotlib.pyplot as plt
from collections import Counter

np.random.seed(42)
def barabasi_albert(n, m, seed=42):
    """
    Build a BA graph via preferential attachment.

    Uses a repeated-node (stubs) list so each new node's
    target selection runs in expected O(m) time:
    sampling uniformly from stubs is equivalent to
    sampling proportionally to degree.

    n : total number of nodes
    m : edges each new node attaches
    """
    rng = np.random.default_rng(seed)
    G = nx.complete_graph(m)          # seed graph K_m

    stubs = [v for v, d in G.degree() for _ in range(d)]

    for new_node in range(m, n):
        G.add_node(new_node)
        targets = set()
        while len(targets) < m:
            t = stubs[rng.integers(len(stubs))]
            if t != new_node:
                targets.add(t)
        for t in targets:
            G.add_edge(new_node, t)
            stubs.extend([new_node, t])

    return G


G_ba = barabasi_albert(n=500, m=2)
print(f'Nodes: {G_ba.number_of_nodes()},  Edges: {G_ba.number_of_edges()}')
fig, ax = plt.subplots(figsize=(9, 7))
deg_map = dict(G_ba.degree())
sizes = [15 * deg_map[v] for v in G_ba.nodes()]
pos = nx.spring_layout(G_ba, seed=42)
nx.draw_networkx(
    G_ba, pos=pos, ax=ax,
    node_size=sizes, node_color='steelblue',
    edge_color='lightgray', alpha=0.8, with_labels=False
)
ax.set_title('Barabasi-Albert Graph  (n=500, m=2)\nNode size proportional to degree', fontsize=13)
ax.axis('off')
plt.tight_layout()
plt.show()
deg_seq = [d for _, d in G_ba.degree()]
counts = Counter(deg_seq)
ks = np.array(sorted(counts))
pk = np.array([counts[k] for k in ks], dtype=float) / len(deg_seq)

# Fit P(k) ~ C * k^{-gamma} on log-log scale
mask = (ks >= 2) & (pk > 0)
gamma_fit, logC = np.polyfit(np.log(ks[mask]), np.log(pk[mask]), 1)
print(f'Estimated gamma = {-gamma_fit:.3f}  (theoretical = 3.0)')

fig, ax = plt.subplots(figsize=(7, 5))
ax.loglog(ks, pk, 'o', color='steelblue', label='Empirical')
k_fit = np.linspace(ks.min(), ks.max(), 200)
ax.loglog(k_fit, np.exp(logC) * k_fit**gamma_fit, 'r--',
          label=f'Power-law fit (gamma={-gamma_fit:.2f}, theory=3.0)')
ax.set_xlabel('Degree k')
ax.set_ylabel('P(k)')
ax.set_title('Degree Distribution (log-log scale)')
ax.legend()
plt.tight_layout()
plt.show()