基于LU分解的阻尼譜修正迭代法在病態線性方程組中的應用

莫春鵬 覃柏英

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

關鍵詞：LU分解;譜修正迭代法;病態矩陣;線性方程組

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

0? ? 引言

線性方程組是一類常見的數學模型，科學技術領域中的許多實際問題都可將其轉化為該模型[1-2].若系數矩陣的條件數較大，則該矩陣是病態的，對應的線性方程組常稱病態線性方程組.對于病態線性方程組的求解，傳統算法常難以保證所求數值解的精度[3-5].因此，研究病態線性方程組的求解，以提高求解的準確性和數值計算的效率，具有重要的意義和價值.

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

因此，本文改進阻尼譜修正迭代法的數值迭代方式，并采用LU分解避免直接求解矩陣的逆矩陣，提出基于矩陣LU分解的新阻尼譜修正迭代法.采用4個經典算例，探討該算法求解病態線性方程組的準確性和效率.

1? ? 阻尼譜修正迭代法

將病態線性方程組表示為：

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

已知[m×n]實系數矩陣[A]和[m×1]實常數項[b].需求解該方程組的[n×1]解向量[X].

為了將系數矩陣統一為實對稱矩陣，可將上述方程組轉化為如下線性方程組：

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

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

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

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

為改善線性系數矩陣[B]的病態性，適當選擇一個阻尼因子[α>0]，通過添加[αX]，可將線性方程組（3）轉化為：

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

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

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

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

在實際應用中，式（5）的迭代法常需轉化為如下的迭代算法求解：

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

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

對式（6）也可構造如下改進迭代算法數值求解而獲得線性方程組（1）和式（3）的解向量[X]：

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

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

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

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

對于DCCV，由式（5），設[Y=UX（k）].有

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

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

其中[W]為結果矩陣.

由于[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）]可根據如下遞推式分別求出[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）.

對于IDCCV，由式（7），設[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）]可根據如下遞推式分別求出[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分解的改進阻尼譜修正迭代法（LUIDCCV）.

3? ? 經典算例

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

此時[X=[0.2， 1.5， 1.6， -2.8]T]為相應線性方程組（1）的真解.由于矩陣[A]非對稱，將其轉化，取[B=ATA]，則其條件數為1.610 1E9.因此.相應線性方程組（3）為病態線性方程組.

算例2 給定如下的系數矩陣[A18×7]和常數? ? ?項[b18×1].

此時[X=[0.2， 2.0， 1.5， -1.6， 4.8， 3.4， -2.1]T]為相應線性方程組（1）的真解.由于矩陣[A]非對稱，將其轉化而取[B=ATA]，則其條件數為2.996 7E5.因此，相應線性方程組（3）為病態線性方程組.

算例3 給定系數矩陣[A]和常數項[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]，此時[X=[1， 1， …， 1]T]為相應線性方程組（1）的真解.由于[A]對稱，不需轉化而直接取[B=A].當[n=10]和[p=5.0×10-3]時，矩陣[A]的條件數為4.000 0E7，因此，相應線性方程組（1）為病態線性方程組.

算例4 給定系數矩陣[A]和常數項[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]，此時[X=[1， 2， …， n]T]為相應線性方程組（1）的真解.由于[A]對稱，不需轉化而直接取[B=A].當[n=8]時，矩陣[A]的條件數為? 1.525 8E10，因此，相應線性方程組（1）為病態線性方程組.

4? ? 算法的性能分析

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

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

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

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

5? ? 高維問題中的應用

為了更充分說明LUDCCV求解病態線性方程組的性能，現將其應用于高維問題，即將LUDCCV應用于算例3和算例4的高維病態線性方程組的求解.

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

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

由表5可知，對于高維弱病態的線性方程組，LUDCCV可獲得其較高精度的數值解，該數值解較高程度滿足線性方程組（1），同時也比較接近線性方程組（1）的真實解.但由表6可知，將LUDCCV應用于高維強病態的線性方程組時，所求的數值解其精度有待提高.

為了提高LUIDCCV求解高維強病態線性方程組的精度，對于線性方程組（3）的常數項，若其各元素之間的數值相差較大，可對常數項[H]先歸一化，再求解線性方程組（3），從而提高LUDCCV求解高維強病態線性方程組的精度，使其數值解更滿足線性方程組（1）和接近其真實解.設[H=[h1， h2， …， hn]T].由此，構造[C=diag（1/h1， 1/h2， …， 1/hn）]，可將線性方程組（3）轉化為：

[CBX=En×1]

其中，[En×1=CH=[1， 1， …， 1]T]，即實現常數項的歸一化.由此，再采用下式進行數值迭代求解：

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

根據上述方法，對算例4的強病態線性方程組進行歸一化處理后，再采用LUIDCCV數值迭代求解強病態線性方程組，結果如表7所示.

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

6? ? 結語

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

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

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

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

參考文獻

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

[2]? ? ?錢進，吳時國，崔若飛，等. 新驛礦區奧陶系灰巖孔隙度預測方法[J]. 桂林理工大學學報，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]? ? ?覃有堂，林賢坤，覃柏英.精細積分阻尼譜修正迭代法在病態線性方程組中的應用[J].廣西科技大學學報，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，FANG 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]? ?王新洲，劉丁酉，張前勇，等. 譜修正迭代法及其在測量數據處理中的應用[J]. 黑龍江工程學院學報，2001（2）：3-6.

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

[16]? ?鄧興升，孫虹虹. 自適應譜修正LU分解法解算高病態法方程[J].大地測量與地球動力學，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

（責任編輯：羅小芬、黎? ?婭）