划分空间有多大——Bell 数与小网络上的穷举最优模块度

# 节点数为 N 的网络共有 Bell 数 B_N 种划分；只考虑等分二分也有 C(N, N/2)/2 种。
# 先打印这些数随 N 的增长，再在一个 10 节点网络上穷举全部 B_10 = 115975 种划分，
# 找到精确的最大模块度划分，并与 Louvain（多个随机种子）的结果比较。
# 输出：N=50 时划分数超过 10^47，穷举不可行；在小网络上 Louvain 达到或接近精确最优。
from math import comb
import networkx as nx
import numpy as np


def bell_numbers(n_max):
    """用 Bell 三角形计算 B_0..B_{n_max}。"""
    row, bells = [1], [1]
    for _ in range(n_max):
        new = [row[-1]]
        for x in row:
            new.append(new[-1] + x)
        row = new
        bells.append(row[0])
    return bells


B = bell_numbers(50)
for N in [6, 10, 20, 50]:
    print("N=%2d  Bell 数 B_N = %.3e   等分二分数 = %.3e" % (N, B[N], comb(N, N // 2) / 2))


def set_partitions(n):
    """以“受限增长串”枚举 {0..n-1} 的全部划分：label[0]=0，label[i] ≤ max(label[:i])+1。"""
    labels = [0] * n

    def rec(i, mx):
        if i == n:
            yield labels
            return
        for c in range(mx + 2):
            labels[i] = c
            yield from rec(i + 1, max(mx, c))

    yield from rec(1, 0)


G = nx.Graph([(0, 1), (0, 2), (1, 2), (2, 3), (3, 4), (3, 5), (4, 5), (5, 6),
              (6, 7), (6, 8), (7, 8), (8, 9), (7, 9), (1, 4), (0, 9)])
n = G.number_of_nodes()
A = nx.to_numpy_array(G, nodelist=range(n))
k = A.sum(1)
m = A.sum() / 2
Bmat = A - np.outer(k, k) / (2 * m)

best_q, best_lab, count = -1.0, None, 0
for lab in set_partitions(n):
    count += 1
    lab_arr = np.array(lab)
    q = Bmat[lab_arr[:, None] == lab_arr[None, :]].sum() / (2 * m)
    if q > best_q + 1e-12:
        best_q, best_lab = q, lab_arr.copy()
print("穷举了 %d 种划分（应等于 B_10 = %d）" % (count, B[10]))
print("精确最大模块度 Q* = %.4f，最优划分标签:" % best_q, best_lab.tolist())

q_louvain = [nx.community.modularity(G, nx.community.louvain_communities(G, seed=s)) for s in range(20)]
print("Louvain 20 个种子：Q 最大 %.4f，最小 %.4f" % (max(q_louvain), min(q_louvain)))

# 作图：左为划分数随 N 的增长（对数纵轴）；右为 10 节点网络的精确最优划分
import matplotlib.pyplot as plt

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))
Ns = np.arange(2, 51)
ax1.semilogy(Ns, [float(B[N]) for N in Ns], "-", label="Bell number B_N (all partitions)")
ax1.semilogy(Ns, [comb(N, N // 2) / 2 for N in Ns], "--", label="balanced bisections")
ax1.set_xlabel("number of nodes N"); ax1.set_ylabel("count (log scale)"); ax1.legend(fontsize=8)
ax1.set_title("Partition space explodes: B_50 = %.1e" % B[50], fontsize=10)
pos = nx.spring_layout(G, seed=1)
nx.draw_networkx_edges(G, pos, ax=ax2, edge_color="0.6", alpha=0.5)
nx.draw_networkx_nodes(G, pos, ax=ax2, nodelist=range(n), node_color=best_lab, cmap="Set2", node_size=300, edgecolors="0.3")
nx.draw_networkx_labels(G, pos, ax=ax2, font_size=9)
ax2.set_title("Exact optimum Q* = %.4f; best Louvain Q = %.4f" % (best_q, max(q_louvain)), fontsize=10); ax2.axis("off")
plt.tight_layout(); plt.show()