复杂网络科学深度实战：从社团发现、模块度与 Louvain 社区检测到 Kuramoto 相振子与 SIR 传播动力学的完整工程链路

复杂网络科学把「谁连着谁」这件事变成可计算的几何对象。我们已经发过《复杂系统科学深度实战》，讲的是元胞自动机、伊辛相变、洛伦兹吸引子这类「涌现」的普适特征；但复杂系统往往把个体当成格子里的点，而真实世界里的脑、基因调控、社交平台、电网、互联网，真正起作用的是连接结构本身。本文是复杂系统的「第二视角」：用纯 NumPy 从零实现图的生成、社团发现、相位同步与传播动力学四条工程主线，全部代码可运行、数值可复现。

一、为什么需要图，而不是向量

机器学习习惯把样本表示成特征向量；但当样本之间天然有依赖（谁关注谁、哪条基因调控哪条、哪个路由连着哪个）时，向量会丢掉最关键的拓扑信息。图的两条黄金性质把真实系统分成两类：

• 小世界（small-world）：聚类系数 C 远高于随机图，而平均路径长度 L 又和随机图一样短。Watts–Strogatz 用「规则环 + 少量随机重连」生成它——重连让远距离被一两条捷径打通，本地三角结构却保留。
• 无标度（scale-free）：度分布近似幂律 P(k) ∼ k^(−γ)，存在极少数度极高的 hub。Barabási–Albert 用「偏好连接」生成它——新节点更愿连已经很火的节点，富者愈富。

二、图的生成与基本拓扑度量

下面一次性实现三种生成模型与三项度量（聚类系数、平均最短路径、度分布幂律指数），全部用邻接集合，无第三方图库。

import numpy as np

def watts_strogatz(N, k, p, rng):
    adj = [set() for _ in range(N)]
    for i in range(N):
        for j in range(1, k//2 + 1):
            t = (i + j) % N
            adj[i].add(t); adj[t].add(i)
    for i in range(N):
        for j in range(1, k//2 + 1):
            if rng.random() < p:
                t = (i + j) % N
                adj[i].discard(t); adj[t].discard(i)
                for _ in range(100):
                    w = rng.integers(N)
                    if w != i and w not in adj[i]:
                        adj[i].add(w); adj[w].add(i); break
    return adj

def barabasi_albert(N, m, rng):
    adj = [set() for _ in range(N)]
    targets = []
    for i in range(m):
        for j in range(i):
            adj[i].add(j); adj[j].add(i); targets.append(i); targets.append(j)
    for i in range(m, N):
        picked = set()
        while len(picked) < m:
            w = targets[rng.integers(len(targets))]
            picked.add(w)
        for w in picked:
            adj[i].add(w); adj[w].add(i)
        targets.extend([i]*m); targets.extend(list(picked))
    return adj

def erdos_renyi(N, p, rng):
    adj = [set() for _ in range(N)]
    for i in range(N):
        for j in range(i+1, N):
            if rng.random() < p:
                adj[i].add(j); adj[j].add(i)
    return adj

def avg_shortest_path(adj):
    N = len(adj); total = 0; cnt = 0
    for s in range(N):
        dist = [-1]*N; dist[s]=0; q=[s]; head=0
        while head < len(q):
            u = q[head]; head += 1
            for v in adj[u]:
                if dist[v] == -1:
                    dist[v] = dist[u]+1; q.append(v)
        for d in dist:
            if d > 0: total += d; cnt += 1
    return total/cnt

def avg_clustering(adj):
    N = len(adj); tot = 0.0
    for i in range(N):
        nb = list(adj[i]); k = len(nb)
        if k < 2: continue
        links = 0; s = set(nb)
        for a in range(k):
            for b in range(a+1, k):
                if nb[b] in s and nb[a] in adj[nb[b]]:
                    links += 1
        tot += 2.0*links/(k*(k-1))
    return tot/N

def degree_gamma(adj):
    degs = np.array([len(a) for a in adj], dtype=float)
    maxd = int(degs.max()); K = np.arange(1, maxd+1)
    P = np.array([np.mean(degs >= k) for k in K])
    mask = P > 0; K = K[mask]; P = P[mask]
    x = np.log(K); y = np.log(P)
    sel = K >= 3
    if sel.sum() >= 3:
        x = x[sel]; y = y[sel]
    slope, _ = np.polyfit(x, y, 1)
    return 1 - slope

rng = np.random.default_rng(0)
N = 400
ws = watts_strogatz(N, 4, 0.1, rng)
ba = barabasi_albert(N, 3, rng)
er = erdos_renyi(N, 6.0/(N-1), rng)
print("WS  clustering=%.3f path=%.3f" % (avg_clustering(ws), avg_shortest_path(ws)))
print("BA  clustering=%.3f path=%.3f gamma=%.2f" % (avg_clustering(ba), avg_shortest_path(ba), degree_gamma(ba)))
print("ER  clustering=%.4f path=%.3f" % (avg_clustering(er), avg_shortest_path(er)))

实测（N=400，固定随机种子）：

模型	聚类系数 C	平均路径 L	度指数 γ
Watts–Strogatz (k=4,p=0.1)	0.393	7.799	—
Barabási–Albert (m=3)	0.058	3.140	2.94
Erdős–Rényi (对照)	0.0167	3.620	—

关键结论：WS 的聚类系数是 ER 对照的约 23 倍（0.393 vs 0.0167），但平均路径只略长——这就是小世界。BA 的 γ=2.94 与理论值 3 几乎吻合，证明偏好连接确实长出幂律尾巴；同时 BA 的聚类系数比 ER 高，说明 hub 周围也保留了局部三角结构。

三、社团发现：模块度与谱聚类

很多网络天然分成内部稠密、外部稀疏的「社团」。度量好坏的标准是模块度（modularity）：

Q = (1 / 2m) · Σ_{i,j} [ A_ij − k_i k_j / (2m) ] · δ(c_i, c_j)

其中 m 是边数，k_i 是节点度，c_i 是节点所属社团。Q 越大，社团划分越「出乎意料地好」。谱聚类把 Q 的近似优化变成线性代数：取对称归一化拉普拉斯 L_sym = I − D^(−1/2) A D^(−1/2) 的第二小特征向量，按符号切分即可。

import numpy as np

def sbm(Nc, pin, pout, rng):
    N = Nc*2
    adj = [set() for _ in range(N)]
    comm = [0]*Nc + [1]*Nc
    for i in range(N):
        for j in range(i+1, N):
            p = pin if comm[i]==comm[j] else pout
            if rng.random() < p:
                adj[i].add(j); adj[j].add(i)
    return adj, comm

def modularity(adj, comm):
    N = len(adj); m = sum(len(a) for a in adj)/2.0
    k = np.array([len(a) for a in adj], dtype=float)
    cidx = {c:i for i,c in enumerate(sorted(set(comm)))}
    labels = np.array([cidx[c] for c in comm])
    Q = 0.0
    for i in range(N):
        for j in adj[i]:
            if labels[j] == labels[i]:
                Q += 1.0 - k[i]*k[j]/(2*m)
    return Q/(2*m)

def spectral_clustering(adj, n_clusters, rng):
    N = len(adj)
    A = np.zeros((N, N))
    for i in range(N):
        for j in adj[i]:
            A[i, j] = 1.0
    deg = A.sum(axis=1)
    Dinv = np.diag(1.0/np.sqrt(deg))
    Lsym = np.eye(N) - Dinv @ A @ Dinv
    w, V = np.linalg.eigh(Lsym)
    idx = np.argsort(w)
    F = V[:, idx[1:1+n_clusters]]
    norm = np.linalg.norm(F, axis=1, keepdims=True); norm[norm[:,0]==0] = 1.0
    F = F/norm
    centers = F[rng.choice(N, n_clusters, replace=False)].copy()
    for _ in range(50):
        d = np.linalg.norm(F[:,None,:] - centers[None,:,:], axis=2)
        lab = d.argmin(axis=1)
        newc = np.array([F[lab==c].mean(axis=0) if np.any(lab==c) else centers[c] for c in range(n_clusters)])
        if np.allclose(newc, centers): centers = newc; break
        centers = newc
    return lab

def agree(a, b):
    best = 0.0
    for flip in [0,1]:
        aa = a if flip == 0 else (1-a)
        best = max(best, np.mean(aa==b), 1-np.mean(aa==b))
    return best

rng = np.random.default_rng(1)
adj, comm = sbm(100, 0.8, 0.1, rng)
Qtrue = modularity(adj, comm)
lab = spectral_clustering(adj, 2, rng)
print("Q_true=%.3f spectral_agree=%.3f" % (Qtrue, agree(lab, np.array(comm))))

实测：真实划分的模块度 Q=0.491，谱聚类恢复的划分与真实社团完全吻合（agreement=1.000）。这说明只要社团间连边足够稀疏（pout=0.1），拉普拉斯谱就能把两个社群干净分开。

四、Louvain：贪心模块化优化

谱聚类要算全图特征分解，O(N³)，上亿节点吃不消。Louvain 是工业界主力：两层迭代——局部移动（把节点搬进让 ΔQ 最大的邻居社团）与聚合（把已得的社团压缩成超级节点再迭代），直到模块度不再上升。注意聚合时必须保留社区内部边的自环权重，否则高层会错误地合并社团。

import numpy as np

def louvain(adj, rng):
    N = len(adj)
    w = [dict() for _ in range(N)]
    for i in range(N):
        for j in adj[i]:
            w[i][j] = w[i].get(j, 0.0) + 1.0
    comm_of_orig = list(range(N))
    while True:
        cur = [dict(d) for d in w]
        m = sum(sum(d.values()) for d in cur) / 2.0
        k = np.array([sum(d.values()) for d in cur], dtype=float)
        comm = list(range(len(cur)))
        tot = k.copy()
        changed = True
        while changed:
            changed = False
            order = rng.permutation(len(cur))
            for i in order:
                ci = comm[i]
                neigh = {}
                for j, wt in cur[i].items():
                    cj = comm[j]; neigh[cj] = neigh.get(cj, 0.0) + wt
                ki = k[i]
                ki_in = neigh.get(ci, 0.0)
                dQ_remove = ki_in - tot[ci]*ki/(2*m)
                bestc = ci; bestgain = 0.0
                for cj, wij in neigh.items():
                    if cj == ci: continue
                    gain = (wij - tot[cj]*ki/(2*m)) - dQ_remove
                    if gain > bestgain + 1e-12:
                        bestgain = gain; bestc = cj
                if bestc != ci:
                    tot[ci] -= ki; tot[bestc] += ki
                    comm[i] = bestc; changed = True
        labels = sorted(set(comm)); newN = len(labels)
        if newN == len(cur):
            break
        uniq = {c:i for i,c in enumerate(labels)}
        comm_of_orig = [uniq[comm[c]] for c in comm_of_orig]
        new_w = [dict() for _ in range(newN)]
        for i in range(len(cur)):
            ci = uniq[comm[i]]
            for j, wt in cur[i].items():
                cj = uniq[comm[j]]
                new_w[ci][cj] = new_w[ci].get(cj, 0.0) + wt
        w = new_w
        if newN <= 1:
            break
    return comm_of_orig

rng = np.random.default_rng(2)
adj, comm = sbm(100, 0.8, 0.1, rng)
lou = louvain(adj, rng)
print("Louvain Q=%.3f n_comm=%d agree=%.3f" % (modularity(adj, lou), len(set(lou)), agree(np.array(lou), np.array(comm))))

实测：Louvain 得到 Q=0.487、2 个社团、与真实划分 agreement=1.000。它和谱聚类的模块度几乎相同（0.487 vs 0.491），但复杂度只有 O(N log N)，是大规模网络的事实标准。

五、标签传播：极简而强大

如果连 Louvain 的 O(N log N) 都嫌贵，还有更朴素的方法。标签传播（Label Propagation）：每轮让每个节点采纳邻居中最流行的标签，随机打破平局，直到收敛。

import numpy as np

def label_propagation(adj, rng, maxiter=100):
    N = len(adj); labels = np.arange(N)
    for it in range(maxiter):
        order = rng.permutation(N); changed = False
        for i in order:
            cs = [labels[j] for j in adj[i]]
            if cs:
                uniq, cnts = np.unique(cs, return_counts=True)
                mx = cnts.max()
                cands = uniq[cnts == mx]
                lab = cands[rng.integers(len(cands))]
                if lab != labels[i]:
                    labels[i] = lab; changed = True
        if not changed:
            break
    return labels

rng = np.random.default_rng(3)
adj, comm = sbm(100, 0.8, 0.1, rng)
lp = label_propagation(adj, rng)
print("LPA n_comm=%d agree=%.3f" % (len(set(lp)), agree(lp, np.array(comm))))

实测：标签传播同样恢复 2 个社团、agreement=1.000，且几乎是线性的。代价是结果依赖初始化与随机种子（可能退化成单一巨社团），所以工程上常跑多次取最稳定的划分。

三种社团发现法对比：

方法	复杂度	模块度	稳定性	适用规模
谱聚类	O(N³)	高	高	中小图
Louvain	O(N log N)	高（近似最优）	高	十亿级
标签传播	O(N)	中	随机	超大规模

六、Kuramoto 模型：相位同步的相变

社团是「静态结构」，同步是「动态秩序」。Kuramoto 模型把每个节点看成一个相位振子 θ_i，全连接耦合：

dθ_i/dt = ω_i + (K / N) · Σ_j sin(θ_j − θ_i)

ω_i 是各自的自然频率（取高斯分布），K 是耦合强度。用序参量 r = |(1/N) Σ_j exp(i θ_j)| 衡量同步程度：r≈0 表示各振子互不相让（无序），r≈1 表示锁相（同步）。临界耦合存在解析值 K_c = 2 / (π · g(0))，其中 g(0) 是频率分布在原点的值（高斯下 g(0)=1/(σ√(2π))）。

import numpy as np

def kuramoto_rss(N, K, sigma, T=200.0, dt=0.05, seed=0):
    rng = np.random.default_rng(seed)
    omega = rng.normal(0, sigma, N)
    theta = rng.uniform(-np.pi, np.pi, N)
    steps = int(T/dt)
    r_list = []
    for s in range(steps):
        dtheta = omega + (K/N)*np.sum(np.sin(theta[None,:] - theta[:,None]), axis=1)
        theta = theta + dt*dtheta
        z = np.mean(np.exp(1j*theta))
        r_list.append(abs(z))
    return np.mean(r_list[-int(20/dt):])

Nk = 200
print("Kuramoto r_ss vs K:")
for K in [0.5, 1.0, 1.6, 2.0, 4.0]:
    print("   K=%.1f r_ss=%.3f" % (K, kuramoto_rss(Nk, K, 1.0, seed=5)))

实测（σ=1，故 K_c ≈ 2·√(2π)/π ≈ 1.60）：

耦合 K	稳态序参量 r
0.5	0.094
1.0	0.128
1.6	0.349
2.0	0.751
4.0	0.968

在 K_c≈1.6 附近，r 从近乎 0 陡然跃升——这是典型的二阶相变。弱耦合下每个振子被自己的自然频率带跑，强耦合下全局涌现出一致的节律。这把「同步」从一个比喻变成了可量化的相变。

七、SIR 传播动力学：网络上的爆发阈值

信息、谣言、病毒都沿边扩散。SIR 模型把节点分成易感 S、感染 I、康复 R：每步感染者以概率 β 传染每个易感邻居，并以概率 γ 康复。在均匀混合人群里爆发条件是 R0 = β⟨k⟩/γ > 1，即 β_c^MF = γ/⟨k⟩。而在真实网络上，阈值由邻接矩阵最大特征值 λ1 决定：β_c = γ / λ1——因为 hub 让路径数暴涨。

import numpy as np

def build_A(adj):
    N = len(adj); A = np.zeros((N, N))
    for i in range(N):
        for j in adj[i]:
            A[i, j] = 1.0
    return A

def power_iter(A, iters=200, seed=0):
    rng = np.random.default_rng(seed)
    N = A.shape[0]; x = rng.normal(size=N); x /= np.linalg.norm(x)
    for _ in range(iters):
        x = A @ x; n = np.linalg.norm(x)
        if n == 0: break
        x /= n
    return float(x @ (A @ x))

def barabasi_albert(N, m, rng):
    adj = [set() for _ in range(N)]
    targets = []
    for i in range(m):
        for j in range(i):
            adj[i].add(j); adj[j].add(i); targets.append(i); targets.append(j)
    for i in range(m, N):
        picked = set()
        while len(picked) < m:
            w = targets[rng.integers(len(targets))]
            picked.add(w)
        for w in picked:
            adj[i].add(w); adj[w].add(i)
        targets.extend([i]*m); targets.extend(list(picked))
    return adj

def sir_sim(adj, beta, gamma, seed_frac=0.01, steps=400, seed=0):
    rng = np.random.default_rng(seed)
    N = len(adj); state = np.zeros(N, dtype=int)
    n0 = max(1, int(seed_frac*N))
    init = rng.choice(N, n0, replace=False); state[init] = 1; ever = state.copy()
    for t in range(steps):
        inf = np.where(state == 1)[0]
        if len(inf) == 0: break
        newI = []; newR = []
        for i in inf:
            for j in adj[i]:
                if state[j] == 0 and rng.random() < beta:
                    newI.append(j)
            if rng.random() < gamma:
                newR.append(i)
        for j in newI:
            if state[j] == 0: state[j] = 1; ever[j] = 1
        for i in newR:
            state[i] = 2
    return ever.sum()/N

rng = np.random.default_rng(6)
ba6 = barabasi_albert(500, 3, rng)
lam1 = power_iter(build_A(ba6), seed=1)
gamma = 0.1
print("BA lambda1=%.2f beta_c(network)=%.4f beta_c(MF)=%.4f" % (lam1, gamma/lam1, gamma/(2*3)))
for beta in [0.004, 0.012, 0.03]:
    fs = np.mean([sir_sim(ba6, beta, gamma, seed=i) for i in range(5)])
    print("   beta=%.3f final_size=%.3f" % (beta, fs))

实测（BA(500,m=3)，γ=0.1）：λ1=12.41，于是 网络阈值 β_c=0.0081，而均匀混合阈值 β_c^MF=γ/⟨k⟩=0.1/6=0.0167。最终感染规模：

感染率 β	是否过网络阈值	是否过 MF 阈值	最终规模
0.004	否	否	0.014（熄灭）
0.012	是	否	0.065（爆发）
0.030	是	是	0.434（大规模爆发）

注意 β=0.012 这一行：已经低于均匀混合阈值，却在无标度网络上仍然爆发。这正是 hub 的威力——少数超级节点把传播路径数放大到远超平均度所能解释的程度，让真实系统比「均匀混合」模型预测的更容易被点燃。

八、中心性：谁是重要的节点

社团回答「谁和谁一伙」，中心性回答「谁最关键」。两种经典定义：

• 特征向量中心性（eigenvector centrality）：x ∝ A x，重要节点的邻居也重要，等价于幂迭代求邻接矩阵主特征向量。
• PageRank：在随机游走视角下加入阻尼因子 d=0.85 与均匀回跳，避免死胡同并凸出权威节点：x = (1−d)/N · 1 + d · Aᵀ x / 出度。

import numpy as np

def build_A(adj):
    N = len(adj); A = np.zeros((N, N))
    for i in range(N):
        for j in adj[i]:
            A[i, j] = 1.0
    return A

def barabasi_albert(N, m, rng):
    adj = [set() for _ in range(N)]
    targets = []
    for i in range(m):
        for j in range(i):
            adj[i].add(j); adj[j].add(i); targets.append(i); targets.append(j)
    for i in range(m, N):
        picked = set()
        while len(picked) < m:
            w = targets[rng.integers(len(targets))]
            picked.add(w)
        for w in picked:
            adj[i].add(w); adj[w].add(i)
        targets.extend([i]*m); targets.extend(list(picked))
    return adj

def eigenvector_centrality(A, iters=300):
    N = A.shape[0]; x = np.ones(N); x /= np.linalg.norm(x)
    for _ in range(iters):
        x = A @ x; n = np.linalg.norm(x)
        if n == 0: break
        x /= n
    return x

def pagerank(A, d=0.85, iters=100):
    N = A.shape[0]; outdeg = A.sum(axis=1)
    M = A.T / np.where(outdeg==0, 1, outdeg)
    x = np.ones(N)/N
    for _ in range(iters):
        x = (1-d)/N + d*(M @ x)
    return x

rng = np.random.default_rng(7)
ba7 = barabasi_albert(300, 3, rng)
A7 = build_A(ba7)
ec = eigenvector_centrality(A7, 300)
pr = pagerank(A7, 0.85, 100)
print("top EC nodes:", np.argsort(ec)[::-1][:5].tolist(), "score range %.3f..%.3f" % (ec.min(), ec.max()))
print("top PR nodes:", np.argsort(pr)[::-1][:5].tolist())

实测：在 BA 网络上，度最高的 hub 同时拿到最高的特征向量中心性与 PageRank——两者高度一致。中心性直接支撑「影响力最大化」（先免疫/先投放给高中心性节点收益最大）与「关键基础设施保护」（高中心性节点是单点故障）。

九、与 AI / 工程映射

网络概念	在 AI / 系统中的对应
社团发现	图聚类、社区检测预训练、用户分群
Kuramoto 同步	联邦平均（FedAvg）的收敛、多智能体共识、模型一致性
SIR 传播	信息扩散、谣言抑制、病毒式营销、缓存预热
中心性	注意力路由、重要性采样、关键节点保护
小世界	小世界图神经网络（GraphSAGE / SGC 的邻居扩张）、高效消息传递

工程边界也要记牢：尺度——聚类系数与平均路径是 O(N²)（需采样估计）；稀疏性——幂律让度差异巨大，均值统计量失效；动态网络——上述都假设静态快照，真实网络随时变；分辨率极限——模块度本身会漏掉小于某尺度的社团（resolution limit），Louvain 不一定找回最自然的划分。

十、总结

本文用纯 NumPy 从零打通了复杂网络的四条工程主线，全部数值可复现：

1. 结构生成：WS 高聚类（0.393，约 ER 的 23 倍）+ 短路径 = 小世界；BA 幂律度分布 γ=2.94 = 无标度。
2. 社团发现：真实模块度 Q=0.491；谱聚类 / Louvain / 标签传播三种方法都把合成双社群恢复得 agreement=1.000，Louvain 以 O(N log N) 胜出。
3. 同步：Kuramoto 在 K_c≈1.6 处发生二阶相变，r 从 0.09 跃升到 0.97。
4. 传播：无标度网络上 β_c=0.0081 远低于均匀混合阈值 0.0167，hub 让系统更易被点燃（β=0.012 即爆发）。
5. 中心性：hub 同时主导特征向量中心性与 PageRank，是影响力与脆弱性的双重焦点。

复杂网络把「结构—动态—功能」连成一条链：社群决定信息怎么分层，同步决定系统怎么协同，传播决定风险怎么扩散，中心性决定该保护谁。它既是复杂系统的第二视角，也是图机器学习、群体智能与系统安全的方法支柱。