2026 年 4 月：从稀疏 Cholesky 分解推导消元树，解锁矩阵分解新路径
稀疏 Cholesky 消元树
2026 年 4 月 9 日，将推导用于计算 `A = LL^T`（其中 `L` 为下三角矩阵，`A` 为稀疏矩阵）的（右视）稀疏 Cholesky 算法的消元树。即便 Cholesky 算法的基本假设（对称性和正定性）不成立，这棵树也是大多数稀疏分解软件的基础。最终，这棵树能告诉我们两件事：1. 即便原矩阵 `A` 中不存在非零元素，矩阵 `L` 中哪些位置会出现非零元素（即“填充”）；2. 分解结果的任务依赖图。传统上，这个概念通常是在稀疏三角求解的背景下提出的，然后将其作为 Cholesky 分解的构建块。而想直接从 Cholesky 分解入手，下面将详细阐述。
消元树的动机
假设首先从下面的纯稠密右视 Cholesky 算法伪代码开始：
// 输入：对称正定矩阵 A，赋值给输出矩阵 L // 输出：L 的下三角部分满足 A = LL^T for (int k = 0; k < n; ++k) { // 1. 分解主元 L[k][k] = sqrt(A[k][k]); // 2. 缩放主元下方的列 for (int i = k + 1; i < n; ++i) { L[i][k] /= L[k][k]; } // 3. 对后续矩阵进行右视秩 1 更新 // 注意，只有这一步会产生填充。 for (int j = k + 1; j < n; ++j) { for (int i = j; i < n; ++i) { // 仅下三角部分 L[i][j] -= L[i][k] * L[j][k]; } } }这种隐含的顺序和循环结构形成了一个任务有向无环图（DAG）。当 `A` 是稀疏矩阵时，事实证明，这个 DAG 以及步骤 (3) 中隐含的填充可以完全由 `A` 的初始非零模式和一个简单的树数据结构（即“列消元树”）确定。用一个简单的 `5x5` 矩阵（为便于阅读，仅显示下三角部分）来说明：
完整的 Cholesky 任务图
对于上述稀疏矩阵，得到：
正如预期的那样，由于 `A` 是稀疏矩阵，这里有些操作根本不需要执行。如果修剪掉这些操作，就会得到一个小得多的图。
注意，这里的列分组目前没有实际用途，但可以将列分组理解为“当完成列组 `i` 中的所有操作后，后续操作不会再引用列 `i`”。
修剪后的（稀疏）Cholesky 任务图
消除所有不必要的操作后，得到以下结果：
这个修剪后的图告诉我们，为了对这个稀疏矩阵进行完整的 Cholesky 分解，需要执行的所有操作。但是，为了确定这个图，不得不进行一次“稠密”分解，然后修剪掉不必要的计算。
然而，事实证明，有一种简单的 `O(n)` 数据结构，结合 `A` 的初始非零模式，可以快速告诉我们最终的填充模式以及修剪后的任务图。这就是列消元树。
列消元树
回顾稠密 Cholesky 分解中的步骤 (3)：
// 3. 对后续矩阵进行右视秩 1 更新 // 注意，只有这一步会产生填充。 for (int j = k + 1; j < n; ++j) { for (int i = j; i < n; ++i) { // 仅下三角部分 L[i][j] -= L[i][k] * L[j][k]; } }这告诉我们一个关键的结构事实：如果 `k < j <= i`，`L[i][k]!=0` 且 `L[j][k]!=0`，那么 `L[i][j]!=0`。
现在，当查看修剪后的任务图时，可以看到一些隐含的列依赖关系，如 `0->1` 和 `0->2`。由于 `0` 有两个“父节点”，这是一个 DAG 模式，而不是树。但是，上述结构规则告诉我们，`0->2` 是一条冗余边，因为 `0->1` 必然意味着 `L[1,0]!=0`，`0->2` 必然意味着 `L[2,0]!=0`，因此必须有 `L[2,1]!=0`，所以存在一条隐含边 `1->2`，`0->2` 是冗余信息，因为无论如何都必须先从 `0->1` 开始。如果从列 DAG 中移除所有冗余边，结果就是一棵树：
消元树告诉我们如何根据 `A` 的初始非零模式和结构规则 `k < j <= i`，`L[i][k]!=0` 且 `L[j][k]!=0` 意味着 `L[i][j]!=0` 来生成 `L` 的填充模式。
将把消元树的计算推迟到本文末尾，现在展示如何使用编码为 `parents` 数组的消元树：
使用消元树
符号分解
// 输入： // n // A_row[i] : A 的原始下三角行模式 // parent[k] : k 在消元树中的父节点，或 -1 // mark[k] : 调用者分配的长度为 n 的临时数组 // // 输出： // L_col[k] : 存在 L[i][k] 的行索引 i >= k // // 临时变量： // mark[k] 会被覆盖。 void symbolic_cholesky( int n, int** A_row, int* A_row_len, int* parent, int* mark, vector_int* L_col ) { for (int k = 0; k < n; ++k) { mark[k] = -1; vector_clear(&L_col[k]); // 包含对角线元素。 vector_push(&L_col[k], k); } for (int i = 0; i < n; ++i) { for (int p = 0; p < A_row_len[i]; ++p) { int k = A_row[i][p]; // 基础非零元素 A[i][k]，k < i while (k != -1 && k < i && mark[k] != i) { mark[k] = i; // 树路径表明 L[i][k] 存在。 vector_push(&L_col[k], i); k = parent[k]; } } } }数值分解
一旦得到了 `L` 的非零模式，使用适当的稀疏矩阵数据结构，根据稠密算法大纲编写数值分解就很直接了。关键步骤是引入填充的步骤，下面是伪 C 代码：
for (int j in L_col[k]) { if (j <= k) continue; for (int i in L_col[k]) { if (i < j) continue; L[i][j] -= L[i][k] * L[j][k]; } }由于预先计算了非零模式，因此无需检查这些数值条目是否存在，它们会被预先分配。
计算消元树
消元树本质上是：
parents[k] = min {j>k s.t. L[j,k]!=0}为了仅使用 `A` 的初始非零元素来计算消元树，而不进行过多的额外工作，需要逐步构建 `parents`。这意味着对于每个 `r`，需要尝试在这棵树中找到通向 `r` 的路径。
因此，如果按顺序遍历 `A` 的行，然后查看列 ID，会得到一个候选路径 `c-> ... -> r`。将 `c->r` 保存在 `ancestors[c]` 中，并沿着祖先节点向上迭代，直到找到一个尚未匹配的列，然后该列的父节点变为 `r`，因为是按顺序遍历行的，所以不可能找到更小的 `r`。
以下是伪 C 代码：
// 计算列消元树。 // // 输入： // n // A_row[r] : A[r,c] 在结构上非零的列 c < r // A_row_len[r] : 行 r 中此类列的数量 // // 输出： // parent[k] : 列 k 在消元树中的父节点， // 如果 k 是根节点，则为 -1 // // 临时变量： // ancestor[k] : 调用者拥有的长度为 n 的临时数组 // // 哨兵值： // -1 表示 “尚未分配父节点/祖先节点” void compute_elimination_tree( int n, int** A_row, int* A_row_len, int* parent, int* ancestor ) { for (int k = 0; k < n; ++k) { parent[k] = -1; ancestor[k] = -1; } for (int r = 0; r < n; ++r) { for (int p = 0; p < A_row_len[r]; ++p) { int c = A_row[r][p]; // A[r,c] 存在，且 c < r // A[r,c] 给出了祖先约束： // // c -> ... -> r // // 从 c 开始，沿着已发现的路径向上遍历。 int v = c; while (v != -1 && v < r) { int next = ancestor[v]; // 记录在处理行 r 时，v 已向上连接到 r。 ancestor[v] = r; if (next == -1) { // 此路径在 v 处停止。由于行是按升序处理的，r 是第一个需要 v 作为祖先的后续行，因此 r 是 parent[v]。 parent[v] = r; break; } // 继续沿着先前发现的路径遍历。 v = next; } } } }总结
消元树的介绍通常涉及一些初步的图论知识，这些知识总让人感觉与线性代数有些脱节。目的不是取代图论，而是让它更贴近底层算法。直接从右视 Cholesky 分解推导出任务 DAG，对该 DAG 进行修剪，然后表明算法步骤 (3) 所隐含的结构规则导致了任务 DAG 的压缩表示，而这个压缩表示恰好是一棵树。