控制中心性——单个节点能控制的子空间维数

# Liu、Slotine 与 Barabási（2012）把“只驱动节点 i 时可控子空间的通用维数”定义为节点 i 的控制中心性。
# 可控子空间就是 Krylov 子空间 span{b, Ab, A^2 b, ...}；这里给非零边赋随机权重（对应“几乎所有权重”），
# 用逐步正交化求其维数，避免直接计算 A^{N-1} 带来的数值病态。
# 对无环有向图（DAG），控制中心性等于节点在层次结构中的层号（Hosoe 定理的特例），本例同时打印两者对照。
# 解读：控制中心性高的节点未必是度最高的节点；它反映的是从该节点出发能“展开”多长的不相交路径与环。
import numpy as np
import networkx as nx

def krylov_dim(A, b, tol=1e-9):
    Q = []
    v = b / np.linalg.norm(b)
    for _ in range(A.shape[0]):
        for q in Q:                              # 两次 Gram–Schmidt，保证正交性
            v = v - (q @ v) * q
        for q in Q:
            v = v - (q @ v) * q
        nv = np.linalg.norm(v)
        if nv < tol:
            break
        v = v / nv
        Q.append(v)
        v = A @ v
    return len(Q)

def control_centrality(G, seed=0):
    nodes = list(G)
    idx = {u: k for k, u in enumerate(nodes)}
    rng = np.random.default_rng(seed)
    A = np.zeros((len(nodes), len(nodes)))
    for u, v in G.edges():
        A[idx[v], idx[u]] = rng.uniform(0.5, 1.5)
    out = {}
    for u in nodes:
        b = np.zeros(len(nodes)); b[idx[u]] = 1.0
        out[u] = krylov_dim(A, b)
    return out

# 一个小的有向无环网络（层次结构）
G = nx.DiGraph([(0, 1), (0, 2), (1, 3), (2, 3), (3, 4), (4, 5), (1, 6), (6, 5), (7, 6)])
cc = control_centrality(G)
# 层号：从节点出发的最长有向路径所含的节点数（DAG 情形，对应 Liu 2012 的层次标记）
longest = {u: 1 + max([0] + [len(p) - 1 for w in nx.descendants(G, u) for p in nx.all_simple_paths(G, u, w)]) for u in G}
for u in sorted(G):
    print("节点 %d：出度 %d，控制中心性 %d，最长下游路径节点数 %d" % (u, G.out_degree(u), cc[u], longest[u]))

# 含环网络：随机有向网络中控制中心性的分布
H = nx.gnm_random_graph(60, 90, directed=True, seed=2)
ccH = control_centrality(H)
vals = np.array(list(ccH.values()))
outdeg = np.array([H.out_degree(u) for u in H])
print("随机有向网络 N=60：控制中心性 最大 %d，中位数 %.0f；与出度的 Spearman 相关 = %.2f"
      % (vals.max(), np.median(vals),
         np.corrcoef(np.argsort(np.argsort(vals)), np.argsort(np.argsort(outdeg)))[0, 1]))

import matplotlib.pyplot as plt
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
layers = {}
for u in sorted(G):
    layers.setdefault(longest[u], []).append(u)
pos = {u: (-L, j - (len(us) - 1) / 2) for L, us in layers.items() for j, u in enumerate(us)}
nx.draw_networkx_edges(G, pos, ax=ax1, edge_color="0.5", arrows=True, arrowsize=12, node_size=350)
nodes_ = nx.draw_networkx_nodes(G, pos, ax=ax1, node_size=350, node_color=[cc[u] for u in G], cmap="YlGnBu", vmin=0, vmax=7, edgecolors="0.3")
nx.draw_networkx_labels(G, pos, ax=ax1, font_size=8)
plt.colorbar(nodes_, ax=ax1, label="control centrality"); ax1.axis("off")
ax1.set_title("DAG: control centrality = layer index")
jit = np.random.default_rng(0).uniform(-0.15, 0.15, (2, len(vals)))
ax2.scatter(outdeg + jit[0], vals + jit[1], s=30, alpha=0.6)
ax2.set_xlabel("out-degree"); ax2.set_ylabel("control centrality (dim of controllable subspace)")
ax2.set_title("Random digraph N = 60: not determined by degree")
plt.tight_layout(); plt.show()