对于大规模稀疏线性方程组$ Ax=b $, 其中$ A \in R^{n\times n}$是非奇异且稀疏的, $ x, b \in R^{n}$为$n$维列向量, 该方程组的求解常使用迭代法, 而在众多的迭代法当中, Krylov子空间方法是当前应用最为广泛的方法. Krylov子空间方法主要基于两类过程来产生Krylov子空间基, 一类是基于Arnoldi过程构造Krylov子空间的正交基的方法, 如GMRES^[1]方法, FOM ^[1]方法等, 这类算法的研究也是目前研究的热点, 如文献[3]基于文献[2]提出一种修正的GMRES方法, 提高了GMRES方法的收敛速度; 另一类是基于Lanczos双正交过程产生的Krylov子空间的算法, 如QMR ^[4], BiCG ^[1]等方法, 文献[5]中提出一种Lanczos双$A$ -正交过程, 数值试验表明该过程能够有效减少迭代过程中的存储量.本文将文献[5]中提到的用Lanczos双$A$ -正交过程来产生的Krylov子空间基应用到QMR方法中, 同时结合文献[2]和[3]中的思想对该方法予以修正, 得到一种基于Lanczos双$A$ -正交的修正的QMR算法, 旨在加快QMR方法收敛速度的同时减少迭代过程中的存储量.数值试验表明新算法在这些方面是优于QMR方法的.

对于两个实向量$x, y\in R^{n}$, 它们的标准内积记为如下形式

$ \langle x, y\rangle=x^{T}y=\sum\limits_{i=1}^n x_{i}y_{i}, $

对于两个给定的向量$\nu_{1}$和$\omega_{1}$, 它们的欧几里德内积满足$\langle \omega_{1}, A \nu_{1}\rangle=1$, 通过以下回代关系产生Lanczos向量$\nu_{j}$和$\omega_{j}$, 以及标量$\delta_{j}$, $\beta_{j}$, $ j=1, 2, \cdots, m, $

$ \delta_{j+1}\nu_{j+1}=A \nu_{j} -\beta_{j}\nu_{j-1}-\alpha_{j}\nu_{j}=s_{j+1}, \\ \beta_{j+1}\omega_{j+1}=A^{T} \omega_{j}-\delta_{j}\omega_{j-1}-\alpha_{j}\omega_{j}=t_{j+1}, $

其中

$ \alpha_{j}=\langle\omega_{j}, A(A\nu_{j})\rangle, \qquad \delta_{j+1}=\sqrt{|\langle\hat{\omega}_{j+1}, A\hat{\nu}_{j+1}\rangle|}, \qquad \beta_{j+1}=\frac{\langle\hat{\omega}_{j+1}, A\hat{\nu}_{j+1}\rangle}{\delta_{j+1}}, $

这些标量的选择决定了产生的向量$\nu_{i}$, $\omega_{i}$是彼此$A$ -正交的, 且两组标量满足以下关系

$ \beta_{j+1}\delta_{j+1}=s^{T}_{j+1}At_{j+1}=\langle s_{j+1}, At_{j+1} \rangle. $

具体的Lanczos双$A$ -正交过程如下:

算法 1  (Lanczos双$A$ -正交过程)

1. 选取两个向量$\nu_{1}$和$\omega_{1}$, 使得$\langle \omega_{1}, A \nu_{1} \rangle=1$;

2. 令$\beta_{1}=\delta_{1}=0$, $\omega_{0}=\nu_{0}=\mathbf{0}$:

3. For $j=1, 2, \cdots.$ Do;

4. $\alpha_{j}=\langle\omega_{j}, A(A\nu_{j})\rangle $;

5. $\hat{\nu}_{j+1}=A\nu_{j}-\alpha_{j}\nu_{j}-\beta_{j}\nu_{j-1}$;

6. $\hat{\omega}_{j+1}=A^{T}\omega_{j}-\alpha_{j}\omega_{j}-\delta_{j}\omega_{j-1}$;

7. $\delta_{j+1}=|\langle \hat{\omega}_{j+1}, A \hat{\nu}_{j+1}\rangle|^{\frac{1}{2}}$;

8. $\beta_{j+1}=\frac{\langle\hat{\omega}_{j+1}, A \hat{\nu}_{j+1}\rangle}{\delta_{j+1}}$;

9. $\nu_{j+1}=\frac{\hat{\nu}_{j+1}}{\delta_{j+1}}$;

10. $\omega_{j+1}=\frac{\hat{\omega}_{j+1}}{\beta_{j+1}}$.

如果算法1进行$m$步, 那么产生的向量是Lanczos双$A$ -正交的, 即$\langle \nu_{i}, A \omega_{j}\rangle=\sigma_{i, j}, ~~~ 1\leq i, j \leq m. $而且, 令$ V_{m}=[\nu_{1}, \nu_{2}, \cdots, \nu_{m}]$, $ W_{m}=[\omega_{1}, \omega_{2}, \cdots, \omega_{m}]$为$n \times m$的矩阵, 那么扩充的三对角矩阵具有如下形式

$ {\bar T_m} = \left( {\begin{array}{*{20}{c}} {{T_m}}\\ {{\delta _{m + 1}}e_m^T} \end{array}} \right), $

其中

$ {T_m} = \left( \begin{array}{l} {\alpha _1}\;\;{\beta _2}\\ {\delta _2}\;\;{\alpha _2}\;\;\;{\beta _3}\\ \;\;\;\;\;\; \ddots \;\;\;\; \ddots \;\;\;\;\; \ddots \;\;\;\;\;\;\;\\ \;\;\;\;\;\;\;\;\;\;\;\;{\delta _{m - 1}}\;\;\;\;{\alpha _{m - 1}}\;\;\;\;\;{\beta _m}\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\delta _m}\;\;\;\;\;\;\;\;{\alpha _m}\\ \; \end{array} \right), $

$T_{m}$中的元素是迭代过程产生的, 其中$\alpha_{1}, \alpha_{2}, \cdots, \alpha_{m}$, $\beta_{1}, \beta_{2}, \cdots, \beta_{m}$是复的, $\delta_{1}, \delta_{2}, \cdots, \delta_{m}$是正的, Lanczos双$A$ -正交过程满足如下关系

$ AV_{m}=V_{m}T_{m}+\delta_{m+1}\nu_{m+1}e^{T}_{m}, \\ A^{T}W_{m}=W_{m}T^{T}_{m}+\beta_{m+1}\omega_{m+1}e^{T}_{m}, \\ W_{m}^{T}AV_{m}=I_{m}, \\ W_{m}^{T}A^{2}V_{m}=T_{m}. $

根据Lanczos双$A$ -正交过程有

$ AV_{m}=V_{m+1}\bar{T}_{m}, $

其中$\bar{T}_{m}$是的三对角矩阵, 如果$\nu_{1}$是由$r_{0}$的倍数定义的, 如$\nu_{1}=\beta r_{0}$, 那么原方程组具有以下形式的近似解

$ x_{m}=x_{0}+V_{m}y, $

该近似解的残差向量可以表示为

$ b-Ax_{m}\;\;=\;\;b-A(x_{0}+V_{m}y) =r_{0}-AV_{m}y\\ \;\;\;\;\;\;\;\;\;\;\;\;\;\;=\beta\nu_{1}-V_{m+1}\bar{T}_{m}y = V_{m+1}(\beta e_{1}-\bar{T}_{m}y), $

残差向量的范数为$\|b-x_{m}\|_{2}=\|V_{m+1}(\beta e_{1}-\bar{T}_{m}y)\|_{2}. $

如果$V_{m+1}$的列向量是正交的, 那么$|b-x_{m}\|_{2}=\|\beta e_{1}-\bar{T}_{m}y\|_{2}$, 但在Lanczos双$A$ -正交过程中$\nu_{i}$不是正交的, 但通过极小化函数$ F(y)=\|\beta e_{1}-\bar{T}_{m}y\|_{2}$得到$y$也是可取的, 从而得到相应的近似解为$ x_{m}=x_{0}+V_{m}y $, 为了和原始的QMR方法加以区别, 将这种算法记为QMRA.具体算法如下:

算法 2  (基于Lanczos双$A$ -正交的QMR算法)

1. 给定$x_{0}$, 计算$r_{0}=b-Ax_{0}$, 令$\beta=\|r_{0}\|_{2}$;

2. 令$\nu_{1}=r_{0}/\beta$, 选取$\omega_{1}$, 使得$\langle\omega_{1}, A\nu_{1}\rangle=1$, 比如$(\omega_{1}=\frac{A\nu_{1}}{\|A\nu_{1}\|_{2}^{2}})$;

3. 执行算法1到$m$步, 得到$ V_{m}=[\nu_{1}, \nu_{2}, \cdots, \nu_{m}]$和$ W_{m}=[\omega_{1}, \omega_{2}, \cdots, \omega_{m}]$以及三对角矩阵$\bar{T}_{m}$;

4. 对$\bar{T}_{m}$进行一系列的Givens变换$\Omega_{m}, \Omega_{m-1}, \cdots, \Omega_{1}$, 得到

$ ({\Omega _m}{\Omega _{m - 1}} \cdots {\Omega _1}){\bar T_m} = \left( {\begin{array}{*{20}{c}} {{R_{m, m}}}\\ 0 \end{array}} \right); $

5. 计算$x_{m}=x_{0}+P_{m}g_{m}$, 其中$P_{m}=V_{m}R^{-1}_{m, m}$, $g_{m}$为

$ g_{m+1, m}=(\Omega_{m}\Omega_{m-1}\cdots\Omega_{1}\beta e_{1})=(\gamma_{1}, \gamma_{2}, \cdots, \gamma_{m+1})^{T}, $

去掉最后一行所得向量;

6. 计算$r_{m}=b-Ax_{m}$, 如果$\frac{\| r_{m}\|}{\| b\|}$达到精度要求, 停止.

将QMRA算法求得的近似解$x_{m}$与$\nu_{m+1}$的线性组合$\tilde{x}_{m}=x_{m}+\theta_{m}\nu_{m+1}$作为新的近似解, 其中$\theta_{m}$是极小化问题$\min\|b-A(x_{m}+\theta_{m}\nu_{m+1})\|$的解, 如果令

$ r_{m}=b-Ax_{m}, \qquad f=A\nu_{m+1}, \qquad \tilde{r}_{m}=b-A \tilde{x}_{m}, $

那么极小化问题$\min\|b-A(x_{m}+\theta_{m}\nu_{m+1})\|$的解为$ \theta_{m}=\frac{f^{*}r_{m}}{\|f\|^{2}}, $新的残差范数为

$ \|\tilde{r}_{m}\|=\frac{\sqrt{\|r_{m}\|^{2}\|f\|^{2}-\|f^{*}r_{m}\|}}{\|f\|}. $

证明见文献[3].于是得到基于Lanczos双$A$ -正交的一种修正的QMR算法(记为MQMRA算法).具体算法如下:

算法 3  (基于Lanczos双$A$ -正交的一种修正的QMR算法)

1. 给定$x_{0}$, 计算$r_{0}=b-Ax_{0}$, 令$\beta=\|r_{0}\|_{2}$;

2. 令$\nu_{1}=r_{0}/\beta$, 选取$\omega_{1}$, 使得$\langle\omega_{1}, A\nu_{1}\rangle=1$, 比如$(\omega_{1}=\frac{A\nu_{1}}{\|A\nu_{1}\|^{2}_{2}});$

3. 执行算法1到m步, 得到$V_{m}=[\nu_{1}, \nu_{2}, \cdots, \nu_{m}]$和$W_{m}=[\omega_{1}, \omega_{2}, \cdots, \omega_{m}]$以及三对角矩阵$\bar{T}_{m}$;

4. 执行算法2得到$x_{m}$, 及$r_{m}=b-Ax_{m}$;

5. 计算$f=A\nu_{m+1}$, $\theta_{m}=\frac{f^{*}r_{m}}{\|f\|^{2}}$;

6. 计算$\tilde{x}_{m}=x_{m}+\theta_{m}\nu_{m+1}$;

7. 计算$\tilde{r}_{m}=b-A\tilde{x}_{m}$, 如果$\frac{\|\tilde{r}_{m}\|}{\| b\|}$达到精度要求, 停止.

下面给出四个数值算例来比较原始的QMR算法, 基于Lanczos双$A$ -正交的QMR算法(QMRA), 以及修正的QMRA算法(MQMRA)的执行结果.为了使数值结果便于重复, 这里具体给出各参数的值, 而不考虑算法的执行细节.选取右端向量$b$使得方程组$ Ax=b $的精确解为$x^{T}= \left( \begin{array}{cccc} 1, 1, \cdots, 1 \end{array} \right) $, 预初始估计值为$x_{0}^{T}= \left( \begin{array}{cccc} 0, 0, \cdots, 0 \end{array} \right), $图中的残差范数均取10为底的对数.

算例 1  选取系数矩阵为Grcar矩阵, 如下:

$ A = {\left( {\begin{array}{*{20}{c}} 1&1&1&1&{}&{}\\ { - 1}&1&1& \ddots&\ddots &{}\\ {}&{ - 1}&1& \ddots&\ddots &1\\ {}&{}& \ddots&\ddots&\ddots &1\\ {}&{}&{}& \ddots&\ddots &1\\ {}&{}&{}&{}&{ - 1}&1 \end{array}} \right)_{n \times n}}. $

当矩阵阶数$n=1500$, 求解精度为$Tol=1e-8$时, 三种方法得到迭代收敛效果如图 1, 由于该矩阵伪谱的凸包中含有零点, 在这种情况下GMRES或其他一些混合算法将导致失败, 但图中可以看出新算法具有满意的收敛效果.

算例 2  本算例中的系数矩阵取自文献[3], 具有如下形式:

$ A = \left( {\begin{array}{*{20}{c}} 1&{}&{}&{}&\alpha \\ {}&2&{}&{}&{}\\ {}&{}& \ddots &{}&{}\\ {}&{}&{}&{n - 1}&{}\\ {}&{}&{}&{}&n \end{array}} \right). $

取求解精度${\rm Tol}=1e-10$, 当矩阵阶数$n=2000, \alpha=1.1$时, $A$接近对称矩阵, 三种方法得到迭代收敛效果如图 2.当矩阵阶数$n=2000, \alpha=20000$时, $A$是典型的非对称矩阵, 三种方法得到迭代收敛效果如图 3.

算例 3  本例中矩阵PDE2961和矩阵CDDE1, 取自矩阵市场(http: //math. Nist. Gov./ Marix Market/), PDE2961条件数9.49e+02, 非零元个数为14585, 阶数2961;矩阵CDDE1条件数4.1e+03, 非零元素个数为4681, 阶数961, 当求解精度为${\rm Tol}=1e-8$时, 三种算法的收敛效果如图 4和图 5.

算例 4^[6]  考虑如下二阶椭圆形偏微分方程

$ \begin{eqnarray*} \left \{ \begin{aligned}{cc} &- \triangle u + 2 p_{1}u_{x}+ 2 p_{2}u_{y} + p_{3}u =f,&u\in S, \\ &u = 0,&u\in \partial S, \end{aligned} \right. \end{eqnarray*} $

其中$S=\{(x, y)\in R^{2}, 0<x, y<1\}$, $\partial S $表示$S$的边界, $f$是定义在$S$上的函数.

用五点中心差分格式离散该微分方程, 把求解区域区$S$划分成网格尺寸为$ h=1/l+1 $的$ (l+2)\times(l+2)$个均匀正方形网格, 网格点为$ s_{i}=ih, t_{j}=jh, 0 \leq i, j \leq l+1 $.把网格点按自然次序排序得$l^{2} \times l^{2}$块三对角矩阵

$ A = {\left( {\begin{array}{*{20}{c}} T&{(\beta - 1){I_l}}&{}&{}&{}\\ { - (\beta + 1){I_l}}&T&{(\beta - 1){I_l}}&{}&{}\\ {}&{ - (\beta + 1){I_l}}& \ddots&\ddots &{}\\ {}&{}& \ddots &T&{(\beta - 1){I_l}}\\ {}&{}&{}&{ - (\beta + 1){I_l}}&T \end{array}} \right)_{{l^2} \times {l^2}}}, $

其中

$ T = {\left( {\begin{array}{*{20}{c}} {4 - \sigma }&{\gamma - 1}&0&{}&{}\\ { - \gamma - 1}&{4 - \sigma }&{\gamma - 1}&{}&{}\\ {}&{ - \gamma - 1}&{4 - \sigma }& \ddots &{}\\ {}&{}& \ddots&\ddots &{\gamma - 1}\\ {}&{}&{}&{ - \gamma - 1}&{4 - \sigma } \end{array}} \right)_{l \times l}}, $

$\gamma=p_{1}h, \beta=p_{2}h, \sigma=p_{3}h^{2}$, 线性系统右端的元素为$h^{2}f(s_{i}, t_{j})$, 为简单起见以下算例右端向量取$b={\rm sum}(A, 2)$, 取$l=50, p_{1}=25, p_{2}=50, p_{3}=30$得到的结果如图 6.

以下两个表格给出了三种算法迭代次数及CPU时间比较, 总体来看, MQMRA算法的运算迭代次数和所用时间都比QMR算法要少, 有的情况下甚至要快很多, 从而加快了收敛速度.

本文将QMR算法中的Lanczos双正交过程用Lanczos双$A$ -正交过程代替, 得到了QMRA算法, 进一步将由该算法得到的近似解与最后一个基向量的线性组合来作为新的近似解, 使新近似解的残差范数满足一个一维极小化问题, 从而得到一种基于Lanczos双$A$ -正交的修正的QMR算法.从数值算例的结果来看, 对于某些大型稀疏矩阵新算法的收敛速度加快、迭代次数减少, 运算效率提高.该方法也可以考虑和预条件技术相结合来进一步提高计算效率.