覃有堂 林賢坤 覃柏英
摘 要:基于譜修正迭代法,引入阻尼因子和結(jié)合精細(xì)積分的矩陣求逆,提出了精細(xì)積分阻尼譜修正迭代法,并將其應(yīng)用于病態(tài)線性方程組的求解.通過4個(gè)經(jīng)典算例,探討矩陣的精細(xì)積分求逆對(duì)阻尼譜修正迭代法求解病態(tài)線性方程組的影響.結(jié)果表明,矩陣的精細(xì)積分求逆可提高阻尼譜修正迭代法求解病態(tài)線性方程組的精度.再通過兩個(gè)經(jīng)典算例,探討了精細(xì)積分阻尼譜修正迭代法在高維病態(tài)線性方程組中的應(yīng)用.結(jié)果表明,該算法也可提高高維病態(tài)線性方程組求解的精度.
關(guān)鍵詞:精細(xì)積分;阻尼因子;譜修正迭代法
中圖分類號(hào):O241.6;O151.2 文獻(xiàn)標(biāo)志碼:A
0 引言
本文考慮如下的病態(tài)線性方程組問題:
AX=b (1)
其中:A——m×n實(shí)系數(shù)矩陣,X——n×1的解向量,b——m×1的實(shí)常數(shù)項(xiàng).
病態(tài)線性方程組在工程設(shè)計(jì)、圖像處理、地震勘探、GNSS衛(wèi)星定位等科學(xué)和工程領(lǐng)域有著重要的應(yīng)用[1-6].由于系數(shù)矩陣的微小擾動(dòng)或數(shù)值計(jì)算過程中存在的舍入誤差,往往引起誤差傳播的難以控制,可能導(dǎo)致該問題的數(shù)值解存在較大擾動(dòng)甚至嚴(yán)重失真[7-9].
目前病態(tài)線性方程組的數(shù)值求解方法很多[10-16],如誤差轉(zhuǎn)移法和增廣方程組法等直接法,殘差校正迭代法和Wilkinson迭代法等迭代法,以及遺傳算法、模擬退火法和混合算法等智能算法.但這些算法都有一定的局限性,如求解的效率或解的精度和穩(wěn)定性等.為了改善矩陣的病態(tài)性,王新洲等[17]提出了譜修正迭代法.該算法在不改變方程的等量關(guān)系基礎(chǔ)上,可改善矩陣的病態(tài)性.該算法雖然理論上收斂于問題的最優(yōu)解,但實(shí)際應(yīng)用中,其收斂效果并不理想,收斂速度較慢.
為了克服譜修正迭代法的缺點(diǎn),實(shí)現(xiàn)高精度和高效率求解病態(tài)線性方程組,引入阻尼因子和結(jié)合矩陣的精細(xì)積分求逆,本文提出精細(xì)積分阻尼譜修正迭代法,將其應(yīng)用于病態(tài)線性方程組的求解,探討該算法求解病態(tài)線性方程組的收斂速度和計(jì)算精度.
設(shè)B=ATA和H=ATb,則B為n×n的實(shí)矩陣,H為n×1的實(shí)列向量,從而
BX=H, (2)
因BT=B,所以實(shí)矩陣B為對(duì)稱矩陣.若B為正定矩陣,則線性方程組(2)存在唯一解.
根據(jù)最小二乘原理,線性方程組(2)的最小二乘解為:
=B-1H, (3)
求解線性方程組(1)的算法稱為最小二乘法,所求的解相應(yīng)稱為最小二乘解.
當(dāng)線性方程組(1)病態(tài)時(shí),其系數(shù)矩陣A的條件數(shù)往往很大,則矩陣B的條件數(shù)也會(huì)很大,此時(shí)采用最小二乘法(3)所求的病態(tài)線性方程組(1)的最小二乘解精度往往較低.
為了改善線性方程組(1)的病態(tài)性,可在式(2)等號(hào)兩邊加上aX,I為n階單位陣,有:
(B+?琢I)X=H+?琢X, (4)
其中?琢稱為阻尼因子.設(shè)N=B+?琢I,則NT=N,因此N為實(shí)對(duì)稱矩陣.設(shè)Y為n×1的任意非零實(shí)向量,則YTNY=YT(B+?琢I)Y=YTBY+?琢YTY≥?琢YTY>0,從而實(shí)矩陣N為正定矩陣.因此線性方程組(4) 也存在唯一解,且該解也是線性方程組(2)的唯一解.
對(duì)于線性方程組(2),要采用如下的阻尼譜修正迭代法進(jìn)行數(shù)值迭代求解:
k=(B+?琢I)-1H+?琢(k-1)(5)
1 阻尼譜修正迭代法的收斂性
對(duì)于對(duì)稱且正定的實(shí)矩陣B,其存在n個(gè)實(shí)特征值,從大到小排序?yàn)椋孔?,?姿2,…,?姿n. 為了敘述的方便,設(shè)N=B+?琢I和P=(B+?琢I)-1,探討阻尼譜修正迭代法的收斂性.
引理1 矩陣N的n個(gè)特征值為?姿1+?琢,?姿2+?琢,…,?姿n+?琢.
證明 設(shè)?姿為矩陣B的任意一特征值,其對(duì)應(yīng)的特征向量為?孜,則B?孜=?姿?孜,從而
N?孜=(B+?琢I)?孜=B?孜+?琢?孜=?姿?孜+?琢?孜=(?姿+?琢)?孜,
即?姿+?琢為N的特征值,從而矩陣N也存在n個(gè)特征值?姿1+?琢,?姿2+?琢,…,?姿n+?琢.
引理2 矩陣P的譜范數(shù)‖P‖2 =(?姿n+?琢)-1.
證明 設(shè)D=diag(?姿1,?姿2,…,?姿n),則D+?琢I=diag(?姿1+?琢,?姿2+?琢,…,?姿n+?琢).因?yàn)榫仃嘊為實(shí)對(duì)稱矩陣,所以存在正交矩陣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個(gè)特征值為(?姿1+?琢)-1,(?姿2+?琢)-1,…,(?姿n+?琢)-1,其最大值(?姿n+?琢)-1.
又因?yàn)榫仃嘠為正交矩陣,從而
‖P‖2 =‖QT(D+?琢I)-1Q‖2 =‖(D+?琢I)-1‖2 =(?姿n+?琢)-1 .
引理3 當(dāng)實(shí)矩陣B正定時(shí),對(duì)任意的阻尼因子?琢滿足?琢>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
因?yàn)榫仃嘊正定,從而其最小特征值?姿n>0.于是?琢>0時(shí),有‖?琢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 當(dāng)矩陣B正定時(shí),對(duì)任意的阻尼因子?琢滿足?琢>0,由n階方陣P構(gòu)造的矩陣冪級(jí)數(shù)I+?琢P+?琢2P 2+…+?琢kP k+…收斂,且其和為(I-?琢P)-1.
證明 因矩陣B正定,則其最小特征值?姿n>0.當(dāng)?琢>0時(shí),有‖?琢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可逆.
對(duì)于任意非負(fù)整數(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 當(dāng)矩陣B正定時(shí),對(duì)任意的阻尼因子?琢,滿足?琢>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知,當(dāng)矩陣B正定時(shí),對(duì)任意的阻尼因子?琢滿足?琢>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不對(duì)稱,取B=ATA,否則取B=A.由此若可能再對(duì)常數(shù)項(xiàng)H歸一化,則可將線性方程組(1)轉(zhuǎn)化線性方程組(2),此時(shí)其實(shí)系數(shù)矩陣B對(duì)稱.由定理1知,對(duì)任意阻尼因子?琢>0,當(dāng)矩陣B正定時(shí),式(5)的阻尼譜修正迭代法收斂于方程組(1)的真解.
2 基于精細(xì)積分的矩陣求逆
為了保證求矩陣N=B+?琢I的逆矩陣P=(B+?琢I)-1的精度,本文提出精確求解逆矩陣P的精細(xì)積分法.為了敘述方便,設(shè)M=-N=-(B+?琢I),則M為負(fù)定矩陣.
設(shè)F(t)=eMsds,則F(t)=eMsds=M-1eMs│=M-1eMt-M-1. M為負(fù)定矩陣,則eMs=0,從而F(t)=-M-1=P.即矩陣N的逆矩陣P=eMsds.
取定一小步長(zhǎng)?駐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,…,則可采用如下迭代式計(jì)算P:
Fm+1=(I+Rm+)Fm,m=0,1,2,…
由上式知,當(dāng)m足夠大時(shí)可得到較精確的P.為了利用上式計(jì)算P,需給定初始值F0和Rm.
對(duì)eMs進(jìn)行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. 又因?yàn)镽m=e3△tM=R,從而有Rm=(I+Tm-1)3=I+3Tm-1+3T+T=I+Tm.可采用如下迭代式計(jì)算Tm:
Tm=3Tm-1+3T+T,m=0,1,…
3 經(jīng)典算例
為了探討矩陣的精細(xì)積分求逆對(duì)阻尼譜修正迭代法求解病態(tài)線性方程組的性能影響,以及精細(xì)積分阻尼譜修正迭代法求解病態(tài)線性方程組的精度與效率,現(xiàn)給出4個(gè)經(jīng)典的算例.
算例1 系數(shù)矩陣A和常數(shù)項(xiàng)b分別如下:
A19×4=, b19×1=
此時(shí)線性方程組(1)的真解為X4 ×1=[0.2,1.5,1.6,-2.8]T.因矩陣A非對(duì)稱,取B=ATA,其條件數(shù)為1.610 1E9,因此線性方程組(1)對(duì)應(yīng)的線性方程組(2)為弱病態(tài)線性方程組.
算例2 系數(shù)矩陣A和常數(shù)項(xiàng)b分別如下:
A18×7=, b18×1=
此時(shí)線性方程組(1)的真解為X7×1=[0.2,2,1.5,-1.6,4.8,3.4,-2.1]T.因矩陣A非對(duì)稱,取B=ATA,其條件數(shù)為2.996 7E5,因此線性方程組(1) 對(duì)應(yīng)的線性方程組(2)為弱病態(tài)線性方程組.
算例3 系數(shù)矩陣A和常數(shù)項(xiàng)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ù)項(xiàng)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對(duì)稱且正定,取B=A,當(dāng)n=10和p=5×10-3時(shí),其條件數(shù)為4.000 0E7,因此方程組(1)對(duì)應(yīng)的方程組(2)為弱病態(tài)線性方程組.
算例4 系數(shù)矩陣為典型的Hilbert病態(tài)矩陣,即系數(shù)矩陣A和常數(shù)項(xiàng)b分別如下:
A=(ai,j)n×n,ai,j=,i,j=1,2,…,n,b=[b1,b2,…,bn]T
若常數(shù)項(xiàng)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對(duì)稱且正定,取B=A,當(dāng)n=8時(shí),其條件數(shù)為1.525 8E10,因此線性方程組(1)對(duì)應(yīng)的線性方程組(2)為強(qiáng)病態(tài)線性方程組.
4 算法的性能分析
若已知病態(tài)線性方程組(1)的真解為X和數(shù)值迭代解為,則有b=AX和=A.由此可設(shè)Eb=-b和Ex=-X,則可定義評(píng)價(jià)數(shù)值迭代解精度的如下3個(gè)殘差:
Eb=‖Eb‖2 ,Ex=‖Ex /n,E∞=‖Ex‖∞ /‖X‖∞
殘差Eb可用于評(píng)價(jià)滿足線性方程組(1)的程度,而殘差Ex實(shí)際是均方誤差(mean squared error,MSE),可評(píng)價(jià)偏離X的程度.殘差E∞實(shí)際是相對(duì)誤差,更能反映的可信程度.
基于軟件MATLAB R2016a,采用LSM,DCCV和PIM-DCCV分別求解上述4個(gè)算例,?琢取值分別為0.28,0.089,4.0×10-14和5.0×10-12,精細(xì)積分求逆矩陣時(shí)△t=1×10-16和m=100,且算例3 的n=10和p=5×10-4,算例4的n=8,T表示算法的迭代次數(shù),結(jié)果分別如表1—表4所示.
由表1—表4知,3種算法都可實(shí)現(xiàn)4個(gè)算例較高精度的求解.但DCCV和PIM-DCCV都優(yōu)于LSM,這說明DCCV有利于病態(tài)線性方程組求解精度的提高.同時(shí),PIM-DCCV也明顯優(yōu)于DCCV,這也說明精細(xì)積分可更精確求解逆矩陣,從而提高算法DCCV的精度.因此,算法PIM-IDCCV可使數(shù)值迭代解更滿足線性方程組(1)和接近其真實(shí)解X.
5 算法應(yīng)用于高維問題
為了更充分說明精細(xì)積分改進(jìn)阻尼譜修正迭代法求解病態(tài)線性方程組的性能,現(xiàn)將其應(yīng)用于高維問題,即將該算法應(yīng)用于算例3和算例4的高維病態(tài)線性方程組的求解.
對(duì)于算例3,分別取維數(shù)n=100,200,500,1 000,2 000,3 000,4 000,參數(shù)p取5×10-6和5×10-4時(shí),矩陣A條件數(shù)分別如表5和表6所示.當(dāng)取bi=ai,j時(shí),阻尼因子?琢=1,當(dāng)取bi=jai,j時(shí),阻尼因子?琢=5.4×10-14,采用PIM-DCCV分別求解,結(jié)果分別如表5和表6所示.
對(duì)于算例4的病態(tài)線性方程組,分別取其維數(shù)n=100,200,500,1 000,2 000,3 000,4 000,系數(shù)矩陣A對(duì)應(yīng)的條件數(shù)分別如表7所示. 當(dāng)常數(shù)項(xiàng)b取bi=ai,j時(shí),?琢=5.0×10-2,當(dāng)常數(shù)項(xiàng)b取bi=jai,j時(shí),?琢=5.0×10-12,采用PIM-DCCV分別求解,結(jié)果分別如表7和表8所示.
采用精細(xì)積分改進(jìn)阻尼譜修正迭代法對(duì)高維病態(tài)線性方程組進(jìn)行數(shù)值求解,由表5—表8可知,對(duì)于高維病態(tài)問題,該算法可獲得較高精度的數(shù)值解,該解也比較滿足方程組(1)和接近其真實(shí)解X.
6 結(jié)束語(yǔ)
為了提高譜修正迭代法求解病態(tài)線性方程組的精度,引入阻尼因子和結(jié)合精細(xì)積分的矩陣求逆,本文提出了精細(xì)積分阻尼譜修正迭代法,并將其應(yīng)用于病態(tài)線性方程組的求解.采用4個(gè)經(jīng)典算例,探討了精細(xì)積分求逆對(duì)阻尼譜修正迭代法求解病態(tài)線性方程組的性能影響.同時(shí),通過兩個(gè)經(jīng)典算例,探討了精細(xì)積分阻尼因子譜修正迭代法在高維病態(tài)線性方程組中的應(yīng)用.由此可得到如下結(jié)論:
1) 阻尼因子的引入可改善系數(shù)矩陣的病態(tài),從而提高算法求解病態(tài)線性方程組的效率與精度;
2) 結(jié)合精細(xì)積分求解逆矩陣,可提高阻尼譜修正迭代法求解病態(tài)線性方程組的精度;
3) 精細(xì)積分阻尼譜修正迭代法,明顯可提高病態(tài)線性方程組的求解精度,且應(yīng)用于高維弱病態(tài)線性方程組,也可獲得較高的精度數(shù)值迭代解.
參考文獻(xiàn)
[1] 熊新, 吳洪濤, 于學(xué)華, 等. 工程車輛油氣懸架參數(shù)化建模與幅頻特性分析[J]. 振動(dòng)、測(cè)試與診斷, 2014,34(5): 926-931.
[2] 曲昕馨, 李桐林, 王飛. 基于數(shù)字圖像分割法的跨孔雷達(dá)走時(shí)層析成像[J]. 吉林大學(xué)學(xué)報(bào)(地球科學(xué)版), 2014,44(4): 1340-1347.
[3] 錢進(jìn), 吳時(shí)國(guó), 崔若飛, 等. 新驛礦區(qū)奧陶系灰?guī)r孔隙度預(yù)測(cè)方法[J]. 桂林理工大學(xué)學(xué)報(bào), 2012,32(1):43-47.
[4] 秦紅磊, 梁敏敏. 基于GNSS的高軌衛(wèi)星定位技術(shù)研究[J]. 空間科學(xué)學(xué)報(bào),2008,28(4): 316-325.
[5] 李曉妮,程濤. 分析微梁尺寸效應(yīng)的傳遞矩陣法[J]. 廣西科技大學(xué)學(xué)報(bào),2017,28(1):91-96.
[6] 余婧妮,向宇,李曉妮,等. 求解梁大變形問題的一種精細(xì)分析方法[J]. 廣西科技大學(xué)學(xué)報(bào),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)線性方程組的自動(dòng)控制步長(zhǎng)法[J]. 太原科技大學(xué)學(xué)報(bào), 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é)報(bào), 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] 王新洲,劉丁酉,張前勇,等. 譜修正迭代法及其在測(cè)量數(shù)據(jù)處理中的應(yīng)用[J]. 黑龍江工程學(xué)院學(xué)報(bào), 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é)科編輯:張玉鳳)