0 引言

在科学和工程领域中, 很多问题最终都可归结为大型线性方程组的求解, 文献[1]提出的基于Galerkin原理的GMRES(m)算法是公认的求解大型非对称线性方程组较为成功的方法之一, 同时也是快速多极边界元法中的求解器。

随着GMRES算法和相关知识被广泛学习, 很多学者也对该算法进行了优化和改进。一种方法是采用预处理技术^[2-4]。刘晓明等^[5]将GMRES方法应用于油藏数值模拟计算问题, 利用矩阵分块技术, 采用块拟消去法(PE)对系数矩阵进行了预处理。于春肖等^[6]提出ILU预处理Householder-GMRES(m)算法, 通过降低系数矩阵的条件数达到加快收敛的目的。关朋燕等^[7]将多项式预处理与加权GMRES(m)相结合, 构造出一种新的方法, 减少了计算量。于春肖等^[8]通过改变重启动参数m的值提出VRP-HGMRES(m)算法, 该算法有效地避免了因重启动参数的选取不当而造成的问题。文献[9]将稀疏矩阵的小特征值对应的特征向量加入Krylov子空间。于春肖等^[10]提出一种截断型IGMRES(m)新算法, 此算法是基于FMM的Krylov子空间投影类方法。郝雪景等^[11]提出基于不完全正交的变参数IGMRES(m)算法。Walker^[12]对GMRES算法进行改进, 提出一种基于Householder变换的GMRES(简称H-GMRES)算法。

本文将利用截断技术^[13-14]给出一种不完全正交的变参数H-IGMRES(m)算法, 使计算量大大减少, 并结合数值算例分析该算法的高效性。

1 H-GMRES(m)算法

在实际计算中, H-GMRES(m)算法构造Hessenberg矩阵元素和Krylov向量时都要用到较长的递推式, 从而导致消耗的计算时间很长, 迭代过程数值也不稳定。这种算法的具体步骤为:

1) 给定重启动参数m和误差上界ε>0, 取初始值x⁽⁰⁾∈Rⁿ, 计算r⁽⁰⁾=b－Ax⁽⁰⁾, β=‖r⁽⁰⁾‖₂。

2) 计算P₁=I－2ωω^Τ, 其中：‖ω‖₂=1；P₁r⁽⁰⁾=βe₁；v₁=P₁e₁。

3) 对于j=1, 2, …, m，令c_j=P^ΤAv_j(当j=1时, 有P=P₁成立), 然后计算

$ {\mathit{\boldsymbol{P}}_{j + 1}} = \mathit{\boldsymbol{I}} - 2{\mathit{\boldsymbol{c}}_j}\mathit{\boldsymbol{c}}_j^{\rm{T}}/\left\| {{\mathit{\boldsymbol{c}}_j}} \right\|_2^2,\mathit{\boldsymbol{P}} = {\mathit{\boldsymbol{P}}_1}{\mathit{\boldsymbol{P}}_2} \cdots {\mathit{\boldsymbol{P}}_{j + 1}}, $

$ {\mathit{\boldsymbol{v}}_{j + 1}} = {\mathit{\boldsymbol{P}}_{j + 1}}{\mathit{\boldsymbol{P}}_j} \cdots {\mathit{\boldsymbol{P}}_1}{\mathit{\boldsymbol{e}}_{j + 1}},{\mathit{\boldsymbol{h}}_j} = {\mathit{\boldsymbol{P}}_{j + 1}}{\mathit{\boldsymbol{c}}_{j}}。$

形成H_m和V_m, 从而变成求解最小二乘问题y_m=arg min‖βe₁－H_my‖₂，y∈R^m。得出原方程组的近似解x^(m)=x⁽⁰⁾+z_m=x⁽⁰⁾+V_my_m。

4) 如果‖r⁽⁰⁾－Az_m‖₂ < ε，则x=x^(m), 迭代停止；否则x⁽⁰⁾=x^(m), m=m+1且转向1)再一次迭代, 直到满足给定的误差上界为止。

2 H-IGMRES(m)算法与收敛性分析 2.1 H-IGMRES(m)算法

当m很大时, H-GMRES(m)算法需要计算和保存所有的{v_i}_i=1^m, 这将需要很大的存储空间。当问题规模很大时, 有可能引起内存不足。基于H-GMRES(m)算法, 利用截断技术, 本文提出H-IGMRES(m)算法, 其基本思想是在采用Householder变换构造Krylov子空间的基向量时, 仅使用部分而非所有前面计算出的向量来构造新的递推式以计算后面的向量, 从而减少算法的计算量, 即仅对前面的至多k个向量v_j0, …, v_j, j₀=max{1, j－k－1}进行正交的Householder变换, 其中k < m, 定义截断指标q₀=j₀=max{1, j－k－1}, 截断比为(m－q₀)/m。具体步骤如下:

1) 初始化:给定重启动参数m和误差上界ε>0, 设置截断指标q₀, 取初始值x⁽⁰⁾∈Rⁿ, 计算

$ {\mathit{\boldsymbol{r}}^{\left( 0 \right)}} = \mathit{\boldsymbol{b}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{x}}^{\left( 0 \right)}},\beta = {\left\| {{\mathit{\boldsymbol{r}}^{\left( 0 \right)}}} \right\|_2}。$

2) 计算P₁=I－2ωω^Τ, 其中‖ω‖₂=1；P₁r⁽⁰⁾=βe₁；v₁=P₁e₁。

3) 对于j=1, 2, …, m有不完全正交化, 令c_j=P^ΤAv_j(当j=1时, 有P=P₁成立), 然后计算

$ {\mathit{\boldsymbol{P}}_{j + 1}} = \mathit{\boldsymbol{I}} - 2{\mathit{\boldsymbol{c}}_j}\mathit{\boldsymbol{c}}_j^{\rm{T}}/\left\| {{\mathit{\boldsymbol{c}}_j}} \right\|_2^2,\mathit{\boldsymbol{P}} = {\mathit{\boldsymbol{P}}_{{q_0}}} \cdots {\mathit{\boldsymbol{P}}_{j + 1}}, $

$ {\mathit{\boldsymbol{v}}_{j + 1}} = \mathit{\boldsymbol{P}}{\mathit{\boldsymbol{e}}_{j + 1}},{\mathit{\boldsymbol{h}}_j} = {\mathit{\boldsymbol{P}}_{j + 1}}{\mathit{\boldsymbol{c}}_j}。$

更新V_j+1与H_j，V_j+1=(V_j, v_j+1)，H_j=(H_j－1, h_j)_(j+1)×j。

H_j为(j+1)×j的带状上Hessenberg矩阵，其非零元素由上式生成。当j=1时第一列省略, 并且有

$ \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{V}}_m} = {\mathit{\boldsymbol{V}}_{m + 1}}{{\mathit{\boldsymbol{\bar H}}}_m}。$

求解最小二乘问题

$ {\mathit{\boldsymbol{y}}_m} = \arg \min{\left\| {\beta {\mathit{\boldsymbol{e}}_1} - {{\mathit{\boldsymbol{\bar H}}}_m}\mathit{\boldsymbol{y}}} \right\|_2},\mathit{\boldsymbol{y}} \in {{\bf{R}}^m}, $

得出原方程组的近似解x^(m)=x⁽⁰⁾+z_m=x⁽⁰⁾+V_my_m。

4) 如果‖r⁽⁰⁾－Az_m‖₂ < ε，则x=x^(m), 迭代停止；否则x⁽⁰⁾=x^(m), m=m+1且转向1)再一次迭代, 直到满足给定的误差上界为止。

2.2 收敛性分析

为研究这种截断型H-IGMRES(m)算法的收敛形态, 这里将H-IGMRES(m)算法和H-GMRES(m)算法结合起来, 以建立H-IGMRES(m)算法的收敛性理论。

定义V_m+1=(v₁, v₂, …, v_m+1), 可以证明, 如果H-IGMRES(m)算法在m步中断且V_m+1列满秩, 则得到的数值解为准确解x^*。另外, 当q₀=1时, H-IGMRES(m)算法与H-GMRES(m)算法相同。为讨论方便, 选取K_m=span{r⁽⁰⁾, Ar⁽⁰⁾, …, A^m－1r⁽⁰⁾}, 将H-IGMRES(m)算法进行m步迭代后的数值解表示为x^(m)(HI), 残余向量表示为r^(m)(HI); 将H-GMRES(m)算法进行m步迭代后的数值解表示为x^(m)(H), 残余向量表示为r^(m)(H)。

定理1  假设V_m+1列满秩，则H-IGMRES(m)算法和H-GMRES(m)算法在K_m上的残值范数满足

$ \left\| {{\mathit{\boldsymbol{r}}^{\left( m \right)}}\left( {\mathit{\boldsymbol{HI}}} \right)} \right\| \le S\left( {{\mathit{\boldsymbol{V}}_{m + 1}}} \right)\left\| {{\mathit{\boldsymbol{r}}^{\left( m \right)}}\left( \mathit{\boldsymbol{H}} \right)} \right\|。$

(1)

其中V_m+1由不完全正交化过程生成, S(V_m+1)=‖V_m+1‖‖V_m+1⁺‖, V_m+1⁺是V_m+1的广义逆矩阵。如果H-IGMRES(m)算法在m步中断, 即h_{m+1, m}=0, 则x^(m)(HI)=x^*。

证明  由r⁽⁰⁾=βV_m+1e₁和AV_m=V_m+1H_m可知, 对应于x^(m)(HI)的残值为

$ {\mathit{\boldsymbol{r}}^{\left( m \right)}}\left( {\mathit{\boldsymbol{HI}}} \right) = \mathit{\boldsymbol{b}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{x}}^{\left( m \right)}}\left( {\mathit{\boldsymbol{HI}}} \right) = {\mathit{\boldsymbol{V}}_{m + 1}}\left( {\beta {\mathit{\boldsymbol{e}}_1} - {{\mathit{\boldsymbol{\bar H}}}_m}{\mathit{\boldsymbol{y}}_m}} \right)。$

因此

$ \left\| {{\mathit{\boldsymbol{r}}^{\left( m \right)}}\left( {\mathit{\boldsymbol{HI}}} \right)} \right\| = \left\| {{\mathit{\boldsymbol{V}}_{m + 1}}\left( {\beta {\mathit{\boldsymbol{e}}_1} - {{\mathit{\boldsymbol{\bar H}}}_m}{\mathit{\boldsymbol{y}}_m}} \right)} \right\|。$

(2)

由r⁽⁰⁾=βV_m+1e₁, 可得V_m+1⁺r⁽⁰⁾=βe₁, 再由AV_m=V_m+1H_m可得H_m=V_m+1⁺AV_m。

因此, 式(2)变成

$ \left\| {{\mathit{\boldsymbol{r}}^{\left( m \right)}}\left( {\mathit{\boldsymbol{HI}}} \right)} \right\| = \left\| {{\mathit{\boldsymbol{V}}_{m + 1}}\left( {{\mathit{\boldsymbol{V}}_{m + 1}}^ + {\mathit{\boldsymbol{r}}^{\left( 0 \right)}} - {\mathit{\boldsymbol{V}}_{m + 1}}^ + \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{V}}_m}{\mathit{\boldsymbol{y}}_m}} \right)} \right\|。$

式$\left\| \boldsymbol{r}_{m}\right\|=\min \limits_{ \boldsymbol{y}_{m} \in \mathbf{R}^{m}}\left\|\beta \boldsymbol{e}_{1}- \boldsymbol{\bar{H}}_{m} \boldsymbol{y}_{m}\right\|$等价于求解最小二乘问题，

$ \mathop {\min }\limits_{{y_m} \in {{\bf{R}}^m}} \left\| {{\mathit{\boldsymbol{V}}_{m + 1}}^ + {\mathit{\boldsymbol{r}}^{\left( 0 \right)}} - {\mathit{\boldsymbol{V}}_{m + 1}}^ + \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{V}}_m}{\mathit{\boldsymbol{y}}_m}} \right\| = \mathop {\min }\limits_{{z_m} \in {K_m}} \left\| {{\mathit{\boldsymbol{V}}_{m + 1}}^ + {\mathit{\boldsymbol{r}}^{\left( 0 \right)}} - {\mathit{\boldsymbol{V}}_{m + 1}}^ + \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{z}}_m}} \right\|。$

因此可得

$ \begin{array}{l} \left\| {{\mathit{\boldsymbol{r}}^{\left( m \right)}}\left( {\mathit{\boldsymbol{HI}}} \right)} \right\| \le \left\| {{\mathit{\boldsymbol{V}}_{m + 1}}} \right\|\left\| {\left( {{\mathit{\boldsymbol{V}}_{m + 1}}^ + {\mathit{\boldsymbol{r}}^{\left( 0 \right)}} - {\mathit{\boldsymbol{V}}_{m + 1}}^ + \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{V}}_m}{\mathit{\boldsymbol{y}}_m}} \right)} \right\| = \\ \left\| {{\mathit{\boldsymbol{V}}_{m + 1}}} \right\|\mathop {\min }\limits_{{z_m} \in {K_m}} \left\| {{\mathit{\boldsymbol{V}}_{m + 1}}^ + {\mathit{\boldsymbol{r}}^{\left( 0 \right)}} - {\mathit{\boldsymbol{V}}_{m + 1}}^ + \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{z}}_m}} \right\| \le \\ \left\| {{\mathit{\boldsymbol{V}}_{m + 1}}} \right\|\left\| {{\mathit{\boldsymbol{V}}_{m + 1}}^ + } \right\|\mathop {\min }\limits_{{z_m} \in {K_m}} \left\| {{\mathit{\boldsymbol{r}}^{\left( 0 \right)}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{z}}_m}} \right\| = \\ S\left( {{\mathit{\boldsymbol{V}}_{m + 1}}} \right)\mathop {\min }\limits_{{z_m} \in {K_m}} \left\| {{\mathit{\boldsymbol{r}}^{\left( 0 \right)}} - \mathit{\boldsymbol{A}}{\mathit{\boldsymbol{z}}_m}} \right\| = S\left( {{\mathit{\boldsymbol{V}}_{m + 1}}} \right)\left\| {{\mathit{\boldsymbol{r}}^{\left( m \right)}}\left( \mathit{\boldsymbol{H}} \right)} \right\|, \end{array} $

所以式(1)成立。

如果h_{m+1, m}=0, 则AV_m=V_m+1H_m变为AV_m=V_mH_m, 其中H_m是去掉H_m最后一行的m×m带状上Hessenberg矩阵, V_m的列张成A的一个不变子空间。所以, H_m的特征值均为A的特征值, 从而H_m非奇异。现在‖r^(m)(HI)‖=‖r⁽⁰⁾－AV_my_m‖=‖r⁽⁰⁾－V_mH_my_m‖=‖V_m(βe₁－H_my_m)‖, 而当y_m=βH_m^－1e₁时, 有$\min \limits_{ \boldsymbol{y}_{m} \in \mathbf{R}^{m}}\left\|\beta \boldsymbol{e}_{1}- \boldsymbol{\bar{H}}_{m} \boldsymbol{y}_{m}\right\|=0$。因此‖r^(m)(HI)‖=0, 所以x^(m)(HI)=x^*。

3 数值算例

例1   考虑两点边值问题$\left\{\begin{array}{l} {p \frac{\mathrm{d}^{2} y}{\mathrm{d} x^{2}}+\frac{\mathrm{d} y}{\mathrm{d} x}=q, 0<q<1} \\ {y_{0}=0, y_{1}=1} \end{array}\right.。$

此问题的精确解为$y=\frac{1-q}{1-{\rm e}^{-\frac{1}{p}}}\left(1-{\rm e}^{-\frac{x}{p}}\right)+q x$。为了求解两点边值问题, 首先将微分法过程离散化, 即将x坐标轴上[0, 1]区间n等分, 并令h=1/n, 由差分法可得

$ Ay = b,A \in {{\bf{R}}^{{{\left( {n - 1} \right)}^2} \times {{\left( {n - 1} \right)}^2}}},y,b \in {{\bf{R}}^{{{\left( {n - 1} \right)}^2}}}。$

如果选取p=0.01, q=0.5, n=40, 初始向量为y⁽⁰⁾=(0, …, 0)^Τ, 误差上界为t=1e-6, 当m=10时, 取截断指标为9, 用H-IGMRES(m)算法求解该线性方程组, 将得到的数值解与精确解进行比较, 如图 1所示。从图 1可看出, 数值解和精确解完全吻合, 证明此算法的可行性。

下面选取n=100, 分别用H-IGMRES(m)算法和H-GMRES(m)算法求解这个线性方程组。当重启动参数不同时, 两种算法所需的迭代次数和得到的残值范数如图 2所示, 两种算法所需的运行时间见表 1。从图 2和表 1可看出, 重启动参数不同时, H-IGMRES(m)算法与H-GMRES(m)算法相比, 前者迭代次数和运行时间较少, 说明新算法提高了计算效率, 并且具有更高的计算精度。

例2   考虑如下二维Possion方程, 其中Ω{(x, y) 0≤x≤1, 0≤y≤1}，u(x, y)为温度分布函数。

$ \left\{ \begin{array}{l} - \left( {{u_{xx}} + {u_{yy}}} \right) = 2{{\rm{ \mathsf{ π} }}^2}\sin {\rm{ \mathsf{ π} }}x\sin {\rm{ \mathsf{ π} }}y,\left( {x,y} \right) \in \mathit{\Omega }\\ u\left( {x,y} \right) = 0,\left( {x,y} \right) \in \partial \mathit{\Omega } \end{array} \right.。$

令f(x, y)=2π²sin πxsin πy, 将求解区域沿x方向和y方向进行等分, 由五点差分法可得Au=b, A∈R^n2×n2, u, b∈R^n2×1。

当n=35时, 选取重启动参数m=20及截断指标q₀=9, 分别用高斯消元法和H-IGMRES(m)算法求解线性方程组, 其结果分别如图 3(a)和3(b)所示。

由图 3可以看出, 两种算法的计算结果完全一致, 这表明H-IGMRES(m)算法是正确的。截断指标可以在1~19中取值, 不同的截断指标对计算速率、计算精度、运行时间的影响如图 4所示。

由图 4可看出, 迭代次数和残值范数随截断指标的增大呈下降趋势, 运行时间也在突降之后逐渐趋于平稳, 说明截断技术在保证计算精度的前提下提高了计算效率。当m=10时, 随着计算规模n的不断增大, 当残值范数满足给定的误差上界时, H-IGMRES(m)算法总的运行时间要比H-GMRES(m)算法少很多, 进一步说明了H-IGMRES(m)算法的优越性, 见表 2。

4 结论

1) 提出了H-IGMRES(m)算法, 对算法的收敛性进行了理论分析。

2) 数值算例表明H-IGMRES(m)算法的计算结果和精确解吻合, 说明该算法的有效性。

3) 分析了截断指标的选取对算法的计算精度和效率的影响, 证明了H-IGMRES(m)算法的计算效率和精度都高于H-GMRES(m)算法。