Skip to content

算法Algorithm

不完全 Cholesky 预条件器

Incomplete Cholesky

在固定稀疏模式内执行近似 Cholesky,逐项追踪丢弃更新,并用精确正定却产生负主元的四阶例子解释失败与移位。

形式陈述 ​

精确 Cholesky 分解可能产生大量填充。若只允许因子保留一部分位置,能否用便宜的两次三角求解构造预条件器?这就是不完全 Cholesky 的目标。

设 A=AT≻0,先固定排列及允许的严格下三角模式 P。IC(0) 取 P={(i,j):i>j, aij 在原结构模式中},不增加原模式之外的因子位置。

为避免手算中反复出现平方根,将近似因子写成

M=L^DL^T,L^ 单位下三角,D=diag(dj).

若所有 dj>0,令 L=L^D1/2 就得到通常的 M=LLT。逐列计算为

dj=ajj−∑k<jl^jk2dk,l^ij=aij−∑k<jl^ikdkl^jkdj((i,j)∈P),

模式外固定为零。每个和只需要取两行已有模式的交集。若某个 dj≤0,应返回非正主元及其位置,而不是继续开方或把它取绝对值。

输入是矩阵、固定模式与排列;设置成功后,输出预条件应用接口:解 L^y=r,再令 w=D−1y,最后解 L^Tz=w。于是 z=M−1r。这是两次三角求解加一次对角缩放;不形成逆矩阵。

直觉

精确符号消元允许剩余邻居补成团。不完全分解在其中插入一条不同规则:某些新耦合即使会出现,也不为它保留因子位置。省掉的并非无用零,而是有意舍弃的数值作用。

保留模式、丢弃更新与失败主元

成功完成后,M 必正定,因为非零 v 满足

vTMv=(L^Tv)TD(L^Tv)>0.

不过这里的“成功”不能由 A≻0 自动推出。精确 Cholesky 的每个主元来自真正的正定 Schur 补;丢弃更新以后,正在分解的矩阵已不再是那个 Schur 补,原证明失去了前提。

保留位置上,递推使 Mij=Aij,包括全部对角位置;模式外则一般不相等。这个逐位置约束不意味着 M 在任何矩阵范数下都是最佳近似,也不保证迭代次数必然减少。

例子与边界

九节点网格中的第一次丢弃 ​

取三乘三内部网格,按行编号 0,…,8,矩阵对角为 4,水平与竖直邻接项为 −1,其余为零。消去节点 0 时,后继邻居为 1,3。精确更新会产生

a31new=0−(−1)(−1)/4=−1/4.

IC(0) 不允许位置 (3,1),于是拒绝这一更新;对角更新仍然保留。因此 d0=4,l^10=l^30=−1/4,接下来的 d1=d3=15/4。这也意味着重构矩阵的 M31=1/4,故 (A−M)31=−1/4,明确显示了近似误差落在哪里。

完整九个主元依次为

4, 154, 5615, 154, 5215, 2507728, 5615, 2507728, 85722507.

均为正,设置成功。对 b=(1,2,3,1,2,3,1,2,3)T、x0=0,比较原残差的 Euclidean 范数:

更新步数 普通 CG IC(0) 预条件 CG
0 6.480741 6.480741
1 3.376736 0.865117
2 1.256596 0.047161
3 0.460952 0.002639
4 0.093822 0.000051

绝对残差阈值 10−4 下,预条件版本四步达标,普通版本需五步。本例还保留一个容易忽视的区别:精确分数算术中,普通 CG 恰好五步终止,而 IC(0) 版本恰好七步才达到严格零残差。预条件打破了一些原有谱重合;“较早达到实际容差”与“更少步精确终止”不是同一句话。

正定输入也能发生精确 breakdown ​

取

A=(13/503/53/513/5003/51−3/53/50−3/51).

写作 A=I+(3/5)C,其中 C2=2I。因此特征值为 1±32/5,全都严格为正。

精确分解的 D 为 (1,16/25,7/16,7/25)。IC(0) 则把位置 (3,1) 固定为零,前三个主元仍相同,但随后

l^32=−3/57/16=−4835,d3=1−925−(4835)2716=−32175.

这是精确负数,与浮点舍入无关。若改为对 A+15I 构造 IC(0),主元成为 (6/5,9/10,4/5,9/20),全部为正。移位修复改变的是预条件器,外层仍求解原来的 Ax=b;它不是把原问题也悄悄替换成移位系统。

推论与应用

固定成功因子给出的 M−1 是线性、自伴、正定算子,可接入标准 PCG。若迭代中重建模式、改变移位或使用输入相关的内层停止规则,就需要重新核对外层算法允许的接口。

设允许模式每列下方有 dj 项,因子存储为 O(n+∑jdj),每次应用同阶。设置需要计算保留位置的内积交集;用右看式邻居对更新及平均常数时间模式查询,可以用 O(∑j(dj+1)2) 作为算术和查询上界。设置不能无条件只写成与非零数成正比。

加入填充、改变排列或使用对角移位,会同时改变稳定性、近似质量和成本。比较时应报告失败状态、设置时间、因子大小、原残差达标步数及总时间。某些具有额外符号和占优结构的矩阵可保证不完全分解存在;一般正定矩阵则必须保留实际主元检查。

参考资料
  • Yousef Saad, Iterative Methods for Sparse Linear Systems, 2nd ed., 2003, §§10.3.1–10.3.2 and §10.8.2:作者教材。固定模式丢弃、残余位置关系,以及一般正定矩阵上 IC 的存在性边界与移位。
关系图谱8 个相邻概念 · 2 类关系

拖动节点调整位置。

显示关系

显示:依赖

  1. 前置三跳
  2. 前置二跳
  3. 前置一跳
  4. 当前条目
  5. 后续一跳
  6. 后续二跳
  7. 后续三跳
文字版关系按与当前条目的最短距离分组
类型化关系