子图同构的问题我最早是在做化学分子结构校验时碰到的给定一个小分子模板要去一个很大的分子库里找有没有包含这个片段的结构。后来在做知识图谱的规则校验、社交网络里的圈子发现时又反复遇到同一类问题。它本质上就是——给定一个模式图和一个目标图判断目标图中是否存在一个子图能和模式图“长得一样”。这个问题看起来简简单单但暴力枚举在节点一多就彻底跑不动所以才需要专门设计算法。Ullmann Algorithm乌尔曼算法就是这方面最经典的算法之一1976年提出至今仍然是理解子图同构和回溯剪枝思想的绝佳教材。本文就带你从零开始用纯Python实现一遍Ullmann算法不依赖networkx、igraph这类现成图数据库所有逻辑都建立在基础的邻接矩阵上。适合刚接触图算法、想搞懂子图同构本质的读者也适合想系统了解经典算法工程化实现的同学。1. 子图同构问题与Ullmann算法的核心思想1.1 先看问题模式图 vs 目标图先把术语说清楚。我们有两个图一个是小图叫模式图记作图A另一个是大图叫目标图记作图B。子图同构的任务就是在B里找到一组节点使得这组节点之间的连接关系完全覆盖A中所有节点之间的连接关系。注意是“覆盖”不是“完全相等”。也就是说A里的边在映射后的目标节点对上必须存在但目标图里多余的那些边不影响结果。举例来说A是一个三角形B是一个正方形含一条对角线那A就可以映射到B中任意一个三角形——B里的对角线是多余的不影响。这类映射关系在数学上叫“单射”每个模式节点映射到不同的目标节点但在边上只要求“保留”不要求“禁止额外”。这个定义很关键因为很多人把子图同构和“诱导子图同构”混在一起。诱导子图同构要求不仅A里的边要存在A里没有的边映射后也必须不存在。Ullmann算法解决的是前者也就是非诱导子图同构。搞清楚这一点后面代码里“最终验证”环节就很好理解了。1.2 为什么还在用Ullmann算法图同构和子图同构问题有一个很有意思的背景图同构问题至今没有多项式复杂度算法子图同构问题更是被证明是NP完全问题。也就是说在最坏情况下我们不可能指望有特别快的解法真正的战场是“剪枝”——在搜索过程中尽早排除不可能的分支。Ullmann算法能在众多算法里经久不衰主要因为三点思路直观它把“找同构”变成“找一个候选映射矩阵”然后不断在矩阵上做布尔运算。整个过程可以用纸笔模拟。剪枝效果好它的细化refinement过程能在搜索早期大幅降低候选空间尤其适合中等规模、结构相对稠密的图。实现门槛低用Python只需要几十行核心代码就能跑通不像VF2等算法需要复杂的可行性判定规则。对于学习者和很多实际场景来说这是很大的优势。后面你会看到这个算法天然分成“构建候选矩阵—细化剪枝—回溯搜索”三个模块模块之间耦合度很低改动起来非常方便。1.3 整体流程三阶段走在写代码之前先把Ullmann算法的整体流程画在脑子里整个过程其实分三步根据节点度数等约束构建一个 d×k 的候选矩阵M。d是模式图A的节点数k是目标图B的节点数。如果M[i][j]1表示模式图节点i可以映射到目标图节点j。对这个矩阵做细化。细化规则是如果M[i][j]1那么对于A中i的每一个邻居i2在B中j的邻居里必须至少存在一个节点j2使得M[i2][j2]1。如果不满足就把M[i][j]改成0。在细化后的矩阵上做回溯搜索。每次给当前模式节点选一个目标节点选完继续细化递归往下走直到所有模式节点都映射完成最后做一次边验证记录结果。这三步里第一步负责“初筛”第二步负责“大规模剪枝”第三步负责“精确枚举”。为什么Ullmann快不是因为最后那步枚举快而是因为第二步已经把绝大多数不可能的分支都砍掉了。2. 算法原理深入解析2.1 邻接矩阵图的“身份证”代码里所有图都用邻接矩阵表示。一个d个节点的图就是一个d×d的0/1矩阵A[i][j]1表示节点i到节点j有一条边。无向图的邻接矩阵是对称的A[i][j]A[j][i]有向图则不一定。做了十几年图算法我现在的习惯是只要不是特别复杂的图结构一律先用邻接矩阵打底因为矩阵运算配合numpy可以写得很简洁而且调试的时候打印出来一目了然。本文为了展示算法本质全部用Python原生list实现不引numpy这样逻辑更透明。要注意的是本文约定邻接矩阵对角线都是0也就是不考虑自环。如果你要处理自环细化和验证条件还得额外判断这个后面在常见问题里我会专门说。2.2 映射矩阵M候选表怎么建Ullmann算法的核心数据结构是映射矩阵M维度是 d×k。M[i][j]1表示“模式图节点i可能映射到目标图节点j”0表示“不可能。”初始矩阵怎么建最简单也最常用的初筛条件是看节点度数如果模式图里节点i的度数大于目标图里节点j的度数那i无论如何也不可能映射到j。原因很直接假如i要映射到j那么i的每条边都要对应到j的一条边如果j的边数比i还少必然有边找不到对应位置。这个条件还可以扩展比如节点如果有标签/颜色属性可以在初始化时加一个判断标签不同直接置0。后面测试部分你会看到这种“标签预剪枝”能带来非常明显的加速效果。代码如下def build_M(A, B): d, k len(A), len(B) deg_A [sum(row) for row in A] deg_B [sum(row) for row in B] M [[0] * k for _ in range(d)] for i in range(d): for j in range(k): if deg_A[i] deg_B[j]: M[i][j] 1 return M这里sum(row)就是对邻接矩阵某一行求和无向图不考虑自环时得到的就是节点度数。如果你要给节点加标签就在循环里多写一个判断条件。2.3 细化Refinement剪掉不可能的组合细化是整个Ullmann算法的灵魂。很多人第一次看论文会被那堆布尔矩阵公式绕晕我把这条规则翻译成人话如果模式图节点i真的能映射到目标图节点j那么i在模式图里的每个邻居节点都必须能在j的目标图邻居里找到对应的落脚点。这个“落脚点”不是随便一个节点是同时满足两个条件的节点第一它在目标图里和j相连第二它在当前候选矩阵M中是模式图那个邻居的候选节点之一。用代码来表达这段逻辑就是遍历M里所有等于1的位置对每个位置做一次邻居检查一旦发现某个邻居找不到合格目标就把这个位置的1改成0。改完后可能又会导致其他位置不满足条件所以要循环反复地检查直到矩阵不再变化。def refine(M, A, B): d, k len(A), len(B) changed True while changed: changed False for i in range(d): for j in range(k): if M[i][j] 1: keep True for i2 in range(d): if A[i][i2] 1: has_neighbor False for j2 in range(k): if B[j][j2] 1 and M[i2][j2] 1: has_neighbor True break if not has_neighbor: keep False break if not keep: M[i][j] 0 changed True return M请注意这里有三层循环最坏情况下是 O(d²k²) 的复杂度但实际操作中每一轮循环都会砍掉一批1循环次数通常很小。这段代码还有几个可以优化的细节如果某一行已经全变成0了那整个映射直接不存在搜索可以提前终止。如果某一列全为0意味着目标图节点j不会成为任何模式节点的映射目标这一列可以直接忽略不再参与后续枚举。细化过程之所以强是因为它反复利用“已经确定不可能的信息”去推导新的不可能。这个思想非常类似数独里从已知格推未知格的逻辑只不过这里的“已知”和“未知”都是在候选矩阵里动态变化的。3. Python实现完整代码逐步拆解3.1 环境与整体骨架实现这个算法只需要Python 3.8以上的标准库不需要pip安装任何第三方包。我平时写这类算法一般就是一个文件搞定方便快速测试和集成。如果你用的是vscode直接新建一个ullmann.py把代码粘进去终端里python ullmann.py就能运行。整个工程的骨架可以分成四个函数build_M(A, B)构建初始候选矩阵refine(M, A, B)细化剪枝dfs(step, cur_M, ...)回溯搜索find_subgraph_isomorphisms(A, B)对外入口整合以上步骤3.2 build_M构造初始候选矩阵这部分在2.2里已经给出了代码。我在实际项目中经常会在这一步根据业务场景加入各种“低成本约束”。比如如果模式图节点i和目标图节点j的颜色、类型、标签不同直接置0如果模式图是带权图甚至可以要求目标边的权重满足上下界如果已知某些节点必须映射到某些固定节点可以直接把M里对应的行和列锁死初始化矩阵做得越“狠”后面回溯搜索的压力就越小。这就像相亲平台先根据年龄、收入、城市筛掉一批后面见面才不会太累。3.3 refine剪枝核心函数这个函数在2.3里已经给出。我写代码时的建议是不要追求一步到位的“极致优化”先把标准版写对再逐层加优化。比如细化过程在“判断邻居存在”时其实可以先用any()生成器压缩代码但我故意写成显式循环目的是让你看清楚每一层到底在干什么。一个我曾经踩过的坑是在细化循环里一边修改M一边继续遍历同一层的M结果导致某轮循环里刚置0的格子又被当作候选来推导。解决方式其实很简单就是像上面的代码一样设置一个changed标志每轮扫描结束后再决定要不要从头再来。这个“扫描直到稳定”的模式在很多图算法里都通用。3.4 回溯搜索与验证候选矩阵准备好了接下来就要枚举实际的映射关系了。回溯搜索的思路是从模式图里挑一个“最不好挑”的节点也就是候选目标最少的节点优先给它分配。这个策略叫MRVMinimum Remaining Values是约束满足问题里的经典启发式。遍历该节点的候选目标列如果该列还没被别的模式节点占用就尝试把模式节点映射到该列。把当前行的候选压缩到一列后重新做一次细化。如果细化后没有任何未分配模式节点的行变成全零就进入下一层递归。所有模式节点都分配完毕时做一次最终的边验证确认模式图里每条边在映射后都存在。这里有个非常关键的细节细化会修改矩阵所以在进入下一层递归时必须复制一份当前矩阵否则这一层改完回溯到上一层再试其他候选列时矩阵已经被“污染”了。这个问题我在实际项目中栽过好几次跟头后面会专门展开讲。最终验证函数可以写得很简短。映射后我们手上有一个长度为d的列表mappingmapping[i]表示模式图节点i对应的目标图节点。只要遍历模式图所有边检查目标图中对应位置是否也是1def verify(mapping, A, B): d len(A) for u in range(d): for v in range(d): if A[u][v] 1 and B[mapping[u]][mapping[v]] ! 1: return False return True注意这里只检查A中为1的边。A中为0的边不检查因为子图同构允许目标图有多余边。3.5 完整代码汇总把上面的函数整合到一起就是一个完整可运行的Ullmann算法def build_M(A, B): d, k len(A), len(B) deg_A [sum(row) for row in A] deg_B [sum(row) for row in B] M [[0] * k for _ in range(d)] for i in range(d): for j in range(k): if deg_A[i] deg_B[j]: M[i][j] 1 return M def refine(M, A, B): d, k len(A), len(B) changed True while changed: changed False for i in range(d): for j in range(k): if M[i][j] 1: keep True for i2 in range(d): if A[i][i2] 1: has_neighbor False for j2 in range(k): if B[j][j2] 1 and M[i2][j2] 1: has_neighbor True break if not has_neighbor: keep False break if not keep: M[i][j] 0 changed True return M def find_subgraph_isomorphisms(A, B): d, k len(A), len(B) if d k: return [] M build_M(A, B) M refine(M, A, B) if any(sum(row) 0 for row in M): return [] order sorted(range(d), keylambda i: sum(M[i])) used [False] * k mapping [None] * d results [] def dfs(step, cur_M): if step d: for u in range(d): for v in range(d): if A[u][v] 1 and B[mapping[u]][mapping[v]] ! 1: return results.append(list(mapping)) return row order[step] for col in range(k): if used[col] or cur_M[row][col] 0: continue next_M [r[:] for r in cur_M] for c in range(k): next_M[row][c] 1 if c col else 0 next_M refine(next_M, A, B) mapping[row] col used[col] True if not any(sum(next_M[u]) 0 for u in range(d) if mapping[u] is None): dfs(step 1, next_M) used[col] False mapping[row] None dfs(0, M) return results整个代码不到70行能胜任任意两个图的子图同构枚举。理解这70行的关键在于搞明白cur_M和next_M的关系cur_M是当前层的候选矩阵next_M是从它派生出的下一层矩阵每一层都持有自己的“世界状态”。4. 测试案例与运行结果分析4.1 案例一三角形模式匹配先来一个最经典的例子。模式图A是三角形3个节点两两相连目标图B是4个节点的图其中前3个节点也构成三角形第4个节点是孤立点A [ [0, 1, 1], [1, 0, 1], [1, 1, 0] ] B [ [0, 1, 1, 0], [1, 0, 1, 0], [1, 1, 0, 0], [0, 0, 0, 0] ] results find_subgraph_isomorphisms(A, B) print(找到, len(results), 个映射) for m in results: print(m)运行结果是“找到6个映射”。为什么是6因为模式图3个节点可以任意排列映射到目标图的3个三角形节点上排列数是3! 6。第4个孤立点永远不会被选中因为它的度数为0而模式图度数至少为2在build_M阶段就被筛掉了。这个结果验证了算法的正确性。4.2 案例二路径模式匹配环图再看一个稍复杂的例子。模式图是长度为2的路径节点0-1-2目标图是4个节点的环每个节点跟相邻两个节点相连A [ [0, 1, 0], [1, 0, 1], [0, 1, 0] ] B [ [0, 1, 0, 1], [1, 0, 1, 0], [0, 1, 0, 1], [1, 0, 1, 0] ] results find_subgraph_isomorphisms(A, B) print(找到, len(results), 个映射) for m in results: print(m)这个结果应该是8个映射。你可以自己验证4环中任意选3个节点形成的子图都是路径一共有4个不同的3点集合每个3点集合里路径的两个端点互换方向又算两种映射所以总共4×28。对比案例一和案例二可以发现Ullmann算法不会重复搜索同一个映射因为有order固定了模式节点的分配顺序而used数组禁止两个模式节点选择同一个目标节点这一组合起来天然避免了重复和非法映射。4.3 案例三加节点标签约束实际项目里很少有不带标签的裸图。给Ullmann算法加上标签约束非常简单只需要在build_M里增加一个判断def build_M_with_labels(A, B, labels_A, labels_B): d, k len(A), len(B) deg_A [sum(row) for row in A] deg_B [sum(row) for row in B] M [[0] * k for _ in range(d)] for i in range(d): for j in range(k): if deg_A[i] deg_B[j] and labels_A[i] labels_B[j]: M[i][j] 1 return M标签可以代表任何业务含义化学元素类型、数据库里的实体类型、社交网络里的用户等级、代码语法树里的节点类型……只要你把节点属性编码成可比较的字符串或数字即可。这在真实场景里往往是收益最大的一步优化因为绝大多数节点对在初始阶段就被属性约束排除了Ullmann的细化反而只需要处理少数难啃的骨头。5. 复杂度分析与优化方向5.1 最坏情况为什么是组合爆炸子图同构是NP完全问题这句话的意思是在最坏情况下算法需要尝试所有可能的注入式映射数量上限是P(k, d) k! / (k - d)!当d和k都等于10时这个数大约是920万。还好Ullmann的细化阶段能砍掉大量候选但碰到下面这种极端情况细化也会失效模式图是空图没有任何边目标图是任意图。此时任何映射都满足边条件算法只能暴力枚举所有注入式映射。模式图和目标图都是完全图。此时度数约束不起作用细化条件也等于空转因为每个目标节点都能找到合适的邻居。遇到这类情况任何基于子图同构的算法都无法回避组合爆炸。所以工程上通常还有一招先判断图是否够“稀疏”如果目标图或者模式图太密集可能需要考虑完全不同的策略。5.2 提升搜索效率的几个方向实测下来Ullmann算法在以下四个方向上的优化收益最大选择度候选最少优先递归前先用sum(M[i])排序优先处理最难分配的行。这个策略几乎没有任何额外开销却能显著减少搜索树宽度。列使用跟踪用used数组禁止同一目标节点被映射两次这是保证“单射”的必要手段同时也能减少无效分支。细化提前终止细化过程中如果发现某一行全为0立即返回失败不再继续计算。这个在递归深层特别管用。矩阵复制时机递归参数用拷贝后的矩阵每层之间互不污染。很多初学者为了省内存试图用回退恢复矩阵结果状态管理一团糟反而更慢更易错。我的经验是小矩阵直接拷贝大矩阵再考虑回退在图的节点数不超过几百时列表拷贝的性能完全不是问题。还可以考虑引入numpy把细化过程的循环改成向量化布尔运算。不过我个人建议先保证逻辑完全正确再考虑性能优化。过早优化是万恶之源这句话在图算法里体现得淋漓尽致。5.3 与VF2等同类算法的选型对比Ullmann算法不是唯一的子图同构算法。网络x里内置的其实主要是VF2算法。这两个算法的定位差异很大对比维度UllmannVF2核心思路候选矩阵细化剪枝候选节点对可行性规则判定数据结构矩阵邻接表和状态向量稀疏图性能一般优秀稠密图性能表现稳定可能退化实现难度低高扩展性容易加标签、权重约束规则复杂扩展成本高简单说如果图比较稀疏、节点数很多优先选VF2如果图相对密集、节点数中等或者你需要灵活加入各种自定义约束标签、权重、类型Ullmann实现起来更顺手。如果你在学习阶段我更推荐先写一遍Ullmann因为它的逻辑直白剪枝思路清晰理解了Ullmann再看VF2会轻松很多。6. 常见问题与排查技巧实录6.1 细化剪枝“剪过头”了曾经有个朋友跑Ullmann发现结果少了。排查后发现问题出在细化条件他要求在模式节点i的每个邻居i2都必须有目标图里的j2与之对应但忘记了i2本身可能还没有被分配目标。这导致细化过程中很多还没被分配的候选被误删了。标准细化条件里M[i2][j2]的意思是“i2在当前候选矩阵里可能与j2匹配”就算i2最终不一定选j2只要存在这个可能性就能保住M[i][j]。这个“可能性”是细化剪枝的底气但也是最容易写错的地方。写的时候务必记住细化只删“必不可能”的候选不删“还没确定”的候选。6.2 回溯状态没有恢复我在3.4里提到过复制矩阵的方式。如果你试图不复制直接在原矩阵上“改一行再改回来”很容易踩雷因为refine不只是改当前行它会网格式地修改整个矩阵的其他行。等你回溯准备试当前行的下一个候选列时矩阵已经处于被上一分支污染过的状态。这里有一个通用的排查方法在所有递归入口处打印当前矩阵的快照如果发现同一层的两次递归进入时矩阵不一致十有八九是状态污染。几个常见的修复方式在递归前复制一份当前矩阵传给下一层这是我推荐的如果担心复制开销可以对refine修改过的位置做一个“撤销日志”回溯时按日志恢复实在不行就不要原地修改始终用不可变结构比如把矩阵转成元组来传递6.3 输出结果重复或漏解结果重复通常是因为没有固定搜索顺序导致同一种映射被不同路径访问到。解决办法就是咱们代码里的order先对模式节点排序搜索顺序固定下来就不会出现多个分支产生同一种映射。漏解的原因一般有三个build_M里初筛条件写得过强比如强行要求deg_A[i] deg_B[j]这会漏掉那些模式节点度数小于目标节点度数的合法映射。记住子图同构允许目标图里有多余边。细化里B[j][j2] 1的方向搞反了把“j的邻居”写成“j2的邻居”导致误删合法候选。最终验证时用了“完全相等”的判断把目标图的多余边也当成错误条件导致漏解。6.4 有向图、无向图与自环的处理Ullmann算法天然支持有向图只要邻接矩阵不要求对称就行。细化和验证条件里A[i][i2]1代表从i到i2有一条有向边B[j][j2]1代表j到j2有一条有向边方向语义完全一致。无向图只是邻接矩阵恰好对称代码不用改。自环的情况则需要小心。如果模式图里i有自环A[i][i]1那么细化时要求目标节点j也必须有自环B[j][j]1同时M[i][j]还得保持为1。我们在2.1约定对角线都是0所以没处理这个细节。万一你要处理带自环的图细化条件和最终验证里都要额外加一行对角判断否则可能误判。7. 个人经验与最终建议最后分享一点实操体会。我最初写Ullmann算法时总觉得“反正最后有验证前面剪枝随便写写就行”结果在中等规模图上跑了快十分钟都没出结果。后来老老实实把细化逻辑写完并且打上标签预剪枝同样的图不到半秒就出结果了。这个反差让我深刻意识到子图同构这类NP完全问题核心根本不是“枚举”而是“剪枝”。如果你今天就想在自己的项目里用上这段代码我给两条建议第一步把build_M里的初始约束尽量做足凡是业务上能确定的“不可能映射”都放进初始矩阵。这一步看似简单收益往往最大。第二步先用小规模测试用例把基本功能跑通再逐渐加大数据量。遇到性能瓶颈优先看细化是否生效在递归入口打印每层候选总数如果某一层候选数下降不明显说明剪枝力度不够。子图同构是一个很深的领域Ullmann算法只是起点。后面你可以继续研究VF2、LAD、GraphQL等更现代的算法也可以把眼光放到大规模图上的近似匹配。但无论走到哪一步回溯加剪枝这个组合拳的底层逻辑都是相通的。把这个Python实现吃透你对图算法、约束满足问题以及对“状态搜索”这类编程范式的理解都会上一个台阶。