1、第 28 卷第 1 期2023 年 2 月新余学院学报JOUNALOFXINYUUNIVESITYVol 28，NO 1Feb 2023求解非厄米特正定线性系统的预处理 ichardson 迭代算法 肖小勇(新余学院数学与计算机学院，江西新余338004)摘要:为了求解非厄米特正定线性系统，引入了一种新的预处理 ichardson 迭代算法(P 迭代算法)。每一次迭代，P 迭代算法只需求解一个带厄米特正定系数矩阵的线性系统。在适当的条件下，分析了 P 迭代矩阵的谱半径，并讨论使上述谱半径取最小值时的最优参数。数值结果表明，不管是否采用实验最优参数，P 迭代算法都是有效的。关键词:非厄米特矩阵;
2、正定;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为未知向量。科学计算中的许多问题及应用工程领域经常会遇
3、见这样的线性方程组，如分子散射、晶格量子色动力学、量子化学、漫射光学层析成像、基于快速傅立叶变换的时变偏微分方程求解等。众所周知，系数矩阵 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(S
4、HSS)迭代法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 年对于求解大型稀疏的非厄米特正
5、定线性系统式(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)
6、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)=P
7、12 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 迭代算法的收敛性分析为了
8、分析 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(
9、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，并 且02P
10、max(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(Pm
11、ax)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=P12SP1
12、2是 斜 厄 米 特 矩 阵，分 别 是 矩 阵 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
13、(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)根据式(1
14、3)、(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 迭代算法收敛。并且，最小化谱半径()的
15、最优参数*可由下式给出*=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的最
16、小特征值，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)是预处理子。
