覃有堂 林賢坤 覃柏英

摘 要：基于譜修正迭代法，引入阻尼因子和結(jié)合精細積分的矩陣求逆，提出了精細積分阻尼譜修正迭代法，并將其應(yīng)用于病態(tài)線性方程組的求解.通過4個經(jīng)典算例，探討矩陣的精細積分求逆對阻尼譜修正迭代法求解病態(tài)線性方程組的影響.結(jié)果表明，矩陣的精細積分求逆可提高阻尼譜修正迭代法求解病態(tài)線性方程組的精度.再通過兩個經(jīng)典算例，探討了精細積分阻尼譜修正迭代法在高維病態(tài)線性方程組中的應(yīng)用.結(jié)果表明，該算法也可提高高維病態(tài)線性方程組求解的精度.

關(guān)鍵詞：精細積分；阻尼因子；譜修正迭代法

中圖分類號：O241.6；O151.2 文獻標志碼：A

0 引言

本文考慮如下的病態(tài)線性方程組問題：

AX=b （1）

其中：A——m×n實系數(shù)矩陣，X——n×1的解向量，b——m×1的實常數(shù)項.

病態(tài)線性方程組在工程設(shè)計、圖像處理、地震勘探、GNSS衛(wèi)星定位等科學(xué)和工程領(lǐng)域有著重要的應(yīng)用[1-6].由于系數(shù)矩陣的微小擾動或數(shù)值計算過程中存在的舍入誤差，往往引起誤差傳播的難以控制，可能導(dǎo)致該問題的數(shù)值解存在較大擾動甚至嚴重失真[7-9].

目前病態(tài)線性方程組的數(shù)值求解方法很多[10-16]，如誤差轉(zhuǎn)移法和增廣方程組法等直接法，殘差校正迭代法和Wilkinson迭代法等迭代法，以及遺傳算法、模擬退火法和混合算法等智能算法.但這些算法都有一定的局限性，如求解的效率或解的精度和穩(wěn)定性等.為了改善矩陣的病態(tài)性，王新洲等[17]提出了譜修正迭代法.該算法在不改變方程的等量關(guān)系基礎(chǔ)上，可改善矩陣的病態(tài)性.該算法雖然理論上收斂于問題的最優(yōu)解，但實際應(yīng)用中，其收斂效果并不理想，收斂速度較慢.

為了克服譜修正迭代法的缺點，實現(xiàn)高精度和高效率求解病態(tài)線性方程組，引入阻尼因子和結(jié)合矩陣的精細積分求逆，本文提出精細積分阻尼譜修正迭代法，將其應(yīng)用于病態(tài)線性方程組的求解，探討該算法求解病態(tài)線性方程組的收斂速度和計算精度.

設(shè)B=ATA和H=ATb，則B為n×n的實矩陣，H為n×1的實列向量，從而

BX=H， （2）

因BT=B，所以實矩陣B為對稱矩陣.若B為正定矩陣，則線性方程組（2）存在唯一解.

根據(jù)最小二乘原理，線性方程組（2）的最小二乘解為：

=B-1H， （3）

求解線性方程組（1）的算法稱為最小二乘法，所求的解相應(yīng)稱為最小二乘解.

當線性方程組（1）病態(tài)時，其系數(shù)矩陣A的條件數(shù)往往很大，則矩陣B的條件數(shù)也會很大，此時采用最小二乘法（3）所求的病態(tài)線性方程組（1）的最小二乘解精度往往較低.

為了改善線性方程組（1）的病態(tài)性，可在式（2）等號兩邊加上aX，I為n階單位陣，有：

（B+？琢I）X=H+？琢X， （4）

其中？琢稱為阻尼因子.設(shè)N=B+？琢I，則NT=N，因此N為實對稱矩陣.設(shè)Y為n×1的任意非零實向量，則YTNY=YT（B+？琢I）Y=YTBY+？琢YTY≥？琢YTY>0，從而實矩陣N為正定矩陣.因此線性方程組（4） 也存在唯一解，且該解也是線性方程組（2）的唯一解.

對于線性方程組（2），要采用如下的阻尼譜修正迭代法進行數(shù)值迭代求解：

k=（B+？琢I）-1H+？琢（k-1）（5）

1 阻尼譜修正迭代法的收斂性

對于對稱且正定的實矩陣B，其存在n個實特征值，從大到小排序為？姿1，？姿2，…，？姿n. 為了敘述的方便，設(shè)N=B+？琢I和P=（B+？琢I）-1，探討阻尼譜修正迭代法的收斂性.

引理1 矩陣N的n個特征值為？姿1+？琢，？姿2+？琢，…，？姿n+？琢.

證明 設(shè)？姿為矩陣B的任意一特征值，其對應(yīng)的特征向量為？孜，則B？孜=？姿？孜，從而

N？孜=（B+？琢I）？孜=B？孜+？琢？孜=？姿？孜+？琢？孜=（？姿+？琢）？孜，

即？姿+？琢為N的特征值，從而矩陣N也存在n個特征值？姿1+？琢，？姿2+？琢，…，？姿n+？琢.

引理2 矩陣P的譜范數(shù)‖P‖2 =（？姿n+？琢）-1.

證明 設(shè)D=diag（？姿1，？姿2，…，？姿n），則D+？琢I=diag（？姿1+？琢，？姿2+？琢，…，？姿n+？琢）.因為矩陣B為實對稱矩陣，所以存在正交矩陣Q使得QTBQ=D.有B+？琢I=Q（D+？琢I）QT.從而：

P=Q（D+？琢I）QT-1=Q-1（D+？琢I）-1（QT）-1=QT（D+？琢I）-1Q .

矩陣P的n個特征值為（？姿1+？琢）-1，（？姿2+？琢）-1，…，（？姿n+？琢）-1，其最大值（？姿n+？琢）-1.

又因為矩陣Q為正交矩陣，從而

‖P‖2 =‖QT（D+？琢I）-1Q‖2 =‖（D+？琢I）-1‖2 =（？姿n+？琢）-1 .

引理3 當實矩陣B正定時，對任意的阻尼因子？琢滿足？琢>0，由n階方陣P.產(chǎn)生的矩陣序列？琢P，？琢2P2，…，？琢kP k，…收斂，且？琢kP k=0.

證明 由矩陣范數(shù)的相容性定義，有：

‖P k‖2 =‖P k-1P‖2≤‖P k-1‖2‖P‖2≤…≤‖P

因為矩陣B正定，從而其最小特征值？姿n>0.于是？琢>0時，有‖？琢P‖2 =？琢/（？姿n+？琢）<1，有：

‖？琢kP k-0‖2 =‖？琢kP k‖2≤‖？琢P=？琢k/（？姿n+？琢）k→0（k→∞）.

從而序列？琢P，？琢2P 2，…，？琢kP k，…收斂，且？琢kP k=0.

引理4 當矩陣B正定時，對任意的阻尼因子？琢滿足？琢>0，由n階方陣P構(gòu)造的矩陣冪級數(shù)I+？琢P+？琢2P 2+…+？琢kP k+…收斂，且其和為（I-？琢P）-1.

證明 因矩陣B正定，則其最小特征值？姿n>0.當？琢>0時，有‖？琢P‖2=？琢/（？姿n+？琢）<1，且：

‖（I-？琢P）X‖2 =‖X-？琢PX‖2 ≥‖X‖2 -？琢‖PX‖2 ≥‖X‖2 -？琢‖P‖2 ‖X‖2 =（1-？琢‖P‖2）‖X‖2 >0

從而（I-？琢P）X≠0，這表明線性方程組（I-？琢P）X=0只有零解，因此矩陣I-？琢P可逆.

對于任意非負整數(shù)k，設(shè)Sk=I+？琢P+？琢2P 2+…+？琢kP k，有：

I-（I+？琢P+？琢2P 2+…+？琢kP k）（I-？琢P）=I-Sk（I-？琢P）=？琢k+1P k+1，

在等式兩邊右乘以矩陣I-？琢P的逆矩陣（I-？琢P）-1，因P（I-？琢P）-1=B-1，則：

（I-？琢P）-1-Sk=？琢k+1P k+1（I-？琢P）-1=？琢k+1P kB-1 ，

則‖（I-？琢P）-1-Sk-0‖2 =‖？琢k+1P kB-1‖2 ≤？琢‖？琢kP k‖2 ‖B-1‖2 =？琢k+1‖B-1‖2 /（？姿n+？琢）k.

因？琢>0和？姿n>0，則？姿n+？琢>？琢>0，有0<？琢/（？姿n+？琢）<1，從而？琢k+1/（？姿n+？琢）k=0.‖（I-？琢P）-1-Sk-0‖2 =0，即（I-？琢P）-1-Sk=0. 亦即：

（I+？琢P+？琢2P 2+…+？琢kP k）=（I-？琢P）-1 .

定理1 當矩陣B正定時，對任意的阻尼因子？琢，滿足？琢>0和任意的初值（0），有阻尼譜修正迭代法k=PH+？琢（k-1）收斂，且k=B-1H.

證明 由k=PH+？琢（k-1）=PH+？琢Pk-1有：

k=PH+？琢PPH+？琢Pk-2=PH+？琢P 2H+？琢2P 2k-2=…=P（I+？琢P+？琢2P2+…+？琢k-1P k-1）H+？琢kP k（0）

由引理3和引理4知，當矩陣B正定時，對任意的阻尼因子？琢滿足？琢>0和任意的初值（0），有：

k=P（I+？琢P+？琢2P 2+…+？琢k-1P k-1）H=P（I-？琢P）-1H=（I-？琢P）-1P-1-1H=（P-1-？琢I）-1H=B-1H

由以上證明知，若系數(shù)矩陣A不對稱，取B=ATA，否則取B=A.由此若可能再對常數(shù)項H歸一化，則可將線性方程組（1）轉(zhuǎn)化線性方程組（2），此時其實系數(shù)矩陣B對稱.由定理1知，對任意阻尼因子？琢>0，當矩陣B正定時，式（5）的阻尼譜修正迭代法收斂于方程組（1）的真解.

2 基于精細積分的矩陣求逆

為了保證求矩陣N=B+？琢I的逆矩陣P=（B+？琢I）-1的精度，本文提出精確求解逆矩陣P的精細積分法.為了敘述方便，設(shè)M=-N=-（B+？琢I），則M為負定矩陣.

設(shè)F（t）=eMsds，則F（t）=eMsds=M-1eMs│=M-1eMt-M-1. M為負定矩陣，則eMs=0，從而F（t）=-M-1=P.即矩陣N的逆矩陣P=eMsds.

取定一小步長？駐t，設(shè)t=3？駐t，m=0，1，…，則P=F（3？駐t）.

F（3？駐t）=eMsds=eMsds+eMsds+eMsds=（I+e3？駐tM+e2 ·3？駐tM）F（3？駐t）

設(shè)Fm=F（3？駐t），Rm=e3？駐 tM，m=0，1，2，…，則可采用如下迭代式計算P：

Fm+1=（I+Rm+）Fm，m=0，1，2，…

由上式知，當m足夠大時可得到較精確的P.為了利用上式計算P，需給定初始值F0和Rm.

對eMs進行Taylor展開，有eMs≈I+Ms+（Ms）2/2+（Ms）3/6+（Ms）4/24，從而

F0=F（？駐t）=eMsds≈？駐t[I+（M？駐t）/2+（M？駐t）2/6+（M？駐t）3/24+（M？駐t）4/120]

設(shè)Rm=I+Tm，m=0，1，….因R0=eM△t≈I+M？駐t+（M？駐t）2/2+（M？駐t）3/6+（M？駐t）4/24，則T0=M？駐t+（M？駐t）2/2+（M？駐t）3/6+（M？駐t）4/24. 又因為Rm=e3△tM=R，從而有Rm=（I+Tm-1）3=I+3Tm-1+3T+T=I+Tm.可采用如下迭代式計算Tm：

Tm=3Tm-1+3T+T，m=0，1，…

3 經(jīng)典算例

為了探討矩陣的精細積分求逆對阻尼譜修正迭代法求解病態(tài)線性方程組的性能影響，以及精細積分阻尼譜修正迭代法求解病態(tài)線性方程組的精度與效率，現(xiàn)給出4個經(jīng)典的算例.

算例1 系數(shù)矩陣A和常數(shù)項b分別如下：

A19×4=， b19×1=

此時線性方程組（1）的真解為X4 ×1=[0.2，1.5，1.6，-2.8]T.因矩陣A非對稱，取B=ATA，其條件數(shù)為1.610 1E9，因此線性方程組（1）對應(yīng)的線性方程組（2）為弱病態(tài)線性方程組.

算例2 系數(shù)矩陣A和常數(shù)項b分別如下：

A18×7=， b18×1=

此時線性方程組（1）的真解為X7×1=[0.2，2，1.5，-1.6，4.8，3.4，-2.1]T.因矩陣A非對稱，取B=ATA，其條件數(shù)為2.996 7E5，因此線性方程組（1） 對應(yīng)的線性方程組（2）為弱病態(tài)線性方程組.

算例3 系數(shù)矩陣A和常數(shù)項b分別如下：

A=（ai，j）n×n，ai，j=1， i≠j，1+p2， i=j， i，j=1，2，…，n，b=[b1，b2，…，bn]T

若常數(shù)項b取bi=ai，j和bi=jai，j，i=1，2，…，n，則線性方程組（1）的真解分別為：Xn×1=[1，1，…，1]T和

Xn×1=[1，2，…，n]T.因矩陣A對稱且正定，取B=A，當n=10和p=5×10-3時，其條件數(shù)為4.000 0E7，因此方程組（1）對應(yīng)的方程組（2）為弱病態(tài)線性方程組.

算例4 系數(shù)矩陣為典型的Hilbert病態(tài)矩陣，即系數(shù)矩陣A和常數(shù)項b分別如下：

A=（ai，j）n×n，ai，j=，i，j=1，2，…，n，b=[b1，b2，…，bn]T

若常數(shù)項b取bi=ai，j和bi=jai，j，其中i=1，2，…，n，則線性方程組（1）的真解分別為Xn×1=[1，1，…，1]T和Xn×1=[1，2，…，n]T.因矩陣A對稱且正定，取B=A，當n=8時，其條件數(shù)為1.525 8E10，因此線性方程組（1）對應(yīng)的線性方程組（2）為強病態(tài)線性方程組.

4 算法的性能分析

若已知病態(tài)線性方程組（1）的真解為X和數(shù)值迭代解為，則有b=AX和=A.由此可設(shè)Eb=-b和Ex=-X，則可定義評價數(shù)值迭代解精度的如下3個殘差：

Eb=‖Eb‖2 ，Ex=‖Ex /n，E∞=‖Ex‖∞ /‖X‖∞

殘差Eb可用于評價滿足線性方程組（1）的程度，而殘差Ex實際是均方誤差（mean squared error，MSE），可評價偏離X的程度.殘差E∞實際是相對誤差，更能反映的可信程度.

基于軟件MATLAB R2016a，采用LSM，DCCV和PIM-DCCV分別求解上述4個算例，？琢取值分別為0.28，0.089，4.0×10-14和5.0×10-12，精細積分求逆矩陣時△t=1×10-16和m=100，且算例3 的n=10和p=5×10-4，算例4的n=8，T表示算法的迭代次數(shù)，結(jié)果分別如表1—表4所示.

由表1—表4知，3種算法都可實現(xiàn)4個算例較高精度的求解.但DCCV和PIM-DCCV都優(yōu)于LSM，這說明DCCV有利于病態(tài)線性方程組求解精度的提高.同時，PIM-DCCV也明顯優(yōu)于DCCV，這也說明精細積分可更精確求解逆矩陣，從而提高算法DCCV的精度.因此，算法PIM-IDCCV可使數(shù)值迭代解更滿足線性方程組（1）和接近其真實解X.

5 算法應(yīng)用于高維問題

為了更充分說明精細積分改進阻尼譜修正迭代法求解病態(tài)線性方程組的性能，現(xiàn)將其應(yīng)用于高維問題，即將該算法應(yīng)用于算例3和算例4的高維病態(tài)線性方程組的求解.

對于算例3，分別取維數(shù)n=100，200，500，1 000，2 000，3 000，4 000，參數(shù)p取5×10-6和5×10-4時，矩陣A條件數(shù)分別如表5和表6所示.當取bi=ai，j時，阻尼因子？琢=1，當取bi=jai，j時，阻尼因子？琢=5.4×10-14，采用PIM-DCCV分別求解，結(jié)果分別如表5和表6所示.

對于算例4的病態(tài)線性方程組，分別取其維數(shù)n=100，200，500，1 000，2 000，3 000，4 000，系數(shù)矩陣A對應(yīng)的條件數(shù)分別如表7所示. 當常數(shù)項b取bi=ai，j時，？琢=5.0×10-2，當常數(shù)項b取bi=jai，j時，？琢=5.0×10-12，采用PIM-DCCV分別求解，結(jié)果分別如表7和表8所示.

采用精細積分改進阻尼譜修正迭代法對高維病態(tài)線性方程組進行數(shù)值求解，由表5—表8可知，對于高維病態(tài)問題，該算法可獲得較高精度的數(shù)值解，該解也比較滿足方程組（1）和接近其真實解X.

6 結(jié)束語

為了提高譜修正迭代法求解病態(tài)線性方程組的精度，引入阻尼因子和結(jié)合精細積分的矩陣求逆，本文提出了精細積分阻尼譜修正迭代法，并將其應(yīng)用于病態(tài)線性方程組的求解.采用4個經(jīng)典算例，探討了精細積分求逆對阻尼譜修正迭代法求解病態(tài)線性方程組的性能影響.同時，通過兩個經(jīng)典算例，探討了精細積分阻尼因子譜修正迭代法在高維病態(tài)線性方程組中的應(yīng)用.由此可得到如下結(jié)論：

1） 阻尼因子的引入可改善系數(shù)矩陣的病態(tài)，從而提高算法求解病態(tài)線性方程組的效率與精度；

2） 結(jié)合精細積分求解逆矩陣，可提高阻尼譜修正迭代法求解病態(tài)線性方程組的精度；

3） 精細積分阻尼譜修正迭代法，明顯可提高病態(tài)線性方程組的求解精度，且應(yīng)用于高維弱病態(tài)線性方程組，也可獲得較高的精度數(shù)值迭代解.

參考文獻

[1] 熊新， 吳洪濤， 于學(xué)華， 等. 工程車輛油氣懸架參數(shù)化建模與幅頻特性分析[J]. 振動、測試與診斷， 2014，34（5）： 926-931.

[2] 曲昕馨， 李桐林， 王飛. 基于數(shù)字圖像分割法的跨孔雷達走時層析成像[J]. 吉林大學(xué)學(xué)報（地球科學(xué)版）， 2014，44（4）： 1340-1347.

[3] 錢進， 吳時國， 崔若飛， 等. 新驛礦區(qū)奧陶系灰?guī)r孔隙度預(yù)測方法[J]. 桂林理工大學(xué)學(xué)報， 2012，32（1）：43-47.

[4] 秦紅磊， 梁敏敏. 基于GNSS的高軌衛(wèi)星定位技術(shù)研究[J]. 空間科學(xué)學(xué)報，2008，28（4）： 316-325.

[5] 李曉妮，程濤. 分析微梁尺寸效應(yīng)的傳遞矩陣法[J]. 廣西科技大學(xué)學(xué)報，2017，28（1）：91-96.

[6] 余婧妮，向宇，李曉妮，等. 求解梁大變形問題的一種精細分析方法[J]. 廣西科技大學(xué)學(xué)報，2014，25（4）：34-39.

[7] DENG X S，YIN L B，PENG S C，et al. An iterative algorithm for solving ill-conditioned linear least squares problems[J].Geodesy and Geodynamics， 2015，6（6）： 453-459.

[8] 李鵬飛， 李鵬舉， 趙建民， 等. 應(yīng)用自適應(yīng)混合遺傳算法求解病態(tài)線性方程組[J]. 科學(xué)技術(shù)與工程， 2010，10（9）： 2098-

2102.

[9] BREZINSKI C ， NOVATI P， REDIVO-ZAGLIA M. A rational Arnoldi approach for ill-conditioned linear systems[J].Journal of

Computational and Applied Mathematics， 2012， 236（8）： 2063-2077.

[10] WU X Y. Error analysis for the predictor-corrector process relating to ill-conditioned linear system of equations[J].Applied Mathematics and Computation， 2007， 186（1）： 530-534.

[11] 李雪芳， 王希云. 求解病態(tài)線性方程組的自動控制步長法[J]. 太原科技大學(xué)學(xué)報， 2012，33（6）： 480-482.

[12] LIU C S. Optimally scaled vector regularization method to solve ill-posed linear problems[J]. Applied Mathematics and Computation， 2012， 218（21）： 10602-10616.

[13] SMOKTUNOWICZ A，SMOKTUNOWICZ A. Iterative refinement techniques for solving block linear systems of equations[J].Applied Numerical Mathematics， 2013， 67（67）： 220-229.

[14] 胡圣榮，戴納新. 病態(tài)線性方程組新解法：增廣方程組法[J]. 華南農(nóng)業(yè)大學(xué)學(xué)報， 2009，30（1）： 119-121.

[15] MATINFAR M， ZAREAMOGHADDAM H， ESLAMI M， et al. GMRES implementations and residual smoothing techniques for solving ill-posed linear systems[J].Computers & Mathematics with Applications， 2012， 63（1）： 1-13.

[16] WU X Y， FANG Y L. Wilkinsons iterative refinement of solution with automatic step-size control for linear system of equations[J].Applied Mathematics and Computation， 2007， 193（2）： 506-513.

[17] 王新洲，劉丁酉，張前勇，等. 譜修正迭代法及其在測量數(shù)據(jù)處理中的應(yīng)用[J]. 黑龍江工程學(xué)院學(xué)報， 2001，15（2）： 3-6.

Iteration method by correcting eigenvalue with damping factor based

on precise integration for ill-conditioned linear system

QIN You-tang1， LIN Xian-kun*2 ， QIN Bo-ying2

（1. Liuzhou No.1 Vocational and Technical School， Liuzhou 545616， China； 2. School of Automobile and

Traffic Engineering， Guangxi University of Science and Technology， Liuzhou 545006， China）

Abstract： To improve the performance of the iteration method by correcting eigenvalue （IMCE）， a damping factor and the precise integration method are introduced to IMCE. Thus an iteration method by correcting characteristic value with damping factor based on the precise integration method （PIM-DIMCE） is presented in the paper. Four classical examples are used to discuss the influence of the precise integration method （PIM） when IMCE solves the ill conditioned linear system （ICLS）. The results show that PIM can improve the accuracy of IMCE for solving ICLS. Meanwhile， two classical examples are used to discuss the performance of PIM-DIMCE for solving the high dimensional ICLS. High computation accuracy is achieved as well when solving the high dimensional ICLS.

Key words： precise integration method； damping factor； iteration method by correcting eigenvalue

（學(xué)科編輯：張玉鳳）