解线性方程组的VRP-GMRES（m）迭代法

白雪婷，杨琴乐

应用科学、物理和工程领域许多问题可建立偏微分方程的数学模型，这些微分方程经过离散后一般都可归结为求解大型稀疏线性方程组Au=b，A∈Rn×n，u，b∈Rn，其中A为非奇异矩阵.这类方程组的求解通常采用迭代法.目前，以Galerkin 原理为基础建立的广义极小残余GMRES（m）迭代法是求解此类方程组最有效的方法之一.因此，近年来GMRES（m）算法及其在其他领域的应用一直是人们研究的热点［1-4］.实际计算表明，当系数矩阵A是良态矩阵时，GMRES（m）算法及一些简单改进后的GMRES（m）算法便可有效地求解这些线性系统［5-6］，但如果系数矩阵A具有很强的病态性，则必须结合一些特定的预处理技术［7-10］.

不管是GMRES（m）算法还是预处理GMRES（m）算法，在实际执行这些算法时，参数m一旦选定，在之后整个迭代过程当中将始终固定不变.所以参数m的选择也是影响算法有效执行的关键因素之一.研究表明，m取值较小时，有可能出现收敛慢甚至不收敛等现象，m选择过大会造成存储空间过多的需求.因此，如何选择一个合适的参数m一直困扰着众多学者，最近有不少学者试图通过变化GMRES（m）算法中的重启动参数m来克服这些困难.BAKE A H 最初通过选择随机的参数m来改善GMRES（m）算法的收敛性，后来又根据两次迭代所得残余向量的范数之比来确定参数m的大小，但这也并不能保证算法的收敛性［11-12］，PEAIRS L 等则通过强化学习方法来决定参数m的取值，该方法进一步表明适当的参数m可以提高GMRES（m）算法的计算效率，但由于强化学习只是一种机器学习方法，因此存在一定局限性［13］.

1 GMRES（m）算法基本理论

1.1 Galerkin 原理

任意给定向量u0∈Rn，令u=u0+z，则方程组Au=b等价于

其中：r0=b-Au0是初始残余向量.给定一个参数m＞0，设Rn中两个m维子空间Km和Lm分别为：

且v1,v2,…,vm和w1,w2,…,wm分别为子空间Km和Lm的基.令求解式（1）的Galerkin 原理为：给定一个m＞0，寻找近似解zm=Vmym∈Km，使(r0-Azm)⊥Lm，即(r0-AVmym)⊥Lm，该 式 也 可 表 示 为(r0-AVmym,ωi)= 0,ym∈Rm.

1.2 Arnoldi 过程及其性质

Arnoldi 过程具体如下（文中若不作特殊说明，‖ · ‖均表示向量的二范数）：

①对m＞0 和初始向量v0∈Rn.计算v1=r0‖r0‖，其中r0=b-Au0；

②对于k= 1,2,…,m，计算hi,k= (Avk,vi)，要 求

③如果hk+1,k= 0，停止；否则vk+1=

Arnoldi 过程有如下重要性质.

定 理1［14］在Arnoldi 过 程 中，令Vm+1=则如下关系成立.

其中：Hm= (hij)m是上海森伯格矩阵，=(0,0,…,1) ∈Rm且I是一个m+ 1 阶单位矩阵.

1.3 GMRES 算法

定理2［15］令A为n×n方阵并假定L=AK，u=u0+z=u0+Vmy，则由Galerkin 原理得到的近似解Zm将使残余向量‖r0-Az‖在Krylov子空间Km中达到最小，反之亦然，并且有

其中：β= ‖r0‖，e1= ( 1,0,…,0)T∈ Rm+1.‖r0-Az‖2等价于在Rm中极小化

这个定理表明，在Km中极小化残余向量GMRES（m）算法基本思想就是固定重启动参数m，通过求解最小二乘问题迭代计算求得zm，使得‖rm‖= ‖r0-Azm‖＜τ，GMRES（m）算法详细步骤如下：

①给定误差上界τ＞0，给定一个合适的m(m＜n)，取u0∈Rn，计 算r0=b-Au0以及β= ‖r0‖；

④计 算zm=Vmym.若‖rm‖= ‖r0-Azm‖＜τ，迭代终止；否则，令um=u0+zm，u0=um，转向步骤②.

在GMRES（m）算法中，并非任意选取m都能使其收敛.事实上，适当地变化m可以有效地改善GMRES（m）的收敛性，且在一定程度上可以避免收敛慢甚至不收敛等现象.

2 VRP-GMRES（m）算法

2.1 VRP-GMRES（m）算法

适当变化重启动参数m，通过迭代使得残余向量‖rm‖= ‖r0-Azm‖＜τ，具体过程为：

Step1：选择一个较小的初始正整数m(m≪n)，给定误差上界τ＞0；

Step3：求解如下最小二乘问题

可得ym，计算近似解zm=Vmym，且令um=u0+zm；

Step4：如 果‖rm‖= ‖r0-Azm‖＜τ，迭 代停止，输出u=um.否则，转向Step5；

Step5：令u0=um，m=m+ 1(m≤n)转至Step2.

实际执行VRP-GMRES（m）算法时，可能出现收敛慢甚至不收敛等现象（虽然这种情况较少发生）.所以，为防止参数m过大而造成存储空间过多的需求，可以先执行若干次VRP-GMRES（m）算法，待m增加到一定程度且残量范数减小到某一适当值时，再固定m循环迭代，即执行GMRES（m）算法.

2.2 收敛分析

定义1［15］设

其中：实数c和s满足，称Gi为平面旋转矩阵，也称为Givens（吉文斯）变换矩阵.

定 理3 在VRP-GMRES（m）算 法 中，令

证明 在VRP-GMRES（m）算法中，令式（5）中Gi∈R(m+1)×(m+1)，

其中：i= 1,2,…,m，经i-1 次Givens 正交变换后对角线及次对角线元素.

令Qm=GmGm-1…G1，对 式（4）中和βe1作正交变换，则可得

其中：Rm是m阶可逆上三角方阵.

解Rmy=gm，得ym.此时

当参数m变为m+ 1 时，由Arnoldi 过 程，有

令式（5）中Gi∈R(m+2)×(m+2)，且

其中：i= 1,2,…,m+ 1，取Qm+1=Gm+1Gm…G1，对式（4）中和βe1作正交变换，则

其中：Rm+1是m+ 1 阶可逆上三角方阵，

其中：e1= ( 1,0,…,0)T∈ Rm+2,gm+1∈ Rm+1.

于是

解Rm+1y=gm+1，得ym+1.此时

由定理3 知，VRP-GMRES（m）算法不仅收敛，而且比GMRES（m）算法计算精度更高.

3 数值实验

本节通过两个不同类型的数值算例说明VRP-GMRES（m）算法的有效性和可行性，分析比较GMRES（m）和VRP-GMRES（m）两种算法的数值结果.在以下算例中均选取u0=作 为 初 始 向 量，‖rm‖=‖r0-Azm‖＜τ= 1e- 8 作为停止迭代标准.算例1 考虑如下一维波动方程

该方程精确解为：u(x,t)= sin( πx)cos( 2πt)+ sin( 2πx)cos( 4πt)，其中：u(x,t)表示振弦，是对特定位置x和特定时间t的波强度的一个测量，如图1（a）所示.用中心差分格式离散后可得

其中：n表示x方向和t方向网格数.用VRPGMRES（m）算法求解式（8），计算结果如图1（b）所示.由图1 知，式（7）的数值解和精确解是一致的，这也说明VRP-GMRES（m）算法的正确性.

分别用GMRES（m）算法和VRP-GMRES（m）算法求解式（8），当n= 10，m= 7时，cond(A)=56.507 9.GMRES（m）算法循环迭代9 次后出现迭代停滞，此时残量范数为1.409 9，但对于VRP-GMRES（m）算法，迭代3 次后残量范数已达到1.423 0e-014.迭代过程各节点绝对误差大小的比较如图2 所示.

图2 迭代过程中各节点绝对误差大小的比较

由图2 可知，当迭代次数相同时，用VRPGMRES（m）算法比用GMRES（m）算法得到的计算结果的绝对误差要小得多，而且随着迭代次数的增加，用VRP-GMRES（m）算法计算产生的绝对误差减少得更快.这说明本文算法不仅具有更高的计算精度，而且能快速收敛.

当参数m分别取不同值时，对两种算法的计算结果进行比较如表1～表4 所示.表1～表4 分别给出两种算法的迭代次数、运行时间、计算精度和收敛速度比较.从表1～表4 可知，文中提出的算法收敛快，稳定性好，计算精度高，运行时间短.参数m的选择对GMRES（m）算法至关重要，参数m选择不当会导致算法失败，而参数m的选择对VRP-GMRES（m）算法的执行基本无任何影响.这表明，GMRES（m）算法对参数m的选择相当敏感，m的微小变化可能会导致GMRES（m）算法出现迭代停滞等现象，而本文提出的VRP-GMRES（m）算法中参数m是变化的，在一定程度上降低了GMRES（m）算法对参数m的敏感度，这也是VRP-GMRES（m）算法的优点之一.

表1 两种算法迭代次数比较

表2 两种算法运行时间比较

表3 两种算法计算精度比较

表4 两种算法收敛性比较

事实上，当网格数n增大时，系数矩阵A的条件数也不断增大，当n=59时，A的条件数已达到2.109 8e+003，方程组（8）呈现严重的病态性.由表1 知，此时即使参数m取值很大，GMRES（m）算法的收敛速度也不是很快，而且在达到给定误差时，GMRES（m）算法的迭代次数是VRP-GMRES（m）算法的11.7 倍，运行时间是其10.7 倍.这表明，随着计算规模地增大，VRP-GMRES（m）算法比GMRES（m）算法更有效，有更广阔的前景.

算例2 考虑如下二维泊松方程

该方程精确解为u(x,y)=x2(x+y2)+ 2，其中：u(x,y)是温度分布函数，温度分布如图3（a）所示.将式（9）所示的求解区域分别沿x方向和y方向n等分，可得

取n=40，m=10，用VRP-GMRES（m）求解式（10）时，各节点处的温度分布如图3（b）所示.由图3 可知，数值解和精确解是一致的.

图3 温度分布（n=40，m=10）

当m=10 时，随着计算规模不断增大，VRP-GMRES（m）算法所需迭代次数变化较小且所需的迭代次数要比GMRES（m）少很多，如图4（a）所示，说明VRP-GMRES（m）算法有更高的计算效率，同时，VRP-GMRES（m）算法还有更高的精度，如图4（b）所示.另外，实际执行VRP-GMRES（m）算法时，虽然每次迭代所耗时可能比GMRES（m）算法长，但是由于VRP-GMRES（m）算法收敛速度很快，所以当残量范数达到给定误差时，如表5 所示，本文提出的VRP-GMRES（m）算法总的运行时间要比GMRES（m）算法少很多.

图4 不同网格下两种算法计算效率和计算精度比较（m=10）

表5 不同网格下两种算法运行时间比较

4 结语

文中提出一种求解线性方程组的VRPGMRES（m）算法，理论证明新算法收敛且计算精度更高.数值试验结果表明本文提出的VRP-GMRES（m）算法不仅有更高的计算效率和计算精度，而且有更好的稳定性.新算法适合用来求解由科学研究和工程问题得到的大型线性方程组.