一类复对称线性系统的不平衡CRI迭代方法

Lopsided CRI Iteration Method for a Class of Complex Symmetric Linear Systems

摘要: 基于实部与虚部组合(CRI)迭代方法，提出了一种不平衡CRI (LCRI)迭代方法，用于求解复对称半正定线性系统。理论上利用谱理论分析了LCRI方法的收敛性，并给出了拟最优参数的表达式，数值上进一步验证了新方法的高效性。

Abstract: Based on the combination of real and imaginary parts (CRI) iteration method, a lopsided CRI (LCRI) iteration method is proposed for solving complex symmetric positive semi-definite linear systems. By using the spectral theory, we not only analyze the convergence property of the LCRI method, but also obtain the quasi-optimal parameter expression. The efficiency of the new method is further verified numerically.

1. 引言

考虑复对称线性系统

$A\text{x}=\text{b},A\in {ℂ}^{n\times n},\text{ x},\text{b}\in {ℂ}^n,$ (1)

其中A是一个大型稀疏复对称矩阵，形式为$A=W+iT$ ，$i=\sqrt{-1}$ 表示虚数单位，矩阵$W,T\in {ℝ}^{n\times n}$ 均为实对称矩阵。系统(1)常出现在科学计算和工程应用中，例如时间依赖的偏微分方程基于快速傅里叶变换(FFT)的数值解问题[1]、光的散射成像问题[2]、结构动力学问题[3]和量子力学问题[4]等。近年来，学者们针对线性系统(1)提出了许多矩阵分裂迭代法。

当矩阵$W,T$ 是对称半正定的，且其中至少有一个是对称正定的时，系统(1)通常被称为复对称正定线性系统。Bai等人[5]基于系数矩阵A的Hermitian和反Hermitian分裂(HSS)，即$A=H+S$ ，其中

$H=\frac{1}{2}\left(A+A^{*}\right)=W,T=\frac{1}{2}\left(A-A^{*}\right)=iT,$

分别为矩阵A的Hermitian部分和反Hermitian部分，率先提出了HSS迭代方法；由于HSS方法在每一步迭代中都需要求解一个平移的反Hermitian方程组，Bai等人[6]为克服此问题，提出了修正的HSS (MHSS)迭代方法；为了加快MHSS方法的收敛速率，Bai等人[7]将预处理技术应用于MHSS方法，建立了预处理的MHSS (PMHSS)迭代方法；随后，Li等人[8]提出了一种不平衡PMHSS (LPMHSS)迭代方法，并在理论和数值上表明了当系数矩阵的实部占主导时，LPMHSS方法的表现优于MHSS和PMHSS方法。

此外，Wang等人[9]在不限制矩阵$W,T$ 至少有一个是对称正定的情况下，即$W,T$ 均为对称半正定矩阵，且满足$\mathcal{N}\left(W\right)\cap \mathcal{N}\left(T\right)=\left\{0\right\}$ ，提出了一种实部与虚部组合(CRI)迭代方法，其格式如下：

方法1 (CRI迭代方法) 给定初始估计向量${\text{x}}^{\left(0\right)}\in {ℂ}^n$ ，直到迭代序列$\left\{{\text{x}}^{\left(k\right)}\right\}\left(k=0,1,2,\cdots \right)$ 收敛，按照如下方式计算新的迭代解${\text{x}}^{\left(k+1\right)}$

$\{\begin{array}{l}\left(\alpha T+W\right){\text{x}}^{\left(k+\frac{1}{2}\right)}=\left(\alpha -i\right)T{\text{x}}^{\left(k\right)}+\text{b},\\ \left(\alpha W+T\right){\text{x}}^{\left(k+1\right)}=\left(\alpha +i\right)W{\text{x}}^{\left(k+\frac{1}{2}\right)}-i\text{b},\end{array}$

其中，$\alpha$ 为给定的正常数。

上述方法可以看作是PMHSS迭代方法的变体，作者证明了CRI方法的理论性质与数值性能均优于PMHSS方法。

为了进一步提高CRI迭代方法的数值性能，本文借鉴文献[8]的不平衡分裂迭代思想，构建了一种用于求解复对称半正定线性系统的不平衡CRI (LCRI)迭代方法。本文结构安排如下：第2节构建LCRI方法的迭代格式；第3节利用谱理论分析LCRI方法的收敛性和拟最优参数；第4节通过数值实验检验LCRI方法的性能；第5节主要总结本文的研究工作。

2. LCRI迭代方法

本节首先给出LCRI迭代方法的迭代格式。线性系统(1)可以写作以下两种等价形式

$\begin{array}{l}W\text{x}=-iT\text{x}+\text{b},\\ \left(\alpha W+T\right)\text{x}=\left(\alpha +i\right)W\text{x}-i\text{b},\end{array}$

通过交替迭代上述两个不动点方程，可以得到LCRI迭代方法。值得一提的是，为了避免$W\text{x}$ 带来额外的工作量，该方法的迭代格式可描述为：

方法2 (LCRI迭代方法) 给定初始估计向量${\text{x}}^{\left(0\right)}\in {ℂ}^n$ ，直到迭代序列$\left\{{\text{x}}^{\left(k\right)}\right\}\left(k=0,1,2,\cdots \right)$ 收敛，按照如下方式计算新的迭代解${\text{x}}^{\left(k+1\right)}$

$\{\begin{array}{l}{\text{x}}^{\left(k+\frac{1}{2}\right)}=-iT{\text{x}}^{\left(k\right)}+\text{b},\\ \left(\alpha W+T\right){\text{x}}^{\left(k+1\right)}=\left(\alpha +i\right){\text{x}}^{\left(k+\frac{1}{2}\right)}-i\text{b},\end{array}$ (2)

其中，$\alpha$ 为给定的正常数。

在LCRI迭代方法的每一步中，我们只需要求解一个线性子系统，由文献[9]可知其系数矩阵$\alpha W+T$ 是实对称正定的，因此可以通过Cholesky分解精确计算，也可以通过共轭梯度法或多重网格法不精确计算。而CRI迭代方法的每一步需要求解系数矩阵分别为$\alpha T+W$ 和$\alpha W+T$ 的两个线性子系统，由此可以断言LCRI迭代方法会比CRI方法在求解线性系统(1)时更加高效。

经过简单的推导，我们可以将LCRI方法的迭代格式(2)重新表述为标准形式

${\text{x}}^{\left(k+1\right)}=L\left(\alpha \right){\text{x}}^{\left(k\right)}\text{+}G\left(\alpha \right)\text{b,}$ (3)

其中

$L\left(\alpha \right)=\left(1-\alpha i\right){\left(\alpha W+T\right)}^{-1}T,G\left(\alpha \right)=\alpha {\left(\alpha W+T\right)}^{-1}.$

显然，$L\left(\alpha \right)$ 是LCRI方法的迭代格式(2)或(3)的迭代矩阵。因此，LCRI迭代方法收敛当且仅当迭代矩阵的谱半径小于1，即$\rho \left(L\left(\alpha \right)\right)<1$ 。此外，若引入矩阵

$B\left(\alpha \right)=\frac{1}{\alpha }\left(\alpha W+T\right),C\left(\alpha \right)=\left(\frac{1}{\alpha }-i\right)T,$

则LCRI迭代方法可以看作是由矩阵分裂

$A=B\left(\alpha \right)-C\left(\alpha \right)$

导出的，且迭代矩阵可以表示为$L\left(\alpha \right)=B{\left(\alpha \right)}^{-1}C\left(\alpha \right)$ 。

3. LCRI迭代方法的收敛性分析

本节从理论上对LCRI迭代方法的收敛性进行分析，我们首先给出一个有用的引理。

引理1 (见[9]) 设$W,T\in {ℝ}^{n\times n}$ 均为对称半正定矩阵，且满足$\mathcal{N}\left(W\right)\cap \mathcal{N}\left(T\right)=\left\{0\right\}$ ，则存在一个非奇异矩阵P使得

$W=P^{\text{T}}\Lambda P,T=P^{\text{T}}\Gamma P,$

其中，$\Lambda =diag\left({\lambda }_1,{\lambda }_2,\cdots ,{\lambda }_n\right)$ 和$\Gamma =diag\left({\gamma }_1,{\gamma }_2,\cdots ,{\gamma }_n\right)$ 满足

${\lambda }_i+{\gamma }_i=1,{\lambda }_i\ge 0,{\gamma }_i\ge 0\left(i=1,2,\cdots ,n\right).$

基于上述引理的结论，下面我们分析LCRI迭代方法的收敛条件。

定理1 设$A=W+iT$ 是一个复对称非奇异矩阵，其中$W,T\in {ℝ}^{n\times n}$ 均为对称半正定矩阵。假设${\gamma }_{\max }$ 是矩阵T的最大特征值，若参数$\alpha$ 满足

$\left(1-2{\gamma }_{\max }\right)\alpha +2{\gamma }_{\max }\left(1-{\gamma }_{\max }\right)>0,$

则用于求解系统(1)的LCRI迭代方法对于任意的初始估计向量${\text{x}}^{\left(0\right)}\in {ℂ}^n$ 都是收敛的。

证明：由于矩阵W和T都是对称半正定的，且满足$\mathcal{N}\left(W\right)\cap \mathcal{N}\left(T\right)=\left\{0\right\}$ ，因此根据引理1，存在一个非奇异矩阵P，使得

$L\left(\alpha \right)=\left(1-\alpha i\right){\left(\alpha W+T\right)}^{-1}T=\left(1-\alpha i\right)P^{-1}{\left(\alpha \Lambda +\Gamma \right)}^{-1}\Gamma P,$

显然，迭代矩阵$L\left(\alpha \right)$ 相似于

$\tilde{L}\left(\alpha \right)=\left(1-\alpha i\right){\left(\alpha \Lambda +\Gamma \right)}^{-1}\Gamma ,$

则有

$\rho \left(L\left(\alpha \right)\right)=\rho \left(\tilde{L}\left(\alpha \right)\right)=\rho \left(\left(1-\alpha i\right){\left(\alpha \Lambda +\Gamma \right)}^{-1}\Gamma \right)=\underset{1\le j\le n}{\max }\left\{\frac{\sqrt{1+{\alpha }^2}{\gamma }_j}{\alpha {\lambda }_j+{\gamma }_j}\right\}=\underset{1\le j\le n}{\max }\left\{\frac{\sqrt{1+{\alpha }^2}{\gamma }_j}{\alpha +\left(1-\alpha \right){\gamma }_j}\right\}.$

为了方便起见，记

$f\left({\gamma }_j\right)=\frac{\sqrt{1+{\alpha }^2}{\gamma }_j}{\alpha +\left(1-\alpha \right){\gamma }_j},$

对函数$f\left({\gamma }_j\right)$ 关于${\gamma }_j$ 求导可得

$\frac{\text{d}f\left({\gamma }_j\right)}{\text{d}{\gamma }_j}=\frac{\alpha \sqrt{1+{\alpha }^2}}{{\left(\alpha +\left(1-\alpha \right){\gamma }_j\right)}^2}>0,$

则函数$f\left({\gamma }_j\right)$ 在$\left[0,+\infty \right)$ 上单调递增，从而

$\underset{1\le j\le n}{\max }f\left({\gamma }_j\right)=f\left({\gamma }_{\max }\right)=\frac{\sqrt{1+{\alpha }^2}{\gamma }_{\max }}{\alpha +\left(1-\alpha \right){\gamma }_{\max }},$

因此，我们有

$\rho \left(L\left(\alpha \right)\right)=f\left({\gamma }_{\max }\right)<1\iff \left(1-2{\gamma }_{\max }\right)\alpha +2{\gamma }_{\max }\left(1-{\gamma }_{\max }\right)>0,$

定理证毕。

由定理1可知，LCRI迭代方法的收敛速度可能取决于两个因素，一个是对称半正定矩阵T的谱，另一个是迭代参数$\alpha$ 的选取。在实际应用中，为了避免LCRI方法选取实验最优参数所造成的额外工作量，我们希望找到最优参数$\alpha$ 的良好估计，以极小化迭代矩阵的谱半径，从而提高LCRI方法的收敛速率。因此，我们在下述推论中给出迭代参数$\alpha$ 的估计式。

推论1 假设定理1的条件成立，设${\gamma }_{\max }$ 是矩阵T的最大特征值，则LCRI迭代方法的拟迭代参数为

${\alpha }^{\text{*}}\text{=arg}\underset{\alpha }{\text{min}}\left\{\rho \left(L\left(\alpha \right)\right)\right\}=\text{arg}\underset{\alpha }{\text{min}}\left\{\frac{\sqrt{1+{\alpha }^2}{\gamma }_{\max }}{\alpha +\left(1-\alpha \right){\gamma }_{\max }}\right\}=\frac{1}{{\gamma }_{\max }}-1,$

此时对应的迭代矩阵的谱半径为

$\rho \left(L\left({\alpha }^{*}\right)\right)=\frac{{\gamma }_{\max }}{\sqrt{2{\gamma }_{\max }^2-2{\gamma }_{\max }+1}}.$

4. 数值实验

本节通过数值实验来检验LCRI迭代方法的数值性能，我们将从迭代步数(IT)和迭代时间(CPU)这两个方面比较PMHSS、LPMHSS、CRI和LCRI迭代方法。在实验中，我们对所有测试方法所含的子系统均使用Cholesky分解进行求解，并在分解之前使用对称近似极小度方法(SYMAMD)进行重排，具体细节可参考文献[10]。

在算法的实现过程中，我们选取初始向量${\text{x}}^{\left(0\right)}$ 为零向量，停止标准均为当前迭代解${\text{x}}^{\left(k\right)}$ 的相对残差满足条件

$\frac{{\Vert \text{b}-A{\text{x}}^{\left(k\right)}\Vert }_2}{{\Vert \text{b}\Vert }_2}\le {10}^{-6},$

或者迭代步数达到$k_{\max }=1000$ 。对于所有测试方法中涉及的迭代参数，我们选取实验最优值${\alpha }_{\exp }$ ，即使得迭代步数最小的值。若最优参数值形成一个区间，则在此区间内选取使得迭代时间最小的值作为最优参数。PMHSS和LPMHSS方法中使用的预处理矩阵V选取为W。所有数值结果均通过MATLAB [版本号9.12.0.1884302 (R2022a)]计算得到，计算机配置：中央处理器1.70 GHz [Intel(R) Core(TM) i5 1240P]，内存15.7 G，操作系统Windows 11。

算例1 (参考[8] [9]) 考虑如下形式的复对称线性系统

$\left[\left(-{\omega }^{2}M+K\right)+i\left(\omega C_V+C_H\right)\right]\text{x}=\text{b},$

其中$\omega$ 为圆周频率，M是惯性矩阵，K是刚度矩阵，$C_V$ 为粘性阻尼矩阵，$C_H$ 为滞后阻尼矩阵。这里取$M=I$ ，$C_V=10M$ ，$C_H=\mu K$ ，其中$\mu$ 是阻尼系数。此外，假设K为二维区域$\Omega =\left[0,1\right]\times \left[0,1\right]$ 上一致网格剖分下对具有齐次Dirichlet边值条件的负拉普拉斯算子的五点中心差分格式近似矩阵，每个方向上网格步长都取$h=1/\left(m+1\right)$ 。因此矩阵$K=I_m\otimes V_m+V_m\otimes I_m$ ，其中$V_m=h^{-2}\text{tridiag}\left(-1,2,-1\right)\in {ℝ}^{m\times m}$ 。实验中，取圆周频率$\omega =0.5$ ，右端向量$b=\left(\text{1}+i\right)A1$ ，其中$1$ 是分量全为1的向量。此外，通过对系统方程两边同时乘以$h^2$ 使问题标准化。

迭代方法的收敛速度在很大程度上取决于迭代矩阵的谱半径，图1描绘了网格数$m=64$ 的情况下，CRI和LCRI迭代方法在不同阻尼系数($\mu =0.1,0.01,0.001$ )下迭代矩阵谱半径的比较。从图1可以看出，当$\mu$ 较小时(即系数矩阵实部占主导)，LCRI方法的迭代矩阵谱半径比CRI方法的迭代矩阵谱半径小很多，但当$\mu$ 较大时(即系数矩阵虚部占主导)，CRI方法的表现更好。

Figure 1. The spectral radius of the iteration matrix $\rho$ versus the parameter $\alpha$

图1. 迭代矩阵的谱半径$\rho$ 与参数$\alpha$ 的关系图

下面我们取阻尼系数$\mu =0.001$ 进行实验。在表1中，我们列出了PMHSS、LPMHSS、CRI和LCRI方法的实验最优参数范围以及LCRI方法的拟最优参数${\alpha }^{*}$ 的值。实验中发现LPMHSS和LCRI方法对于参数$\alpha$ 的取值不太敏感，这里我们人为地将参数$\alpha$ 的上界设为1000，而PMHSS和CRI方法的实验最优参数取值范围则相对较小，且随着网格数m的增大而逐渐变窄，这说明LPMHSS和LCRI方法在参数选取上更加灵活。此外，LCRI方法的拟最优参数${\alpha }^{*}$ 的值均在其实验最优参数取值范围内，这表明由推论1所得到的拟最优参数${\alpha }^{*}$ 是LCRI方法迭代参数的一个不错的选择。

Table 1. The ranges of experimental optimal parameter or the values of quasi-optimal parameter of the tested methods

表1. 测试方法的实验最优参数取值范围或拟最优参数值

方法	参数	$m=64$	$m=128$	$m=256$	$m=\text{512}$
PMHSS	${\alpha }_{\exp }$	[0.75,1.32]	[0.76,1.31]	[0.76,1.30]	[0.76,1.30]
LPMHSS	${\alpha }_{\exp }$	[0.96,1000]	[1.10,1000]	[1.42,1000]	[0.62,1000]
CRI	${\alpha }_{\exp }$	[0.63,1.61]	[0.63,1.58]	[0.64,1.56]	[0.66,1.53]
LCRI	${\alpha }_{\exp }$	[0.96,1000]	[1.10,1000]	[1.42,1000]	[0.62,1000]
LCRI	${\alpha }^{*}$	107.95	119.49	122.83	123.71

所有测试方法对应于不同网格大小的数值结果如表2所示，包括实验最优参数、迭代步数和迭代时间。可以观察到，PMHSS方法的迭代步数不受网格数m的影响，LPMHSS、CRI和LCRI方法的迭代步数均随网格数m的增大而逐渐减小。无论是从迭代步数还是迭代时间来看，LPMHSS、CRI和LCRI方法均对PMHSS方法有着显著的改进效果。LCRI方法和LPMHSS方法的迭代步数相同，但前者的迭代时间小于后者。事实上，LCRI方法可以看作将LPMHSS方法的预处理矩阵V取为W，且相较于LPMHSS方法少求解一个系数矩阵为W的线性子系统，从而降低了计算成本。此外，当网格数$m=64,128,\text{256}$ 时，LCRI方法的迭代步数和迭代时间均优于CRI方法。当网格数$m=512$ 时，这两种方法的迭代步数相同，但LCRI方法在迭代时间上的优势更为明显，这表明LCRI方法能够更高效地处理大规模问题。

Table 2. The numerical results of the tested methods with experimental optimal parameter values

表2. 测试方法利用实验最优参数所得的数值结果

方法		$m=64$	$m=128$	$m=256$	$m=\text{512}$
PMHSS	${\alpha }_{\exp }$	0.99	1.15	1.01	0.76
	IT	34	34	34	34
	CPU	0.0847	0.1961	0.6949	3.0806
LPMHSS	${\alpha }_{\exp }$	940	630	420	130
	IT	6	5	4	4
	CPU	0.0219	0.0457	0.0919	0.3787
CRI	${\alpha }_{\exp }$	1.17	0.80	1.02	0.66
	IT	7	6	5	4
	CPU	0.0321	0.0533	0.1443	0.4407
LCRI	${\alpha }_{\exp }$	130	690	70	60
	IT	6	5	4	4
	CPU	0.0179	0.0332	0.0668	0.2008

因此，当系数矩阵的实部W相对于虚部T占主导且问题规模较大时，我们倾向于使用LCRI迭代方法来求解复对称线性系统(1)。

5. 总结

本文基于CRI迭代方法，利用不平衡分裂迭代思想，构建了一种不平衡CRI (LCRI)迭代方法来求解复对称半正定线性系统。理论分析表明，当迭代参数$\alpha$ 满足一定条件时，LCRI方法可以收敛到线性系统(1)的唯一精确解。并通过极小化迭代矩阵的谱半径，给出了拟最优参数${\alpha }^{*}$ 的表达式。数值结果表明，当线性系统(1)的系数矩阵实部占主导且问题规模较大时，LCRI方法具有高效性与稳健性。