AMD近似最小度排序：稀疏矩阵分解前的关键预处理

做仿真的和数值计算的工程师，打交道最多的矩阵之一就是“稀疏对称矩阵”。刚度矩阵、质量矩阵、电路方程里的节点导纳矩阵，几乎天天见。这类矩阵规模可以到几十万甚至上千万阶，但每行非零元往往只有几十个。直接做 Cholesky 或 LU 分解之前，几乎所有人都会先做一步重排序，而这里面最常被推荐、也最经典的一个算法，就是 AMD——Approximate Minimum Degree Ordering Algorithm，翻译过来叫近似最小度排序算法。

一句话讲明白 AMD 解决什么问题：它通过重新排列矩阵的行列顺序，让稀疏矩阵在分解过程中少产生填充元（fill-in），也就是让分解因子里的非零元素尽量少、内存尽量省、计算尽量快。这个需求听起来很朴素，但背后牵扯到图论、贪心算法和超大规模数据结构的工程实现。而且排序选得好不好，能把一次分解的耗时和峰值内存拉开一个数量级，真不是开玩笑。

下面我从问题本身讲起，把这个算法掰开揉碎了聊，同时把实际使用中容易踩的坑一并整理出来。内容按“能直接看懂、能直接上手”的标准来组织，适合刚接触数值线性代数的学生，也适合被大型矩阵折磨已久、想回头把预处理环节彻底搞明白的工程师。

1. 为什么排序这么重要：填充元问题

1.1 先搞明白“填充元”是什么

假设我们有对称正定矩阵 A，要做 Cholesky 分解。一个稀疏矩阵在分解之前，很多位置本来就是零，比如行 1 和行 3 之间没有直接连接，按零存储完全没有负担。但在消元过程中，原本的零位置会被动地变成非零值。这些新出现的非零值，就叫填充元。

举个具体的例子。一个 4×4 矩阵，非零模式如下：

x x 0 x x x x 0 0 x x x x 0 x x

自然顺序是节点 1、2、3、4 依次消去。消去节点 1 的时候，节点 2 和节点 4 会被“拉在一起”，因为它们都和节点 1 相邻，于是位置 (2,4) 就产生了一个新的非零元。这个位置原来明明是零，结果因为消元，硬生生被填上了。接下来消去节点 2 的时候，节点 4 的邻居又多了一个节点 3，位置上又可能产生新的填充。一步传一步，最后矩阵往往变得比原来稠密得多。

可以把它想象成聚会上的熟人网络：你拉走一个人，他认识的两拨人从此就得互相认识。人越早拉走，他的熟人越少，新认识的组合就越少。反过来，如果先把社交达人拉走，他带出来的新关系能把整个场面搅成一锅粥。填充元的传播就是这样一种“连锁反应”。

1.2 填充元多起来会发生什么连锁反应

填充元不是数字上的零变非零这么简单，它有三层影响：

第一是内存。分解过程中产生的 L 因子必须把所有非零位置存下来，包括填充元。原来的稀疏矩阵可能只有 1% 的非零密度，但分解因子可能达到 20% 甚至更高。千万阶的矩阵，每多 1% 的密度，可能就是几个 GB 的内存差距。

第二是计算量。数值分解的每一轮消元，都要对因子中的所有非零位置做乘减运算。填充元多，意味着你要花大量时间去算那些“本来不存在”的数值。对于 Cholesky、LU 这类分解，填充元的数目和浮点运算量几乎是线性挂钩的。

第三是缓存和并行效率。稀疏矩阵本来就因为访存不连续而难以发挥 CPU 性能，如果填充元把因子结构撑得很散，缓存命中率会更低，超节点也难以形成。很多时候你会发现矩阵明明不算大，但分解就是慢得离谱，八成就是填充元太多、数据结构太乱。

1.3 重排序到底在做什么

所谓排序（ordering），就是找到一个置换矩阵 P，把原来的方程 Ax = b 改写成 (PAPᵀ)(Px) = Pb。A 和 PAPᵀ 本质上代表同一个线性系统，只是把行和列都重新换了个顺序。

为什么换顺序就能减少填充元？因为 Cholesky 分解的消元顺序就是矩阵的行列顺序。你的矩阵第 1 行先消，就相当于图里节点 1 先被删掉。而哪个节点先删、哪个后删，直接决定了填充元产生的路线图。排序算法做的工作，就是在 n! 种排列里找一个能让填充元尽量少的顺序。

当然 n! 是天文数字，不可能穷举。实际用的都是一类启发式搜索算法，AMD 就是其中成绩最好、应用最广的代表。

2. AMD 的核心逻辑：度、贪心与近似

2.1 最小度排序的直觉：谁是社恐先删谁

把对称矩阵的非零模式看成一个无向图 G = (V, E)。矩阵有多少行，图就有多少个节点；矩阵里 (i, j) 位置非零，图里节点 i 和节点 j 之间就有一条边。矩阵做 Cholesky 分解，等价于一个节点一个节点地把图“消掉”。

每次消去一个节点 v，v 的所有邻居都要两两连边。如果 v 的邻居少，这个过程中新增的边就少；如果 v 的邻居多，新增的边会成倍增长。所以最直接的贪心策略就是：每一步都选择当前“度数”（邻居数量）最小的节点消去。这种策略被称为最小度（Minimum Degree）排序。

持续选度数最小的节点，直观上就是在反复问：剩下的人里谁认识的人最少，那就先让谁退场。认识的人少，退场时引出的新关系最少。这个策略很朴素，但效果出奇得好，长期霸占稀疏矩阵排序的主流位置。

2.2 精确计算度数的代价为什么扛不住

如果只是“选度数最小的节点”，那算法本身并不复杂。麻烦在于：每次消去一个节点，图的边要更新，所有邻居节点的度数都要重新计算。而且被消掉的节点还可能形成一个中间结构，度数的定义也会变得麻烦。

朴素实现下，每次消去都要扫一遍全部未消节点，算一遍各自的邻接情况，整体复杂度会接近三次方的量级。矩阵规模上万的时候还能忍，上百万简直是一场灾难。所以要真正落地到大规模场景，必须把度数计算的路径压缩到极致，哪怕牺牲一点点排序质量的确定性。

2.3 AMD 的“近似”到底近似在哪儿

AMD 全名的关键就在“Approximate”。它做的最小度并不是“每一步都精确更新全图所有节点度数再选最优”，而是允许部分节点的度数信息落后于图的真实状态。核心技巧可以概括为三个层面。

第一个层面是度数的延迟更新。图结构每次消元之后都发生了变化，但 AMD 不会立刻把每个节点的度数全部刷新。它只在某些关键节点附近做局部更新，让旧度数继续用一段时间。度数的“不新鲜”会带来一定的排序质量损失，但换来的是更新成本的大幅下降。

第二个层面是吸收（absorbance）。图里经常出现这样的情况：某个节点 w 的未消邻居，完全包含在另一个节点 v 的邻域里面。这时候 w 的度数对消元顺序基本没有独立影响，不管 v 选不选，w 的边都会随着 v 的消去看不见。AMD 会把 w 直接“吸收”进 v 的势力范围，后续更新度数时根本不用处理 w，进一步减少操作量。

第三个层面是超节点和元素的压缩表示。消掉 v 之后，v 的邻居两两相连，形式上形成了一个小团。如果真把这些新增边一条条地写进邻接表，坏的代价依旧很大。AMD 的做法是把被消节点当成一个“元素”记录下来，谁跟它相邻就记一条引用，而不是把两两之间的边全部建出来。这样图的更新从“O(邻居数²)”降成“O(邻居数)”量级。

这三个机制叠加在一起，让排序算法的整体复杂度从立方级别降到了近乎线性，而排序结果与精确最小度相比，填充元数目只是略有增加，工程上完全可接受。

2.4 AMD 的完整流程长什么样

用大白话描述一遍典型实现流程：

1. 读入矩阵的稀疏结构，把非零模式转换成图的邻接表。
2. 如果需要，进行对称化处理，确保行列结构一致。
3. 初始化度数组，给每个未消节点记录邻居个数。
4. 进入贪心循环：从当前未消节点里挑出度数最小的那个，把它确定为下一号消元顺序。
5. 将这个节点标记为已消去，更新它邻居之间的连接关系。
6. 利用吸收规则合并可吸收的节点，更新受影响节点的近似度。
7. 重复 4-6 步，直到所有节点都排完。
8. 输出排列，完成从原矩阵到重排后矩阵的对应关系。

这个流程写出来只有几行，但工程实现里每一步都藏着大量细节。我见过不少自己照着论文写排序代码的项目，写出来能用，但一旦矩阵规模上去就玩不转，原因通常都出在数据结构和度数更新的设计上。

3. 工程实现必须抠的三个细节

3.1 度数的“桶结构”比你想的重要

最小度算法需要反复做两件事：取当前度数最小的节点，以及修改某个节点的度数。如果每次取最小都遍历一遍所有节点，那就回到了 O(n²) 的老路。

经典实现里最常用的数据结构是按度数分桶的链表。具体说就是开一个数组，下标代表度数，每个下标挂一个链表，链表里串着当前度数等于这个值的所有节点。取最小度节点时，只要从度数 0 的桶开始往前找，找到第一个非空桶，取走桶头节点就行。更新节点度数时，把它从原桶摘下，挂到新桶的头部。这些操作都是 O(1) 的。

我第一次用这个结构是重写一个排序实验，起初图省事用了优先队列，结果节点的度数频繁变化，堆里的陈旧信息越来越多，逻辑变得相当痛苦。换成桶结构以后，不仅速度上来了，连代码思路都清爽很多。从这以后我对“数据结构决定算法可落地性”这句话有了切身体会。

3.2 填充边不值得一条条写进图里

前面提到超节点和元素，这个设计我认为是整个 AMD 实现里最精彩的部分。如果不做压缩，每消去一个高精度节点，它的邻居集合两两相连，就是 k 条节点的完全图，要新增 k(k-1)/2 条边。一个万度节点直接让邻接表爆炸。

AMD 的处理方法绕开了物化所有边：已消去节点 v 对应一个元素，记录“当前还活着的邻居有哪些”。两个存活节点之间有没有边，不需要真正在邻接表里画出来，只需要判断它们是否同时属于某个已消元素的邻居集合。度数更新的核心，也变成了维护“哪些活节点挂着哪些元素”的统计。

这有点像数据库里的“物化视图”和普通视图的区别：把每条边都显式存下来，维护起来又慢又贵；存一个抽象的关联结构，查询时再算，反而省。AMD 牺牲了一点边查询的直接性，换取整体更新的高效，这笔账算得非常划算。

3.3 内存预分配和复用直接决定成败

大型排序算法跑起来以后，最怕的是遍布代码的动态内存分配。邻接表会增长，节点会被移来移去，如果每步都 new 一个数组，程序会大量时间耗在申请和释放上，而且内存碎片化会让数据访问模式变得更差。

稍微成熟点的实现都会做两件事：一是预先分配足够容纳原始边和估计填充边的内存池；二是把节点移出桶、移入桶的操作设计成移动指针而非复制数据。这样排序过程中几乎没有真正的内存申请，全部是已有内存的重复利用。

我自己调试过一个上百万节点的模型，用简陋实现跑排序，光分配内存就占掉整个排序时间的四成。后来改成预分配与复用方案之后，那部分开销几乎没有存在感。如果自己写排序相关代码，开局就先做内存池设计，后面能省下大量排查时间。

4. 什么时候用 AMD：场景选择和替代者

4.1 对称正定矩阵：AMD 的主场

AMD 最理想的应用对象是结构对称、数值对称正定的矩阵，典型代表就是有限元方法形成的刚度矩阵。这类矩阵的填充元会天然对称，用 AMD 做 Cholesky 分解前的重排序，效果非常显著。

实际项目里我经常做的一件事是：拿到网格模型后直接告诉中游的同学，先对刚度矩阵做一次符号排序，再看非零元的统计。很多时候只做这一步，分解所需内存就掉了一个数量级。这个排序和数值计算本身无关，纯粹是给分解器一个更好的顺序，但收益往往比升级硬件还明显。

4.2 非对称矩阵：不要照搬 AMD，要转向 COLAMD

如果你的矩阵结构不对称，直接用 AMD 做行和列的同步排列会有点拧巴。因为 LU 分解的填充元结构不对称，行排列和列排列各自影响不同。这时候应该使用列近似最小度排序（COLAMD），它对列进行单独排序，让 LU 因子尽量稀疏。

判断标准也很简单：你用的分解器吃的是对称结构还是非对称结构。Cholesky 对应的必须是对称排序，LU 一般配合列排序。很多工程库的默认参数会替你选好，但如果你自己控制预处理流程，千万别拿 AMD 的排列去硬套 LU 分解。

4.3 超节点并行：AMD 不一定总是最佳答案

AMD 在串行 Cholesky 里表现接近天花板，但它的消元顺序很“深”，适合按单条依赖链跑。在大规模并行直接法里，嵌套剖分（Nested Dissection，简称 ND）往往更吃香，因为它先把大规模矩阵分成互不相干的子块，各子块可以独立消除，并行度高得多。

实际工程中不少求解器会做混合策略：先用嵌套剖分把矩阵分到若干子域，再在每个子域内部用 AMD 做细粒度排序。这种“顶层剖分、底层最小度”的组合几乎成了大型稀疏直接法的标准配方。如果你发现算法换了之后并行效率上不去，可以考虑是不是排序策略根本没换。

5. 主流排序算法怎么选：一张表理出头绪

从实测经验来说，中小规模矩阵（一万阶到几十万阶），AMD 几乎无脑好用；上百万元素且并行资源充分的场景，优先尝试嵌套剖分与 AMD 的混合排序。RCM 更多是给“带状结构明显、想靠带宽压缓存”的场合，比如一些隐式差分网格的排序，直接用它做预处理成本不高。

很多成熟的直接法求解器内部已经把这几种算法都集成好了，你只需要了解它们的差异，好判断自己的场景是否需要手动干预默认策略。

6. 常见问题与排查技巧实录

6.1 排序之后填充元依旧很大

不少人第一个想到的就是排序算法不行，但更常见的原因是送入排序的非零结构不干净。

排查的第一步是检查有没有显式零。矩阵里如果存在“数值是零但存储位置非零”的项，排序会把它们当作真实边看待，导致图和实际数值结构不对称。可以用矩阵查看器或者脚本扫一遍，把所有绝对值等于零的项剔除掉再排序，往往填充元统计立刻大幅下降。

第二步是确认结构对称性。如果矩阵的数值不对称，但你打算做 Cholesky 分解，排序前得用 A+Aᵀ 的结构来排序。有人直接拿 A 的单边结构去做对称排序，结果因子填充元凭空多出一大截。

第三步是检查编号和映射有没有错位。做过网格重编号的项目，特别注意行列映射是否对齐。最典型的故障就是行用了新编号、列用了旧编号，所有排序努力全部白费，填充还是不降。

6.2 分解仍然慢，问题不在排序而在后续

排序质量很好，但 Cholesky 还是很慢，这时候要看看分解器有没有真正利用排序的“超节点结构”。AMD 排完序后，大量连续节点会被压缩成超节点，分解器应当用稠密分块算法对超节点做 BLAS3 级的计算。如果你的分解器只是做逐元素稀疏追赶，那么再好的排序也发挥不出来。

另外，检查符号分解是否与数值分解分开进行。一般流程是先做一次符号分解，把 L 因子的结构确定下来，分配好内存，再做数值分解。如果每次分解都从头分析结构，排序带来的收益就会被重复分析的开销抵消掉大部分。

6.3 峰值内存超出预期怎么办

遇到内存爆炸，先要做的是收集符号因子的统计量，也就是排序之后的 L 因子总非零数目，再估一下按双精度存储需要多少字节。对比机器的实际内存就知道瓶颈在哪。

如果差距不大，多半是稀疏结构里存在一个特别高密的子矩阵，导致局部填充急剧增多。这时候可以先检查“是不是有一小撮节点度数异常”，很多情况下是网格的单元连接设置出错，把本不该相连的节点连成了一个大团。

如果差距很大，考虑换用递归终止更早的嵌套剖分或者 AMD 与 ND 的混合方案，同时把分解器的内存缓冲参数调小。不要一上来就改排序算法，先定位是不是数据本身有问题，能省掉大量瞎试时间。

6.4 不同平台排序结果不一致

有些团队在 Windows 和 Linux 两个平台上跑同一套代码，发现排序结果有差异。这种现象多半来自哈希表、指针地址或并行归并时的不确定性。排序结果只要都能有效降低填充元，细微差异其实无伤大雅。

但如果项目要求可复现性，可以规定统一使用单线程排序、统一数据装载顺序，并且在排序结果写盘时把轮次和序号也记录下来。这样每次重跑结果就完全一致，遇到问题也能快速回溯。

7. 实操心得：一次内存救场的经历

有一回我处理一个五十万自由度的结构力学模型，刚度矩阵本身只有不到三千万个非零元，感觉规模不大，结果直接上稀疏 Cholesky 分解，峰值内存冲到八十多个 GB，机器直接卡死。后来把矩阵结构导出来，用 AMD 做了一次排序，符号分解的结果让我记忆犹新：L 因子非零元从预期的一亿两千万掉到一千八百万左右，峰值内存降到二十多 GB，整个分解时间也缩短到原来的三分之一。

这次经历让我彻底扭转了对预处理环节的态度。以前总觉得排序是可有可无的“加速技巧”，后来意识到它其实是稀疏直接法里决定生死的一步。矩阵结构乱七八糟的时候，硬件配置再高也扛不住指数式增长的填充元。

从那以后我的项目流程里固定加了一环：所有大规模稀疏矩阵，在进入分解器之前先做符号阶段的排序和质检，看一眼填充统计再往下走。这一步只要花几分钟，却能在后面省下几小时甚至几天的运行时间。

最后还想提醒一句：AMD 再强也不是银弹。遇到对称正定、万级到百万级、希望尽快出结果的场景，优先试它；遇到非对称矩阵，改用 COLAMD；遇到超大规模并行场景，考虑嵌套剖分打底、AMD 收尾。排序算法的选择没有绝对正确答案，只有结合你的矩阵形态、求解器和机器资源去综合判断，才能把性能榨到极限。