精細積分阻尼譜修正迭代法在病態線性方程組中的應用

覃有堂 林賢坤 覃柏英

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

關鍵詞：精細積分；阻尼因子；譜修正迭代法

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

0 引言

本文考慮如下的病態線性方程組問題：

AX=b （1）

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

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

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

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

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

BX=H， （2）

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

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

=B-1H， （3）

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

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

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

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

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

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

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

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

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

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

證明 設？姿為矩陣B的任意一特征值，其對應的特征向量為？孜，則B？孜=？姿？孜，從而

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

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

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

證明 設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.產生的矩陣序列？琢P，？琢2P2，…，？琢kP k，…收斂，且？琢kP k=0.

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

‖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構造的矩陣冪級數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可逆.

對于任意非負整數k，設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

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

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

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

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

取定一小步長？駐t，設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）

設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]

設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 經典算例

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

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

A19×4=， b19×1=

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

算例2 系數矩陣A和常數項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，其條件數為2.996 7E5，因此線性方程組（1） 對應的線性方程組（2）為弱病態線性方程組.

算例3 系數矩陣A和常數項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

若常數項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時，其條件數為4.000 0E7，因此方程組（1）對應的方程組（2）為弱病態線性方程組.

算例4 系數矩陣為典型的Hilbert病態矩陣，即系數矩陣A和常數項b分別如下：

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

若常數項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時，其條件數為1.525 8E10，因此線性方程組（1）對應的線性方程組（2）為強病態線性方程組.

4 算法的性能分析

若已知病態線性方程組（1）的真解為X和數值迭代解為，則有b=AX和=A.由此可設Eb=-b和Ex=-X，則可定義評價數值迭代解精度的如下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表示算法的迭代次數，結果分別如表1—表4所示.

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

5 算法應用于高維問題

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

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

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

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

6 結束語

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

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

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

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

參考文獻

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

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

[3] 錢進， 吳時國， 崔若飛， 等. 新驛礦區奧陶系灰巖孔隙度預測方法[J]. 桂林理工大學學報， 2012，32（1）：43-47.

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

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

[6] 余婧妮，向宇，李曉妮，等. 求解梁大變形問題的一種精細分析方法[J]. 廣西科技大學學報，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] 李鵬飛， 李鵬舉， 趙建民， 等. 應用自適應混合遺傳算法求解病態線性方程組[J]. 科學技術與工程， 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] 李雪芳， 王希云. 求解病態線性方程組的自動控制步長法[J]. 太原科技大學學報， 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] 胡圣榮，戴納新. 病態線性方程組新解法：增廣方程組法[J]. 華南農業大學學報， 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] 王新洲，劉丁酉，張前勇，等. 譜修正迭代法及其在測量數據處理中的應用[J]. 黑龍江工程學院學報， 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

（學科編輯：張玉鳳）