一種求解病態(tài)復(fù)線性方程組的混合算法

陳鳳坤，雷秀仁

(華南理工大學(xué) 數(shù)學(xué)學(xué)院，廣東 廣州 510640)

一種求解病態(tài)復(fù)線性方程組的混合算法

陳鳳坤，雷秀仁

(華南理工大學(xué) 數(shù)學(xué)學(xué)院，廣東 廣州 510640)

病態(tài)復(fù)線性方程的求解是現(xiàn)代應(yīng)用數(shù)學(xué)和很多工程應(yīng)用面臨的難題，用一般算法進(jìn)行求解時(shí)，得到的誤差較大，因此在一些高精度的工程應(yīng)用上，其結(jié)果往往不是特別理想。而隨著科技的發(fā)展，現(xiàn)代很多工程應(yīng)用對(duì)數(shù)據(jù)具有越來越高的精度要求(尤其是國(guó)家航天航空)，因此一個(gè)能求解病態(tài)復(fù)線性方程組的高精度算法是很有必要的。從病態(tài)復(fù)線性方程組求解的特點(diǎn)出發(fā)，對(duì)模擬退火法進(jìn)行改進(jìn)，并將其全局的收斂能力與雙共軛梯度法的高精度求解能力結(jié)合起來，提出了一種BCG-SA混合算法。數(shù)據(jù)實(shí)驗(yàn)表明，模擬退火法能對(duì)雙共軛梯度法求出的解進(jìn)行微調(diào)動(dòng)，幫助雙共軛梯度法在概率意義上跳出局部極小值點(diǎn)，從而提高求解精度。

病態(tài)復(fù)線性方程組；模擬退火算法；雙共軛梯度法；混合算法；希爾伯特矩陣

0 引 言

求解復(fù)線性方程組是很多工程實(shí)際應(yīng)用常遇到的問題，當(dāng)方程組的規(guī)模較大或者系數(shù)矩陣的條件數(shù)很大時(shí)，系數(shù)矩陣易呈現(xiàn)病態(tài)特性。當(dāng)前求解病態(tài)方程組效果較好的有預(yù)處理ILUCG[1]、主元加權(quán)迭代法[2]、新Jacobi迭代法[3]、粒子群算法[4]等，但用它們求解病態(tài)的復(fù)線性方程組的效果不是很理想。目前求解病態(tài)復(fù)線性方程組效果較好的有雙共軛梯度法，但其存在不收斂或者收斂速度慢的潛在問題，且在搜索的過程中可能陷入局部極小值。

現(xiàn)代很多高科技的工程應(yīng)用對(duì)數(shù)據(jù)具有越來越高的精度要求(尤其是國(guó)家航天航空)，因此一個(gè)能求解病態(tài)復(fù)線性方程組的高精度算法是很有必要的。為此，在分析研究雙共軛梯度算法和模擬退火法的基礎(chǔ)上，根據(jù)復(fù)線性方程組求解的特點(diǎn)對(duì)模擬退火法進(jìn)行了適當(dāng)改進(jìn)，將雙共軛梯度法高精度求解的特性和模擬退火法的全局搜索能力有機(jī)結(jié)合，提出了BCG-SA混合算法。實(shí)驗(yàn)結(jié)果表明，該算法雖然增加了迭代次數(shù)，但明顯提高了計(jì)算精度，對(duì)于一些有較高精度要求的工程應(yīng)用具有重要的現(xiàn)實(shí)意義。

1 雙共軛梯度法(BCG)原理

雙共軛梯度法是C.Lanczos[5]提出的一種用于求解一般復(fù)線性方程組的方法，具有計(jì)算量少、收斂速度快等優(yōu)點(diǎn)。在求解規(guī)模較大或者條件數(shù)較大的復(fù)線性方程組時(shí)往往效果較好，是目前求解病態(tài)方程組較為常用的方法之一。

記(·，·)為通常復(fù)內(nèi)積。若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：計(jì)算

并置k=k+1。

Step3：若‖xk+1-xk‖2<ε或k>kmax，則算法終止，轉(zhuǎn)到Step4,否則轉(zhuǎn)入Step2繼續(xù)計(jì)算。

Step4：輸出數(shù)值解xk+1和迭代次數(shù)k。

2 模擬退火法(SA)原理

模擬退火算法是由N.Metropolis[7]等于1953年提出的，它是基于Monte-Carlo迭代求解策略的一種隨機(jī)尋優(yōu)算法，其出發(fā)點(diǎn)是基于物理中固體物質(zhì)的退火過程與一般組合優(yōu)化問題之間的相似性。模擬退火算法從某一較高初溫出發(fā)，伴隨溫度參數(shù)的不斷下降，結(jié)合概率突跳特性在解空間中隨機(jī)尋找目標(biāo)函數(shù)的全局最優(yōu)解，是一種通用的優(yōu)化算法，具有概率性地跳出局部最優(yōu)解,并最終趨于全局最優(yōu)的特性。1982年Kinkpatrick等將內(nèi)能E模擬為目標(biāo)函數(shù)f，溫度T演化成控制參數(shù)t，首次將模擬退火算法用于解決組合優(yōu)化問題。

通過變分原理將式(3)轉(zhuǎn)化為求解如下無約束條件的多元非線性方程組。

(4)

則模擬退火法應(yīng)用于求解病態(tài)線性方程組的算法(算法2)[8-9]如下：

Step2：給x一個(gè)擾動(dòng)x'=x+Δx，其中Δx是整個(gè)實(shí)軸上服從N(0，σ2)的隨機(jī)變量，或者限制當(dāng)前解的一個(gè)鄰域Δx=scale×(rand-0.5)，其中rand是0到1之間的隨機(jī)數(shù)，計(jì)算相應(yīng)的目標(biāo)函數(shù)值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)

雙共軛梯度法是一種收斂速度快，且精度較高的算法，在實(shí)踐中表明雙共軛梯度法的效果非常好，但由于其在搜索過程中可能陷入局部極小值，所以在一些對(duì)于數(shù)據(jù)精度有較高要求的實(shí)際應(yīng)用中，需對(duì)雙共軛梯度法進(jìn)行一些改進(jìn)，以幫助它跳過局部極小值點(diǎn)。

模擬退火算法[10]具有良好的全局優(yōu)化性，雖然它在解決病態(tài)方程組方面效果不是特別好，但它能夠以隨機(jī)搜索技術(shù)從概率意義上找出目標(biāo)函數(shù)的全局最小值點(diǎn)。

當(dāng)雙共軛梯度算法達(dá)到較高精度后會(huì)出現(xiàn)收斂速度慢或者不再繼續(xù)收斂的情況，這表明算法陷入了局部極小值，而模擬退火在搜索過程引入了隨機(jī)因素，以一定的概率來接受一個(gè)比當(dāng)前解要差的解，因此有可能跳出這個(gè)局部最優(yōu)解。所以為了得到更高精度的解，可以在雙共軛梯度算法陷入局部極小值時(shí)，引入模擬退火法來幫助其跳出局部極小值點(diǎn)。由此受到啟發(fā)，將雙共軛梯度法和模擬退火法各自的優(yōu)點(diǎn)結(jié)合起來，提出了混合算法。具體的混合算法(算法3)[11-12]如下：

Step1：輸入矩陣A、右端向量b、誤差精度ε、ε0(ε0>ε),最大迭代次數(shù)kmax，SA算法最大調(diào)用次數(shù)Nmax(不宜設(shè)置過大，可取6～10)，置算法迭代次數(shù)k=0，SA算法調(diào)用次數(shù)N=0。

Step3：計(jì)算

并置k=k+1。

Step4：若‖xk+1-xk‖2<ε或k>kmax或N>Nmax，則算法終止，轉(zhuǎn)到Step5；若‖xk+1-xk‖<ε0，調(diào)用SA算法，置N=N+1,轉(zhuǎn)入Step2，否則轉(zhuǎn)入Step3繼續(xù)計(jì)算。

Step5：輸出數(shù)值解xk+1和迭代次數(shù)k。

4 實(shí)驗(yàn)及分析

為了說明算法BCG-SA的有效性和優(yōu)越性，進(jìn)行了數(shù)據(jù)模擬仿真實(shí)驗(yàn)。實(shí)驗(yàn)所用到的矩陣均從經(jīng)典的病態(tài)矩陣演變而來，并且已經(jīng)證明演變后的矩陣一方面保留了原矩陣的病態(tài)性，另一方面擁有復(fù)線性，可以代表病態(tài)復(fù)線性方程組。

將算法BCG-SA與在求解病態(tài)方程組上性能優(yōu)越的共軛梯度法(CG)、預(yù)處理共軛梯度法(IC-CG)、雙共軛梯度法(BCG)、預(yù)處理雙共軛梯度法(IC-CBCG)、模擬退火法(SA)進(jìn)行比較。數(shù)值實(shí)驗(yàn)結(jié)果表明，BCG-SA是有效的且在精度上更優(yōu)越。需要指出的是進(jìn)行數(shù)值實(shí)驗(yàn)時(shí)，也將BCG-SA混合算法與蟻群算法、粒子群算法、神經(jīng)網(wǎng)絡(luò)算法進(jìn)行了比較，但這些群智能算法的實(shí)驗(yàn)效果并不好，隨著矩陣階數(shù)的增加，誤差太大，甚至出現(xiàn)了解的失真。

算例1：考查Hilbert演變矩陣,對(duì)于具有病態(tài)方程組代表性的Hilbert矩陣B，其中

(5)

令演變矩陣A=B+B*j(j是復(fù)數(shù)單位，j2=-1),并取方程組右端項(xiàng)b為:

(6)

算例2:考查Pascal演變矩陣，對(duì)于具有病態(tài)方程組代表性的Pascal矩陣C[13-14]，其中

(7)

(8)

當(dāng)矩陣的階數(shù)n=20時(shí)，Hilbert演變矩陣的條件數(shù)為1.453 609×1018，Pascal演變矩陣的條件數(shù)為8.077 406×1017，已是高度病態(tài)的矩陣。實(shí)驗(yàn)時(shí)均取初始解x0為零向量，初始溫度Tmax=1 000，Lmax=n(n為矩陣A的階)，最大的試探次數(shù)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，實(shí)驗(yàn)結(jié)果如圖1和圖2所示(其中IC-CG在求解算例2時(shí)失真)。

圖1 Hilbert演變矩陣在不同階數(shù)下算法求解誤差的比較

圖2 Pascal演變矩陣在不同階數(shù)下算法求解誤差的比較

數(shù)據(jù)實(shí)驗(yàn)表明：

(1)與BCG-SA混合算相比，蟻群算法、粒子群算法等智能算法在求解高階病態(tài)方程組上效果并不理想。

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

(3)在算例2中，對(duì)CG、BCG進(jìn)行預(yù)處理的效果并不理想。

(4)相對(duì)于其他算法，混合算法能得到較高精度的解，從而表明引入模擬退火法對(duì)于雙共軛梯度法的改進(jìn)是有效的，它能夠幫助雙共軛梯度法在概率意義上跳過局部極小值點(diǎn)，從而得到較高精度的解。

5 結(jié)束語

為了提高病態(tài)復(fù)線性方程組的求解精度，在研究分析病態(tài)復(fù)線性方程的特性和求解方法的基礎(chǔ)上，提出了一種高精度混合算法。該算法是在對(duì)模擬退火法進(jìn)行適當(dāng)改進(jìn)的基礎(chǔ)上，將雙共軛梯度法高精度求解的能力和模擬退火法的全局搜索能力有機(jī)結(jié)合。實(shí)驗(yàn)結(jié)果表明，提出的混合算法，雖然增加了計(jì)算的迭代次數(shù)，但明顯提高了計(jì)算精度，對(duì)于一些對(duì)數(shù)據(jù)有較高精度要求的工程應(yīng)用具有重要的現(xiàn)實(shí)意義。下一步的工作是進(jìn)一步優(yōu)化模擬退火算法的參數(shù)。

[1] 于春肖，苑潤(rùn)浩，穆運(yùn)峰.新預(yù)處理ILUCG法求解稀疏病態(tài)線性方程組[J].數(shù)值計(jì)算與計(jì)算機(jī)應(yīng)用，2014，35(1)：21-27.

[2] 唐 麗,李鵬飛.主元加權(quán)迭代法求解病態(tài)線性方程組[J].科學(xué)技術(shù)與工程，2012，12(2)：381-383.

[3] 孔祥強(qiáng).病態(tài)線性方程組新的Jacobi迭代解法[J].大學(xué)數(shù)學(xué)，2013,29(5)：50-54.

[4] 賀天宇,李國(guó)望.病態(tài)線性方程組的粒子群算法[J].科技資訊，2012(8)：15.

[5] 張永杰,孫 泰.大型復(fù)線性方程組預(yù)處理雙共軛梯度法[J].計(jì)算機(jī)工程與應(yīng)用,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] 鄭洲順，黃光輝，楊曉輝.求解病態(tài)線性方程組的混合算法[J].貴州工業(yè)大學(xué)學(xué)報(bào):自然科學(xué)版，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].數(shù)學(xué)的實(shí)踐與認(rèn)識(shí)，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 網(wǎng)絡(luò)出版時(shí)間：2017-03-13

國(guó)家基金數(shù)學(xué)天元基金(B13-B5071130)；國(guó)家教育部高校博士點(diǎn)基金(B13-C7070170)

陳鳳坤(1990-)，男，碩士研究生，研究方向?yàn)槿褐悄芩惴?；雷秀仁，副教授，通信作者，研究方向?yàn)榭茖W(xué)與工程計(jì)算。

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