DOI:10.12677/AAM.2018.711164,PDF,HTML,XML,被引量下载: 1,603浏览: 3,413国家自然科学基金支持
作者:张 婷,张乃敏：温州大学数学系，浙江 温州
关键词:复对称线性系统；MPMHSS迭代法；动量项；收敛；Complex Symmetric Linear System；MPMHSS Iterative Method；Momentum Term

文章引用：张婷, 张乃敏. 求解复对称线性系统的MPMHSS迭代法[J]. 应用数学进展, 2018, 7(11): 1405-1410. https://doi.org/10.12677/AAM.2018.711164

1. 引言

复对称线性系统问题为 $Ax=b$， $A\in C^{n\times n}$， $x$， $b\in C^n$(1)

其中 $A\in C^{n\times n}$是一个形式为 $A=W+iT$的复对称矩阵， $W$和 $T$分别是实对称正定阵和实对称半正定阵， $i=\sqrt{-1}$为虚数单位。事实上，如果令且 $b=p+iq$，复对称线性系统(1)可以转变为实的块线性系统：

$C{\rm K}\equiv \left(\begin{array}{cc}W& -T\\ T& W\end{array}\right)\left(\begin{array}{c}y\\ z\end{array}\right)=\left(\begin{array}{c}p\\ q\end{array}\right)\equiv g$， $C\in R^{2n\times 2n}$， $g\in R^{2n}$， $y,z,p,q\in R^n$(2)

这一类复对称线性系统广泛用于解决科学计算和工程应用领域的问题，如计算流体力学 [1] 等。为了有效地解决复对称线性系统(1)和(2)，文章 [2] 通过添加预处理子提出PMHSS迭代法，从而提高MHSS迭代法的收敛速率。由于复对称线性系统(2)可以被视为一种特殊的鞍点问题，因此本文结合文献 [3] 提出的动量项加速求解鞍点问题的想法，构造带动量项加速的PMHSS迭代法(MPMHSS)，以此更有效地求解复对称线性系统问题(2)。

2. MPMHSS迭代法

首先将文章 [2] 提出的PMHSS迭代法转化为块实值形式，然后使用类似于 [3] 中的动量项加速求解鞍点问题的方法，建立MPMHSS迭代法。

算法1 (PMHSS迭代法)：令 $x^{\left(0\right)}\in C^n$是一个任意的初始向量， $\alpha$是被给定的正参数。对于 $k=0,1,2\cdots$，根据下面的步骤计算迭代 $x^{k+1}$，直到迭代序列 ${\left\{x^{\left(k\right)}\right\}}_0^{\infty }$收敛：

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

其中， $V\in R^{n\times n}$是一个给定的对称正定阵。

根据文章 [4] 可知，PMHSS迭代法对于(2)中的系数矩阵 $C$而言，其预处理子为

$M=\left(1+\frac{1}{\alpha }\right)Q\left(\begin{array}{cc}\alpha W+T& 0\\ 0& \alpha W+T\end{array}\right)=\frac{\alpha +1}{2\alpha }\left(\begin{array}{cc}\alpha W+T& -\left(\alpha W+T\right)\\ \alpha W+T& \alpha W+T\end{array}\right)$，其中 $Q=\frac{1}{2}\left(\begin{array}{cc}I& -I\\ I& I\end{array}\right)$

通过对系数矩阵 $C$分裂如下 $C=M-N$，可知PMHSS迭代的迭代矩阵为 $H=M^{-1}N=I-M^{-1}C$。通过详细计算可得

$H=I-M^{-1}C=\frac{\alpha }{\alpha +1}\{\begin{array}{l}\begin{array}{cc}{\left(\alpha W+T\right)}^{-1}\left(\alpha W+\frac{1}{\alpha }T\right)& {\left(\alpha W+T\right)}^{-1}\left(T-W\right)\end{array}\\ \begin{array}{cc}{\left(\alpha W+T\right)}^{-1}\left(W-T\right)& {\left(\alpha W+T\right)}^{-1}\left(W+T\right)\end{array}\end{array}$

其中 $\alpha V+T,\alpha V+W,\alpha W+T$都为对称正定阵。注意到PMHSS迭代法可写为

$\{\begin{array}{c}{y}^{\left(k+1\right)}=\frac{\alpha }{\alpha +1}{\left(\alpha W+T\right)}^{-1}\left[\left(\alpha W+\frac{1}{\alpha }T\right)y^{\left(k\right)}+\left(T-W\right)z^{\left(k\right)}+p+q\right]\\ z^{\left(k+1\right)}=\frac{\alpha }{\alpha +1}{\left(\alpha W+T\right)}^{-1}\left[\left(W-T\right)y^{\left(k\right)}+\left(\alpha W+\frac{1}{\alpha }T\right)z^{\left(k\right)}-p+q\right]\end{array}$

下面我们对算法1添加动量项加速，得到如下算法2：

算法2 (MPMHSS迭代法)：令 ${\left(y^{{\left(0\right)}^T},z^{{\left(0\right)}^T}\right)}^T\in R^{2n}$是一个任意的初始向量，其中 $y^{\left(0\right)}\in R^n,z^{\left(0\right)}\in R^n$， $\mu$为动量项参数， $\alpha$是给定的正参数， ${(\cdot )}^T$定义了一个向量或矩阵的转置。

对于 $k=0$时，

$\{\begin{array}{c}{y}^{\left(1\right)}=\frac{\alpha }{\alpha +1}{\left(\alpha W+T\right)}^{-1}\left[\left(\alpha W+\frac{1}{\alpha }T\right)y^{\left(0\right)}+\left(T-W\right)z^{\left(0\right)}+p+q\right]\\ z^{\left(1\right)}=\frac{\alpha }{\alpha +1}{\left(\alpha W+T\right)}^{-1}\left[\left(W-T\right)y^{\left(0\right)}+\left(\alpha W+\frac{1}{\alpha }T\right)z^{\left(0\right)}-p+q\right]\end{array}$

对于 $k=1,2\cdots$时，

$\{\begin{array}{c}{y}^{\left(k+1\right)}=\frac{\alpha }{\alpha +1}{\left(\alpha W+T\right)}^{-1}\left[\left(\alpha W+\frac{1}{\alpha }T\right)y^{\left(k\right)}+\left(T-W\right)z^{\left(k\right)}+p+q\right]+\mu \left(y^{\left(k\right)}-y^{\left(k-1\right)}\right)\\ z^{\left(k+1\right)}=\frac{\alpha }{\alpha +1}{\left(\alpha W+T\right)}^{-1}\left[\left(W-T\right)y^{\left(k\right)}+\left(\alpha W+\frac{1}{\alpha }T\right)z^{\left(k\right)}-p+q\right]+\mu \left(z^{\left(k\right)}-z^{\left(k-1\right)}\right)\end{array}$

若令 ${{\rm K}}^{\left(k\right)}={\left(y^{{\left(k\right)}^T},z^{{\left(k\right)}^T}\right)}^T$，则算法2按如下格式计算 ${{\rm K}}^{\left(k+1\right)}$：

${{\rm K}}^{\left(k+1\right)}=H{{\rm K}}^{\left(k\right)}+M^{-1}b+\mu \left({{\rm K}}^{\left(k\right)}-{{\rm K}}^{\left(k-1\right)}\right)\left(k>0\right)$(3)

因此，算法2也可以写为：

$\{\begin{array}{l}{{\rm K}}^{\left(1\right)}=H{{\rm K}}^{\left(0\right)}+M^{\text{-}1}g\\ {{\rm K}}^{\left(k+1\right)}=H{{\rm K}}^{\left(k\right)}+M^{-1}g+\mu \left({{\rm K}}^{\left(k\right)}-{{\rm K}}^{\left(k-1\right)}\right)\left(k=1,2\cdots \right)\end{array}$(4)

令 $\Phi =\left(\begin{array}{cc}\begin{array}{l}C-\mu M\\ -I\end{array}& \begin{array}{l}\mu M\\ I\end{array}\end{array}\right)$， $\Psi =\left(\begin{array}{cc}M& 0\\ 0& I\end{array}\right)$， ${\xi }^{\left(k\right)}=\left(\begin{array}{l}{{\rm K}}^{\left(k\right)}\\ {{\rm K}}^{\left(k-1\right)}\end{array}\right)$， $F=\left(\begin{array}{c}g\\ 0\end{array}\right)$，则易知线性系统(2)等价于

$\Phi \xi =F\iff \left(\begin{array}{cc}\begin{array}{l}C-\mu M\\ -I\end{array}& \begin{array}{l}\mu M\\ I\end{array}\end{array}\right)\left(\begin{array}{c}{\rm K}\\ {\rm K}\end{array}\right)=\left(\begin{array}{c}g\\ 0\end{array}\right)$(5)

其中 $M$是算法1的预处理子。

对于线性系统(5)，对系数矩阵 $$进行分裂： $\Phi =\Psi -\left(\Psi -\Phi \right)$，则算法2对应的迭代矩阵为

$\Gamma =I-{\Psi }^{-1}\Phi =I-{\left(\begin{array}{cc}M& 0\\ 0& I\end{array}\right)}^{-1}\left(\begin{array}{cc}\begin{array}{l}C-\mu M\\ -I\end{array}& \begin{array}{l}\mu M\\ I\end{array}\end{array}\right)=\left(\begin{array}{cc}\begin{array}{l}\mu I+H\\ I\end{array}& \begin{array}{l}-\mu I\\ 0\end{array}\end{array}\right)$(6)

从而(4)也可以被写为

$\{\begin{array}{l}\begin{array}{cc}{{\rm K}}^{\left(1\right)}=H{{\rm K}}^{\left(0\right)}+M^{\text{-}1}g,& {\xi }^{\left(1\right)}=\left({{\rm K}}^{{\left(1\right)}^T},{{\rm K}}^{{\left(2\right)}^T}\right)\end{array}\\ \begin{array}{cc}{\xi }^{\left(k+1\right)}=\Gamma {\xi }^{\left(k\right)}+r,& k=1,2\cdots ,\end{array}\end{array}$，其中 $r=\left(\begin{array}{l}{M}^{-1}g\\ 0\end{array}\right)$(7)

由于线性系统(2)等价于线性系统(5)，因此 $\Phi$是非奇异阵当且仅当 $C$是非奇异阵，而且我们可以通过(7)注意到MPMHSS迭代收敛的充要条件是谱半径 $\rho \left(\Gamma \right)<1$。

3. MPMHSS迭代法的收敛性分析

本节主要讨论MPMHSS迭代法收敛的收敛理论和使迭代矩阵的谱半径最小的最优参数 ${\alpha }_{opt}$和 ${\mu }_{opt}$，文中 $\sigma (\cdot )$为谱集。

引理1 [5] ：二次方程 $x^2-dx+c=0$的两根的模都小于1的充要条件是 $\left|c\right|<1$和 $\left|d-\bar{d}c\right|<1-{\left|c\right|}^2$，其中 $\bar{d}$是 $d$的共轭数。

引理2 [3] ：对于PMHSS迭代法的迭代阵 $H$和MPMHSS迭代法的迭代阵 $\Gamma$，

1) 若 $\mu =0$，则 $H$的特征根和 $\Gamma$的特征根的关系是： $\sigma \left(\Gamma \right)=\sigma \left(H\right)\cup \left\{0\right\}$；

2) 若 $\mu \ne 0$，则对任意的 $\epsilon \in \sigma \left(H\right)$，都有 $\lambda \in \sigma \left(\Gamma \right)$，使 ${\lambda }^2-\left(\mu +\epsilon \right)\lambda +\mu =0$(8)

相反地，对于任意的 $\lambda \in \sigma \left(\Gamma \right)$，都存在 $\epsilon \in \sigma \left(H\right)$，使得 ${\lambda }^2-\left(\mu +\epsilon \right)\lambda +\mu =0$成立。

从引理2中，我们可以看出 $\Gamma$是非奇异阵当且仅当 $\mu \ne 0$。

当取定PMHSS方法的最优参数(见文献 [2] ) ${\alpha }_{opt}$后，下述定理给出了最小化的最优动量项参数。

定理3： $W$是对称正定阵， $T$是对称半正定阵。对于MSBTS迭代法， $\rho \left(\Gamma \right)<1$成立的充要条件是迭代阵 $H$的特征值 $\epsilon =a+ib\left(a,b\in R\right)$满足点 $\left(a,b\right)$位于椭圆内部：

$\frac{{\left(x+\mu \right)}^2}{{\left(1+\mu \right)}^2}+\frac{y^2}{{\left(1-\mu \right)}^2}=1,-1<\mu <1$(9)

此外，当取定PMHSS方法的最优参数 ${\alpha }_{opt}=\sqrt{{\eta }_{\min }{\eta }_{\max }}$， ${\eta }_{\min }$和 ${\eta }_{\max }$分别是 $V^{-1}W$的最小特征值和最大特征值时，若 $H$的所有特征值满足 $-3<{\epsilon }_1\le {\epsilon }_2\le \cdots \le {\epsilon }_s<1\left(s\le 2n\right)$，那么最小化 $\rho \left(\Gamma \right)$的最优动量项参数为 ${\mu }_{opt}=\max \left\{{\left(1-\sqrt{1-{\epsilon }_s}\right)}^2,{\left(1-\sqrt{1-{\epsilon }_1}\right)}^2\right\}$，且 $\underset{\mu }{\min }\rho \left(\Gamma \right)=\sqrt{{\mu }_{opt}}$。

证明：由引理1和引理2可知，令 $d=\mu +\epsilon$， $c=\mu$，则

$\rho \left(\Gamma \right)<1\iff \left|d-\bar{d}c\right|<1-{\left|c\right|}^2\iff \left|\mu +\epsilon -\mu \left(\mu +\bar{\epsilon }\right)\right|<1-{\mu }^2\iff \frac{{\left(a+\mu \right)}^2}{{\left(1+\mu \right)}^2}+\frac{b^2}{{\left(1-\mu \right)}^2}=1,\left(-1<\mu <1\right)$

此外，由引理1可知 $\rho \left(\Gamma \right)<1\iff \{\begin{array}{l}\left|\mu \right|<1\\ \left|\mu +\epsilon \right|<1+\mu \end{array}\iff \{\begin{array}{l}\epsilon <1\\ \epsilon >-3\end{array}$，则 $H$的所有特征值满足 $\epsilon \in \left(-3,1\right)$。

接下来证明最优动量项为 ${\mu }_{opt}=\max \left\{{\left(1-\sqrt{1-{\epsilon }_s}\right)}^2,{\left(1-\sqrt{1-{\epsilon }_1}\right)}^2\right\}$，且 $\underset{\mu }{\min }\rho \left(\Gamma \right)=\sqrt{{\mu }_{opt}}$。

由引理2可知，迭代阵 $\Gamma$的特征值是：

${\lambda }_{j1}=\frac{\mu +{\epsilon }_j+\sqrt{{\left(\mu +{\epsilon }_j\right)}^2-4\mu }}{2}$， ${\lambda }_{j2}=\frac{\mu +{\epsilon }_j-\sqrt{{\left(\mu +{\epsilon }_j\right)}^2-4\mu }}{2}$

其中 $-1<\mu <1,j=1,2,\cdots ,s$。

1) 当 $-1<\mu <1$且 ${\left(\mu +{\epsilon }_j\right)}^2-4\mu \le 0$时 ${\left(1-\sqrt{1-{\epsilon }_j}\right)}^2\le \mu <1$，则 $\max \{\left|{\lambda }_{j1}\right|,\left|{\lambda }_{j2}\right|\}=\sqrt{\mu }$。

2) 当 $-1<\mu <1$且 ${\left(\mu +{\epsilon }_j\right)}^2-4\mu >0$时， $-1<\mu <{\left(1-\sqrt{1-{\epsilon }_j}\right)}^2$且 ${\lambda }_{j1}$和 ${\lambda }_{j2}$是两个不同实根。

定义：当 $\mu +{\epsilon }_j>0$时， $f_{j1}:=f_{j1}\left(\mu \right)=\frac{\mu +{\epsilon }_j+\sqrt{{\left(\mu +{\epsilon }_j\right)}^2-4\mu }}{2}={\lambda }_{j1}$；

当 $\mu +{\epsilon }_j<0$时， $f_{j2}:=f_{j2}\left(\mu \right)=\frac{-\left(\mu +{\epsilon }_j\right)+\sqrt{{\left(\mu +{\epsilon }_j\right)}^2-4\mu }}{2}=-{\lambda }_{j2}$。

容易看出：

$f_{j1}\left(\mu \right)=\frac{1}{2}\left[\left(\mu +{\epsilon }_j\right)+\sqrt{{\left(\mu +{\epsilon }_j\right)}^2-4\mu }\right]>\sqrt{\mu }$， $f_{j2}\left(\mu \right)=\frac{1}{2}\left[-\left(\mu +{\epsilon }_j\right)+\sqrt{{\left(\mu +{\epsilon }_j\right)}^2-4\mu }\right]>\sqrt{\mu }$

由于 $\frac{\text{d}{f}_{j1}\left(\mu \right)}{\text{d}\mu }=\frac{\sqrt{{\left(\mu +{\epsilon }_j\right)}^2-4\mu }+\left(\mu +{\epsilon }_j\right)-2}{2\sqrt{{\left(\mu +{\epsilon }_j\right)}^2-4\mu }}<0$， $\frac{\text{d}{f}_{j2}\left(\mu \right)}{\text{d}\mu }=\frac{1}{2}\left(-1-\frac{2-\mu -{\epsilon }_j}{\sqrt{{\left(\mu +{\epsilon }_j\right)}^2-4\mu }}\right)<0$

所以 $f_{j1}\left(\mu \right)$和 $f_{j2}\left(\mu \right)$都关于 $\mu$单调递减 $\left(j=1,2,\cdots ,s\right)$。因此为使 $\rho \left(\Gamma \right)$最小化，动量项参数 $\mu$应最大化。

定义： $l\left(\mu \right)=\{\begin{array}{cc}0,& -1<\mu <0\\ \sqrt{\mu },& 0\le \mu <1\end{array}$，因此对于任意的 $\mu$，有 $\rho \left(\Gamma \right)=\max \left\{l\left(\mu \right),f_{s1}\left(\mu \right),f_{12}\left(\mu \right)\right\}$。

当 $\mu =0$时，

$f_{s1}\left(0\right)=\frac{{\epsilon }_s+\sqrt{{\epsilon }_s{}^2-0}}{2}=\{\begin{array}{cc}{\epsilon }_s& \left({\epsilon }_s>0\right)\\ 0& \left({\epsilon }_s\le 0\right)\end{array}$， $f_{12}\left(0\right)=\frac{-{\epsilon }_1+\sqrt{{\epsilon }_1{}^2-0}}{2}=\{\begin{array}{l}\begin{array}{cc}-{\epsilon }_1& \left({\epsilon }_1<0\right)\end{array}\\ \begin{array}{cc}0& \left({\epsilon }_1\ge 0\right)\end{array}\end{array}$

可以看出 $f_{s1}\left(0\right)$和 $f_{12}\left(0\right)$都是非负的，因此我们可以得到 $\max \left\{f_{s1}\left(0\right),f_{12}\left(0\right)\right\}\ge 0$。

由 $f_{s1}\left(\mu \right)$和 $f_{12}\left(\mu \right)$关于 $\mu$单调递减， $l\left(\mu \right)$关于 $\mu$单调递增， $\max \left\{f_{s1}\left(0\right),f_{12}\left(0\right)\right\}\ge 0$且 $l\left(0\right)=0$可知：最优动量项参数 ${\mu }_{opt}\in \left[0,1\right)$，且 $\max \left\{f_{s1}\left(\mu \right),f_{12}\left(\mu \right)\right\}=\sqrt{\mu }$。

接下来分别求解 $f_{s1}\left(\mu \right)=\sqrt{\mu }$和 $f_{12}\left(\mu \right)=\sqrt{\mu }$，可得 $\mu ={\left(1\pm \sqrt{1-{\epsilon }_s}\right)}^2$和 $\mu ={\left(1\pm \sqrt{1-{\epsilon }_1}\right)}^2$。又因 ${\left(1+\sqrt{1-{\epsilon }_j}\right)}^2>1\left(j=1,2,\cdots \cdots ,s\right)$与 $-1<\mu <1$矛盾，因此

${\mu }_{opt}=\max \left\{{\left(1-\sqrt{1-{\epsilon }_s}\right)}^2,{\left(1-\sqrt{1-{\epsilon }_1}\right)}^2\right\}$， $\underset{\mu }{\min }\rho \left(\Gamma \right)=\sqrt{{\mu }_{opt}}$

4. 数值实验

这一节，我们通过数值实验来说明MPMHSS迭代法求解复对称线性系统的有效性。所有数值实验将利用Matlab解决，在Intel (R) Celeron (R) CPUN3060@1.60GHz，内存为4.0 GB的个人电脑上进行。在下面的实验中，取零向量 ${\left({\left(y_0\right)}^T,{\left(z_0\right)}^T\right)}^T$为初始向量，选择右端向量 $\left(p^T,q^T\right)\in R^{n+n}$使线性系统的精确解为 ${\left({\left(y^{\ast }\right)}^T,{\left(z^{\ast }\right)}^T\right)}^T\in R^{n+n}$。其中IT表示迭代步数，CPU表示运行时间，RES表示残差比，停机准则为

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

算例1 [2] ：复对称线性系统(1)的形式为

$\left(K+\frac{3-\sqrt{3}}{\tau }I\right)+i\left(K+\frac{3+\sqrt{3}}{\tau }I\right)x=b$(10)

其中 $\tau$为时间步长， $K$是五点中心差分阵，近似于负的拉普拉斯算符 $L=-\Delta$，并且是在方形域 $\left[0,1\right]\times \left[0,1\right]$内，基于均匀网格 $m\times m$的细分， $h={\left(m+1\right)}^{-1}$。(10)中的矩阵 $K\in R^{n\times n}$拥有张量积形式 $K=I\otimes V_m+V_m\otimes I$，其中 $V_m=h^{-2}tridiag\left(-1,2,-1\right)\in R^{m\times m}$且 $n=m^2$。

在实验中，我们令 $\tau =h$，并在等式两边同时乘 $h^2$正规化系数矩阵和右端向量，详细请参考 [6] 。

通过对 $m$的不同取值，我们得到了不同大小的系数矩阵。在表1中，我们列出了PMHSS迭代法和MPMHSS迭代法分别对应的最优参数、IT、CPU和残差比(RES)。通过比较表1中的结果，可以看出MPMHSS迭代法的迭代步数和CPU消耗都比PMHSS迭代的小，因此MPMHSS迭代的收敛速度更快，对求解复对称线性系统更有效。

致谢

本研究获国家自然科学基金(项目编号61572018)和浙江省自然科学基金(项目编号LY15A010016)的资助。

参考文献

[1]	Elman, H.C., Silvester, D.J. and Wathen, A.J. (1997) Iterative Methods for Problems in Computational Fluid Dynamics. In: Iterative Methods in Scientific Computing, Springer-Verlag, Singapore, 271-327.
[2]	Bai, Z.-Z., Benzi, M. and Chen, F. (2011) On Preconditioned MHSS Iteration Methods for Complex Symmetric Linear Systems. Numerical Algorithms, 56, 297-317. https://doi.org/10.1007/s11075-010-9441-6
[3]	Zhang, N.M. (2015) On Parameter Acceleration Methods for Saddle Point Problems. Journal of Computational and Applied Mathematics, 288, 169-179. https://doi.org/10.1016/j.cam.2015.04.028
[4]	Bai, Z.-Z. (2013) Rotated Block Triangular Preconditioned Based on PMHSS. Science China Mathematics, 56, 2523-2538. https://doi.org/10.1007/s11425-013-4695-9
[5]	Miller, J.J.H. (1971) On the Location of Zeros of Certain Classes of Polynomials with Applications to Numerical Analysis. IMA Journal of Applied Mathematics, 8, 397-406. https://doi.org/10.1093/imamat/8.3.397
[6]	Axelsson, O. and Kucherov, A. (2000) Real Valued Iterative Methods for Solving Complex Symmetric Linear Systems. Numerical Linear Algebra with Applications, 7, 197-218. https://doi.org/10.1002/1099-1506(200005)7:4<197::AID-NLA194>3.0.CO;2-S