精细积分阻尼谱修正迭代法在病态线性方程组中的应用

覃有堂 林贤坤 覃柏英

摘 要：基于谱修正迭代法，引入阻尼因子和结合精细积分的矩阵求逆，提出了精细积分阻尼谱修正迭代法，并将其应用于病态线性方程组的求解.通过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

（学科编辑：张玉凤）