求解正定线性方程组的PPS迭代方法

刘仲云，梅聪辉

(长沙理工大学 数学与统计学院，湖南 长沙，410114)

文中考虑迭代求解大型稀疏线性方程组

Ax=b

(1)

其中，A∈Cn×n为正定矩阵；x,b∈Cn。

对于任意的矩阵，A∈Cn×n很自然地拥有一种 Hermitian/skew-Hermitian分裂(HSS)

A=H+S

(2)

这里

(3)

基于上述分裂和古典的交替方向法，BAI等[1]提出了如下求解非 Hermitian 正定线性方程组 (1) 的 HSS 迭代方法。

HSS 迭代方法：给定初始向量x(0)，计算

(4)

直到由式(4)产生的迭代序列{x(k)}收敛，其中，α>0是给定的常数，k=0,1,2，…。

文献[1]中指出 HSS 迭代方法无条件收敛到方程组 (1) 的唯一解，并且它的收缩因子的上界只与系数矩阵 Hermitian 部分的谱有关，与 skew-Hermitian 部分的谱以及迭代涉及的矩阵的特征向量无关。

HSS 迭代方法求解方程组 (1) 时，每一步迭代都需要解2个子方程，一个系数矩阵为αI+H，另一个为αI+S。系数矩阵为αI+H的子方程可以被直接法 (如：Cholesky 分解) 有效求解，而系数矩阵为αI+S的子方程可以被迭代法 (如：Krylov 子空间方法) 求解。另外，当方程组 (1) 的系数矩阵是 Hermitian 时，它的收缩因子的上界与共轭梯度法的收缩因子上界相同。因此，HSS 这类迭代方法在古典分裂迭代法和子空间投影法之间架起了联系的桥梁。

为了加快 HSS 迭代方法的收敛速度，BENZI[4]提出了一种广义的 HSS (GHSS) 迭代方法，它主要是将H分裂成两个 Hermitian 半正定矩阵：H=U+V，其中，V是一种简单矩阵 (如：对角阵)，然后，将V分配给 skew-Hermitian 部分，这样就得到一个新的分裂：A=U+(S+V)，进而得出如下 GHSS 迭代方法。

GHSS 迭代方法：给定初始向量x(0)，计算

(5)

直到由式(5)产生的迭代序列{x(k)}收敛，其中，α>0是给定的常数，k=0,1,2，…。

文献[4]中指出 GHSS 迭代方法收缩因子的上界为

‖(αI-U)(αI+U)-1‖2‖(αI-S-V)(αI+S+V)-1‖2，

如果U和V至少有一个是正定的，那么上式中的两个矩阵范数至少有一个严格小于1，另一个至多是1。而 HSS 迭代方法收缩因子的上界为

‖(αI-H)(αI+H)-1‖2‖(αI-S)(αI+S)-1‖2，

其中，‖(αI-H)(αI+H)-1‖2<1，‖(αI-S)(αI+S)-1‖2=1。虽然这还不足以说明GHSS迭代方法的收敛速度一定快于HSS迭代方法，但是文献[4]用大量的数值实验验证了这一论断。

对于求解大型稀疏线性方程组 (1)，HSS 迭代方法及其变型可能会破坏原系数矩阵的稀疏性，进而增加存储单元，使得求解所需的存储量和计算量增大。

受GHSS迭代思想的启发，文中考虑一种式(1)系数矩阵的新的分裂positive-definite/positive-definite分裂(PPS)：A=P1+P2，其中，P1和P2能保留原系数矩阵的稀疏结构，且均为正定矩阵，于是便得到如下的 PPS 迭代方法。

PPS迭代方法：给定初始向量x(0)，计算

(6)

直到由式(5)产生的迭代序列{x(k)}收敛，其中，β>0是给定的常数，k=0,1,2，…。

显然，该方法较好地保留了原系数矩阵的稀疏性，同时，收敛速度比 HSS 迭代方法更快。

1 收敛性分析

将式(6)方程组中的两步合并，便可得到一步定常迭代

x(k+1)=H(β)x(k)+G(β)b

(7)

其中，

(8)

H(β)为 PPS 迭代的迭代矩阵。

为证 PPS 迭代方法的收敛性，先给出以下引理。

引理2.1[3]令

L(β)=(βI-P)(βI+P)-1，

如果P∈Cn×n是一个正定阵，那么

有了上面这个引理，接下来分析 PPS 迭代方法的收敛性。

定理2.2令A∈Cn×n是一个正定阵，分裂A=P1+P2中的P1和P2也均为正定阵，那么对任意的β>0，迭代式(6)无条件收敛到方程组 (1) 的唯一解。

证明:式(7)中的迭代矩阵H(β)与

相似。显然，H(β)的谱半径

根据引理2.1以及P1和P2均为正定阵，可以得出：不等式右边的两个矩阵范数均小于1，因此，对于任意的β>0均有ρ(H(β))<1，证毕。

下面给出β*取值的一种粗糙估计。

。

于是，由文献[6]中的引理3可知：

ρ(H(β))=ρ(L1(β)L2(β))≤φ(β)

2 PPS迭代方法在对流扩散问题中的应用

本部分讨论PPS 迭代方法在对流扩散问题中的应用[7]，并与 HSS 迭代方法进行了比较，进而说明 PPS 迭代方法对于求解线性方程组 (1) 的有效性。数值实验在Matalb(R2018b) 环境中实现，并且初始向量选取为零向量，终止条件为‖r(k)‖2/‖r(0)‖2≤tol，其中,r(k)表示第k次迭代的残差向量，tol设置为10-6。

在本应用中，求解的线性方程组是来自对(0,1)×(0,1)区域上带有狄利克雷边界条件的二维对流扩散方程

-(uxx+uyy)+η(x,y)ux+ζ(x,y)uy=f(x,y)

(9)

的有限差分离散，其中,系数η和ζ分别表示沿x和y方向的速度分量，且均为非负常数[8]。这里考虑一个均匀的n×n的网格，并且采用标准的五点中心差分近似拉普拉斯算子uxx+uyy，用一阶中心差分近似ux和uy，这样得到线性方程组

Au=v

(10)

其中,u是有限维空间中的解向量且按照自然排序(u11,u21,…,un1,u12,…,unn)。

式(10)的系数矩阵A是一个分块三对角矩阵，它的第k行分别包含了次对角、主对角和超对角块，阶数均为n，

Ak,k-1=dIn,Ak,k=tridiag(b,a,c),Ak,k+1=eIn

(11)

(12)

在第(j,k)个网格点，式(10)右端向量满足vjk=h2fjk，其中，fjk=f(jh,kh)。

当η=ζ时，即b=d，c=e方程组 (10) 的系数矩阵A可以写成A=In⊗B+B⊗In，其中，B=tridiag(b,a/2,c)，定义f(x,y)=ex+y，于是得到下列数值结果，见表1。

表1 HSS和PPS的IT和CPUTable 1 IT and CUP of HSS and PPS

表1中，展示了HSS和PPS迭代方法求解方程组(10)所需的IT和CPU对于不同取值的η和h，数值结果表明:无论是从收敛速度还是计算复杂度的角度考虑，PPS迭代方法求解方程组(10) 的效果都比HSS迭代方法要好，其中,无论η=2还是η=10，随着h的减小(n的增大)，HSS迭代方法求解 (10) 的迭代步数始终超过 PPS 迭代方法的两倍，并且在求解时间上，HSS迭代方法也接近甚至超过PPS 迭代方法的两倍。