形式陈述
给定一个对称正定稀疏矩阵 公理库 稀疏矩阵表示与运算 Sparse matrix 只存非零项及其位置,以稀疏模式组织矩阵运算、图结构、重排序与填充成本。 ,能否在计算任何平方根之前,就知道它的 Cholesky 因子 公理库 Cholesky 分解 Cholesky factorization · Cholesky decomposition 把 Hermitian 正定矩阵唯一分解为正对角下三角因子与其共轭转置的乘积。 需要哪些存储位置?这就是符号分解 的问题。
先固定对称排列 A ~ = Q A Q T ,以后编号 0 , … , n − 1 都指排列后的位置。用无向图 公理库 有限简单无向图 Graph · Finite simple undirected graph · 图 由有限顶点集与无序二元顶点子集组成的边集所确定的简单无向图。 表示原始非对角结构:位置 ( i , j ) 被声明为可能非零,就连边 i j 。对角位置始终保留。这里预测的是结构模式,不依赖恰好发生的数值相消。
按编号消去顶点 j 。设当前尚未消去的邻居集合为 N j ,则执行:
记录因子第 j 列的行集合 { j } ∪ N j 。
将 N j 中每对不同顶点连接,使其成为一个团;新添的边称为填充。
删除顶点 j 及其关联边,进入下一列。
运行结束,所有记录的列集合就是 L 的符号模式。输入还应指定原始结构零的解释及排列;输出至少包括列指针、每列行号和填充计数。数值分解再把实际值填入这些已分配的位置。
直觉
Schur 补 公理库 Schur 补与静态凝聚 Schur complement · Static condensation 把内部变量的精确响应折入界面矩阵与右端,推导 Schur 补、恢复公式及最小能量性质,并逐项凝聚五节点链。 说明,消去 j 会对剩余系数执行
a u v ← a u v − a u j a j v / a j j . 只要 u , v 都与 j 耦合,乘积就有可能非零。图算法把“有可能”记成一条边,不提前猜测它是否会与原项抵消。这使同一模式的多次数值分解能够共享一次符号分析。
图片加载失败 邻居补团与完整列模式 一条跨越多步的路径判据
对 i > j ,位置 L i j 是结构非零,当且仅当原图中存在从 j 到 i 的路径,其所有内部顶点编号都小于 j 。直接原边对应没有内部顶点的路径。
证明可随消元归纳。消去 k 时,新增边 u v 来自 u − k − v ;若两条已有边已经分别代表经过更小编号点的路径,把它们接起来就得到内部点不超过 k 的通路,再去掉重复顶点即可。反方向,给定一条内部点都将先被消去的路径,逐个消去这些内部点会把路径缩短,最终留下端点边。
这个判据解释了为何填充不只是“原矩阵两步路径”的一次统计。早先产生的边还会参与以后消元,把长路径的作用继续传递。
例子与边界
同一六节点图,两份完整账本
原边为六环 01 , 12 , 23 , 34 , 45 , 50 加弦 14 ,共七条边。自然顺序 0 , 1 , 2 , 3 , 4 , 5 给出:
当前顶点
尚存邻居
本轮新增边
0
1 , 5
15
1
2 , 4 , 5
24 , 25
2
3 , 4 , 5
35
3
4 , 5
无
4
5
无
5
空
无
所以逐列行号是
( 0 , 1 , 5 ) , ( 1 , 2 , 4 , 5 ) , ( 2 , 3 , 4 , 5 ) , ( 3 , 4 , 5 ) , ( 4 , 5 ) , ( 5 ) . 按列接起来有 17 项;从零开始的 CSC 列指针为 ( 0 , 3 , 7 , 11 , 14 , 16 , 17 ) 。原矩阵下三角含对角有 6 + 7 = 13 项,四条填充恰好使它增加到 17 。
改用原顶点顺序 π = ( 0 , 2 , 5 , 3 , 1 , 4 ) 。先消去 0 添边 15 ,再消去 2 添边 13 ;随后 5 的剩余邻居是 1 , 4 ,二者原本已经相连,3 的剩余邻居也相同,均无新增边。最后消去 1 , 4 。这次只有两条填充,共 15 个因子位置。
注意排列后列号已经改变。按新位置记录的列集合是
( 0 , 2 , 4 ) , ( 1 , 3 , 4 ) , ( 2 , 4 , 5 ) , ( 3 , 4 , 5 ) , ( 4 , 5 ) , ( 5 ) . 第一列中的行 2 , 4 分别对应原顶点 5 , 1 。如果一边用原顶点号画图,一边把它直接当 CSC 行号,就会把正确的消元过程存成错误矩阵。
结构存在,数值也可以恰为零
取正定矩阵
A = ( 1 1 1 1 2 1 1 1 2 ) . 原图是三角形,所以符号预测包含 L 21 。但消去第零列后,剩余非对角项为 1 − 1 ⋅ 1 = 0 ,实际因子为
L = ( 1 0 0 1 1 0 1 0 1 ) . 符号位置 L 21 可以存一个数值零。这是安全的结构上界,并不构成预测失败。反过来,不能因为这一次相消,就在后续系数变化时继续省略该槽位。
推论与应用
记 d j = | N j | 。因子的结构存储恰为
nnz struct ( L ) = n + ∑ j d j . 一次稀疏右看式数值消元要更新 d j × d j 的对称子块,因此工作尺度为 Θ ( ∑ j ( d j + 1 ) 2 ) ;若只数下三角乘减更新,次数是 ∑ j d j ( d j + 1 ) / 2 ,另有对角处理和缩放。上例两种顺序的平方和分别为 55 与 41 ,对称更新次数为 19 与 13 。非零数相近的两个因子,也可能因大列分布不同而有不同算术量。
朴素补团实现枚举这些邻居对。在哈希邻接集合的平均常数操作假设下,时间为 O ( n + m + ∑ j d j 2 ) ,保存图和列模式需要 O ( n + m + fill ) 空间;使用平衡集合则还需计入查找插入的对数因子。更好的符号算法利用消元树 公理库 消元树与稀疏列依赖 Elimination tree · Cholesky elimination tree 从因子每列首个下方结构位置构造父指针,证明非零必指向祖先,并用两条独立子链解释依赖、并行与祖先判据的单向性。 压缩依赖,不能把它们的复杂度直接赋给这份朴素实现。
排列的价值是改变未来邻居集,而不只是让图看起来整齐。嵌套剖分 公理库 嵌套剖分消元顺序 Nested dissection ordering 递归把二维网格分隔器排在子域之后,用完整填充账本与逐层前沿计数推导存储和算术界,并说明三维为何更昂贵。 把分隔器排到最后;不完全分解则主动拒绝保存某些填充,二者分别改变精确消元的次序与数值近似的规则。
本页依赖正定问题不需数值选主元。若一般 LU 为稳定性临时交换行列,实际消元顺序可能偏离预先计划,不能照搬固定排列下的精确模式结论。
参考资料
Yousef Saad, CSCI 8314, “Sparse Direct Methods”, slides 8-2–8-18:官方课程讲义 。符号与数值阶段、列模式传播和消元依赖。
Manpreet S. Khaira, Gary L. Miller and Thomas J. Sheffler, Nested Dissection: A Survey and Comparison of Various Nested Dissection Algorithms , CMU-CS-92-106R, 1992, §2 and Lemma 3.3:作者机构报告 。消元图与经过先消去节点的路径刻画。