在知识图谱、分子检索和程序分析这些场景里“子图同构”是个绕不开的问题给定一个小模式图和一个大目标图找出大图中所有与小图同构的子结构。1976年Ullmann提出的Ullmann Algorithm用一张候选矩阵加剪枝回溯在今天依然是理解这类问题的绝佳起点也很适合用Python快速实现。我第一次真正吃透这个算法是在做分子结构检索的时候。当时数据集里只有几百个分子用排列组合的暴力搜索一秒钟都等不出来。问题不在代码跑得慢而是指数爆炸本身就无解。后来换用Ullmann的候选矩阵加refine剪枝同样规模的问题毫秒级出结果。这篇就把完整的Python实现、推导过程以及我踩过的几个坑一次说清楚。1. 子图同构到底在求解什么定义、难度与应用1.1 从“图同构”到“子图同构”的一步之遥先看最基础的问题。图同构是问两个图是不是“结构完全相同”存在一个顶点之间的一一映射让A图中任意两点有边当且仅当B图中对应两点也有边。这个问题有意思的地方在于它至今没有找到多项式时间算法也没被证明是NP完全问题属于那类“卡在中间”的怪问题。子图同构比图同构更实际也更难。它的定义可以写成这样给定模式图 P(Vp,Ep) 和目标图 G(Vg,Eg)寻找一个单射 f:Vp→Vg使得对P中的每条边(u,v)都有边(f(u),f(v))存在于G中。也就是说我们要在大图中找到一个小图的“完整副本”允许大图里还有其他顶点和边。这里还有个关键分支容易搞混诱导子图同构induced subgraph isomorphism和非诱导子图同构。上面的定义只要求模式图的边在目标图里都存在这叫非诱导匹配相当于“我只关心这组关系在不在不关心映射后的顶点之间还有没有额外联系”。诱导子图同构则额外要求模式图中不存在的边映射后的顶点之间也不能有边。用一个例子说明。给定模式图是一条三节点路径 [A-B-C]在目标图 [a-b-c-d] 这个四节点链上[A→b, B→c, C→d] 是合法匹配因为它们之间恰好有两条边。再看 [A→a, B→b, C→c]也合法。但如果目标图在b和d之间加一条边那么非诱导匹配的 [A→b, B→c, C→d] 仍然合法因为A-B对应b-c有边B-C对应c-d有边额外多出的b-d边不影响。诱导匹配就会拒绝它因为模式图里A和C没有边映射后的b和d却连在了一起。两者的选择取决于实际需求。化学分子检索通常用非诱导因为官能团匹配时片段上的两个原子映射到分子里两个原子之后它们之间还可能通过其他原子间接相连这是允许的。而社区发现、模式识别中的很多场景则要求映射后的子图“形状严丝合缝”必须用诱导匹配。1.2 NP完全的背面是大量的可剪枝空间子图同构是经典的NP完全问题。这意味着不存在多项式时间的通用算法最坏情况下模式图有n个顶点、目标图有m个顶点时需要考虑的组合数是排列数 P(m,n)m!/(m-n)!。当 n10、m50 时这个数字已经超过 3.7×10^16暴力枚举就是灾难。但NP完全说的是“最坏情况”现实场景里的图普遍稀疏还有节点标签、度数分布这些额外信息可以利用。好的算法可以通过剪枝把大量不可能的分支在很早的层级切断。Ullmann算法做的就是这件事先用一个候选矩阵粗略框定可能的映射再在回溯搜索过程中反复收紧这个矩阵。这也是我觉得Ullmann算法值得认真实现一遍的原因。它把“子图同构”这个抽象问题拆成了三个非常具体的工程问题候选矩阵怎么建、剪枝规则怎么写、回溯怎么搜。这三个问题搞清楚你再看VF2、TurboIso这些更现代的算法会轻松很多因为核心思路是一脉相承的。1.3 Ullmann算法的历史定位1976年J.R. Ullmann在论文《An Algorithm for Subgraph Isomorphism》里提出了这个算法。那还是计算机科学的早期论文里的符号和证明在今天读起来相当费劲但算法思想却异常朴素且有效。它的奠基性体现在“候选矩阵约束传播”这个框架上。后来的VF2算法用状态空间表示替代矩阵本质上还是“先缩小搜索空间再回溯”TurboIso利用候选区域和邻域关系也是在变着法子剪枝。可以说Ullmann算法是理解整个子图同构算法家族的钥匙。2. 手工推演一遍Ullmann的候选矩阵与剪枝逻辑2.1 候选矩阵M是怎么来的Ullmann算法的核心数据结构是一张布尔矩阵 M形状是 (n,m)。其中 n 是模式图顶点数m 是目标图顶点数。M[i][j]True 表示“模式顶点 i 有可能映射到目标顶点 j”。整个算法就是不断猜测、剪枝、验证这张矩阵里有哪些 True 能保留到最后。初始矩阵的构建并不复杂两条过滤规则标签过滤如果模式顶点 i 和目标顶点 j 有节点标签比如分子的原子类型那么标签不同则直接置 False。度数过滤如果模式顶点 i 的度数大于目标顶点 j 的度数那么 M[i][j]False。度数过滤的逻辑值得细说。这里的度数是邻接矩阵对应行中 True 的个数。一个顶点要映射到另一个顶点前提是目标顶点有足够多的邻居来容纳模式顶点的邻居。因为子图同构允许目标顶点有多余的邻居所以条件是“小于等于”而不是“等于”。以一条3节点路径作为模式图为例模式顶点0的度数是1那么目标图中任何度数大于等于1的顶点都可能是它的映射目标。如果有人把条件写成了“度数相等”就会漏掉大量解。在实际代码里我会用A[i].sum() B[j].sum()而不是。2.2 refine的核心邻居存在性检查初始矩阵只是粗筛真正让Ullmann算法区别于暴力搜索的是refine这个剪枝过程。refine的规则可以概括为一句话如果 M[i][j]True那么对模式图中 i 的每一个邻居 k目标图里 j 的邻居集合中必须至少有一个顶点 l 使得 M[k][l]True。为什么因为如果 i 映射到了 j那么 i 的邻居们必须映射到 j 的邻居们身上。假如某个邻居 k 的所有候选顶点都不在 j 的邻居集合里说明 k 无论怎么映射都无法和 j 形成边那 M[i][j]True 就是一个不可能的条件必须剪掉。这个检查是一个迭代过程。删掉一个候选之后可能引发连锁反应比如某个行只剩下一个候选删除后又导致其他行的候选不满足邻居条件。所以refine要用一个changed标志循环进行直到矩阵稳定。不过要特别提醒一点这个检查是必要条件不是充分条件。它只保证了“每个邻居都有退路”但没有保证这些退路互不相同也没有保证邻居之间的边最终真的存在。所以refine只能剪掉一部分不可能剩下的要靠回溯搜索和最终验证。这也是Ullmann算法的精髓——剪枝是安全的永远不会误删正确解。2.3 手工推演三角形模式在5节点图中的匹配过程理论讲多了容易虚直接手推一个例子。模式图 P 是三角形三个顶点两两相连。目标图 G 有5个顶点边是0-1、1-2、2-0、2-3、3-4。也就是说G中有一个三角形0-1-2三角形的一个顶点2上又连了一条尾巴3-4。P的三个顶点度数都是2。G中个顶点度数是0度2、1度2、2度3、3度2、4度1。按度数过滤P中每个顶点只能映射到0、1、2、3这四个顶点。初始候选矩阵如下True记作1False记作0M01234p011110p111110p211110第一次全局refine几乎不会剪掉任何候选因为每个候选顶点的邻居集合里都能找到下一行的候选顶点。比如 M[0][3]True表示p0可能映射到G的顶点3。顶点3的邻居是2和4而不论p1还是p2候选集合里都包含2所以这一步暂时保留。真正的剪枝发生在回溯赋值之后。假如先让p0映射到G的顶点0把M第0行改成只有M[0][0]True重新进入refine先看p1行。M[1][0]表示p1映射到顶点0但是p1的邻居p0已经固定在顶点0了而p0不在顶点0自己的邻居集合里——顶点0的邻居是1和2不包含0。所以M[1][0]被剪掉。同理M[1][3]也被剪掉因为顶点3的邻居是2和4不包含顶点0。p1这一行只剩下1和2两个候选。再看p2行。M[2][0]和M[2][3]也由于同样的原因被删掉p2的候选只剩1和2。接着DFS进入p1行尝试p1→1。再进入refinep2行的候选进一步收缩M[2][1]被剪掉因为顶点1的邻居是0和2而p1已经映射到1顶点1的邻居里找不到p1对应的候选。最后p2只能映射到顶点2。此时完整的映射是p0→0、p1→1、p2→2验证三条边都存在记录一个解。然后回退尝试p1→2经过类似推理得到p2→1映射p0→0、p1→2、p2→1也通过验证这是三角形的另一个方向。继续回退到第一层换p0的其他候选顶点最终能找出全部6个映射对应G中三角形0-1-2的全部对称映射。从这个推演里能直观看到Ullmann算法的节奏矩阵粗筛-回溯-局部refine-验证每一层递归都在用邻居关系快速收紧候选集。3. Python实现数据结构与完整代码3.1 统一输入接口邻接矩阵与networkx图为了让代码足够通用我会让主函数同时支持两种输入二维list或numpy数组形式的邻接矩阵以及networkx的Graph对象。这样一个工具函数make_adjacency_matrix就能把各种输入统一成numpy的bool矩阵。import numpy as np import networkx as nx def make_adjacency_matrix(graph, nodesNone, directedFalse): 将networkx图或二维list转换为邻接矩阵(bool)。 if isinstance(graph, (list, tuple, np.ndarray)): return np.asarray(graph, dtypebool) # networkx.Graph n graph.number_of_nodes() if nodes is None: nodes list(graph.nodes()) idx {nd: i for i, nd in enumerate(nodes)} A np.zeros((n, n), dtypebool) for u, v in graph.edges(): if u in idx and v in idx: A[idx[u], idx[v]] True if not directed: A[idx[v], idx[u]] True return A默认是无向图所以每条边会在邻接矩阵里对称地置两次。如果处理有向图把directedTrue传入即可。注意这里用bool矩阵而不是0/1整数矩阵后面做np.any、逻辑与运算时语义更干净。3.2 initial_matrix标签和度数双过滤初始候选矩阵的代码写起来非常直白。标签过滤是可选的传入labels_p和labels_g两个列表时才会生效。def initial_matrix(A, B, labels_pNone, labels_gNone): n, m A.shape[0], B.shape[0] M np.zeros((n, m), dtypebool) degree_p A.sum(axis1) degree_g B.sum(axis1) for i in range(n): for j in range(m): if labels_p is not None and labels_p[i] ! labels_g[j]: continue if degree_p[i] degree_g[j]: M[i, j] True return M这一步写完后矩阵里的True数量通常已经比全排列少很多。比如在上面的三角形例子里12个候选只剩12个因为4度顶点被排除了在带标签的大图里标签过滤往往能砍掉90%以上。3.3 refine剪枝函数的实现refine的核心逻辑就是把2.2节的规则翻译成代码。这里我用最直观的三层循环写便于理解后面第5章再给出优化思路。def refine(M, A, B): n, m A.shape[0], B.shape[0] changed True while changed: changed False for i in range(n): for j in range(m): if not M[i, j]: continue neighbor_i np.nonzero(A[i])[0] for k in neighbor_i: # 目标图j的邻居里是否存在模式顶点k的一个候选 if not np.any(B[j] M[k]): M[i, j] False changed True break return Mnp.nonzero(A[i])[0]拿出模式顶点i的所有邻居下标。B[j] M[k]是逻辑与取目标图j的邻居集合和模式顶点k的候选集合的交集。只要交集为空就说明这个候选不满足邻居存在性直接删掉。注意外层while changed。候选删除后其他候选的邻居集合可能因此失去交集形成连锁反应。这个循环保证了refine达到不动点才返回。3.4 回溯枚举与最终验证完整代码候选矩阵建好、refine写好后剩下的就是DFS回溯。这里我固定在每一层处理一行顺序从第0行到第n-1行。每行选择一个尚未被占用的列复制矩阵并局部refine如果后续某行候选为空就提前终止。def validate_mapping(mapping, A, B, inducedTrue): n A.shape[0] for i in range(n): for k in range(n): if i k: continue if A[i, k]: if not B[mapping[i], mapping[k]]: return False elif induced: if B[mapping[i], mapping[k]]: return False return True def ullmann(pattern_graph, target_graph, labels_pNone, labels_gNone, inducedTrue, node_list_pNone, node_list_gNone): A make_adjacency_matrix(pattern_graph, nodesnode_list_p) B make_adjacency_matrix(target_graph, nodesnode_list_g) n, m A.shape[0], B.shape[0] if n m: return [] M initial_matrix(A, B, labels_p, labels_g) M refine(M, A, B) if any(not M[i].any() for i in range(n)): return [] solutions [] assignment [] # assignment[i] 表示模式顶点 i 映射到了哪个目标顶点 def dfs(M_local, row): if any(not M_local[i].any() for i in range(n)): return if row n: if validate_mapping(assignment, A, B, inducedinduced): solutions.append(assignment.copy()) return used set(assignment) for j in range(m): if not M_local[row, j]: continue if j in used: continue M_next M_local.copy() M_next[row, :] False M_next[row, j] True M_next refine(M_next, A, B) if not M_next[row, j]: continue assignment.append(j) dfs(M_next, row 1) assignment.pop() dfs(M, 0) return solutions几个值得说明的代码细节M_next M_local.copy()是numpy的拷贝浅拷贝和深拷贝在这里效果一致因为元素都是布尔值。used set(assignment)保证了单射性避免两个模式顶点映射到同一个目标顶点。dfs开头对每一行检查not M_local[i].any()这一步不能省。因为在递归过程中refine可能把前面已经赋值的行重新剪成空尤其当选择了一个不可能的分支时提前发现并返回能节省大量无用递归。3.5 先用小例子把代码跑起来用一个简单例子验证正确性模式图是一条3节点路径目标图是一条5节点路径0-1-2-3-4。pattern np.array([ [0, 1, 0], [1, 0, 1], [0, 1, 0] ], dtypebool) target np.array([ [0, 1, 0, 0, 0], [1, 0, 1, 0, 0], [0, 1, 0, 1, 0], [0, 0, 1, 0, 1], [0, 0, 0, 1, 0] ], dtypebool) sols ullmann(pattern, target, inducedFalse) print(sols)理论上5节点链里能找到6个长度为2的路径子图。输出大致是这样[[0, 1, 2], [2, 1, 0], [1, 2, 3], [3, 2, 1], [2, 3, 4], [4, 3, 2]]每个列表表示模式顶点0、1、2分别映射到目标图的哪个顶点。这个结果能帮你快速确认实现没有方向性问题。4. 几个容易踩的坑映射方向、诱导要求与矩阵状态管理4.1 行列对应模式顶点永远是行M[i][j] 里i永远来自模式图j永远来自目标图。这个约定一旦搞反refine的邻居检查就会语义错乱而且很难通过单元测试发现因为小例子上可能碰巧也能跑出结果。我的建议是在代码里显式写上n A.shape[0]作为模式图顶点数m B.shape[0]作为目标图顶点数并写上注释。一旦出现维度不一致第一时间检查是不是赋值的时候把行列写反了。4.2 诱导子图同构 vs 非诱导需求要先想清楚这是我在实际项目里踩过的最深的一个坑。接口里默认inducedTrue也就是默认要求诱导子图同构。原因是我在多数知识图谱和程序分析场景里需要“严丝合缝”的匹配额外边往往是错误匹配。但很多人手动调用时并不清楚这个语义。比如在化学里查找官能团通常应该传inducedFalse。如果误把诱导当成非诱导会漏掉大量真实存在的匹配反过来如果把非诱导当成诱导会出现大量假阳性。一个稳妥的办法拿到需求后先问一句“映射之后的顶点之间如果存在模式里没有的边算不算匹配成功”。这比在代码里纠结语义重要得多。4.3 每次递归都要复制矩阵但不能无脑复制refine是原地修改矩阵的。DFS回溯时需要回到修改前的状态处理方法有两类拷贝法每次递归前M_local.copy()修改后传入下一层。优点是不会污染上一层状态缺点是每层多一次O(n×m)的拷贝开销。备份回滚法refine过程中记录所有被删除的候选递归结束后恢复。优点是省拷贝缺点是代码复杂度明显上升容易出错。我的实现选择拷贝法因为它和算法的递归结构完美对应调试时也能直观看到每一层的矩阵快照。当 n×m 在几千到几万这个量级时拷贝开销完全可接受。只有当矩阵规模更大、且性能敏感时才值得去做in-place回滚。5. 让算法跑得更快MRV启发式与numpy向量化5.1 MRV从约束最强的行开始分配固定顺序搜索的一个问题是如果特别“松”的行排在最前面会产生大量无效分支。经典的做法是用MRV启发式也就是每次选择尚未分配、且候选数最少的行进行扩展。这就好比填数独时先填候选数最少的格子。Ullmann的MRV版本只需要改动dfs里的选行逻辑def dfs_mrv(M_local, assigned_rows, mapping): if any(not M_local[i].any() for i in range(n)): return if len(assigned_rows) n: if validate_mapping(mapping, A, B, inducedinduced): solutions.append(mapping.copy()) return remaining [i for i in range(n) if i not in assigned_rows] row min(remaining, keylambda i: int(M_local[i].sum())) used set(mapping.values()) for j in range(m): if not M_local[row, j] or j in used: continue M_next M_local.copy() M_next[row, :] False M_next[row, j] True M_next refine(M_next, A, B) mapping[row] j assigned_rows.add(row) dfs_mrv(M_next, assigned_rows, mapping) assigned_rows.remove(row) del mapping[row]注意这里mapping改成了字典因为行的处理顺序不再固定无法用列表的assignment[i]来对应行号。输出时需要用[mapping[i] for i in range(n)]恢复顺序。MRV的收益在模式图顶点度数很不均匀时特别明显。比如模式图是一个“星形”中心顶点度数很高、叶子顶点度数都是1直接从中心开始匹配几分钟能解决的问题可能几秒就完成。5.2 向量化refine的提速思路基础版refine里三层循环的写法适合理解但在 m 比较大的图上会成为性能瓶颈。一个直接的向量化思路是用矩阵乘法替代内层的“交集判断”。核心公式是对于模式顶点 k 和第 j 个目标顶点(B M[k].astype(int))[j] 的值就是目标图顶点 j 的邻居集合与 M[k] 候选集合的交集大小。只要这个值大于0就表示存在至少一个共同的候选。def refine_vec(M, A, B): n, m A.shape[0], B.shape[0] while True: neighbor_ok np.zeros((n, m), dtypebool) for k in range(n): neighbor_ok[k] (B M[k].astype(int)) 0 # ok[i, j] True 当且仅当 i 的所有邻居 k 都有 neighbor_ok[k, j] True ai_sum A.sum(axis1, keepdimsTrue) ok np.zeros((n, m), dtypebool) for i in range(n): ok[i] (A[i] neighbor_ok) ai_sum[i] if np.array_equal(ok, M): return M M ok这个版本把原来最内层的大循环换成了矩阵乘法在numpy底层用BLAS加速。对于稠密矩阵和大图性能提升非常显著。需要注意的是这个版本是“同步更新”的每次迭代用上一轮完整的M计算新的ok和基础版的异步逐元素更新略有不同但剪枝的逻辑语义是一致的。5.3 静态重排让度数大的顶点先出场MRV是动态选择行而静态重排是在运行前重新排列模式图顶点顺序。思路是把模式图的顶点按度数从大到小排序让约束最强的顶点先被分配。实现起来更简单而且效果通常和MRV不相上下degree_p A.sum(axis1) order np.argsort(-degree_p) A_perm A[order][:, order] sols_perm ullmann(A_perm, B, ...) solutions [] for sol in sols_perm: mapping [0] * n for new_i, old_i in enumerate(order): mapping[old_i] sol[new_i] solutions.append(mapping)静态重排的另一个好处是它和任何版本的refine都兼容是一个完全独立的预处理步骤。如果你的项目里模式图相对固定、目标图经常变化甚至可以把排序后的模式图缓存起来。6. 实际场景分子片段匹配上的快速验证6.1 模拟分子片段查询用标签数组过滤Ullmann算法最有名的应用之一就是化学分子中的子结构检索。我这里用一个简化的例子来演示不涉及化合价和立体化学只看原子类型和连接关系。假设模式图是一个“C-O-H”片段表示一个羟基连着一个碳。目标图是一个更大的假想分子。用networkx构造pattern nx.Graph() pattern.add_node(c, atomC) pattern.add_node(o, atomO) pattern.add_node(h, atomH) pattern.add_edge(c, o) pattern.add_edge(o, h) target nx.Graph() target.add_node(a, atomC) target.add_node(b, atomC) target.add_node(c, atomO) target.add_node(d, atomH) target.add_node(e, atomH) target.add_edge(a, b) target.add_edge(b, c) target.add_edge(c, d) target.add_edge(c, e)这里目标图相当于一个乙二醇片段C-C-O连着两个H原子。传入显式节点列表保证内部编号和输出命名一致p_nodes list(pattern.nodes()) g_nodes list(target.nodes()) labels_p [pattern.nodes[nd][atom] for nd in p_nodes] labels_g [target.nodes[nd][atom] for nd in g_nodes] sols ullmann(pattern, target, labels_plabels_p, labels_glabels_g, inducedFalse, node_list_pp_nodes, node_list_gg_nodes) for sol in sols: print([g_nodes[i] for i in sol])标签过滤会在initial_matrix阶段就把大量标签不匹配的候选挑出去例如模式里的O不会去匹配目标里的C。度数过滤再砍掉一批。最终输出的两个解应该分别是[b, c, d]和[b, c, e]对应两个氢原子位置。我在做这类查询时经常配合后处理检查化合价和环信息。图同构只能保证拓扑关系一致真实世界的约束还要靠业务规则叠加。6.2 Ullmann算法的适用边界与替代方案纯Python的Ullmann实现配上numpy和MRV优化处理模式图顶点数不超过10、目标图顶点数不超过几百的规模通常毫无压力。我自己在实验里匹配一个6顶点的模式图、300顶点的目标图经常是几十毫秒到一两秒就结束。模式图顶点数超过15、目标图超过几千时就要认真考虑换工具了。常见的替代方案networkx自带的GraphMatcher基于VF2算法接口成熟处理大规模稀疏图表现更好igraph的subisomorphic_lad基于LAD算法对带标签的图非常快如果是数据库里的子图查询Neo4j等图数据库内置了索引和匹配优化比在应用层自己写要可靠得多。这并不意味着Ullmann过时了。恰恰相反当你需要对匹配过程做深度定制——比如自定义剪枝规则、限定时长、返回所有匹配而不只是判断是否存在——自己实现的Ullmann框架远比黑盒库灵活。最后再分享两个小技巧。第一调试Ullmann时在refine之后打印一下候选矩阵的非零个数如果某次递归后数量没有变化说明剪枝没有生效优先检查refine里的邻居集合逻辑第二如果目标图非常大先用节点标签做一次倒排索引快速排除不可能映射到某个模式顶点的目标顶点集合再调用initial_matrix能省下不少无效的np.any计算。根据我的实际使用经验这两招对性能的影响往往比后面换更快的算法还明显。