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

覃有堂1,林贤坤*2,覃柏英2

(1.柳州市第一职业技术学校,广西 柳州 545616;2.广西科技大学 汽车与交通学院,广西 柳州 545006)

摘 要:基于谱修正迭代法,引入阻尼因子和结合精细积分的矩阵求逆,提出了精细积分阻尼谱修正迭代法,并将其应用于病态线性方程组的求解.通过4个经典算例,探讨矩阵的精细积分求逆对阻尼谱修正迭代法求解病态线性方程组的影响.结果表明,矩阵的精细积分求逆可提高阻尼谱修正迭代法求解病态线性方程组的精度.再通过两个经典算例,探讨了精细积分阻尼谱修正迭代法在高维病态线性方程组中的应用.结果表明,该算法也可提高高维病态线性方程组求解的精度.

关键词:精细积分;阻尼因子;谱修正迭代法

0 引言

本文考虑如下的病态线性方程组问题:

其中: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的实列向量,从而

因BT=B,所以实矩阵B为对称矩阵.若B为正定矩阵,则线性方程组(2)存在唯一解.

根据最小二乘原理,线性方程组(2)的最小二乘解为:

求解线性方程组(1)的算法称为最小二乘法,所求的解相应称为最小二乘解.

当线性方程组(1)病态时,其系数矩阵A的条件数往往很大,则矩阵B的条件数也会很大,此时采用最小二乘法(3)所求的病态线性方程组(1)的最小二乘解精度往往较低.

为了改善线性方程组(1)的病态性,可在式(2)等号两边加上aX,I为n阶单位阵,有∶

其中α称为阻尼因子.设N=B+αI,则NT=N,因此N为实对称矩阵.设Y为n×1的任意非零实向量,则YTNY=YT(B+αI)Y=YTBY+αYTY≥αYTY>0,从而实矩阵N为正定矩阵.因此线性方程组(4)也存在唯一解,且该解也是线性方程组(2)的唯一解.

对于线性方程组(2),要采用如下的阻尼谱修正迭代法进行数值迭代求解:

1 阻尼谱修正迭代法的收敛性

对于对称且正定的实矩阵B,其存在n个实特征值,从大到小排序为λ1,λ2,…,λn.为了叙述的方便,设N=B+αI和 P=(B+αI)-1,探讨阻尼谱修正迭代法的收敛性.

引理1 矩阵N的n个特征值为 λ1+α,λ2+α,…,λn+α.

证明 设λ为矩阵B的任意一特征值,其对应的特征向量为ξ,则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 的 n 个特征值为(λ1+α)-1,(λ2+α)-1,…,(λn+α)-1,其最大值(λn+α)-1.

又因为矩阵Q为正交矩阵,从而

引理3 当实矩阵B正定时,对任意的阻尼因子α满足α>0,由n阶方阵P.产生的矩阵序列αP,α2P2,…,αkPk,…收敛,

证明 由矩阵范数的相容性定义,有:

因为矩阵 B 正定,从而其最小特征值 λn>0.于是 α>0 时,有‖αP‖2=α/(λn+α)<1,有:

从而序列 αP,α2P2,…,αkPk,…收敛,

引理4 当矩阵B正定时,对任意的阻尼因子α满足α>0,由n阶方阵P构造的矩阵幂级数I+αP+α2P2+…+αkPk+…收敛,且其和为(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+α2P2+…+αkPk,有:

在等式两边右乘以矩阵 I-αP的逆矩阵(I-αP)-1,因 P(I-αP)-1=B-1,则:

定理1 当矩阵B正定时,对任意的阻尼因子α,满足α>0和任意的初值,有阻尼谱修正迭代法=P[H+α]收敛,且=B-1H.

证明 由有:

由引理3和引理4知,当矩阵B正定时,对任意的阻尼因子α满足α>0和任意的初值,有:

由以上证明知,若系数矩阵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为负定矩阵.

为负定矩阵,则即矩阵N的逆矩阵

取定一小步长 Δt,设 t=3mΔt,m=0,1,…,则

设 Fm=F(3mΔt),Rm=,m=0,1,2,…,则可采用如下迭代式计算 P:

由上式知,当m足够大时可得到较精确的P.为了利用上式计算P,需给定初始值F0和Rm.对 eMs进行 Taylor展开,有 eMs≈I+Ms+(Ms)2/2+(Ms)3/6+(Ms)4/24,从而代式计算Tm

3 经典算例

为了探讨矩阵的精细积分求逆对阻尼谱修正迭代法求解病态线性方程组的性能影响,以及精细积分阻尼谱修正迭代法求解病态线性方程组的精度与效率,现给出4个经典的算例.

算例1系数矩阵A和常数项b分别如下:

此时线性方程组(1)的真解为 X4×1=[0.2,1.5,1.6,-2.8]T.因矩阵 A 非对称,取 B=ATA,其条件数为 1.6101E9,因此线性方程组(1)对应的线性方程组(2)为弱病态线性方程组.

算例2系数矩阵A和常数项b分别如下:

此时线性方程组(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分别如下:

若常数项b取则线性方程组(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分别如下:

若常数项b取其中 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 和.由此可设,则可定义评价数值迭代解精度的如下3个残差:

残差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 3种算法求解算例1的结果
Tab.1 The results of three algorithms for the first example

表2 3种算法求解算例2的结果
Tab.2 The results of three algorithms for the second example

表3 3种算法求解算例3的结果
Tab.3 The results of three algorithms for the third example

表4 3种算法求解算例4的结果
Tab.4 The results of three algorithms for the fourth example

由表1—表4知,3种算法都可实现4个算例较高精度的求解.但DCCV和PIM-DCCV都优于LSM,这说明DCCV有利于病态线性方程组求解精度的提高.同时,PIM-DCCV也明显优于DCCV,这也说明精细积分可更精确求解逆矩阵,从而提高算法DCCV的精度.因此,算法PIM-IDCCV可使数值迭代解Xˆ更满足线性方程组(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所示.当取时,阻尼因子 α=1,当取时,阻尼因子 α=5.4×10-14,采用PIM-DCCV分别求解,结果分别如表5和表6所示.

表5 当bi=ai,j时,算例 3 对应的结果
Tab.5 The results of the third example when bi=ai,j

表6 当bi=jai,j时,算例 3对应的结果
Tab.6 The results of the third example when bi=jai,j

对于算例 4 的病态线性方程组,分别取其维数 n=100,200,500,1 000,2 000,3 000,4 000,系数矩阵 A对应的条件数分别如表7所示.当常数项b取时,α=5.0×10-2,当常数项 b 取时,α=5.0×10-12,采用PIM-DCCV分别求解,结果分别如表7和表8所示.

表7 当bi=ai,j时,算例 4 对应的结果
Tab.7 The results of the fourth example when bi=ai,j

表8 当bi=jai,j时,算例 4对应的结果
Tab.8 The results of the fourth example when bi=jai,j

采用精细积分改进阻尼谱修正迭代法对高维病态线性方程组进行数值求解,由表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.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.

[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

中图分类号:O241.6;O151.2

文献标志码:A

(学科编辑:张玉凤)

文章编号:2095-7335(2017)03-0001-08

DOI:10.16375/j.cnki.cn45-1395/t.2017.03.001

收稿日期:2017-03-19

基金项目:广西自然科学基金(2012GXNSFAA053208);广西教育厅科研项目(200103YB105)资助.

* 通信作者:林贤坤,博士,教授,研究方向:计算流体力学,E-mail: linxl0209@163.com.