线性代数方程组中的预处理共轭梯度法

线性代数方程组迭代法预处理共轭梯度法孟庆彬November26,20201引入•线性代数方程组的解法•直接法：高斯消去法，分解法•迭代法：•古典迭代法：Jacobi，Gauss-Seidel，SOR，SSOR•现代迭代法：（投影方法，子空间法）•正交化的误差投影型Krylov：FOM，IOM，DIOM•对称情形误差投影型Krylov：Lanczos，CG，PCG•正交化的残量投影型Krylov：GMRES，GCR•双正交化投影型Krylov方法：BiCG，CGS…..1引入•共轭梯度法（1950S）•对于固定的k，在Rn空间中取k个线性无关的量p1,p2,...,p𝑘使𝑥𝑘满足：𝜑(𝑥𝑘)=min𝜑(𝑥)𝑠𝑝𝑎𝑛(𝑝1,𝑝2,...,𝑝𝑘)当k=n时，𝑥𝑛一定是方程Ax=b的精确解，CG的出发点•数值稳定性不好•舍入误差必然存在•条件数较大时收敛慢•预处理共轭梯度法（1970S）•预处理系数矩阵cond(𝐴)=||𝐴||•||𝐴−1||2知识回顾•CG算法残差向量搜索步长搜索方向3算法原理•PCG算法•矩阵A的对称正定的•矩阵M相似于A•原始残差𝑟𝑗=𝑏−𝐴𝑥𝑗•预处理残差𝑧𝑗=𝑀−1𝑟𝑗•特殊情况，M为单位矩阵•问题是预处理矩阵M应该如何选取呢？•M是对称正定•M与A具有相似性•方程Mz𝑘=𝑟𝑘•𝐴的特征值比较集中，即cond(𝐴)≈1或尽可能小3算法原理•PCG算法•取M的LU分解M=LLT•L−1AL−𝑇u=L−1𝒃•x=L−1u•𝑝𝑗=𝐿𝑇𝑝𝑗•𝑢𝑗=𝐿𝑇𝑥𝑗•𝑟𝑗=𝐿𝑇𝑧𝑗=𝐿−1𝑟𝑗•𝐴=𝐿−1𝐴𝐿−𝑇•𝑟𝑗,𝑧𝑗)=(𝑟𝑗,𝐿−𝑇𝐿−1𝑟𝑗)=(𝐿−1𝑟𝑗,𝐿−1𝑟𝑗)=(𝑟𝑗,𝑟𝑗•𝐴𝑝𝑗,𝑝𝑗)=(𝐴𝐿−𝑇𝑝𝑗,𝐿−𝑇𝑝𝑗)=(𝐿−1𝐴𝐿−𝑇𝑝𝑗,𝑝𝑗)=(𝐴𝑝𝑗,𝑝𝑗4预处理方法•预处理方法•取预优矩阵（预处理矩阵）为A的一个小带宽部分（如三对角或对角线部分）•矩阵分裂，尤其是线性稳定迭代中的矩阵A的分裂构造预处理矩阵•通过A的各种近似分解得到预处理矩阵（如不完全分解）•通过矩阵A的多项式构造预处理矩阵•子结构，区域分裂，EBE预处理途径等等4预处理方法•预处理方法•对角预处理法：若A是严格对角占优的矩阵，则可以选择𝑀=𝑑𝑖𝑎𝑔(𝑎11,𝑎22,...,𝑎𝑛𝑛)为预处理矩阵；若A是分块对角阵，则可以选择𝑀=(𝐴11,𝐴22,...,𝐴𝑛𝑛)为预处理矩阵但是这样简单的选择，往往不能带来很好的加速收敛效果。4预处理方法•预处理方法•分裂方法：•将A分裂为：A＝P－Q；•通过矩阵P，Q来构造预处理矩阵M，一般取线性稳定迭代法中相应分裂•如：Jacobi分裂（P＝diag(A)）•Gauss-Seidel分裂（P＝diag(A)–L，L是A的严格下三角）•SSOR分裂：𝐴=𝐷−𝐶−𝐿𝐶𝐿𝑇，𝐶𝐿为严格下三角4预处理方法•预处理方法•不完全分解预处理方法：主要有不完全LU分解、不完全Cholesky分解以及相应的块不完全分解等•设𝐿𝐿𝑇是A的一个近似分解𝐴=𝐿𝐿𝑇+𝑅，L是下三角矩阵，R称为剩余矩阵，M=𝐿𝐿𝑇为预处理矩阵。•此方法与CG迭代方法相结合，就形成了不完全Cholesky分解预处理共轭梯度法(ICCG)，由于参数𝜔可以调整，以达到最好的收敛效果，所以它是目前人们使用最多的方法之一，而且从这个方法上改进出的方法也很多。5高效实现•PCG计算量•主要工作量是𝑀−1和向量的乘法，以及一次矩阵A和向量的乘法•A和M是相似的•第一种高效实现不完全分解•M=A-R•𝑀−1𝐴𝑣=𝑀−1(𝑀+𝑅)𝑣=𝑣+𝑀−1𝑅𝑣•R中的非零元素远远少于A。R需要显式存储•第二种高效实现分裂•A是对称的•𝐴=𝐷0−𝐸−E𝑇，-E严格下三角，D0是对角线矩阵•M=(𝐷−𝐸)𝐷−1(𝐷−E𝑇，D不完全等于D0•SSOR-PCG谢谢！