陈凤坤，雷秀仁

(华南理工大学 数学学院，广东 广州 510640)

一种求解病态复线性方程组的混合算法

陈凤坤，雷秀仁

(华南理工大学 数学学院，广东 广州 510640)

病态复线性方程的求解是现代应用数学和很多工程应用面临的难题，用一般算法进行求解时，得到的误差较大，因此在一些高精度的工程应用上，其结果往往不是特别理想。而随着科技的发展，现代很多工程应用对数据具有越来越高的精度要求(尤其是国家航天航空)，因此一个能求解病态复线性方程组的高精度算法是很有必要的。从病态复线性方程组求解的特点出发，对模拟退火法进行改进，并将其全局的收敛能力与双共轭梯度法的高精度求解能力结合起来，提出了一种BCG-SA混合算法。数据实验表明，模拟退火法能对双共轭梯度法求出的解进行微调动，帮助双共轭梯度法在概率意义上跳出局部极小值点，从而提高求解精度。

病态复线性方程组；模拟退火算法；双共轭梯度法；混合算法；希尔伯特矩阵

0 引 言

求解复线性方程组是很多工程实际应用常遇到的问题，当方程组的规模较大或者系数矩阵的条件数很大时，系数矩阵易呈现病态特性。当前求解病态方程组效果较好的有预处理ILUCG[1]、主元加权迭代法[2]、新Jacobi迭代法[3]、粒子群算法[4]等，但用它们求解病态的复线性方程组的效果不是很理想。目前求解病态复线性方程组效果较好的有双共轭梯度法，但其存在不收敛或者收敛速度慢的潜在问题，且在搜索的过程中可能陷入局部极小值。

现代很多高科技的工程应用对数据具有越来越高的精度要求(尤其是国家航天航空)，因此一个能求解病态复线性方程组的高精度算法是很有必要的。为此，在分析研究双共轭梯度算法和模拟退火法的基础上，根据复线性方程组求解的特点对模拟退火法进行了适当改进，将双共轭梯度法高精度求解的特性和模拟退火法的全局搜索能力有机结合，提出了BCG-SA混合算法。实验结果表明，该算法虽然增加了迭代次数，但明显提高了计算精度，对于一些有较高精度要求的工程应用具有重要的现实意义。

1 双共轭梯度法(BCG)原理

双共轭梯度法是C.Lanczos[5]提出的一种用于求解一般复线性方程组的方法，具有计算量少、收敛速度快等优点。在求解规模较大或者条件数较大的复线性方程组时往往效果较好，是目前求解病态方程组较为常用的方法之一。

记(·，·)为通常复内积。若x=(x1,x2,…,xn)T，y=(y1,y2,…,yn)T，xi,yi∈C,i=1,2,…,n，则

(1)

并记

(2)

考虑如下形式的线性方程组：

Ax=b

(3)

其中，A∈Cn×n，x,b∈Cn。

则具体双共轭梯度算法[6](算法1)如下：

Step2：计算

并置k=k+1。

Step3：若‖xk+1-xk‖2<ε或k>kmax，则算法终止，转到Step4,否则转入Step2继续计算。

Step4：输出数值解xk+1和迭代次数k。

2 模拟退火法(SA)原理

模拟退火算法是由N.Metropolis[7]等于1953年提出的，它是基于Monte-Carlo迭代求解策略的一种随机寻优算法，其出发点是基于物理中固体物质的退火过程与一般组合优化问题之间的相似性。模拟退火算法从某一较高初温出发，伴随温度参数的不断下降，结合概率突跳特性在解空间中随机寻找目标函数的全局最优解，是一种通用的优化算法，具有概率性地跳出局部最优解,并最终趋于全局最优的特性。1982年Kinkpatrick等将内能E模拟为目标函数f，温度T演化成控制参数t，首次将模拟退火算法用于解决组合优化问题。

通过变分原理将式(3)转化为求解如下无约束条件的多元非线性方程组。

(4)

则模拟退火法应用于求解病态线性方程组的算法(算法2)[8-9]如下：

Step2：给x一个扰动x'=x+Δx，其中Δx是整个实轴上服从N(0，σ2)的随机变量，或者限制当前解的一个邻域Δx=scale×(rand-0.5)，其中rand是0到1之间的随机数，计算相应的目标函数值E(x')，得到ΔE=E(x')-E(x)。

Step3：若ΔE≤0，则令x=x'；否则若exp(-ΔE/t)>rand，令x=x',i=i+1。

Step4：若i

Step5：若t

3 混合算法(BCG-SA)

双共轭梯度法是一种收敛速度快，且精度较高的算法，在实践中表明双共轭梯度法的效果非常好，但由于其在搜索过程中可能陷入局部极小值，所以在一些对于数据精度有较高要求的实际应用中，需对双共轭梯度法进行一些改进，以帮助它跳过局部极小值点。

模拟退火算法[10]具有良好的全局优化性，虽然它在解决病态方程组方面效果不是特别好，但它能够以随机搜索技术从概率意义上找出目标函数的全局最小值点。

当双共轭梯度算法达到较高精度后会出现收敛速度慢或者不再继续收敛的情况，这表明算法陷入了局部极小值，而模拟退火在搜索过程引入了随机因素，以一定的概率来接受一个比当前解要差的解，因此有可能跳出这个局部最优解。所以为了得到更高精度的解，可以在双共轭梯度算法陷入局部极小值时，引入模拟退火法来帮助其跳出局部极小值点。由此受到启发，将双共轭梯度法和模拟退火法各自的优点结合起来，提出了混合算法。具体的混合算法(算法3)[11-12]如下：

Step1：输入矩阵A、右端向量b、误差精度ε、ε0(ε0>ε),最大迭代次数kmax，SA算法最大调用次数Nmax(不宜设置过大，可取6～10)，置算法迭代次数k=0，SA算法调用次数N=0。

Step3：计算

并置k=k+1。

Step4：若‖xk+1-xk‖2<ε或k>kmax或N>Nmax，则算法终止，转到Step5；若‖xk+1-xk‖<ε0，调用SA算法，置N=N+1,转入Step2，否则转入Step3继续计算。

Step5：输出数值解xk+1和迭代次数k。

4 实验及分析

为了说明算法BCG-SA的有效性和优越性，进行了数据模拟仿真实验。实验所用到的矩阵均从经典的病态矩阵演变而来，并且已经证明演变后的矩阵一方面保留了原矩阵的病态性，另一方面拥有复线性，可以代表病态复线性方程组。

将算法BCG-SA与在求解病态方程组上性能优越的共轭梯度法(CG)、预处理共轭梯度法(IC-CG)、双共轭梯度法(BCG)、预处理双共轭梯度法(IC-CBCG)、模拟退火法(SA)进行比较。数值实验结果表明，BCG-SA是有效的且在精度上更优越。需要指出的是进行数值实验时，也将BCG-SA混合算法与蚁群算法、粒子群算法、神经网络算法进行了比较，但这些群智能算法的实验效果并不好，随着矩阵阶数的增加，误差太大，甚至出现了解的失真。

算例1：考查Hilbert演变矩阵,对于具有病态方程组代表性的Hilbert矩阵B，其中

(5)

令演变矩阵A=B+B*j(j是复数单位，j2=-1),并取方程组右端项b为:

(6)

算例2:考查Pascal演变矩阵，对于具有病态方程组代表性的Pascal矩阵C[13-14]，其中

(7)

(8)

当矩阵的阶数n=20时，Hilbert演变矩阵的条件数为1.453 609×1018，Pascal演变矩阵的条件数为8.077 406×1017，已是高度病态的矩阵。实验时均取初始解x0为零向量，初始温度Tmax=1 000，Lmax=n(n为矩阵A的阶)，最大的试探次数Nmax=8，最低温度Tmin=1×10-4，邻域范围scale=1×10-2，温度下降因子α=0.9。并在算例1取误差精度ε=1×10-10，ε0=1×10-12；算例2取误差精度ε=1×10-11，ε0=1×10-12，实验结果如图1和图2所示(其中IC-CG在求解算例2时失真)。

图1 Hilbert演变矩阵在不同阶数下算法求解误差的比较

图2 Pascal演变矩阵在不同阶数下算法求解误差的比较

数据实验表明：

(1)与BCG-SA混合算相比，蚁群算法、粒子群算法等智能算法在求解高阶病态方程组上效果并不理想。

(2)相比于算例1，各类算法在求解算例2时能得到更高的精度。

(3)在算例2中，对CG、BCG进行预处理的效果并不理想。

(4)相对于其他算法，混合算法能得到较高精度的解，从而表明引入模拟退火法对于双共轭梯度法的改进是有效的，它能够帮助双共轭梯度法在概率意义上跳过局部极小值点，从而得到较高精度的解。

5 结束语

为了提高病态复线性方程组的求解精度，在研究分析病态复线性方程的特性和求解方法的基础上，提出了一种高精度混合算法。该算法是在对模拟退火法进行适当改进的基础上，将双共轭梯度法高精度求解的能力和模拟退火法的全局搜索能力有机结合。实验结果表明，提出的混合算法，虽然增加了计算的迭代次数，但明显提高了计算精度，对于一些对数据有较高精度要求的工程应用具有重要的现实意义。下一步的工作是进一步优化模拟退火算法的参数。

[1] 于春肖，苑润浩，穆运峰.新预处理ILUCG法求解稀疏病态线性方程组[J].数值计算与计算机应用，2014，35(1)：21-27.

[2] 唐 丽,李鹏飞.主元加权迭代法求解病态线性方程组[J].科学技术与工程，2012，12(2)：381-383.

[3] 孔祥强.病态线性方程组新的Jacobi迭代解法[J].大学数学，2013,29(5)：50-54.

[4] 贺天宇,李国望.病态线性方程组的粒子群算法[J].科技资讯，2012(8)：15.

[5] 张永杰,孙 泰.大型复线性方程组预处理双共轭梯度法[J].计算机工程与应用,2007，43(36)：19-20.

[6]GarciaN.ParallelpowerflowsolutionsusingabigconjugategradientalgorithmandaNewtonmethod:aGPU-basedapproach[C]//Power&energysocietygeneralmeeting.[s.l.]:[s.n.],2010:1-4.

[7]SteinbrunnM,MoerkotteG,KemperA.Heuristicandrandomizedoptimizationforthejoinorderingproblem[J].TheVLDBJournal,1997,6(3):191-208.

[8]KeikhaMM.Improvedsimulatedannealingusingmomentumterms[C]//Secondinternationalconferenceonintelligentsystems,modellingandsimulation.[s.l.]:[s.n.],2011:44-48.

[9]WangShengli,ZuoXingquan,LiuXueqing,etal.Solvingdynamicdoublerowlayoutproblemviacombiningsimulatedannealingandmathematicalprograming[J].AppliedSoftComputing,2015,37:303-310.

[10]LinZhenhai.Missionplanningforelectromagneticenvironmentmonitorssatellitebasedonsimulatedannealingalgorithm[C]//28thCanadianconferenceonelectricalandcomputerengineering.[s.l.]:IEEE,2015:530-535.

[11] 郑洲顺，黄光辉，杨晓辉.求解病态线性方程组的混合算法[J].贵州工业大学学报:自然科学版，2008,37(3)：12-15.

[12]BellioR,CeschiaS,GasperoLD,etal.Featuer-basedtuningofsimulatedannealingappliedtothecurriculum-basedcoursetimetablingproblem[J].Computers&OperationResearch,2016,65:83-92.

[13]SrisangngamP,ChivapreechaS,DejhanK.Evenorderbi-quaddigitafilterusingpascalmatrix[C]//Internationalsymposiumoncommunicationsandinformationtechnologies.[s.l.]:[s.n.],2008:327-330.

[14] 杨胜良.Pascal三角形与Pascal矩阵[J].数学的实践与认识，2003，33(2)：96-100.

A Hybrid Algorithm of Ill-conditioned Complex Linear Equations

CHEN Feng-kun，LEI Xiu-ren

(School of Mathematics，South China University of Technology，Guangzhou 510640，China)

Solving ill-conditioned complex linear equations is difficult in modern applied mathematics and many engineering application,and it is easy to produce significant error and bad result for general algorithms in some high-accuracy application with a usual algorithm.With the development of science and technology,there is more and more restrictions on data accuracy for modern industry especially space flight and aviation.Therefore,it is necessary and impending to find a high accuracy algorithm for ill-conditioned complex linear equations.According to the characteristics of ill-conditioned complex linear equations,a hybrid algorithm of BCG-SA improved with simulated annealing algorithm has been proposed with the advantages of global convergence and high precision solution for bi-conjugate gradient algorithm.The experimental results show that the hybrid algorithm has promoted the precision of solution for bi-conjugate gradient algorithm which can jump out of the neighborhoods of local minimum points in the sense of probability.

ill-conditioned complex linear equations;simulated annealing algorithm;bi-conjugate gradient algorithm;hybrid algorithm;Hilbert matrix

2016-06-02

2016-09-09 网络出版时间：2017-03-13

国家基金数学天元基金(B13-B5071130)；国家教育部高校博士点基金(B13-C7070170)

陈凤坤(1990-)，男，硕士研究生，研究方向为群智能算法；雷秀仁，副教授，通信作者，研究方向为科学与工程计算。

http://kns.cnki.net/kcms/detail/61.1450.tp.20170313.1545.036.html

O24

A

1673-629X(2017)05-0016-04

10.3969/j.issn.1673-629X.2017.05.004