特征向量、Katz、PageRank 与 HITS 的幂迭代实现

# 用同一个 6 节点有向小图说明谱类中心性的差别（约定 A[i, j] = 1 表示边 i -> j）：
#   特征向量中心性在有向无环部分会退化为 0；Katz 加常数项 β 解决这一问题；
#   PageRank 再按出度分摊并加随机跳转；HITS 同时给出权威值（authority）与枢纽值（hub）。
# 每种方法都用幂迭代手工实现，并与 NetworkX 核对。
# 输出：λ1 = 1（来自环 3→4→5→3）；不在环上的节点 0、1、2 特征向量中心性都为 0，
#       而在 Katz 与 PageRank 中为正，节点 0（无入边）最小；
#       HITS 中节点 0 是最强的“枢纽”，被多方指向的节点 3 是最强的“权威”。
import networkx as nx
import numpy as np

edges = [(0, 1), (0, 2), (0, 3), (1, 3), (2, 3), (3, 4), (4, 5), (5, 3), (1, 2)]
D = nx.DiGraph(edges)
nodes = sorted(D)
A = nx.to_numpy_array(D, nodelist=nodes)       # A[i, j] = 1: i -> j
N = len(nodes)


def power_iter(M, x0, n_iter=1000):
    x = x0.copy()
    for _ in range(n_iter):
        y = M @ x
        x = y / np.linalg.norm(y)
    return x


# 特征向量中心性（按入边累积）：λ x = A^T x
lam = max(abs(np.linalg.eigvals(A)))
x_eig = power_iter(A.T, np.ones(N))
print("λ1 = %.3f；特征向量中心性:" % lam, np.round(x_eig / x_eig.sum(), 3))

# Katz：x = α A^T x + β 1  =>  x = β (I - α A^T)^{-1} 1，要求 α < 1/λ1
alpha = 0.5 / lam
x_katz = np.linalg.solve(np.eye(N) - alpha * A.T, np.ones(N))
ref = nx.katz_centrality_numpy(D, alpha=alpha, beta=1.0)
print("Katz（α=%.3f）:" % alpha, np.round(x_katz / np.linalg.norm(x_katz), 3),
      " NetworkX 一致:", np.allclose(x_katz / np.linalg.norm(x_katz), [ref[v] for v in nodes]))

# PageRank：x = d A^T D_out^{-1} x + (1-d)/N；本图没有出度为 0 的节点
d = 0.85
kout = A.sum(axis=1)
P = (A / kout[:, None]).T                      # 列随机矩阵
x_pr = np.ones(N) / N
for _ in range(200):
    x_pr = d * P @ x_pr + (1 - d) / N
ref = nx.pagerank(D, alpha=d)
print("PageRank:", np.round(x_pr, 3), " NetworkX 一致:", np.allclose(x_pr, [ref[v] for v in nodes], atol=1e-6))

# HITS：a = A^T h, h = A a  =>  a ∝ A^T A a, h ∝ A A^T h
a = power_iter(A.T @ A, np.ones(N)); h = power_iter(A @ A.T, np.ones(N))
h_ref, a_ref = nx.hits(D, max_iter=1000)
print("authority:", np.round(a / a.sum(), 3), " NetworkX 一致:", np.allclose(a / a.sum(), [a_ref[v] for v in nodes], atol=1e-6))
print("hub      :", np.round(h / h.sum(), 3), " NetworkX 一致:", np.allclose(h / h.sum(), [h_ref[v] for v in nodes], atol=1e-6))

# 作图：5 种谱类中心性（各自归一化到和为 1）的分组柱状图
import matplotlib.pyplot as plt

scores = {"eigenvector": x_eig / x_eig.sum(), "Katz": x_katz / x_katz.sum(), "PageRank": x_pr,
          "authority": a / a.sum(), "hub": h / h.sum()}
fig, ax = plt.subplots(figsize=(6.4, 4.2))
w = 0.16
for j, (name, v) in enumerate(scores.items()):
    ax.bar(np.arange(N) + (j - 2) * w, v, width=w, label=name, hatch=["", "//", "", "..", "xx"][j],
           edgecolor="white" if j % 2 == 0 else "0.3")
ax.set_xticks(range(N)); ax.set_xticklabels(["node %d" % v for v in nodes])
ax.set_ylabel("score (sums to 1)"); ax.legend(fontsize=8, ncol=2)
ax.set_title("Nodes 0-2 vanish under eigenvector centrality only", fontsize=10)
plt.tight_layout(); plt.show()