求解線性逆問題的譜共軛梯度法*

王華軍，王　碩，曹義超

(桂林電子科技大學(xué)數(shù)學(xué)與計(jì)算科學(xué)學(xué)院，廣西桂林　541004)

?

求解線性逆問題的譜共軛梯度法*

王華軍，王碩**，曹義超

(桂林電子科技大學(xué)數(shù)學(xué)與計(jì)算科學(xué)學(xué)院，廣西桂林541004)

(School of Mathematics and Computational Science,Guilin University of Electronic Technology,Guilin,Guangxi,541004,China)

摘要：針對線性逆問題，把原問題的算子方程轉(zhuǎn)化為帶有Tikhonov正則項(xiàng)的無約束優(yōu)化問題，提出一個(gè)求解線性逆問題的新譜共軛梯度法，并證明算法的全局收斂性.數(shù)值結(jié)果表明，新算法是有效的.

關(guān)鍵詞：線性逆問題譜共軛梯度法Tikhonov正則項(xiàng)

0　引言

【研究意義】在理論研究和工程實(shí)踐中,許多問題的解決需要根據(jù)觀測數(shù)據(jù)來恢復(fù)或重構(gòu)物理介質(zhì)的結(jié)構(gòu)信息,通常稱這類問題為逆問題.逆問題分為線性逆問題和非線性逆問題,其中線性逆問題廣泛出現(xiàn)在地球物理學(xué)、圖像處理、地下水文學(xué)、全球海洋模擬、腫瘤檢測、非破壞性檢測等領(lǐng)域中[1].【前人研究進(jìn)展】對于線性逆問題的研究，已有不少數(shù)值解法，文獻(xiàn)[2]提出迭代Tikhonov正則化方法求解該問題，取得了較好的數(shù)值效果.文獻(xiàn)[3]提出了Levenberg-Marquardt方法求解該問題，且在數(shù)值試驗(yàn)中，用偏差原則作為停止準(zhǔn)則，但該方法所需CPU時(shí)間較長.文獻(xiàn)[4]提出并討論了一類NCPs逆問題.【本研究切入點(diǎn)】譜共軛梯度法是一類非常重要的優(yōu)化方法，他是將譜梯度法和共軛梯度法的思想相結(jié)合的一種方法，具有迭代簡單、易于編程以及存儲要求低等優(yōu)點(diǎn)，所以一般適用于求解大規(guī)模無約束優(yōu)化問題.近年來，文獻(xiàn)[5-7]用譜共軛梯度法求解線性和非線性優(yōu)化問題，得到很好的數(shù)值效果.因此，本文希望把譜共軛梯度法應(yīng)用到線性逆問題中.【擬解決的關(guān)鍵問題】針對線性逆問題，提出一個(gè)新的譜共軛梯度法，在理論上證明新提共軛梯度法的下降性和全局收斂性，并把該方法應(yīng)用到卷積方程反演和圖像處理的反演中.

1　新的譜共軛梯度法

在許多科學(xué)研究和工程應(yīng)用中，很多問題可以轉(zhuǎn)化為下面形式的算子方程

Hx=u.

(1)

我們研究的問題是已知H和u反解x，這里H表示線性算子，u為實(shí)驗(yàn)測得的數(shù)據(jù)，因而不可避免地帶有一定的噪聲[1]，即有

ue=u+η，

(2)

這里，η表示與u維數(shù)相同的噪聲。即問題(1)轉(zhuǎn)化為求解

Hx=ue.

(3)

為了便于求解，將問題(3)等價(jià)轉(zhuǎn)化為求解下面的最小二乘問題:

min f(x)=‖Hx-ue‖2.

(4)

由于H是一個(gè)條件數(shù)很大的病態(tài)矩陣，所以問題(3)是一個(gè)不適定問題，因此沒有唯一解.為了克服問題(3)的不適定性，本文在式(4)的基礎(chǔ)上通過增加一個(gè)Tikhonov正則項(xiàng)‖x‖2，將不適定問題轉(zhuǎn)化為適定問題，即求解下面的無約束優(yōu)化問題：

min f(x)=‖Hx-u‖2+α‖x‖2,

(5)

其中，α為正則化參數(shù)，而且該參數(shù)的選取對問題的求解發(fā)揮重要作用.

對于無約束優(yōu)化問題:

min {f(x)|x∈Rn}，

(6)

其中f(x)為式(5)的目標(biāo)函數(shù).共軛梯度法是求解問題(5)的有效方法[8-10]，其標(biāo)準(zhǔn)迭代形式為

xk+1=xk+αkdk，

(7)

其中αk是Wolfe搜索產(chǎn)生的步長，dk為搜索方向.dk定義為

(8)

其中g(shù)k=f(xk)，βk為方向調(diào)控參數(shù).不同的βk對應(yīng)不同的共軛梯度法，其中著名的βk計(jì)算公式[11-12]有

(9)

上述由參數(shù)βk構(gòu)成的方法分別稱為CD方法和LS方法.

Birgin等[13]結(jié)合譜梯度的思想，提出了下面形式的譜共軛梯度法：

(10)

并取

(11)

(12)

其中yk-1=gk-gk-1，sk-1=xk-xk-1.基于迭代格(10)，本文給出新的βk和θk(文獻(xiàn)[14])分別為

(13)

θk=

(14)

稱式(10),(13),(14)為新的譜共軛梯度法(SCG).以下結(jié)合Wolfe線搜索建立新的譜共軛梯度算法，并分析算法的收斂性質(zhì).

算法1

步驟1給定初值x1∈Rn,0<δ<1/2<σ<1,d1∶=-g1,ε≥0.如果‖g1‖≤ε,則停止迭代.

步驟2由Wolfe線搜索準(zhǔn)則計(jì)算步長αk，即αk滿足：

(15)

(16)

步驟3計(jì)算xk+1=xk+αkdk.若‖gk+1‖≤ε，則停止迭代.

步驟4由式(13)計(jì)算βk+1，由式(14)計(jì)算θk+1，由dk+1=-θk+1gk+1+βk+1dk計(jì)算dk+1.

步驟5k∶=k+1,轉(zhuǎn)步驟2.

以下均假設(shè)‖gk‖≠0，否則算法因已找到穩(wěn)定點(diǎn)而停止.

引理1設(shè){gk}和{dk}為算法1生成的序列，則

(17)

(18)

則對于k的情形就有

(19)

從而由數(shù)學(xué)歸納法知引理1成立.

引理1說明算法1產(chǎn)生的搜索方向是下降方向.

引理2如果步長因子αk滿足Wolfe線搜索條件,則有

(20)

證明由式(13)可得

(21)

(22)

(23)

所以，引理2得證.

2　算法的全局收斂性

為了證明算法1的全局收斂性，作如下假設(shè)[14-16]：

(H1)目標(biāo)函數(shù)f(x)在其水平集Ω={x∈Rn|f(x)≤f(x1)}上有界.

(H2)設(shè)f(x)的梯度g(x)在Ω上Lipschitz連續(xù),即存在L>0，使

‖g(y)-g(x)‖≤L‖y-x‖， ?x,y∈Ω.

(24)

引理3假設(shè)(H1)，(H2)成立，{gk}和{dk}為算法 1 生成的序列，則

(25)

證明由式(7)和式(16)及假設(shè)(H2)，可得

(26)

從而得到

(27)

由Wolfe搜索準(zhǔn)則和dk的下降性及假設(shè)(H1)知，{fk}為單調(diào)遞減的收斂數(shù)列，結(jié)合式(27)可得

(28)

對式(28)左右兩端分別從k=1,2,…求和，再利用{fk}的收斂性，可得式(25)成立,則引理3得證.

定理1假設(shè)(H1)，(H2)成立，{gk}為算法1生成的序列，則

(29)

證明若定理1不成立，則必存在常數(shù)r>0，使得對任意的k≥1，有

‖gk‖≥r.

(30)

式(10)移項(xiàng)，得dk+θkgk=βkdk-1，等式兩端取模平方，移項(xiàng)并利用式(20)，可得

(31)

(32)

由式(32)的遞推，并利用d1=-g1，可得

(33)

結(jié)合式(30)，可得

(34)

對式(34)兩端分別求和可得

(35)

顯然，式(35)與式(25)矛盾，所以定理1成立.

3　數(shù)值實(shí)驗(yàn)

算法測試的環(huán)境為Matlab 2013a，Windows 7操作系統(tǒng)，Intel(R)Core(TM) i3-2370M CPU @6.40 GHz.用Matlab的rand函數(shù)在真實(shí)解上添加一定水平的噪聲.參數(shù)δ=0.025，σ=0.9，α=0.014.為了評價(jià)不同方法的求解效果，引入真實(shí)解和數(shù)值解的均方根誤差(RMSE)的指標(biāo)[17-18]來衡量：

(36)

其中，x和xk分別是真實(shí)解和數(shù)值解，顯然RMSE越小，數(shù)值解對應(yīng)的算法越好.停止準(zhǔn)則是

‖gk‖≤10-5.

(37)

例1卷積方程的反演

考察在位勢理論中有重要應(yīng)用的卷積方程[1]：

(38)

其中u(t)與位勢方程邊界條件有關(guān),弱奇異性核H(t,t′)=1-ln(4sin2(t-t′)/2).現(xiàn)將式(38)記為算子形式：

(39)

用配置法離散式(39)得到離散的算子方程H(t)x(t)=u(t).

假設(shè)觀測數(shù)據(jù)是在真實(shí)數(shù)據(jù)上添加一定噪聲水平的高斯白噪聲[19]，則有

ue=u+e*rand(size(u))，

(40)

其中，e表示噪聲水平，rand(·)表示維數(shù)與u相同的高斯白噪聲，即有

H(t)x(t)=ue(t).

(41)

問題(41)等價(jià)于求解下列無約束優(yōu)化問題:

min f(x(t))=‖H(t)x(t)-ue(t)‖2.

(42)

為了克服問題(41)的不適定性，引入Tikhonov正則項(xiàng)，即該問題轉(zhuǎn)化為求解下述問題:

min f(x(t))=‖H(t)x(t)-ue(t)‖2+

γ‖x(t)‖2，

(43)

這里，α表示正則化參數(shù).

選擇核函數(shù)H(t)和真實(shí)解xtrue(x):

xtrue=exp (3sin t)，

(44)

(45)

在實(shí)驗(yàn)中，離散點(diǎn)設(shè)置為256個(gè)單位，即矩陣H(t)∈R256×256，真實(shí)解xture(w,t)∈R256，觀測數(shù)據(jù)ue(t)的噪聲水平設(shè)置為e=0.01，現(xiàn)用觀測數(shù)據(jù)ue(t)來反演x(t).為了更公平的測試算法的速度，每個(gè)實(shí)驗(yàn)重復(fù)做10次，結(jié)果的平均值見表1.從表1可以看出本文提出的譜共軛梯度法(SCG)的REMS值相對較小,且求解過程所需的迭代次數(shù)和CPU時(shí)間也是最小的，說明SCG方法所需更少的計(jì)算代價(jià).

從圖1～4可以看出SCG方法具有一定的優(yōu)勢.從圖5的效果圖也可以得出SCG方法是有效的.

表14種算法的數(shù)值結(jié)果

Table 1Numerical results of four algorithms

方法Method迭代次數(shù)Iterationtimes時(shí)間TimeREMSLandweber[20]81.0450.0044TSVD[21]121.20120.0141TV[21]62.0040.0251SCG50.14360.0019

圖1Landweber法求解例1的數(shù)值結(jié)果

Fig.1Numerical results of sample 1 obtained from Landweber method

圖2TSVD法求解例1的數(shù)值結(jié)果

Fig.2Numerical results of sample 1 obtained from TSVD method

圖3TV法求解例1的數(shù)值結(jié)果

Fig.3Numerical results of sample 1 obtained from TV method

圖4 SCG法求解例1的數(shù)值結(jié)果

Fig.4Numerical results of sample 1 obtained from SCG method

圖5例1效果圖

Fig.5Effects diagram for sample 1

例2圖像處理的反演問題

考察簡單的一維情形圖像識別問題,如方程：

(46)

其中x(t′)代表光源強(qiáng)度，u(t)代表觀測的圖像，通常是一個(gè)模糊核的連續(xù)圖像；h(t-t′)表示成像的過程，通常是一個(gè)點(diǎn)擴(kuò)展函數(shù)(PSF).問題給定光源函數(shù)x(t′)和點(diǎn)擴(kuò)展函數(shù)h(t-t′)來求觀測圖像u(t).與此相關(guān)的反問題：給定點(diǎn)擴(kuò)展函數(shù)h(t-t′)和觀測圖像u(t)來確定光源強(qiáng)度x(t′).

由有限差分法離散式(39)得到離散的算子方程

H(t)x(t)=u(t)，

(47)

假設(shè)觀測數(shù)據(jù)為在真實(shí)數(shù)據(jù)上添加imnoise函數(shù)產(chǎn)生的高斯噪聲，則有

ue=u+e*imnoise(size(u))，

(48)

其中，e表示噪聲水平，imnoise(·)表示維數(shù)與u相同的高斯噪聲，即有

H(t)x(t)=ue(t).

(49)

問題(49)等價(jià)于求解下列無約束優(yōu)化問題:

min f(x(t))=‖H(t)x(t)-ue(t)‖2.

(50)

為了克服問題(49)的不適定性，引入Tikhonov正則項(xiàng)，即轉(zhuǎn)該問題化為求解下述問題:

min f(x(t))=‖H(t)x(t)-ue(t)‖2+

γ‖x(t)‖2，

(51)

這里，γ表示正則化參數(shù).

選擇核函數(shù)H(t)和真實(shí)解xtrue(x):

xtrue=sin4(2πt),

(52)

(53)

在實(shí)驗(yàn)中，離散點(diǎn)設(shè)置為512個(gè)單位，即矩陣H(t)∈R512×512，真實(shí)解xture(w,t)∈R512，觀測數(shù)據(jù)ue(t)的噪聲水平設(shè)置為e=0.02，現(xiàn)用觀測數(shù)據(jù)ue(t)來反演x(t).

為了更公平的測試算法的速度，每個(gè)實(shí)驗(yàn)重復(fù)做10次(表2).從表2可以看出,SCG方法的RMSE值是0.012，在4個(gè)算法中相對較小，而且求解過程所需的迭代次數(shù)和CPU時(shí)間都相對較小，因此SCG方法所需求解成本相對較低.

從圖6～9可以看出SCG方法的數(shù)值解比較理想的逼近真實(shí)解.從圖10的效果圖也可以得出,SCG方法比其他3種算法更有效.

表24種算法的數(shù)值結(jié)果

Table 2Numerical results of four algorithms

方法Method迭代次數(shù)Iterationtimes時(shí)間TimeREMSLangweber[20]91.440.066TSVD[21]193.950.121TV[21]183.630.152SCG60.530.012

圖6Landweber法求解例2的數(shù)值結(jié)果

Fig.6Numerical results of sample 2 obtained from Landweber method

圖7TSVD法求解例2的數(shù)值結(jié)果

Fig.7Numerical results of sample 2 obtained from TSVD method

圖8TV法求解例2的數(shù)值結(jié)果

Fig.8Numerical results of sample 2 obtained from TV method

圖9SCG法求解例2的數(shù)值結(jié)果

Fig.9Numerical results of sample 2 obtained from SCG method

圖10例2效果圖

Fig.10Effects diagram for sample 2

4　結(jié)論

本文提出一個(gè)求解線性逆問題的新譜共軛梯度法，證明了該算法的全局收斂性,并利用該方法對卷積方程的反演問題和圖像處理的反演問題進(jìn)行數(shù)值實(shí)驗(yàn)，發(fā)現(xiàn)該方法比Landweber，TSVD和TV方法更有效.

參考文獻(xiàn)：

[1]肖庭延,于慎根,王彥飛.反問題的數(shù)值解法[M].北京:科學(xué)出版社,2003:166-174. XIAO T Y,YU S G,WANG Y F.Numerical Solution of Inverse Problem[M].Beijing:Science Press,2003:166-174.

[2]LIU W,WU C S.A predictor-corrector iterated Tik- honov regularization for linear ill-posed inverse problems[J].Applied Mathematics and Computation,2013,221:802-818.

[3]JIN Q N.On a regularized Levenberg-Marquardt method for solving nonlinear inverse problems[J].Numerische Mathematik,2010,115(2):229-259.

[4]ZHANG J Z,JIAN J B,TANG C M.Inverse problems and solution methods for a class of nonlinear complementarity problems[J].Computational Optimization and Applications,2011,49(2):271-297.

[5]江羨珍,馬國棟,簡金寶.Wolfe線搜索下一個(gè)新的全局收斂共軛梯度法[J].工程數(shù)學(xué)學(xué)報(bào),2011,28(6):779-786. JIANG X Z,MA G D,JIAN J B.A new globally convergent conjugate gradient method with Wolfe line search[J].Chinese Journal of Engineering Mathematics,2011,28(6):779-786.

[6]江羨珍,韓麟,簡金寶.Wolfe線搜索下一個(gè)全局收斂的混合共軛梯度法[J].計(jì)算數(shù)學(xué),2012,34(1):103-112. JIANG X Z,HAN L,JIAN J B.A globally convergent mixed conjugate gradient method with Wolfe line search[J].Mathematica Numerica Sinica,2012,34(1):103-112.

[7]JIAN J B,HAN L,JIANG X Z.A hybrid conjugate gradient method with descent property for unconstrained optimization[J].Applied Mathematical Modelling,2015,39(3/4):1281-1290.

[8]簡金寶,尹江華,江羨珍.一個(gè)充分下降的有效共軛梯度法[J].計(jì)算數(shù)學(xué),2015,37(4):415-424. JIAN J B,YIN J H,JIANG X Z.An efficient conjugate gradient method with sufficient descent property[J].Mathematica Numerica Sinica,2015,37(4):415-424.

[9]韓麟,簡金寶,江羨珍.Armijo線搜索下一個(gè)雜交共軛梯度法及其強(qiáng)收斂性[J].玉林師范學(xué)院學(xué)報(bào):自然科學(xué),2013,34(2):17-21. HAN L,JIAN J B,JIANG X Z.A mixed conjugate gradient method with Armijo linear search and its strong convergence[J].Journal of Yulin Normal University:Natural Science,2013,34(2):17-21.

[10]JIANG X Z,JIAN J B.Two modified nonlinear conjugate gradient methods with disturbance factors for unconstrained optimization[J].Nonlinear Dynamics,2014,77(1/2):387-397.

[11]POLAK E,RIBIRE G.Note sur la convergence de méthodes de directions conjuguées[J].Revue Fran?aise D'informatique et de Recherche Opérationnelle,1969,16:35-43.

[12]POLYAK B T.The conjugate gradient method in extreme problems[J].Ussr Computational Mathematics and Mathematical Physics,1969,9(4):94-112.

[13]BIRGIN E G,MARTNEZ J M.A spectral conjugate gradient method for unconstrained optimization[J].Applied Mathematics and Optimization,2001,43(2):117-128.

(Continue on page 427)(Continue from page 421)

[14]JIANG X Z,JIAN J B.A sufficient descent Dai-Yuan type nonlinear conjugate gradient method for unconstrained optimization problems[J].Nonlinear Dynamics,2013,72(1/2):101-112.

[15]江羨珍,簡金寶,馬國棟.具有充分下降性的兩個(gè)共軛梯度法[J].數(shù)學(xué)學(xué)報(bào),2014,57(2):365-372. JIANG X Z,JIAN J B,MA G D.Two conjugate gradient methods with sufficient descent property[J].Acta Mathematica Sinica,2014,57(2):365-372.

[16]馬國棟,簡金寶,江羨珍.一個(gè)具有下降性的改進(jìn) Fletcher-Reeves共軛梯度法[J].應(yīng)用數(shù)學(xué)學(xué)報(bào),2015,38(1):89-97. MA G D,JIAN J B,JIANG X Z.An improved Fletcher-Reeves conjugate gradient method with descent property[J].Acta Mathematicae Applicatae Sinica,2015,38(1):89-97.

[17]ZHANG B X,ZHU Z B,LI S A.A modified spectral conjugate gradient projection algorithm for total variation image restoration[J].Applied Mathematics Letters,2014,27:26-35.

[18]ZHANG B X,ZHU Z B,WANG S.A simple primal-dual method for total variation image restoration[J].Journal of Visual Communication and Image Representation,2016,38:814-823.

[19]ZHANG B X,ZHU Z B.A modified quasi-Newton diagonal update algorithm for total variation denoising problems and nonlinear monotone equations with applications in compressive sensing[J].Numerical Linear Algebra with Applications,2015,22(3):500-522.

[20]曹莉,王偉.求解非線性不適定問題的連續(xù)Landweber型正則化方法[J].黑龍江大學(xué)自然科學(xué)學(xué)報(bào),2014,31(2):165-170. CAO L,WANG W.Continuous Landweber type regularization method for solving nonlinear ill-posed problem[J].Journal of Natural Science of Heilongjiang University,2014,31(2):165-170.

[21]VOGEL C R.Computational Methods for Inverse Pro- blems[M].Philadelphia:Society for Industrial and Applied Mathematics,2011:1-11.

(責(zé)任編輯：尹闖)

A Spectral Conjugate Gradient Method for Solving Linear Inverse Problems

WANG Huajun,WANG Shuo,CAO Yichao

Key words：linear inverse problems,spectral conjugate gradient method,Tikhonov regularization

Abstract：In this study,a new spectral conjugate gradient method is presented to solve linear inverse problems,which are transferred into the linear unconstrained optimization with Tikhonov regularization.The global convergence of the proposed scheme is analyzed.The new algorithm is compared with Landweber,TSVD and TV methods,and numerical results illustrate the efficiency of this method.

收稿日期：2016-03-15

作者簡介：王華軍(1990-),男,碩士研究生,主要從事優(yōu)化算法在反問題中的應(yīng)用研究。

中圖分類號：O224

文獻(xiàn)標(biāo)識碼：A

文章編號：1005-9164(2016)05-0416-06

修回日期：2016-04-15

*國家自然科學(xué)基金項(xiàng)目(11361018),廣西自然科學(xué)基金項(xiàng)目(2014GXNSFFA118001),廣西高校中青年教師基礎(chǔ)能力提升項(xiàng)目(ky2016YB167)和桂林電子科技大學(xué)研究生教育創(chuàng)新計(jì)劃項(xiàng)目(2016YJCX46)資助。

**通信作者：王碩(1987-),女,講師,主要從事最優(yōu)化方法及其應(yīng)用研究,E-mail:optimization_zhu@163.com。

廣西科學(xué)Guangxi Sciences 2016,23(5):416～421,427

網(wǎng)絡(luò)優(yōu)先數(shù)字出版時(shí)間：2016-11-21【DOI】10.13656/j.cnki.gxkx.20161121.013

網(wǎng)絡(luò)優(yōu)先數(shù)字出版地址：http://www.cnki.net/kcms/detail/45.1206.G3.20161121.1546.026.html