带收缩的MINRES种子投影方法_朱景福.pdf

第 33 卷第 1 期广东石油化工学院学报Vol 33No 12023 年 2 月Journal of Guangdong University of Petrochemical TechnologyFebruary 2023带收缩的 MINES 种子投影方法朱景福，李欣，林靖杰(广东石油化工学院 理学院，广东 茂名 525000)摘要:研究实对称线性方程组 AX=B 的数值解法。种子投影方法是求解线性方程组的一种常见方法，但是当系数矩阵为对称矩阵时，种子投影方法的有效性往往会降低。考虑把收缩技术和种子投影方法相结合，提出求解实对称线性方程组的收缩 MINES 种子投影方法，并分析算法的残量。数值实验结果表明新方法是有效的。关键词:实对称线性方程组;收缩技术;MINES 种子投影方法;Krylov 子空间中图分类号:O242 2文献标识码:A文章编号:2095 2562(2023)01 0076 04考虑实对称线性方程组 AX=B，式中:An n是非奇异实对称矩阵，X，Bn p。对 AX=B 数值解法的理论研究、算法设计、软件研发及工程应用等问题一直以来是计算科学研究的热点问题。求解 AX=B 的数值方法中常见的有以下两类:块方法和种子投影方法。但是，块方法1 3 通常需要大量的计算量和存储量，有时比逐次求解 AX=B 的 p 个方程组速度更慢。块方法在求解多右端对称线性方程组时的有效性会降低。种子投影法首先是选择一个方程组作为种子方程组，对该种子方程组使用某种 Krylov 子空间方法进行求解，同时产生一个 Krylov 子空间 K，然后利用正交投影，把非种子方程组的残量投影到 K 上并取极值，以获得更好的近似。当种子方程组收敛到精度要求时，在剩下的还没有收敛的非种子方程组中，再重新选择一个新的种子方程组，继续循环，直到得到所有方程组的解。朱文跃和顾桂定提出了 GMES 种子投影方法4，5。刘皞等提出求解多右端线性方程组的块种子投影方法6 等。在实践中数学工作者们不断地尝试改进已有算法。自 1995 年 Morgan 提出了带特征向量的重新开始GMES 方法7，8 以来，收缩技术越来越得到人们的重视，Chapman 和 Saad 进一步发展了 Morgan 提出的收缩技术，提出了收缩增广子空间技术9，2000 年 Saad 等提出了收缩共轭梯度法10。2009 年李欣等提出了求解非对称线性方程组的循环收缩 QM 方法11 等。考虑在求解实对称线性方程组 AX=B 时，采用 MINES 种子投影方法6，同时在 Lanczos 过程利用增广子空间技术，向 Krylov 子空间加入少量模较小的特征值所对应的特征向量进行收缩，本文拟提出两种格式的求解方程组 AX=B 的新方法;最后对新方法做数值试验，验证新方法的有效性。1预备知识考虑到求解方程组 AX=B 的新方法需要用到 Lanczos 过程6，下面先回顾求解对称线性方程组 Ax=b 的 Lanczos 过程，其中 A 是 n n 实对称矩阵，x，bn。设 x0是方程组 Ax=b 的初始估计，并记 r0=b Ax0，=r0，v1=r0/。应用 Lanczos 过程可产生 Krylov 子空间 km(A，r0)=span(r0，Ar0，Am 1r0)的一组标准正交基 v1，v2，vm，记 Vm=(v1，v2，vm)，则 VTmAVm是对称三对角阵，记为收稿日期:2022 09 30;修回日期:2022 10 20基金项目:广东石油化工学院人才引进启动项目(2018rc44，2018rc45)作者简介:朱景福(1970)，男，黑龙江克山人，博士，教授，主要从事计算机应用研究。通信作者:李欣(1968)，女，黑龙江讷河人，硕士，教授，主要从事高等数学研究。Hm=12002230330m00mm则 AVm=VmHm+m+1vm+1eTm，其中 em=(0，0，1)Tm，记?Hm=Hmm+1eTm，em=(0，0，1)Tm，则AVm=VmHm+m+1vm+1eTm=Vm+1?Hm2带收缩的 MINES 种子投影算法种子的选取有多种方法。一个理想的种子方程组应在其收敛后能使每一个非种子方程组得到一个较好的初始近似。这样，当某个非种子方程组作为种子方程组，在下一次迭代时可能会有更好的收敛速度。为此，选取种子方程组的一个方法是，把方程组残量的线性组合作为种子系统的残量。另一个改进的方式是把当前残量模最大的方程组作为种子系统，见参考文献 4、6。本文考虑选取方程组 AX=B 的初始估计后，计算残量?r(k)，k=1，2，p，依次选择第 k 个方程组为种子方程组。2 1收缩 MINES 种子投影方法下面利用收缩技术对 Lanczos 过程产生的子空间添加矩阵 A 的 k 个近似特征向量以扩大子空间维数。具体做法是，执行 m 步 Lanczos，得到向量组 Vm=v1，v2，vm 和矩阵 Hm，再求 Hm的特征值和特征向量，选取 k 个模最小特征值所对应的特征向量 g1，g2，gk，求对应的 itz 向量 yi=Vmgi，(i=1，2，k)，作为 A 的近似特征向量。把 Y=Vm g1，g2，gk=y1，y2，yk 添加到对称 Lanczos 过程中，建立Krylov 子空间 K=span rs，Ars，Am 1rs，y1，y2，yk。记 W=v1，v2，vm，y1，y2，yk，则AW=AVm，A(y1，y2，yk)=Vm+1?Hm，A(y1，y2，yk)记 H 是(m+k+1)(m+k)上 Hessenberg 阵，其前 m 阶顺序主子矩阵是 Lanczos 过程产生的对称三对角阵，其后 m+1 行后 k 列是对应 y1，y2，yk的特征值，则AW=Vm+1，y1，y2，yk H=v1，v2，vm，vm+1，y1，y2，yk H记 Q=(v1，v2，vm，vm+1，y1，y2，yk)，则有 AW=Q H通过前面的分析，首先提出一种简单格式的收缩 MINES 种子投影方法(DMSeed1)。算法 2 1收缩 MINES 种子投影方法(DMSeed1)(1)给出线性方程组 AX=B 的初始估计 X=(x1，x2，xp)，计算残量 =B AX=(r1，r2，rp);(2)l=1 p(求解所有方程组的解)(3)令 s=l(即选取当前第 l 个方程组作为种子方程组)，计算 rs=b(s)Ax(s)，=r(s)2，q1=r(s)/;(4)执行 m 步 Lanczos(A，q1)过程，得到标准正交向量组 Vm=v1，v2，vm 和对称三对角矩阵 Hm;(5)建立增广子空间。求 Hm的特征值和特征向量，构造 Y=y1，y2，yk，构造 W=(v1，v2，vm，y1，y2，yk)，Q=(v1，v2，vm，vm+1，q1，q2，ql)，以及 H(见前面分析);(6)求解最小二乘问题 d(s)=argmmdme1 Hd2并计算?x(s)=xs+Wd(s)，?r(s)=b(s)A?x(s);(7)对 j=l+1，p(即对所有非种子方程组)，计算 d(j)=armindmQTr(j)Hd2，?x(j)=x(j)+Wd(j)，?r(j)=b(j)A?x(j)(或?r(j)=r(j)Q Hd(j);若?r(j)2，删除该已收敛的方程组;(8)若?r(j)2，删除该种子方程组，转步骤(2);77第 1 期朱景福等:带收缩的 MINES 种子投影方法(9)End l。2 2收缩 MINES 种子投影方法 2建立增广子空间时也可通过 W 张成的子空间来寻找近似特征向量，根据前面推导，求解如下特征值WTA2Wg=WTAWg取其 k 个模较小的特征值所对应的特征向量 g1，g2，gk，则 yi=Wgi(i=1，2，k)为所求的近似特征向量。这种选择较小特征值所对应的特征向量的办法实际上是参考文献 8 中的方法的简化形式。提出求解线性方程组 AX=B 的另外一个格式的收缩 MINES 种子投影方法(DMSeed2)。算法 2 2收缩 MINES 种子投影方法(DMSeed2)(1)给出线性方程组 AX=B 的初始估计 X=(x1，x2，xp)及 Y=y1，y2，yk，计算残量 =B AX=(r1，r2，rp);(2)令 l=1 p(求解所有方程组的解)(3)令 s=l(即取第 l 个方程组作为种子方程组)，计算 r(s)=b(s)Ax(s)，=r(s)2，q1=r(s)/;(4)执行 m 步 Lanczos(A，q1)过程，得到标准正交向量组 Vm=v1，v2，vm，vm+1;(5)建立增广子空间。构造 W=(v1，v2，vm，y1，y2，yk)，Q=(v1，v2，vm，vm+1，y1，y2，yl)，满足(3)式，以及 H(见前面分析);(6)求解最小二乘问题 d(s)=argmindme1 Hd2并计算?x(s)=x(s)+Wd(s)，?r(s)=b(s)A?x(s);(7)对 j=l+1，p(即对所有非种子方程组)计算 d(j)=argmindmQTr(j)Hd2，?x(j)=x(j)+Wd(j)，?r(j)=b(j)A?x(j)(?or r(j)=r(j)Q Hd(j);若r(j)(删除该已收敛的方程组);(8)若?r(s)2，删除该种子方程组，转步骤(2);(9)求解广义特征值问题解特征值问题 WTA2Wg=WTAWg，取其 k 个绝对值最小特征值所对应的特征向量 g1，g2，gk，计算 yi=Wgi(i=1，2，k)，令 Y=y2，y2，yk，p=p 1 转步骤(1);(10)End l。3数值实验数值实验均在主频 1 8GHz 的酷睿 i7CPU、内存 8GB、Win10 企业版操作系统的笔记本电脑上采用Matlab 编程实现。用两个新方法(算法 2 1(DMSeed1)、算法 2 2(DMSeed2)与极小向后扰动块 Lanczos方法(BMINBACK)、Minres 种子投影方法做对比实验求解方程组 AX=B。例 1:已知 n=100 阶实对称矩阵A=1 615600015221 600061 6220000000221 660001 6221 600061622矩阵 B 中元素为 B(i，j)=cos(5 cos(i 2 (j 1)/128)，算法中 Lanczos 执行步数为 8，增广子空间增加列数为 4，初始估计为全 0 矩阵，精度要求为 106。收缩种子投影方法(算法 DMSeed1 和算法 DMSeed2)与极小向后扰动块方法(BMINBACK)，Minres 种子投影方法做对比数值实验，数值结果见表 1。87广东石油化工学院学报2023 年表 1四种算法数值对比实验B 列数DMSeed1时间/s残量范数DMSeed2时间/s残量范数BMINBACK时间/s残量范数Minres 种子投影时间/s残量范数4001018 8379E 08000658 8347E 080065586221E 060 00541 6373E 078002441 3113E 07001581 3410E 070061412605E 050 01274 2078E 0712004291 5113E 07002061 6376E 070030111474E 050 01854 2274E 0716006891 6465E 07003092 4762E 070083728896E 050 02581 0319E 0620009211 8801E 07005212 6624E 070091619034E 050 03841 3598E 06注:1E 08 表示 1 108图 1 反映了比较算法 DMSeed1、DMSeed2 和Minres 算法的残量范数随着右端项列数增加时不断变化的情况。从表 1 可见，算法 DMSeed1、DMSeed2在残量性能上要优于 BMINBACK 算法和 Minres 算法。从时间性能上看，算法 DMSeed1、DMSeed2 的时间性能居于 BMINBACK 方法和 Minres 种子投影方法之