基于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

（責任编辑：罗小芬、黎   娅）