稀疏Cholesky消除树

稀疏 Cholesky 消除树 2026 年 4 月 9 日在这里，我推导了（右视）稀疏 Cholesky 算法的消除树，用于计算下三角 L 和稀疏矩阵 A 的 A = LL^T 。即使 Cholesky 的基本假设（对称和定）不适用，该树也构成了大多数稀疏分解软件的基础。最终，这棵树告诉我们两件事：1. 非零值出现在矩阵 L 中，即使原始 A 中不存在（即“填充”）；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. 尾随矩阵的右视Rank-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 的初始非零模式和结构规则生成 L 的填充模式： k < j <= i ， L[i][k]!=0 和 L[j][k]!=0 意味着 L[i][j]!=0 。我将把消除树的计算推迟到本文的末尾，并展示如何使用编码为父数组的消除树： 使用消除树符号分解 // 输入： // n // A_row[i] ： A 的原始下三角行模式 // Parent[k] ： k 的消除树父级，或 -1 // mark[k] ： 调用者分配的长度为 n 的暂存数组 // // 输出： // L_col[k] ： 行索引 i >= k ，其中 i >= k L[i][k] 存在 // // Scratch: // 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 ;矢量_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 ]; }