基于模型函數(shù)與L-曲線的正則化參數(shù)選取方法

胡 彬，徐會(huì)林，王澤文，喻建華

(東華理工大學(xué)理學(xué)院，江西南昌330013;2.河南理工大學(xué)數(shù)學(xué)與信息科學(xué)學(xué)院，河南焦作454000)

0 引言

反問題研究已是計(jì)算數(shù)學(xué)及應(yīng)用數(shù)學(xué)領(lǐng)域研究的熱點(diǎn)問題之一.反問題一般是不適定的［1］，其不適定性在數(shù)值計(jì)算上表現(xiàn)為解不連續(xù)依賴測(cè)量數(shù)據(jù).即使在實(shí)際測(cè)量過程中測(cè)量誤差非常小，也會(huì)引起解的巨大波動(dòng).目前求解不適定問題最具有普遍性、完備性的是Tikhonov正則化方法，但該方法的有效性取決于選取到合適的正則化參數(shù).目前比較有代表性的正則化參數(shù)選取策略有:Morozov偏差原理、廣義交叉檢驗(yàn)、L-曲線(L-curve)等準(zhǔn)則.當(dāng)誤差水平已知或可估計(jì)時(shí)，Morozov偏差原理是最常采用的正則化參數(shù)選取策略之一，但是當(dāng)不適定算子方程右端項(xiàng)的誤差水平δ未知時(shí)，Morozov偏差原理就失效了.為了克服Morozov偏差原理需要已知誤差水平的局限性，P.C.Hansen等提出基于L-曲線的正則化參數(shù)選取方法［2-3］，可以在誤差水平未知的情形下找到近似最佳的正則化參數(shù).

在利用Morozov偏差原理確定正則化參數(shù)時(shí)需要借助牛頓迭代，而牛頓迭代的計(jì)算量比較大，為了克服這一缺點(diǎn)，文獻(xiàn)［4-7］中提出了1種新的確定正則化參數(shù)的方法—模型函數(shù)法.它是將Tikhonov泛函定義為關(guān)于α的函數(shù)，再用1個(gè)簡(jiǎn)單的具有顯示表達(dá)式的模型函數(shù)來近似F(α)，通過簡(jiǎn)單迭代確定正則化參數(shù).

本文在上述2種方法的基礎(chǔ)上研究了基于模型函數(shù)的修正L-曲線準(zhǔn)則，證明了所得序列點(diǎn)是局部收斂的，給出了相應(yīng)的迭代算法.通過算例驗(yàn)證了該準(zhǔn)則的有效性，同時(shí)也指出了目前還存在的不足和今后進(jìn)一步研究的方向.

1 模型函數(shù)法

線性反問題一般可歸結(jié)為解第1類不適定算子方程:

其中K是Hilbert空間X到Y(jié)上的有界線性算子.由于觀測(cè)數(shù)據(jù)存在誤差，所以一般把方程(1)改寫為

其中δ為誤差水平，‖yδ-y‖≤δ.Tikhonov正則化方法是將求不適定方程(2)的解轉(zhuǎn)為求Tikhonov泛函

的最優(yōu)解，其中α＞0是正則化參數(shù).對(duì)于固定的正則化參數(shù)α，記最優(yōu)函數(shù)為

引理1［4］設(shè)K:X→Y是有界線性算子，X，Y均為Hilbert空間.?α ＞0，(3)式的唯一解x(α)是無窮次可微的，且g=dnx(α)dαn∈X可由

遞推求解得到.

定理1［4-5］最優(yōu)函數(shù)F(α)在(0，+∞)內(nèi)是無限次可微的，且?α ＞0，F(xiàn)(α)滿足

模型函數(shù)的方法就是在αk附近構(gòu)造F(α)的局部近似函數(shù)Fk(α)，使其滿足微分方程(4)，即

若假設(shè)‖Kx(α)‖2≈ Tk‖x(α)‖2，代入(5)式并注意到 F'(α)= ‖x(α)‖2，解得

即得到雙曲模型函數(shù):

由于 Kx(α)≈ yδ，故可設(shè) ‖Kx(α)‖2≈‖yδ‖2-Tk，代入(5)式得

解微分方程(6)，得Fk(α)=Tk+Ckα，其中Ck、Tk是待定參數(shù)且由方程組

確定，即

于是，得到更為簡(jiǎn)潔的雙曲模型函數(shù):

文獻(xiàn)［8］研究了非精確數(shù)據(jù)下的線性模型函數(shù)選取正則化參數(shù).雖然線性模型函數(shù)具有計(jì)算簡(jiǎn)單、收斂性較好等優(yōu)點(diǎn)，但是卻不能將它與L-曲線準(zhǔn)則相結(jié)合而得到選取正則化參數(shù)的算法.文獻(xiàn)［9］將雙曲模型函數(shù)m1(α)應(yīng)用于L-曲線選取準(zhǔn)則獲得選取正則化參數(shù)的新算法.本文是前述研究的繼續(xù)和深入，基于雙曲模型函數(shù)m2(α)研究修正的L-曲線選取準(zhǔn)則，從而獲得選取正則化參數(shù)選取的新方法.

2 修正的L-曲線準(zhǔn)則

L-曲線準(zhǔn)則是殘差‖Kx(α)-yδ‖2與正則化解‖x(α)‖2在一組正則化參數(shù)下所構(gòu)成的圖像，也就是由(‖Kx(α)-yδ‖2，‖x(α)‖2)所構(gòu)成的平面曲線.最優(yōu)的正則化參數(shù)α出現(xiàn)在曲線的拐點(diǎn)處. 通常轉(zhuǎn)化為對(duì)應(yīng)的(2log‖Kx(α)-yδ‖，2log‖x(α)‖)曲線，因?yàn)榍€形狀如字母L(見圖1)，故稱為L(zhǎng)-曲線準(zhǔn)則.從圖像上很容易找到曲線的拐點(diǎn)，即最優(yōu)正則化參數(shù)α在對(duì)應(yīng)曲線的“角點(diǎn)”出現(xiàn).記u(α)=2log‖Kx(α)-yδ‖，v(α)=2log‖x(α)‖，則L曲線上各點(diǎn)的曲率公式［10-11］為

曲率最大的點(diǎn)對(duì)應(yīng)正則化參數(shù)即為所需正則化參數(shù).

圖1 L-曲線示意圖

由文獻(xiàn)［12］知，若L-曲線在點(diǎn)α=α*處取到最大曲率，且在該點(diǎn)處曲線的斜率為 -1/μ，則下列泛函:

在α=α*處取得極小值.因此，選取正則化參數(shù)的L-曲線準(zhǔn)則等價(jià)為求泛函(8)的極小值點(diǎn).修正的L-曲線準(zhǔn)則就是通過求(8)的極小值來獲得合適的正則化參數(shù).由于是α的非線性隱式函數(shù)，不便于計(jì)算.本文利用模型函數(shù)的方法將(8)式顯化，從而簡(jiǎn)化計(jì)算，提高計(jì)算效率.

3 基于模型函數(shù)的修正L-曲線算法

對(duì)于任意給定的α＞0，最優(yōu)函數(shù)F(α)為

且F'(α)= ‖x(α)‖2，則(8)式可改寫成F(α)的形式，即

模型函數(shù)m(α)是F(α)的局部近似，則(9)式的局部近似為

本文主要研究了雙曲模型函數(shù)m2(α)，代入(10)式得 ω(α)=(T+2C/α)(-C/α2)μ.

對(duì) ω(α)求導(dǎo)，得

令 ω'(α)=0，解得 α=-C(1+2μ)(μT).

對(duì)于正則化參數(shù) αk，求得對(duì)應(yīng)的正則化解x(αk)，再根據(jù)方程組(7)可確定參數(shù) Ck，Tk，則

由(11)式及Ck＜0，Tk＞0知αk+1＞0且它是唯一的.該算法比文獻(xiàn)［9］中的更為簡(jiǎn)單，因?yàn)槲墨I(xiàn)［9］給出的是關(guān)于α的2次方程，選取其中較大的作為正則化參數(shù)αk+1.

綜合上面分析，得到基于雙曲模型函數(shù)m2(α)與修正的L-曲線正則化參數(shù)選取策略的新算法:

算法1 給定 ε ＞ 0，yδ，K，μ ＞ 0;

Step 1 給定1個(gè)初始值α0＞α*，置k=0;

Step 2 解正則化方程αkx+K*Kx=Kyδ;

Step 3 求出 Ck，Tk，αk+1;

Step 5 停止迭代，輸出正則化參數(shù)αk+1.

定理2 如果

則在算法1中函數(shù)ωk(α)在點(diǎn)αk處是局部嚴(yán)格單調(diào)遞減的.

證算法1中產(chǎn)生的函數(shù)ωk(α)為

ωk(α)=(Tk+2Ck/α)(-Ck/α2)μ，則 ωk'(α)=(-Ck)μ(-2Ck+(αTk+2Ck)·(-2μ))/α2μ+2.取 α= αk，把(7)式中解得 Ck，Tk的表達(dá)式代入得

由已知條件得ω'k(αk)＜0，即函數(shù)ωk(α)在點(diǎn)αk處是局部嚴(yán)格單調(diào)遞減的.

定理3 給定初始值α0滿足(12)式，則算法1產(chǎn)生嚴(yán)格單調(diào)遞減序列{αk}且收斂.

證由定理2 的證明過程知，

令 φ(α)=-2Ck+(αTk+2Ck)(-2μ)，顯然φ(αk)＜ 0，φ'(αk)＜ 0.所以函數(shù)φ(α)是嚴(yán)格單調(diào)遞減的.由算法1知φ(αk+1)=0，所以αk＞ αk+1.又?k均有αk＞0，故收斂性成立.

4 數(shù)值算例

例1 求解第1類Fredholm積分方程［13］:

其中，核 K(s，t)=1 ［1+100(t-s)2］，

本文用等距節(jié)點(diǎn)復(fù)化梯形公式來離散第1類Fredholm 積分方程(13)，得線性方程組Ax－=y－，其中把積分核K(s，t)離散成矩陣Am×n，x(s)離散成n維列向量x－.對(duì)方程組右端加入隨機(jī)擾動(dòng)為

其中r為Matlab中的隨機(jī)函數(shù).取不同的誤差值δ，求出對(duì)應(yīng)不同誤差值δ的正則化參數(shù)，并把不同正則化參數(shù)求出的正則化解對(duì)比(見圖2～4)，其中實(shí)線表示無擾動(dòng)下的精確解，星形線表示不同正則化參數(shù)下的數(shù)值近似解.

情形1 當(dāng)δ=5.8378e-4，α=3.7320e-4，正則化解的相對(duì)誤差ρ=0.0678.正則化解與真解如圖2和圖3所示，其中圖2是該情形下未作正則化處理的計(jì)算所得解(即當(dāng)α=0時(shí)的最小二乘解)，圖3為算法1計(jì)算所得解.

情形2 當(dāng)δ=8.5401e-8，α=6.6470 e-4，正則化解的相對(duì)誤差ρ=0.0500，算法1計(jì)算結(jié)果如圖4所示.

圖2 α=0

圖3 α=3.7320e-4

圖4 α=6.6470e-4

例2 求解第1類Fredholm積分方程［14-15］:

x(s)=a1e-c1(s-t1)2+a2e-c2(s-t2)2，a1=2，a2=1，c1=6，c2=2，t1=0.8，t2=-0.5，K(s，t)=(cos(s)+cos(t))2(sin(u)u)，u= πsin(s).

情形1 當(dāng)δ=0.0022，α=0.0470時(shí)，正則化解的相對(duì)誤差ρ=0.1746.正則化解與真解如圖5和圖6所示，其中圖5是未作正則化處理的計(jì)算所得解，圖6為算法1計(jì)算所得解.

圖5 α=0

圖6 α=0.0470

情形2 當(dāng)δ=2.2268e-5，α=0.0466時(shí)，正則化解的相對(duì)誤差ρ=0.1744，算法1計(jì)算結(jié)果如圖7所示.

圖7 α=0.0466

5 結(jié)論

本文基于模型函數(shù)方法研究了正則化參數(shù)選取的修正L-曲線準(zhǔn)則，使得計(jì)算上更加簡(jiǎn)單.從算例的模擬結(jié)果可以看出，基于雙曲模型函數(shù)m2(α)的所得修正后的L-曲線準(zhǔn)則是有效的.但是，本文只證明了初始正則化參數(shù)選取滿足一定前提條件下，算法1產(chǎn)生的序列是局部收斂的.對(duì)于是否存在全局收斂性的算法［16］，還有μ值選取原則等問題還需進(jìn)一步研究.

［1］劉繼軍.不適定問題的正則化方法及應(yīng)用［M］.北京:科學(xué)出版社，2005.

［2］Hansen P C，O’Leary D P.The use of the L-curve in the regularization of discrete ill-posed problems［J］.SIAM J Sci Comput，1993，14(6):1487-1503.

［3］ Hansen P C.Analysis of discrete ill-posed problems bymeans of the L-curve ［J］.SIAM Review，1992，34(4):561-580.

［4］Xie Jianli，Zou Jun.An improvedmodel functionmethod for choosing regularization parameters in linear inverse problems［J］.Inverse Problems，2002，18(5):631-643.

［5］王澤文，徐定華.線性不適定問題中選取Tikhonov正則化參數(shù)的線性模型函數(shù)方法［J］.工程數(shù)學(xué)學(xué)報(bào)，2013，30(3):451-466.

［6］Wang Zewen，Liu Jijun.Newmodel functionmethods for determining regularization parameters in linear inverse problems［J］.Applied Numerical Mathematics，2009，59(10):2489-2506.

［7］Wang Zewen.Multi-parameter Tiknonov regularization andmodel function approach to the damped Morozov principle for choosing regularization parameters［J］.Journal of Computational andApplied Mathematics.2012，236(7):1815-1832.

［8］胡彬，夏赟，喻建華.算子非精確條件下確定正則化參數(shù)的一種方法［J］，江西師范大學(xué)學(xué)報(bào):自然科學(xué)版，2014，38(1):65-69.

［9］Heng Yi，Lu Shuai，MhamdiA，et al.Model functions in themodified L-curvemethod-case study:the heat flux reconstruction in pool boiling［J］.Inverse Problems，2010，26(5):1-13.

［10］張立濤，李兆霞，張宇峰，等.結(jié)構(gòu)識(shí)別計(jì)算中基于L-曲線的模型確認(rèn)方法研究［J］.運(yùn)動(dòng)與沖擊刊，2011，30(11):36-41.

［11］王宏志，趙爽，胡艷君.基于L-曲線正則化的MAP超分辨率圖像復(fù)原［J］.吉林大學(xué)學(xué)報(bào):理學(xué)版，2008，46(2):275-278.

［12］ Reginska T.A regularization parameter in discrete illposed problems［J］.SIAM J Sci Comput，1996，17(3):740-749.

［13］樊樹芳，馬青華，王彥飛.算子及觀測(cè)數(shù)據(jù)都非精確情況下一種新的正則化參數(shù)選擇方法［J］.北京師范大學(xué)學(xué)報(bào):自然科學(xué)版，2006，42(1):25-31.

［14］王彥飛.反問題的計(jì)算方法及應(yīng)用［M］.北京:高等教育出版社，2007.

［15］郭文彬.奇異值分解及其廣義逆理論中的應(yīng)用［M］.北京:中國(guó)科學(xué)院研究生院，2003.

［16］高煒，朱林立，梁立.基于圖正則化模型的本體映射算法［J］.西南大學(xué)學(xué)報(bào):自然科學(xué)版，2012，34(3):118-121.