求解非厄米特正定线性系统的...ichardson迭代算法_肖小勇.pdf

第 28 卷第 1 期2023 年 2 月新余学院学报JOUNALOFXINYUUNIVESITYVol 28，NO 1Feb 2023求解非厄米特正定线性系统的预处理 ichardson 迭代算法 肖小勇(新余学院数学与计算机学院，江西新余338004)摘要:为了求解非厄米特正定线性系统，引入了一种新的预处理 ichardson 迭代算法(P 迭代算法)。每一次迭代，P 迭代算法只需求解一个带厄米特正定系数矩阵的线性系统。在适当的条件下，分析了 P 迭代矩阵的谱半径，并讨论使上述谱半径取最小值时的最优参数。数值结果表明，不管是否采用实验最优参数，P 迭代算法都是有效的。关键词:非厄米特矩阵;正定;HSS 迭代;谱半径;收敛性分析中图分类号:O241 6;O242 2文献标识码:A文章编号:2095 3054(2023)01 0021 08收稿日期:2022 07 03基金项目:国家自然科学基金项目“求解几类大规模线性系统的随机算法及其应用研究”(12061048)。作者简介:肖小勇(1980 )，男，湖南茶陵人，教授，博士，主要从事数值线性代数、概率极限理论研究。本文考虑如下形式的大型稀疏线性系统的迭代解法Ax=b(1)其中 A nn是非厄米特正定的(即 A 的厄米特部分是正定的)。在本文中，右端的常数项 bn是已知的，x n为未知向量。科学计算中的许多问题及应用工程领域经常会遇见这样的线性方程组，如分子散射、晶格量子色动力学、量子化学、漫射光学层析成像、基于快速傅立叶变换的时变偏微分方程求解等。众所周知，系数矩阵 A 可以改写成厄米特和斜厄米特分裂(HSS)形式:A=H+S，其中 H=A+A*2，S=A A*2。这里 A*表示矩阵 A 的共轭转置。基于这种分裂，Bai 等人首先提出了 HSS 迭代算法1。从那时起，HSS 方法立即引起了极大关注，研究人员提出了许多变形算法。如潘春平等提出的外推 HSS 迭代算法2、温瑞萍等提出的衍生多分裂方法3、张建华提出的若干迭代算法及其预处理4、张迎春等提出的若干预处理迭代算法5。最近，Li 和 Wu 提出了一种单步 HSS(SHSS)迭代法6(I+H)x(k+1)=(I S)x(k)+b，其中 0 是一个给定的常数，I 是单位矩阵。数值结果表明，在一定条件下 SHSS 迭代算法优于HSS 迭代算法。在 SHSS 算法的基础上，Zeng 和 Ma提出了一种参数化的 SHSS 迭代算法来求解复对称的线性方程组7。Wu 等建立了如下的非交替 PHSS(NPHSS)迭代算法8(P+H)x(k+1)=(P S)x(k)+b，其中 P nn是给定的厄米特正定矩阵。其理论和数值结果表明，当系数矩阵的厄米特部分占主导的情况下，NPHSS 迭代算法比 PHSS 迭代算法更有效。22新余学院学报2023 年对于求解大型稀疏的非厄米特正定线性系统式(1)，本文构造了一种新的单步迭代算法，它可以看作一种预处理 ichardson 迭代算法(P 迭代算法)，该算法已被用于增广线性系统的求解9。本文分析了 P 迭代算法的收敛性质，包括收敛条件、迭代矩阵的谱半径和拟最优参数，并给出了新的预处理矩阵的性质。此外，给出了几个数值实验来说明 P 迭代算法的有效性。1P 迭代算法对于任何厄米特正定矩阵 P nn和常数 0，受到 Px=Px (Ax b)=(P A)x+b 的启发，可以构造如下的 P 迭代算法:给定初始值 x(0)n，对于 k=0，1，2，3，直到迭代序列 x(k)收敛，按如下方案计算 x(k+1):Px(k+1)=(P A)x(k)+b(2)这里 0 是一个给定的正常数，P nn是一个指定的厄米特正定矩阵。由于式(2)左边的系数矩阵 P nn是厄米特正定矩阵，因此可以采用Cholesky 分解或共轭梯度法从式(2)求出 x(k+1)。另一方面，式(2)等价于x(k+1)=(I P1A)x(k)+P1b，它可以看作一种预处理 ichardson 迭代，因为它可以改写成 x(k+1)=x(k)P1(Ax(k)b)。此外，它也可以看作是不动点迭代的外推加速，对应的矩阵分裂为 A=F G，这里 F=P。接下来将式(2)改写为x(k+1)=M(，P)x(k)+N(，P)b，其中M(，P)=I P1A=I P1(H+S)=P12 I P12HP12 P12SP12 P12(3)并且 N(，P)=P1。令F(，P)=1P 和 G(，P)=1P A，可以得到A=F(，P)G(，P)，M(，P)=F(，P)1G(，P)。因此，分裂矩阵 F(，P)=1P 可以作为 A nn的预处理矩阵，称为 P 的预处理子，记为 P=1P。在 Krylov 子空间算法中应用预处理子 P的每一步迭代，都需要求解如下形式的广义残差方程序列:Pz=1Pz=r(4)其中 r，z n分别为当前残差向量和广义残差向量。由于 P 是厄米特正定矩阵，因此，上述的线性系统(4)可以用 Cholesky 分解或共轭梯度法求解。2P 迭代算法的收敛性分析为了分析 P 迭代算法的收敛性，特引入如下定义:Pmin=minjsp(HP)j，Pmax=maxjsp(HP)j，Pmax=maxjsp(iSP)|j|(5)其中 i=1 是虚数单位，XP=P12XP12，X=H，S，并且 sp(Y)表示矩阵 Y 的谱集。定理 1设 A nn是一个正定矩阵，H=A+A*2和 S=A A*2分别是它的厄米特和斜厄米特部分，设 P 是厄米特正定矩阵，对于 P 迭代算法，其谱半径满足(M(，P)max(1 Pmin)2，(1 Pmax)2+2(Pmax)2=:P()(6)如果0 min2Pmin(Pmin)2+(Pmax)2，2Pmax(Pmax)2+(Pmax)2(7)那么 P()1，即 P 迭代算法收敛。证明:既然 HP是厄米特正定的，而 SP是斜厄米特的，那么根据(M(，P)=(I P12HP12 P12SP12)，由(3)可以很容易得到(6)。经过简单讨论和直第 1 期肖小勇:求解非厄米特正定线性系统的预处理 ichardson 迭代算法23接计算，可以得到(P()2=(1 )2+2(Pmax)2，0 2Pmin+Pmax(1 )2+2(Pmax)2，2Pmin+Pmax(8)因此，从(8)式可以得出如下结论:如果 PminPmax(Pmax)2，并且0 2Pmin(Pmin)2+(Pmax)2;或者如果PminPmax(Pmax)2，并 且02Pmax(Pmax)2+(Pmax)2，那么 P()1 成立，而这些条件可以简化为不等式(7)。证毕。推论 1 对于 P 迭代算法，如果(Pmin)2+2(Pmax)2 PminPmax(9)则最小化谱半径(M(，P)的上界 P()的最优参数 P*由下式给出P*=Pmin(Pmin)2+(Pmax)2(10)相应地P(P*)=Pmax(Pmin)2+(Pmax)2(11)证明:通过对(8)式的分析，谱半径(M(，P)达到最小值P(P*)=Pmax(Pmin)2+(Pmax)2，P*=Pmin(Pmin)2+(Pmax)22Pmin+PmaxPmax(Pmax)2+(Pmax)2，=Pmax(Pmax)2+(Pmax)22Pmin+Pmax然而不等式Pmax(Pmax)2+(Pmax)22Pmin+Pmax任何时候都不成立，因为 Pmax(Pmin+Pmax)2(Pmax)2。由此，根据条件(9)可以得到式(10)和式(11)。定理 2设 A nn是一个正定矩阵，H=A+A*2和 S=A A*2分别是它的厄米特和斜厄米特部分。设 P 是厄米特正定矩阵，那么预处理矩阵1PA 是正稳定，即 1PA 的每个特征值的实部为正数。证明:因为1PA=P1A=P12(P12AP12)P12，所以1PA 与 P12AP12=HP+SP相似。这里 HP=P12HP12是厄米特正定矩阵，SP=P12SP12是 斜 厄 米 特 矩 阵，分 别 是 矩 阵 AP=P12AP12的厄米特部分和斜厄米特部分。由于 HP=P12HP12是正定的，因此，P12AP12是一个正定矩阵，这意味着矩阵 1PA 也是正定的，故为正稳定的。证毕。注 1:根据推论 1，若(9)成立，则理论上使得P(P*)最小的最优厄米特正定矩阵 P*可由下式给出P*=argminPPmax(Pmin)2+(Pmax)2=argminPPmaxPmin。因为 f(x):=x1+x2是严格单调递增函数，通常在实践中，可以用厄米特正定矩阵 H 代替矩阵 P。实际上，令 P=H，P 迭代算法就变成Hx(k+1)=(1 )P S x(k)+b(12)接下来将(12)改写为x(k+1)=M()x(k)+N()b，其中M()=(1 )I H1S，并且 N()=H1。令F()=1H，G()=1 H S，那么A=F()G()，M()=F()1G()。因此，分裂矩阵 F()=1H 可以作为 Ann的预处理子，记为 PH=1H。引理 1设 H 为厄米特正定矩阵，S 为斜厄米特矩阵。那么 H1S 的每个特征值要么为 0，要么为纯虚数。特别地，S 的每个特征值要么是 0，要么是纯24新余学院学报2023 年虚数。证明:设 为 H1S 对应于特征向量 0 的特征值，即 S=H，那么*S=*H(13)并且*S*=*H，从而*S=*H(14)根据式(13)、(14)，可得(+*)*H=0，这意味着 +*=0 因为 H 是厄米特正定的。证毕。定理 3设 A nn是一个正定矩阵，H=A+A*2和 S=A A*2分别是它的厄米特和斜厄米特部分，那么预处理矩阵 1PHA 的每个特征值的实部为正数 。证明:因为1PHA=H1A=H1(H+S)=(I+H1S)，然后利用引理 1 就完成了证明。证毕。由于谱半径(M()=(1 )I H1S，运用引理 1，可以得到如下结果。定理 4对于采用 P=H 的 P 迭代算法，谱半径(M()=(1 )2+22(H1S)=:()(15)如果0 21+2(H1S)，那么()1，即P 迭代算法收敛。并且，最小化谱半径()的最优参数*可由下式给出*=11+2(H1S)(0，1)，相应地(*)=(H1S)1+2(H1S)(16)证明:这里只证明(16)式，其余的结论可以从定理 1 和推论 1 直接得到。通过求导计算，可以得到 2()=2 1+2(H1S)2。因此，很 明 显，谱 半 径(M()在*=11+2(H1S)取最小值。将*代入()，就可以得到()的最小值。证毕。在定理 4 的条件下，若(H1S)1，则有(M(*)=(*)12。根据 SHSS 迭代法中的推论 2 16，其谱半径(I+H)1(I S)上限 的最小值为*=max2min+2max=maxmin1+max()min2(17)其中 min为矩阵H的最小特征值，max为矩阵S 的最大奇异值。由(16)和(17)式，利用(H1S)|H1S|2|H1|S|2=maxmin，可知(*)*。因此，分别采用最优参数*时，P=H 的 P 迭代算法比 SHSS 迭代算法更有效。3数值实验及结果本节将采用 P=H 的数值实验来测试 P 迭代算法的有效性，并与 SHSS6 和 NPHSS8 迭代算法比较，考察迭代步数(记为 IT)和 CPU 时间(以秒为单位，记为 CPU)。这里 NPHSS 迭代算法中使用的预处理子 P 与 WU 等8 的 论 文 相 同，即 P=diag(a11，a22，ann)，其中 a11，a22，ann为矩阵 A 主对角线元素。所有 实 验 均 在 MATLAB(版 本 8 1 0 604，2013a)中进行。电脑的运行环境为英特尔()核(TM)CPU 1 8 2GHz 和 4 00GB 的内存。在计算中，所有的迭代算法都是从 0 向量开始，停止标准为|b Ax(k)|2|b|2 106，其中 x(k)是当前的近似值。接下来，求解非厄米特正定线性系统(1)的数值算法均采用预处理GMES 算法，F()或 F(，P)是预处理子。