基于LU分解的阻尼譜修正迭代法在病態(tài)線性方程組中的應(yīng)用

莫春鵬 覃柏英

摘? 要：基于阻尼譜修正迭代法，結(jié)合矩陣LU分解和新數(shù)值迭代方式，提出了基于矩陣LU分解的阻尼譜修正迭代法，將其應(yīng)用于病態(tài)線性方程組的求解.采用經(jīng)典算例，探討矩陣LU分解和新數(shù)值迭代方式對(duì)阻尼譜修正迭代法求解病態(tài)線性方程組的性能影響.結(jié)果表明，矩陣LU分解和新數(shù)值迭代方式都可提高阻尼譜修正迭代法求解病態(tài)線性方程組的精度，且提出的算法可提高高維病態(tài)線性方程組求解的精度.

關(guān)鍵詞：LU分解;譜修正迭代法;病態(tài)矩陣;線性方程組

中圖分類號(hào)：O151.2;O241.6? ? ? ? ? ? ?DOI：10.16375/j.cnki.cn45-1395/t.2021.03.019

0? ? 引言

線性方程組是一類常見的數(shù)學(xué)模型，科學(xué)技術(shù)領(lǐng)域中的許多實(shí)際問題都可將其轉(zhuǎn)化為該模型[1-2].若系數(shù)矩陣的條件數(shù)較大，則該矩陣是病態(tài)的，對(duì)應(yīng)的線性方程組常稱病態(tài)線性方程組.對(duì)于病態(tài)線性方程組的求解，傳統(tǒng)算法常難以保證所求數(shù)值解的精度[3-5].因此，研究病態(tài)線性方程組的求解，以提高求解的準(zhǔn)確性和數(shù)值計(jì)算的效率，具有重要的意義和價(jià)值.

目前，研究者提出了多種數(shù)值算法求解病態(tài)線性方程組[6-15].王新洲等[14]提出了譜修正迭代法.該算法借助系數(shù)矩陣病態(tài)性的改善，使病態(tài)線性方程組求解的準(zhǔn)確性得以提高，但該算法的數(shù)值計(jì)算效率較低.由此，鄧興升等[16]提出了阻尼譜修正迭代法.該算法雖可提高數(shù)值計(jì)算的效率，但能否收斂至病態(tài)線性方程組的真解，與逆矩陣求解和數(shù)值迭代等方式有極大關(guān)系.

因此，本文改進(jìn)阻尼譜修正迭代法的數(shù)值迭代方式，并采用LU分解避免直接求解矩陣的逆矩陣，提出基于矩陣LU分解的新阻尼譜修正迭代法.采用4個(gè)經(jīng)典算例，探討該算法求解病態(tài)線性方程組的準(zhǔn)確性和效率.

1? ? 阻尼譜修正迭代法

將病態(tài)線性方程組表示為：

[AX=b]? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （1）

已知[m×n]實(shí)系數(shù)矩陣[A]和[m×1]實(shí)常數(shù)項(xiàng)[b].需求解該方程組的[n×1]解向量[X].

為了將系數(shù)矩陣統(tǒng)一為實(shí)對(duì)稱矩陣，可將上述方程組轉(zhuǎn)化為如下線性方程組：

[ATAX=ATb]? ? ? ? ? ? ? ? ? ? ? ? ? ?（2）

設(shè)[B=ATA]和[H=ATb]，則[B]為[n×n]的實(shí)對(duì)稱矩陣，[H]為[n×1]的實(shí)列向量，從而

[BX=H]? ? ? ? ? ? ? ? ? ? ? ? ? ? ? （3）

若采用最小二乘法（LSM）求解可獲得上述線性方程組的最小二乘解[X=B-1H].

為改善線性系數(shù)矩陣[B]的病態(tài)性，適當(dāng)選擇一個(gè)阻尼因子[α>0]，通過添加[αX]，可將線性方程組（3）轉(zhuǎn)化為：

[（B+αI）X=H+αX]? ? ? ? ? ? ? ? ? ? ? ? （4）

其中[I]為n階單位矩陣.

該線性方程組可采用如下迭代法數(shù)值求解而獲得線性方程組（1）和式（3）的解向量[X]：

[（B+αI）X（k）=H+αX（k-1）]? ? ? ? ? ? ? （5）

在實(shí)際應(yīng)用中，式（5）的迭代法常需轉(zhuǎn)化為如下的迭代算法求解：

[X（k）=（B+αI）-1H+αX（k-1）]? ? ? ? ? ? （6）

該算法稱為阻尼譜修正迭代法（DCCV）.

對(duì)式（6）也可構(gòu)造如下改進(jìn)迭代算法數(shù)值求解而獲得線性方程組（1）和式（3）的解向量[X]：

[R（k）=H-BX（k）D（k）=（B+αI）-1R（k）X（k+1）=X（k）+D（k）]? ? ? ? ? ? ? ?（7）

其中：[R]為誤差矩陣，[D]為對(duì)角矩陣.該算法稱為改進(jìn)阻尼譜修正迭代法（IDCCV）.

2? ? 基于LU分解的阻尼譜修正迭代法

由于數(shù)值計(jì)算存在的舍入誤差，矩陣的求逆常給病態(tài)線性方程組的準(zhǔn)確求解帶來較大影響.因此，為了避免式（6）和式（7）中對(duì)矩陣[B]+[αI]的直接求逆，可對(duì)其進(jìn)行LU分解產(chǎn)生一個(gè)下三角形矩陣[L]和一個(gè)上三角形矩陣[U]以及一個(gè)置換矩陣[P]，使之滿足[LU=P（B+αI）].

對(duì)于DCCV，由式（5），設(shè)[Y=UX（k）].有

[LUX（k）=PH+αPX（k-1）]

[LY=PH+αPX（k-1）=W]

其中[W]為結(jié)果矩陣.

由于[L]為下三角形矩陣，即[L]可逆，有[Y=L-1W].利用中間解[Y]，由[UX（k）=Y]，且[U]為上三角形矩陣，即[U]可逆，有[X（k）=U-1Y].

因[Li， i=1，Y1=W1]以及[X（k）n=Yn/Un， n]，因此，中間解[Y]和解[X（k）]可根據(jù)如下遞推式分別求出[Yi（i=2， …， n）]和[X（k）i（i=n-1， …， 1）]：

[Yi=Wi-j=1i-1Li， jYj]

[X（k）i=Yi/Ui， i-j=i+1nUi， jX（k）j/Ui， i]

該求解線性方程組（1）和式（3）的算法，稱為基于LU分解的阻尼譜修正迭代法（LUDCCV）.

對(duì)于IDCCV，由式（7），設(shè)[Z=UD（k）]，有

[LUD（k）=PR（k），LZ=PR（k）=V]

由于[L]為下三角形矩陣，即[L]可逆.從而[Z=L-1V]，[V]為列向量.利用中間解[Z]，由[UD（k）=Z]，且[U]為上三角形矩陣，即[U]可逆，從而[D（k）=U-1Z].

因[Li， i=1]，[Z1=V1]以及[D（k）n=Zn/Un， n]，因此，中間解[Z]和解[D（k）]可根據(jù)如下遞推式分別求出[Zi（i=2， …， n）]和[D（k）i（i=n-1， …， 1）]：

[Zi=Vi-j=1i-1Li， jZj]

[D（k）i=Zi/Ui， i-j=i+1nUi， jD（k）j/Ui， i]

該求解線性方程組（1）和式（3）的算法，稱為基于LU分解的改進(jìn)阻尼譜修正迭代法（LUIDCCV）.

3? ? 經(jīng)典算例

算例1? 給定線性方程組（1）如下的系數(shù)矩陣[A19×4]和常數(shù)項(xiàng)[b19×1]：

此時(shí)[X=[0.2， 1.5， 1.6， -2.8]T]為相應(yīng)線性方程組（1）的真解.由于矩陣[A]非對(duì)稱，將其轉(zhuǎn)化，取[B=ATA]，則其條件數(shù)為1.610 1E9.因此.相應(yīng)線性方程組（3）為病態(tài)線性方程組.

算例2 給定如下的系數(shù)矩陣[A18×7]和常數(shù)? ? ?項(xiàng)[b18×1].

此時(shí)[X=[0.2， 2.0， 1.5， -1.6， 4.8， 3.4， -2.1]T]為相應(yīng)線性方程組（1）的真解.由于矩陣[A]非對(duì)稱，將其轉(zhuǎn)化而取[B=ATA]，則其條件數(shù)為2.996 7E5.因此，相應(yīng)線性方程組（3）為病態(tài)線性方程組.

算例3 給定系數(shù)矩陣[A]和常數(shù)項(xiàng)[b]：

[A=（ai， j）n×n，ai， j=1， i≠j， 1+p2， i=j.]

[b=[b1， b2， …， bn]T，bi=j=1nai， j]

其中：[i， j=1， 2 ， …， n]，此時(shí)[X=[1， 1， …， 1]T]為相應(yīng)線性方程組（1）的真解.由于[A]對(duì)稱，不需轉(zhuǎn)化而直接取[B=A].當(dāng)[n=10]和[p=5.0×10-3]時(shí)，矩陣[A]的條件數(shù)為4.000 0E7，因此，相應(yīng)線性方程組（1）為病態(tài)線性方程組.

算例4 給定系數(shù)矩陣[A]和常數(shù)項(xiàng)[b]：

[A=（ai， j）n×n，ai， j=1/（i+j-1）]

[b=[b1， b2， …， bn]T，bi=j=1njai， j]

其中：[i， j=1， 2 ， …， n]，此時(shí)[X=[1， 2， …， n]T]為相應(yīng)線性方程組（1）的真解.由于[A]對(duì)稱，不需轉(zhuǎn)化而直接取[B=A].當(dāng)[n=8]時(shí)，矩陣[A]的條件數(shù)為? 1.525 8E10，因此，相應(yīng)線性方程組（1）為病態(tài)線性方程組.

4? ? 算法的性能分析

利用算法LSM、DCCV、IDCCV、LUDCCV和LUIDCCV求解病態(tài)線性方程組（1）的數(shù)值解為[X]，可計(jì)算[Eb=AX]和[Ex=X-X]，從而可定義為絕對(duì)誤差的殘差[Eb=║Eb║2]以評(píng)價(jià)數(shù)值解[X]對(duì)方程組（1）的滿足程度，以及可定義為均方誤差的殘差

[Ex=║Ex║22/n]和相對(duì)誤差的殘差[E∞=║Ex║∞/X∞]以評(píng)價(jià)數(shù)值解[X]偏離真實(shí)解[X]的程度.因此，利用殘差[Eb、Ex、][E∞]可評(píng)價(jià)上述5種算法數(shù)值求解病態(tài)線性方程組的準(zhǔn)確性，以及比較與分析上述5種算法數(shù)值求解病態(tài)線性方程組的性能.

現(xiàn)采用5種算法LSM、DCCV、IDCCV、LUDCCV和LUIDCCV分別數(shù)值求解上述4個(gè)算例，結(jié)果分別如表1—表4所示.其中[α]的值在4個(gè)算例中分別為[0.280、0.089、4.000×10-14]和[5.000×10-12]，算法的迭代次數(shù)用F表示.

由表1—表4可知，5種算法都可實(shí)現(xiàn)上述4個(gè)算例的病態(tài)線性方程組較高精度求解，且LUIDCCV、IDCCV、LUDCCV和DCCV都優(yōu)于LSM.同時(shí)，LUDCCV和LUIDCCV分別優(yōu)于DCCV和IDCCV.這說明當(dāng)線性方程組病態(tài)時(shí)，采用LU分解避免系數(shù)矩陣直接求逆的LUDCCV和LUIDCCV優(yōu)于系數(shù)矩陣直接求逆的DCCV和IDCCV，因此，LU分解可減弱矩陣求逆數(shù)值計(jì)算過程中舍入誤差存在的影響，從而提高算法求解病態(tài)線性方程組數(shù)值解的精度.IDCCV和LUIDCCV分別優(yōu)于DCCV和LUDCCV，這也說明采用新數(shù)值迭代方式，也可提高算法求解病態(tài)線性方程組數(shù)值解的精度.因此，迭代算法LUIDCCV相比他4種迭代算法最優(yōu)，可明顯提高阻尼譜修正迭代法求解病態(tài)線性方程組數(shù)值解的精度，使所求數(shù)值解[X]更高程度滿足線性方程組（1）和接近線性方程組（1）的真實(shí)解[X].

5? ? 高維問題中的應(yīng)用

為了更充分說明LUDCCV求解病態(tài)線性方程組的性能，現(xiàn)將其應(yīng)用于高維問題，即將LUDCCV應(yīng)用于算例3和算例4的高維病態(tài)線性方程組的求解.

對(duì)于算例3的病態(tài)線性方程組，分別取其維數(shù)[n=100、200、500、1 000、2 000、3 000、4 000]，參數(shù)[p]取[5.0×10-6]時(shí)，矩陣[A]的條件數(shù)[κ]分別如? 表5所示，阻尼因子[α] 取1.000.采用LUIDCCV分別求解，殘差[Eb、Ex、E∞]的值分別如表5所示.

對(duì)于算例4的病態(tài)線性方程組，分別取其維數(shù)[n=100、200、500、1 000、2 000、3 000、4 000]，矩陣A的條件數(shù)[κ]分別如表6所示，阻尼因子[α]取[5.000×10-12]，殘差[Eb、Ex、E∞]的值分別如表6所示.

由表5可知，對(duì)于高維弱病態(tài)的線性方程組，LUDCCV可獲得其較高精度的數(shù)值解，該數(shù)值解較高程度滿足線性方程組（1），同時(shí)也比較接近線性方程組（1）的真實(shí)解.但由表6可知，將LUDCCV應(yīng)用于高維強(qiáng)病態(tài)的線性方程組時(shí)，所求的數(shù)值解其精度有待提高.

為了提高LUIDCCV求解高維強(qiáng)病態(tài)線性方程組的精度，對(duì)于線性方程組（3）的常數(shù)項(xiàng)，若其各元素之間的數(shù)值相差較大，可對(duì)常數(shù)項(xiàng)[H]先歸一化，再求解線性方程組（3），從而提高LUDCCV求解高維強(qiáng)病態(tài)線性方程組的精度，使其數(shù)值解更滿足線性方程組（1）和接近其真實(shí)解.設(shè)[H=[h1， h2， …， hn]T].由此，構(gòu)造[C=diag（1/h1， 1/h2， …， 1/hn）]，可將線性方程組（3）轉(zhuǎn)化為：

[CBX=En×1]

其中，[En×1=CH=[1， 1， …， 1]T]，即實(shí)現(xiàn)常數(shù)項(xiàng)的歸一化.由此，再采用下式進(jìn)行數(shù)值迭代求解：

[（CB+αI）X（k）=En×1+αX（k-1）]

根據(jù)上述方法，對(duì)算例4的強(qiáng)病態(tài)線性方程組進(jìn)行歸一化處理后，再采用LUIDCCV數(shù)值迭代求解強(qiáng)病態(tài)線性方程組，結(jié)果如表7所示.

比較表6與表7的結(jié)果可知，對(duì)常數(shù)項(xiàng)[H]先歸一化后再采用LUIDCCV求解線性方程組（3），高維強(qiáng)病態(tài)線性方程組求解的精度提高明顯，從而LUIDCCV求解病態(tài)線性方程組時(shí)可使其數(shù)值解更滿足線性方程組（1），接近其真實(shí)解[X].

6? ? 結(jié)語

為提高譜修正迭代法的計(jì)算效率，使阻尼譜修正迭代法收斂至病態(tài)線性方程組的真解并提高其收斂速度，實(shí)現(xiàn)病態(tài)線性方程組高精度和高效率的求解，本文改進(jìn)了阻尼譜修正迭代法的數(shù)值迭代方式，并采用矩陣的LU分解避免直接求解系數(shù)矩陣的逆矩陣，由此提出基于矩陣LU分解的新阻尼譜修正迭代法.采用4個(gè)經(jīng)典算例，將基于LU分解的改進(jìn)阻尼譜修正迭代法應(yīng)用于高維病態(tài)線性方程組，探討了矩陣的LU分解和新的迭代方式對(duì)阻尼譜修正迭代法求解病態(tài)線性方程組的性能影響.可得到如下結(jié)論：

1）結(jié)合矩陣的LU分解或新的迭代方式都可提高阻尼譜修正迭代法求解病態(tài)線性方程組的精度，但新的迭代方式對(duì)解精度的提高效果更明顯，且同時(shí)結(jié)合矩陣的LU分解和新迭代方式的阻尼譜修正迭代法精度最高;

2）相比于其他4種算法，結(jié)合矩陣LU分解和新迭代方式的阻尼譜修正迭代法求解病態(tài)線性方程組的精度可明顯提高，且將其應(yīng)用于高維弱病態(tài)線性方程組求解，其數(shù)值解的精度也較高;

3）若常數(shù)項(xiàng)[H]的各元素間的數(shù)值相差較大且非零，采用歸一化法對(duì)常數(shù)項(xiàng)[H]進(jìn)行處理后再求解，可提高基于矩陣LU分解的新阻尼譜修正迭代法求解病態(tài)線性方程組的精度.

參考文獻(xiàn)

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

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

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

[4]? ? ?REICHEL? L，RODRIGUEZ? G. Old and new parameter choice rules for discrete ill-posed problems[J].Numerical Algorithms，2013，63（1）：65-87.

[5]? ? ?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.

[6]? ? ?覃有堂，林賢坤，覃柏英.精細(xì)積分阻尼譜修正迭代法在病態(tài)線性方程組中的應(yīng)用[J].廣西科技大學(xué)學(xué)報(bào)，2017，28（3）：1-8.

[7]? ? ?SMOKTUNOWICZ Alicja，SMOKTUNOWICZ Agata. Iterative refinement techniques for solving block linear systems of equations[J]. Applied Numerical Mathematics，2013，67：220-229.

[8]? ? ?NEUMAN? A，REICHEL? L，SADOK? H. Implementations of range restricted iterative methods for linear discrete ill-posed problems[J]. Linear Algebra and Its Applications，2012，436（10）：3974-3990.

[9]? ? ?VILOCHE BAZ?N F S，CUNHA M C C，BORGES L S. Extension of GKB-FP algorithm to large-scale general-form Tikhonov regularization[J]. Numerical Linear Algebra with Applications，2014，21（3）：316-339.

[10]? ?WU X Y. An effective predictor-corrector process for large scale linear system of equations[J]. Applied Mathematics and Computation，2006，180（1）：160-166.

[11]? ?WU X Y. Note on predictor-corrector process for ill-conditioned linear system of equations[J]. Applied Mathematics and Computation，2006，183（1）：596-602.

[12]? ?WU X Y，F(xiàn)ANG Y L. Wilkinson's iterative refinement of solution with automatic step-size control for linear system of equations[J]. Applied Mathematics and Computation，2007，193（2）：506-513.

[13]? ?RILEY J D. Solving systems of linear equations with a positive definite，symmetric，but possibly ill-conditioned matrix[J].Mathematical Tables and other Aids to Computation，1955，9（51）：96-101.

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

[15]? ?劉麗華.解二維Poisson方程離散化線性方程組的新型二次PEk方法[J].廣西科技大學(xué)學(xué)報(bào)，2016，27（2）：100-103.

[16]? ?鄧興升，孫虹虹. 自適應(yīng)譜修正LU分解法解算高病態(tài)法方程[J].大地測量與地球動(dòng)力學(xué)，2014，34（6）：135-139.

Iteration method by correcting characteristic value with damping

factor based on LU decomposition for ill-conditioned system of

linear equations

MO Chunpeng， QIN Boying*

（College of Science， Guangxi University of Science and Technology， Liuzhou 545006， China）

Abstract：Based on the iteration method by correcting characteristic value with damping factor，? ? ? ?combined with LU matrix decomposition and a new numerical iteration method， this paper proposes an improved iteration method by correcting characteristic value with damping factor based on LU matrix? decomposition， which is applied to solve ill-conditioned system of linear equations.Some classical? ? ? examples are used to investigate the influence of LU matrix decomposition and the new numerical? ? ? ? iteration method on the performance of the algorithm for solving ill-conditioned system of linear? ? ?equations.The results show that both the LU matrix decomposition and the new numerical iteration? ?method can improve the accuracy of the algorithm for solving ill-conditioned system of linear? ? ? ? ?equations， and the proposed algorithm can also improve the accuracy of solving high-dimensional? ? ? ?ill-conditioned system of linear equations.

Key words： LU decomposition; iteration method by correcting characteristic value; ill-conditioned? ? matrix; system of linear equations

（責(zé)任編輯：羅小芬、黎? ?婭）