形式陈述
预条件共轭梯度法(PCG)求解 A x = b ,其中 A = A ∗ ≻ 0 。除矩阵—向量乘 v ↦ A v 外,算法需要一个固定的线性接口
apply ( v ) = M − 1 v , M = M ∗ ≻ 0. 这里的Hermitian 正定 公理库 正定与半正定矩阵 Positive definite matrix · Positive semidefinite matrix · PSD matrix 由二次能量严格为正或非负定义的实对称与复 Hermitian 矩阵。 条件同时约束原矩阵与预条件器。实现可以只提供 apply,不显式存储 M 或其逆;所提供的映射必须在所有迭代中保持同一个线性、自伴、正定算子。M 近似 A 是为了改善收敛,M − 1 v 易于计算则是获得实际加速的前提。
给定 b , x 0 ,先计算真实残差 r 0 = b − A x 0 。若它已满足停止准则,就直接返回 x 0 ;否则初始化
z 0 = M − 1 r 0 , ρ 0 = r 0 ∗ z 0 , p 0 = z 0 . 第 k 步先计算
v k = A p k , α k = ρ k p k ∗ v k , 再更新
x k + 1 = x k + α k p k , r k + 1 = r k − α k v k . 若新近似通过真实残差检查,则结束;若继续迭代,则计算
z k + 1 = M − 1 r k + 1 , ρ k + 1 = r k + 1 ∗ z k + 1 , β k = ρ k + 1 ρ k , p k + 1 = z k + 1 + β k p k . 精确算术中,只要尚未得到解,ρ k > 0 、p k ∗ A p k > 0 ,步长和继续迭代时的 β k 都是正实数;复数数据必须使用共轭转置。M = I 时,这些公式逐项退化为普通共轭梯度法 公理库 共轭梯度法 Conjugate gradient method · CG method 在 Hermitian 正定系统的 Krylov 子空间中最小化能量误差,以三项递推获得共轭方向和短存储迭代。 。两处标量递推都使用 ρ k = r k ∗ M − 1 r k ,不能只把方向中的 r k 换成 z k ,却保留普通 CG 的 r k ∗ r k 。
接口还应接受非负的绝对、相对容差和工作预算,并返回原变量 x 、真实残差、迭代次数、A 与 apply 的调用次数及退出状态。本文采用的成功准则是
‖ b − A x ‖ 2 ≤ atol + rtol ‖ b ‖ 2 . 这一定义对 b = 0 仍有意义;若此时也取 atol = 0 ,要求的便是零残差。内部的 ρ k 不是这个停止量。
直觉
预条件 公理库 预条件 Preconditioning · Preconditioner 用易应用的近似逆改变等价线性系统的尺度与 Krylov 几何,并计入构造、存储和每步应用成本。 改变 CG 观察问题的坐标。取分解 M = C C ∗ ,令
B = C − 1 A C − ∗ , c = C − 1 b , y = C ∗ x , C − ∗ = ( C ∗ ) − 1 . C 非奇异,且对非零 u 有
u ∗ B u = ( C − ∗ u ) ∗ A ( C − ∗ u ) > 0 , 所以 B 仍是 Hermitian 正定矩阵,可以对 B y = c 应用普通 CG。以下推导使用普通 CG 已有的正交性与最小化定理,证明它如何成为上述 PCG;并不另行重证普通 CG 定理。
记变换坐标中的残差和搜索方向为 s k , q k ,映回原变量时定义
s k = c − B y k = C − 1 r k , p k = C − ∗ q k . 普通 CG 的初始方向 q 0 = s 0 因而给出
p 0 = C − ∗ C − 1 r 0 = M − 1 r 0 = z 0 . 它的步长分子和分母分别变成
s k ∗ s k = r k ∗ C − ∗ C − 1 r k = r k ∗ z k = ρ k , q k ∗ B q k = p k ∗ A p k . 因此 α k = ( s k ∗ s k ) / ( q k ∗ B q k ) 正好给出 PCG 的步长。对普通更新 y k + 1 = y k + α k q k 左乘 C − ∗ ,便得到原变量的 x 更新;对 s k + 1 = s k − α k B q k 左乘 C ,便得到 r k + 1 = r k − α k A p k 。剩下两条普通 CG 递推为
β k = s k + 1 ∗ s k + 1 s k ∗ s k = ρ k + 1 ρ k , q k + 1 = s k + 1 + β k q k . 第二式左乘 C − ∗ 后成为 p k + 1 = M − 1 r k + 1 + β k p k ,完整解释了 apply 和方向更新的来源。C 只用于证明,执行算法不必构造或调用它。
同一映射也确定了正确的内积几何 公理库 内积空间 Inner product space 带正定对称双线性形式或正定 Hermitian 半双线性形式的向量空间。 。普通 CG 的 s i ∗ s j = 0 与 q i ∗ B q j = 0 分别给出
r i ∗ M − 1 r j = 0 , p i ∗ A p j = 0 ( i ≠ j ) . 原残差在 M − 1 内积中正交,搜索方向仍在原来的 A 内积中共轭。通常并没有 r i ∗ r j = 0 ;把普通 CG 的 Euclidean 残差正交性直接搬到原变量,会遗漏整个坐标变化。
令 T = M − 1 A 。它通常不是 Euclidean Hermitian 矩阵,因为 T ∗ = A M − 1 ;但它在 M 内积中满足
⟨ u , T v ⟩ M = u ∗ A v = ⟨ T u , v ⟩ M , ⟨ u , T u ⟩ M = u ∗ A u > 0. 此外 C ∗ T C − ∗ = B ,故 T 与正定矩阵 B 相似,具有相同的正特征值。PCG 的几何来自这种带权自伴性,不能由“左乘以后看起来更好解”替代。
例子与边界
取一个只需二阶块求解和一次缩放的预条件器:
C = ( 1 0 0 1 2 0 0 0 4 ) , M = C C T = ( 1 1 0 1 5 0 0 0 16 ) . 令
B = ( 2 1 0 1 2 1 0 1 2 ) , A = C B C T = ( 2 4 0 4 14 8 0 8 32 ) , b = ( 1 1 0 ) , x 0 = 0. B 的特征值为 2 − 2 , 2 , 2 + 2 ,所以 A , M 均正定。这里 M ≠ A ,且 ( A M ) 12 = 22 、( M A ) 12 = 18 ,两者并不交换;因而此例保留了左预条件算子不对称的真实情形。对向量 r = ( r ( 1 ) , r ( 2 ) , r ( 3 ) ) T ,apply 明确为
M − 1 r = ( ( 5 r ( 1 ) − r ( 2 ) ) / 4 ( − r ( 1 ) + r ( 2 ) ) / 4 r ( 3 ) / 16 ) . 以下均为精确分数,表中三元组表示列向量。先列出每个时刻的解、残差和预条件残差:
k
x k
r k
z k
0
( 0 , 0 , 0 )
( 1 , 1 , 0 )
( 1 , 0 , 0 )
1
( 1 / 2 , 0 , 0 )
( 0 , − 1 , 0 )
( 1 / 4 , − 1 / 4 , 0 )
2
( 5 / 6 , − 1 / 6 , 0 )
( 0 , 0 , 4 / 3 )
( 0 , 0 , 1 / 12 )
3
( 1 , − 1 / 4 , 1 / 16 )
( 0 , 0 , 0 )
无须计算
产生这些状态的三个搜索方向及其矩阵乘积为
p 0 = ( 1 , 0 , 0 ) T , A p 0 = ( 2 , 4 , 0 ) T , p 1 = ( 1 / 2 , − 1 / 4 , 0 ) T , A p 1 = ( 0 , − 3 / 2 , − 2 ) T , p 2 = ( 2 / 9 , − 1 / 9 , 1 / 12 ) T , A p 2 = ( 0 , 0 , 16 / 9 ) T . 标量更新与误差检查单独列出。最后一行的 ρ 3 = 0 是由零残差得出的值,不需要再执行 apply。
k
ρ k
α k
β k
‖ r k ‖ 2 2
‖ x ∗ − x k ‖ A 2
0
1
1 / 2
1 / 4
2
3 / 4
1
1 / 4
2 / 3
4 / 9
1
1 / 4
2
1 / 9
3 / 4
停止
16 / 9
1 / 12
3
0
—
—
0
0
例如第二步的分母为 p 1 T A p 1 = 3 / 8 ,所以 α 1 = ( 1 / 4 ) / ( 3 / 8 ) = 2 / 3 ;它把 ( 0 , − 1 , 0 ) T 更新为 ( 0 , 0 , 4 / 3 ) T 。第三步的分母为 4 / 27 ,所以 α 2 = ( 1 / 9 ) / ( 4 / 27 ) = 3 / 4 ,直接得到 A x 3 = b 。这完整执行了三步递推,而非仅展示初始与最终状态。
这个历史同时区分了三种量。首先 r 0 T r 1 = − 1 ,但 r 0 T M − 1 r 1 = 0 ,正交性确实依赖权重。其次,第二步的原残差范数由 1 增到 4 / 3 ,能量误差平方却由 1 / 4 降到 1 / 12 。最后,ρ k 与能量误差平方的数值也不相同;本例中二者恰好都下降,不能据此把预条件残差当作能量误差。
谱改善也可以精确检查:
κ 2 ( B ) = 2 + 2 2 − 2 = 3 + 2 2 < 6. 原矩阵的坐标 Rayleigh 商给出 λ max ( A ) ≥ 32 、λ min ( A ) ≤ 2 ,从而 κ 2 ( A ) ≥ 16 。廉价块预条件确实缩小了谱跨度;它没有把 B 变成单位阵,也没有消除所有耦合。
固定线性正定的 apply 是上述推导的边界。若每一步改变内层求解容差,或使用依赖输入的非线性近似过程,就不再有同一个 C 和同一个 B ;普通 ρ k + 1 / ρ k 公式因而失去这里给出的依据。即使每次返回的向量都看似合理,也需要适用于该变化方式的算法分析。
推论与应用
令 e k = x ∗ − x k 。坐标映射保持能量误差,因为
( C ∗ e k ) ∗ B ( C ∗ e k ) = e k ∗ A e k . 又因 C − ∗ B = T C − ∗ 且 C − ∗ s 0 = z 0 ,有
C − ∗ K k ( B , s 0 ) = K k ( T , z 0 ) . 将普通 CG 的最小化定理沿这两个等式映回,就得到 PCG 在原变量中的精确结论:
x k ∈ x 0 + K k ( M − 1 A , z 0 ) , ‖ x ∗ − x k ‖ A = min x ∈ x 0 + K k ( M − 1 A , z 0 ) ‖ x ∗ − x ‖ A . 因此能量误差随嵌套的Krylov 子空间 公理库 Krylov 子空间 Krylov subspace · Arnoldi iteration · Lanczos iteration 从初始向量反复施加矩阵,生成只依赖矩阵—向量乘法的递增子空间及其 Arnoldi、Lanczos 正交基。 单调不增。若 s 0 关于 B 的最小多项式次数为 d ,精确算术中至多 d ≤ n 步终止;浮点运算中不保证第 n 步得到零残差。
普通 CG 的谱界同样映回为
‖ e k ‖ A ≤ 2 ( κ − 1 κ + 1 ) k ‖ e 0 ‖ A , κ = κ 2 ( B ) = λ max ( B ) λ min ( B ) . 这里使用对称正定 B 的条件数 公理库 线性方程组的条件数与扰动 Conditioning of linear systems · Matrix condition number 把一般问题条件性具体化为可逆线性系统的右端、系数矩阵与联合扰动界。 。虽然 T = M − 1 A 与 B 的特征值相同,非酉相似一般不保持奇异值,因此不能把上式中的 κ 换成 Euclidean 奇异值条件数 κ 2 ( T ) 。该界只使用谱端点;特征值聚集与初始误差的分量还会影响实际迭代数。
对大型稀疏 SPD 系统,标准 PCG 每个继续迭代的周期需要一次 A 乘法、一次 apply、两个主要内积和 O ( n ) 向量更新;除算子及预条件结构外,只需常数个长向量。若完成 k ≥ 1 次更新并在最后一次更新后直接停止,apply 的调用数是初始一次加中间 k − 1 次,共 k 次。初始残差计算和另行重算真实残差所需的 A 乘法也必须记录;已知 x 0 = 0 时可省略初始的零向量乘法。
把设置、应用和检查分开,总成本可写为
C setup + N A C A + N M C M + O ( k n ) + C checks , 其中 N A 包含初始及验证残差的矩阵乘法,N M 是实际 apply 次数,C checks 计入额外范数与通信归约。预条件器的因子、填充或层级存储也属于成本。比较方法时,应同时报告设置时间与总求解时间;多个右端能摊销设置,而一次短求解未必能回本。
有限精度会使递推残差偏离 b − A x k ,并使精确正交性逐渐丢失。可以用递推残差筛选候选停止点,再以重算的真实残差决定成功;长迭代也可定期检查漂移。如果把重算残差替换进迭代并继续,应重新设置 z = M − 1 r 、p = z 、ρ = r ∗ z ,作为一次重启,而不是保留旧方向并继续声称完整共轭性。
非零残差下的 ρ ≤ 0 、p ∗ A p ≤ 0 或非有限输出,与精确固定 SPD 合约不相容,实际实现应停止并报告数值故障或接口失败;达到预算上限则返回未收敛状态及真实残差。初始零残差是正常成功路径,须在任何 apply 或步长除法之前处理。这样,理论中的有限步终止、用户要求的原系统精度和机器上的实际退出原因才有各自明确的含义。
参考资料