求解大規(guī)模線性離散不適定問題的RRArnoldi-Fractional Tikhonov正則化算法

張 慧

(南京航空航天大學 金城學院 基礎(chǔ)部，南京 211156)

?

求解大規(guī)模線性離散不適定問題的RRArnoldi-Fractional Tikhonov正則化算法

張 慧

(南京航空航天大學 金城學院 基礎(chǔ)部，南京 211156)

不適定問題廣泛出現(xiàn)在地球物理、自動控制等多種領(lǐng)域.正則化方法是求解此類問題近似解的有效算法.將Fractional Tikhonov正則化算法應(yīng)用于投影算法，提出了求解大規(guī)模線性離散不適定問題的Arnoldi-Fractional Tikhonov正則化算法.進一步提出限制值域的Arnoldi-Fractional Tikhonov正則化算法.并針對經(jīng)典算例，進行了數(shù)值試驗和比較.數(shù)值試驗結(jié)果表明了新算法是有效且具有優(yōu)勢的.

不適定問題；正則化方法；Fractional正則化方法；Fractional Tikhonov正則化算法

隨著科技日新月異的發(fā)展，不適定問題[1]廣泛出現(xiàn)在圖像恢復(fù)、數(shù)學物理等科學領(lǐng)域[2-3]，實際上，上述問題轉(zhuǎn)化為線性離散不適定問題，即求解

Ax=btrue+e，A∈Rn×n

(1)

其中b=btrue+e是實際觀測到的帶有噪音的右端項，btrue是理想狀態(tài)下的輸出結(jié)果，e是觀測過程中無法消除的噪音數(shù)據(jù).而真正需要求解的線性系統(tǒng)為：

Axtrue=btrue，A∈Rn×n

(2)

因此，求解的難度就在于如何由已知的觀測系統(tǒng)式(1)來求解未知的xtrue.

對于不適定問題求解的有意義研究開始于20世紀60年代，由前蘇聯(lián)A.N.Tikhonov[4]院士及其領(lǐng)導的工作小組在從事反問題理論研究中提出的，即當今最為著名和廣為使用的Tikhonov正則化方法，它也一直是學者們研究的熱點.2011年，Hochstenbach[5]和Reichel在此基礎(chǔ)上進一步提出了Fractional Tikhonov正則化方法.本文在此基礎(chǔ)上提出基于大規(guī)模系數(shù)矩陣A的投影Fractional Tikhonov正則化方法.

1 Fractional Tikhonov正則化方法

首先，考慮這樣一個極小化問題：

(3)

其中α>0.若0<α<1，則矩陣W可以借助AAT的廣義逆來表示.參數(shù)α的介入，會使得式(3)式的解xλ,α更接近線性系統(tǒng)式(1)的精確解.我們把上述方法稱作Fractional Tikhonov正則化方法或Weighed Tikhonov正則化方法.在經(jīng)典的Tikhonov正則化方法中α=1，即矩陣W=I.

為了求極小化問題(3)的解，我們可以轉(zhuǎn)化為求解如下線性系統(tǒng)：

(4)

任意給定一個λ>0,α>0，就能求出一個正則解xλ,α，解的精確度與λ,α的選取息息相關(guān).將系數(shù)矩陣A進行奇異值分解后即可轉(zhuǎn)化為求解

(5)

其中

(6)

那么正則化解表示如下：

(7)

其中x=Vy.Fractional Tikhonov正則化方法的濾波因子為

(8)

假設(shè)已知右端輸出項b中噪音部分e滿足‖e‖≤ε，那么可以利用廣義交叉準則確定正則化參數(shù)λ.令α>0，是事先定好的一個數(shù)值，定義

δ=ηε

(9)

這里取η>1，它是一個獨立于ε的參量，第六部分的數(shù)值試驗中均取η=1.1.要求正則化解xλ,α滿足

‖b-Axλ,α‖2=‖δ‖2

(10)

將式(7)代入式(10)的左端，有

(11)

將式(8)代入式(10)右端(1-σjφtikh,W(σ))中，有

其中μ=λ-1，那么求解式(10)相當于求解下式的零點，

(12)

我們采用牛頓法求解上式，得到μ值進而確定λ值，代入式(8)求得正則化解.

2 Arnoldi Fractional Tikhonov 正則化方法

Fractional Tikhonov正則化方法建立在對系數(shù)矩陣A進行奇異值分解的基礎(chǔ)上，然而對于大規(guī)模矩陣，它的奇異值分解需要非常大的計算量，因此我們選擇先將大規(guī)模問題投影到維數(shù)較低的Krylov子空間上，對小規(guī)模問題進行求解.Lewis和Reichel在2008年提出Arnoldi Tikhonov方法[6]，詳細地介紹了該類方法.

2.1 Arnoldi Tikhonov方法

1)選初始向量x0，計算殘量r0=b-Ax0，β=‖r0‖,ν1=r0/‖r0‖;

2)for j=1,2,…,m

ωj=Aνj,

for i=1,2,…,j

hij=(ωj,νi),

ωj=ωj-hijνi,

end for

hj+1,j=‖ωj‖,

if hj+1,j=0

m=j，exit；

else

νj+1=wj/hj+1,j;

end if

end for

4)計算近似正則解xλ=x0+Vmyλ.

2.2 Arnoldi Fractional Tikhonov方法

1)選初始向量x0，計算殘量r0=b-Ax0，β=‖r0‖,ν1=r0/‖r0‖;

2)for j=1,2,…,m

ωj=Aνj,

for i=1,2,…,j

hij=(ωj,νi),

ωj=ωj-hijνi,

end for

hj+1,j=‖ωj‖,

if hj+1,j=0

m=j，exit；

else

νj+1=wj/hj+1,j;

end if

end for

4)計算近似正則解xλ=x0+Vmyλ.

3 限制值域的Arnoldi Fractional Tikhonov正則化方法

Arnoldi類算法是建立在Krylov子空間span{b,Ab,A2b,…,Am-1b}的.這個子空間的最大缺點在于其包含了右端向量中的噪聲成分，這意味著由b,Ab,A2b,…,Am-1b線性組合得到的解很有可能包含了大量的噪聲成分.于是改進Krylov子空間如下

span{Ab,A2b,…,Amb}

這個子空間的優(yōu)點在于，其噪聲成分被乘上了矩陣A，而這樣做，可以對噪聲中的高頻部分起到一種阻尼和光滑效果.這就是所謂的限制值域的Arnoldi類算法，記作RRArnoldi方法.

(13)

RRArnoldi算法簡單描述如下：

3.1 RRArnoldi過程

1)選初始向量ν1=Ab/‖Ab‖;

2)for j=1,2,…,m

ωj=Aνj,

for i=1,2,…,j

hij=(ωj,νi),

ωj=ωj-hijνi,

end for

hj+1,j=‖ωj‖,

if hj+1,j=0

m=j，exit；

else

νj+1=wj/hj+1,j;

end if

end for

(14)

其中

Vm=[ν1,ν2,…,νm], Vm+1=[Vm,νm+1]

(15)

xj=Vjyj, yj∈Rj

那么

(16)

那么通過RRArnoldi過程將原問題轉(zhuǎn)化為極小化問題

(17)

假設(shè)yj是式(17)的極小范數(shù)最小二乘解，那么原問題的解為xj=Vjyj,y∈Rj.

4 數(shù)值試驗

文中提出的所有改進算法都進行了數(shù)值試驗，并與已有的算法進行了數(shù)值比較，以顯示其優(yōu)越性.數(shù)值試驗從matlab(2010b)工具箱“regularization tool”[7]中得到系數(shù)矩陣A，右端向量b和精確解x.算例中的α的選取與文獻[8]中的取法類似.

例 求解方程Ax=b，其中系數(shù)矩陣A分別取baart、foxgood、gravity、i_laplace(2)、phillips、shaw等6類矩陣，矩陣階數(shù)N取1 000、5 000，投影到的子空間維數(shù)m=10,20，取參數(shù)α=0.2，分別取右端擾動項ν=1%,5%,10%，其中ν=‖e‖/‖b‖.比較RRArnoldi Fractional Tikhonov正則化方法(簡記為RRAFT)與RRArnoldi Tikhonov(簡記為RRAT) 正則化方法的計算解與精確解相對誤差范數(shù)結(jié)果如表1、2所示.

表1 N=1 000，m=10，α=0.2 RRAFT和RRAT解的相對誤差比較

表2 N=5 000，m=20，α=0.2 RRAFT和RRAT解的相對誤差比較

根據(jù)表1、2的結(jié)果可知，RRAFT的計算結(jié)果精度比RRAT的計算結(jié)果高，參數(shù)α的選取會直接影響到求解精度，數(shù)值結(jié)果顯示，通過參數(shù)α的調(diào)節(jié)，計算結(jié)果精度提升較明顯.

6 結(jié)語

不適定問題廣泛出現(xiàn)在不同領(lǐng)域.它的求解難度在于系數(shù)矩陣的維數(shù)過大和右端項的噪音干擾.Arnoldi Fractional Tikhonov正則化方法和RRArnoldi Fractional Tikhonov正則化方法在一定程度上減輕了求解的難度.數(shù)值試驗表明，新方法在求解的精度上具有一定的優(yōu)勢.

[1] 黃光遠，劉小軍.數(shù)學物理反問題[M].濟南：山東科學技術(shù)出版社，1993.

[2] KIRSCH A.An introduction to the mathematical theory of inverse problems[M].NewYork:Springer- Verlag，1996.

[3] 王彥飛.反演問題的計算方法及其應(yīng)用[M].北京：高等教育出版社，2007.

[4] TIKHONOV A N.Solution of incorrectly formulated problems and the regularization method[J].Soviet.Math.Dokl.，1963，4:1035-1038.

[5] HOCHSTENBACH M E，REICHEL L.Fractional tikhonov regularization for linear discrete ill-posed problems[J].BIT，2011,51:197-215.

[6] LEWIS B,REICHEL L.Arnoldi-Tikhonov regularization methods[J].J.Comput.Appl.Math.，2009，226：92-102.

[7] 肖庭延，于滲根，王彥飛.反問題數(shù)值解法[M].北京：科學出版社，2003.

[8] MORIKUNI K，REICHEL L,HAYAMI K.FGMRES for linear discrete ill-posed problems[J].Appl.Numer.Math.，2014,75：175-187.

[責任編輯 王新奇]

An Arnoldi-Fractional Tikhonov Regularization Algorithm forSolving Large Scale Linear Discrete Ill Posed Problems

ZHANG Hui

(Department of Basic, Jincheng College, Nanjing University of Aeronautics and Astronautics, Nanjing 211156, China)

The ill posed problems are widely presented in many fields such as geophysics, automatic control and so on. The regularization method is an effective method to solve the approximate solution of this kind of problem. In this paper, the Fractional Tikhonov regularization algorithm is applied to the projection algorithm. Firstly, an Arnoldi-Fractional Tikhonov regularization algorithm is proposed for solving large-scale linear discrete ill posed problems. Secondly, an Arnoldi-Fractional Tikhonov regularization algorithm of restricted range is put forward. And numerical experiments and comparisons are carried out for the classical examples. The numerical experiment results show that the new algorithm is effective and has more advantages.

ill posed problem; regularization method; Fractional regularization method; Fractional Tikhonov regularization algorithm

1008-5564(2016)04-0017-05

2016-02-09

張 慧(1990—)，女，安徽馬鞍山人，南京航空航天大學金城學院基礎(chǔ)部助教，碩士，主要從事數(shù)值代數(shù)研究.

O411

A