劉仲云，李莉

(長(zhǎng)沙理工大學(xué) 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，湖南 長(zhǎng)沙，410114)

文中考慮大型稀疏線性方程組

Ax=b

(1)

的迭代解法，其中x，b∈Rn，A∈Rn×n為含Toeplitz塊的全正塊三對(duì)角隨機(jī)Toeplitz矩陣，即A有如下結(jié)構(gòu)：

其中:Ti(i=-1，0，1)為Toeplitz矩陣[1]。

一個(gè)非奇異矩陣稱為全正矩陣[2]，如果它所有的子式都是非負(fù)的。

一個(gè)矩陣是一個(gè)隨機(jī)矩陣，如果它是非負(fù)矩陣且每一行的行和都為1。

這類線性方程組頻繁出現(xiàn)在三次均勻B樣條曲面擬合及二階微分方程數(shù)值解等問題中。

1 精化迭代法

求解(1)式的精化迭代法的迭代格式為

x(k+1)=x(k)+d(k)

(2)

Ad(k)=r(k)，r(k)=b-Ax(k)

(3)

其中:x(0)為初始迭代向量，k=0，1，…

當(dāng)d(k)為(3)的精確解時(shí)，迭代格式(2-3)退化為單步迭代精化。不難看出，該迭代格式可以提高解x(k)的精度[3]。另外，倘若用高精度方法求解(3)，則迭代格式(2-3)為迭代精化[4-5]，把該迭代格式稱為精化迭代法。

引理1 假設(shè)不考慮舍入誤差，精確計(jì)算殘差r(k)，d(k)為(3)的精確解，則x(k+1)為(1)的精確解。

證明 由(3)式知Ad(k)=r(k)，r(k)=b-Ax(k)，從而Ax(k+1)=Ax(k)+Ad(k)=b。這意味著x(k+1)是(1)的精確解，證畢。

由于A的正定性，因此，用廣義極小殘差法(GMRES)非精確求解(3)，這樣就得到精化迭代法的基本框架。

設(shè)d(k，l)為用GMRES方法求解(3)得到的第l個(gè)近似解，其中，初始迭代向量d(k，0)=0，終止檢驗(yàn)公式為

(4)

其中：τk是終止判別。

倘若d(k，Lk)滿足終止檢驗(yàn)公式(4)，從而

x(k+1)=x(k)+d(k，Lk)

(5)

通常，為了確保計(jì)算效率，(5)的終止檢驗(yàn)公式為

‖r(k+1)‖≥θ‖r(k)‖，θ∈[0，1]。

(6)

為方便算法實(shí)現(xiàn)，將精化迭代法每次迭代的具體步驟總結(jié)為下面的流程圖(見下圖1)。

輸入初始迭代向量x(0),終止判別τk和精度θ。 Step 1.計(jì)算殘差r(k); Step 2.利用GMERES方法求解Ad(k)=r(k),尋找滿足終止檢驗(yàn)公式(4)的d(k,Lk); Step 3.更新x(k+1)=x(k)+d(k,Lk); Step 4.若滿足終止檢驗(yàn)公式(6),則輸出近似解x(k+1),計(jì)算結(jié)束;否則,令 x(k)=x(k+1),k=k+1,轉(zhuǎn)Step 1。

圖1精化迭代法每次迭代的具體步驟
Fig.1Specific steps of each iteration of the refined iterative method

顯然，對(duì)于不同的τk值，d(k，Lk)不同，精化迭代法的收斂速度也不同。當(dāng)τk等于0時(shí)，可以把算法1得到的解認(rèn)為是公式(1)的精確解。因此，當(dāng)選取的τk值較小時(shí)，迭代序列有較好的收斂速度。

定理2 令x*為線性方程組(1)的精確解。對(duì)第k步迭代，用GMRES方法求解(3)的初始迭代向量均為d(k，0)=0，相應(yīng)的迭代序列為{d(k，l)}Lkl=0且d(k，Lk)滿足終止檢驗(yàn)公式(4)，則迭代序列{x(k)}k=0滿足

‖x(k+1)-x*‖≤τkκ(A)‖x(k)-x*‖，

其中:κ(A)=‖A‖·‖A-1‖。

特別地，當(dāng)

ηmax=τmaxκ(A)<1，

迭代序列{x(k)}k=1收斂到x*，其中

證明 由(4)式知

x(k+1)-x*=x(k)-x*+d(k，Lk)=
A-1[Ad(k，Lk)-rk]，

從而

‖x(k+1)-x*‖≤‖A-1‖‖Ad(k，Lk)-rk‖≤τk‖A-1‖‖rk‖≤
τk‖A-1‖·‖A‖·‖x(k)-x*‖=τkκ(A)‖x(k)-x*‖≤ηmax‖x(k)-x*‖。

當(dāng)ηmax<1時(shí)，迭代序列{x(k)}k=1收斂，證畢。

2 數(shù)值實(shí)驗(yàn)

下面給出精化迭代法應(yīng)用于曲面擬合的兩個(gè)例子。采用三次均勻B樣條曲面插值相應(yīng)的數(shù)據(jù)點(diǎn)，利用精化迭代法求相應(yīng)的控制點(diǎn)，從而得到插值給定點(diǎn)集的曲面。第k步迭代得到的三次均勻B樣條曲面的擬合誤差計(jì)算式為

右端向量b等于數(shù)據(jù)點(diǎn){Pij}(i=0，1，…，m；j=0，1，…，n)，系數(shù)矩陣均為A=B1?B2，符號(hào)?表示Kronecker積，其中B1∈Rm×m，B2∈Rn×n為配置矩陣，形如

例1 利用三次均勻B樣條曲面逼近peak函數(shù)

中的點(diǎn)

例2 利用三次均勻B樣條曲面逼近

中的點(diǎn)

利用三次均勻B樣條曲面迭代逼近例1、2中的點(diǎn)。所有的數(shù)值實(shí)驗(yàn)是在Matlab7.4.0.287(R2014a)環(huán)境中實(shí)現(xiàn)。在數(shù)值實(shí)驗(yàn)中，取τk=10-2，ε=10-8。為了比較，還測(cè)試了文獻(xiàn)[6]中加權(quán)漸進(jìn)迭代逼近法(WPIA)。表1、2給出了系數(shù)矩陣的階數(shù)不同時(shí)，三次均勻B樣條曲面的精化迭代法和WPIA逼近例1、2中的點(diǎn)所需迭代次數(shù)和計(jì)算時(shí)間。

表1 精化迭代法與加權(quán)漸進(jìn)迭代逼近法求解例1

Table 1 Comparisons between the refined iterative method and WPIA for Case 1

矩陣階數(shù)/n加權(quán)漸進(jìn)迭代逼近迭代次數(shù)/次Time/s精化迭代法迭代次數(shù)/次Time/s相對(duì)效率/%900591.24E-01323.93E-0268.31 600616.85E-01292.12E-0169.02 500621.71E+00294.94E-0171.13 600623.67E+00227.40E-0179.84 900636.63E+00201.29E+0080.5

表2 精化迭代法與加權(quán)漸進(jìn)迭代逼近法求解例2

Table 2 Comparisons between the refined iterative method and WPIA for Case 2

矩陣階數(shù)/n加權(quán)漸進(jìn)迭代逼近迭代次數(shù)/次Time/s精化迭代法迭代次數(shù)/次Time/s相對(duì)效率/%900511.32E-01313.92E-0270.21 600546.69E-01211.55E-0176.92 500561.55E+00213.28E-0178.93 600563.25E+00206.67E-0179.54 900565.81E+00201.37E+0076.5

在表1、2中，n為系數(shù)矩陣A的階數(shù)；Time為單獨(dú)計(jì)算100次的平均計(jì)算時(shí)間，相對(duì)效率計(jì)算公式為

其中：tw為加權(quán)漸進(jìn)迭代逼近法的計(jì)算時(shí)間；tr為精化迭代法的計(jì)算時(shí)間。

由表1、2可看出，在迭代次數(shù)和計(jì)算時(shí)間上，精化迭代法都要優(yōu)于加權(quán)漸進(jìn)迭代逼近法，且矩陣規(guī)模越大，優(yōu)勢(shì)就愈明顯。