鄭 寧，殷俊鋒

（同濟(jì)大學(xué)數(shù)學(xué)科學(xué)學(xué)院，上海 200092）

線性離散不適定問題來源于許多科學(xué)計算和工程應(yīng)用領(lǐng)域［1］。例如第一類Fredholm積分方程數(shù)值離散后的問題；例如圖像去模糊復(fù)原問題，已知觀測到的模糊帶噪聲圖像和模糊算子，需要盡可能復(fù)原得到清晰圖像。也就是說，在數(shù)學(xué)上考慮一個不相容的線性方程組

其中A∈Rm×n是一個大型稀疏的矩陣，表示離散的模糊算子，可以是方陣，也可以是長方形矩陣。右端項b∈Rm×1表示觀測得到的數(shù)據(jù)，即帶有噪聲的模糊圖像。式（1）和通常的線性方程組不同之處在于，矩陣A是長方形矩陣，且往往行數(shù)m大于列數(shù)n，因此式（1）是一個超定的線性系統(tǒng)。從數(shù)值角度上分析，矩陣A來源于帶積分項的模糊算子的離散，往往大型稀疏，且奇異值聚集在原點附近。這樣的矩陣往往帶有壞條件數(shù)，甚至接近奇異矩陣。觀測的數(shù)據(jù)b往往帶有觀測誤差。

式中：為未知的沒有噪聲的向量，且?為未知的清晰的真實圖像；向量e表示某種誤差或噪聲，其中每個分量都獨立服從相同的隨機(jī)分布。因此，線性系統(tǒng)（1）是不相容的，只能求解其對應(yīng)的最小二乘解。

令A(yù)?為矩陣A的摩爾-彭若斯廣義逆，則不適定問題（1）的最小二乘解為

然而，由于不適定問題中矩陣A病態(tài)，式（1）的解對右端項b中的誤差e非常敏感，由噪聲引起的分量A?e會將所需要的信息x掩蓋，使得解A?b沒有意義。為了解決此問題，需要正則化方法。常用的正則化方法分為2類：當(dāng)噪聲e的取值范圍已知時，可以利用偏差原理（discrepancy principle）方法；當(dāng)噪聲e未知時，可以利用Tikhonov正則化方法、廣義交叉校驗準(zhǔn)則和L曲線準(zhǔn)則等確定正則化參數(shù)。本文采用Tikhonov正則化方法，即最小化帶有懲罰項的目標(biāo)函數(shù)

式中：L為一個離散微分算子；μ為正則化參數(shù)，μ＞0。

在實際應(yīng)用中，數(shù)據(jù)往往是在一定的取值范圍內(nèi)選取的。例如，圖像的像素值是非負(fù)的；灰度圖像的像素值則在［0，255］之間。因此，考慮一類帶盒子約束的不適定問題［2］，即

其中，l∈Rn×1和u∈Rn×1分別是下界向量和上界向量。

針對大規(guī)模帶約束不適定問題（4），Bertero和Boccacci［3］提出一 類 投 影Landweber方法。Hanke等［4］考慮變量無上界的特殊情形，利用顯性變量替換將原問題轉(zhuǎn)換成一類非線性方程組，然后采用擬牛頓法進(jìn)行求解。Nagy和Strakos［5］將最速下降法與投影技巧結(jié)合，在目標(biāo)函數(shù)下降的同時，迭代解滿足約束條件。Rojas和Steihaug［6］從優(yōu)化算法出發(fā)，提出懲罰函數(shù)法，使得數(shù)值解序列在可行區(qū)域內(nèi)部收斂到最優(yōu)點。Calvetti等［2］將正則化方法和Krylow子空間迭代法相結(jié)合，提出一類內(nèi)外迭代法。文獻(xiàn)［6-7］在此基礎(chǔ)上引入積極集策略，外層迭代快速更新積極集，而內(nèi)層迭代則在子空間通過迭代求解一系列僅在非積極集上有約束的問題。數(shù)值實驗驗證此算法的有效性，并且數(shù)值效率高于前人提出的算法，所需的計算復(fù)雜度和CPU時間都更少。但此方法的缺點在于缺乏理論收斂性，尤其是無法保證目標(biāo)函數(shù)值隨著迭代步數(shù)逐漸單調(diào)下降。

針對內(nèi)層迭代中需要求解的無約束最小二乘問題，近年來諸多學(xué)者提出高性能隨機(jī)迭代方法。Strohmer和Vershynin［8］提出隨機(jī)Kaczmarz方法在求解相容系統(tǒng)上的收斂理論。Leventhal和Levis［9］提出求解不相容系統(tǒng)，即最小二乘問題的隨機(jī)坐標(biāo)下降方法，并給出當(dāng)系數(shù)矩陣列滿秩時的收斂理論。隨著新的有效概率準(zhǔn)則的構(gòu)造，更多高性能的隨機(jī)迭代方法相繼被提出［9-11］，在理論和數(shù)值實驗中均展示出比傳統(tǒng)確定性方法更高效的計算性能。

本文提出一類基于積極集策略的內(nèi)外迭代法。該算法在內(nèi)層迭代上利用隨機(jī)坐標(biāo)下降方法求解無約束最小二乘問題，并利用Armijo下降準(zhǔn)則對迭代步長進(jìn)行選擇，這樣就可以保證目標(biāo)函數(shù)值隨著迭代步數(shù)的增加而單調(diào)下降；同時在外層迭代上更新積極集和對應(yīng)的非積極集，并采用投影算子，將不在可行域中的數(shù)值解分量投影到可行域邊界上。數(shù)值實驗驗證此算法的高效性。在偏差準(zhǔn)則的收斂條件下，新的積極集內(nèi)外迭代法所利用的計算量、迭代步數(shù)和CPU時間都比前人提出的算法更少。

1 基于積極集策略的內(nèi)外迭代法

簡單介紹求解問題（4）的2種內(nèi)外迭代法。第1種方法是基于在內(nèi)迭代中將數(shù)值解投影到盒子約束可行域Ω={x∈Rn：l≤x≤u}中。第2種方法是基于積極集策略，通過在每次外迭代中更新積極集，從而在內(nèi)迭代中只需要求解一個規(guī)模較小的子問題。

容易驗證，問題（4）中的目標(biāo)函數(shù)等價于

如果不考慮對變量的盒子約束，極小化f(x)在數(shù)學(xué)上等價于求解一個正規(guī)方程組

令為此方程組的解，則通常。定義一個正交投影算子PΩ，則為

其中j=1，2，3…，n，這樣可以保證x∈Ω。通過反復(fù)迭代，使得數(shù)值解對應(yīng)的目標(biāo)函數(shù)值逐漸下降，從而得到所需要的解。

投影內(nèi)外迭代算法［3］求解帶盒子約束的問題（4）如下：

算法1［3］投影內(nèi)外迭代法。①選取一個迭代向量x0。②計算。③開始迭代k=0，1，2，…，直到收斂。④通過迭代法求得如下問題的數(shù)值解w k：

Morigi等［7］在投影內(nèi)外迭代算法2的基礎(chǔ)上引入積極集策略。令拉格朗日算子λk=?，其中x k是第k步的數(shù)值解，是x k對應(yīng)的殘量。根據(jù)問題（4）的Karush-Kuhn-Tucker條件可知，x*是問題（4）的最優(yōu)解當(dāng)且僅當(dāng)其滿足

對于x k，定義積極集為B(x k)={j：x kj=lj，λkj≥0}∪{j：x kj=uj，λkj≤0}，其對應(yīng)的補(bǔ)集自由變量為F(x k)={1，2，...，n}B(x k)?；诜e極集算法的思想是：如果已知x*的積極集，則約束最優(yōu)化問題（4）可以轉(zhuǎn)化成無約束優(yōu)化問題使得l F≤x F≤u F，其中是矩陣?中列指標(biāo)屬于集合F的子矩陣。由于x*未知，所以先選取初始迭代向量不斷更新積極集，直到收斂為止。

基于積極集的投影內(nèi)外迭代算法［5］求解帶盒子約束的問題（4）如下。

算法2［5］基于積極集的投影內(nèi)外迭代法。①選取一個迭代初始向量x0。②計算r0=??。③開始迭代k=0，1，2，…，直到收斂。④計算拉格朗日算子λk=?。⑤更新積極集B(x k)和自由變量集F(x k)。⑥通過迭代法求得如下問題的數(shù)值解。⑦令x kB，并投影到可行區(qū)域中，即⑧計算。⑨結(jié)束迭代。

在算法2中，外層迭代快速更新積極集，而內(nèi)層迭代則在子空間通過迭代求解一系列僅在非積極集上有約束的問題。和算法1類似，內(nèi)層迭代可以采用CGLS算法快速求解。數(shù)值實驗［5］表明，基于積極集的算法收斂速度更快，且所需的迭代步數(shù)和CPU時間都更少。

算法2的缺點在于理論上無法保證目標(biāo)函數(shù)值隨著迭代單調(diào)下降。因此，在第2節(jié)中構(gòu)造新的基于積極集的算法。

2 結(jié)合Armijo策略的新積極集算法

將Armijo充分下降準(zhǔn)則和積極集思想相結(jié)合，同時，在內(nèi)層無約束最小二乘問題的求解上采用隨機(jī)坐標(biāo)下降策略。

算法3基于積極集和Armijo條件的隨機(jī)投影內(nèi)外迭代法。①選取一個迭代初始向量x0。②計算。③開始迭代k=0，1，2，...，直到收斂。④計算拉格朗日算子。⑤更新積極集B(x k)和自由變量F(x k)。⑥通過隨機(jī)坐標(biāo)下降迭代法求得如下問題的數(shù)值解w k：。⑦令x kB，其中m是滿足如下Armijo充分下降準(zhǔn)則的最小非負(fù)整數(shù)：

并投影到可行區(qū)間中

3 數(shù)值實驗

用圖像去模糊復(fù)原的數(shù)值實驗來驗證基于積極集和Armijo條件的投影內(nèi)外迭代法求解帶盒子約束的不適定問題的可行性和有效性。

在數(shù)值實驗中，采用的收斂準(zhǔn)則為偏差原理。另外，選取二階微分算子為正則化光滑子

在圖1中分別畫出圖像復(fù)原算例“satellite”的原始圖像、模糊帶噪聲圖像、算法2復(fù)原圖像、算法3復(fù)原圖像。圖2、圖3分別為相對誤差和相對殘量隨著迭代步數(shù)變化的曲線。從圖2可以看出，2類積極集內(nèi)外迭代算法隨著迭代步數(shù)的增加，相對誤差都先下降后增加。最優(yōu)外迭代步數(shù)在5步附近。從圖3可以看出，相比于算法2，算法3的相對殘量下降更快，振蕩更小，這說明隨機(jī)坐標(biāo)下降方法結(jié)合Armijo下降準(zhǔn)則可以使計算更快更平穩(wěn)地收斂，所需要的計算量更少。

圖1 圖像復(fù)原算例“satellite”Fig.1 Restored image of‘statellite’

圖2 圖像復(fù)原算例“satellite”相對誤差隨迭代步數(shù)變化Fig.2 Relative error of restored image of‘satel?lite’versus the number of iterative steps

4 結(jié)論和展望

考慮帶盒子約束的不適定問題，構(gòu)造一類基于積極集策略的內(nèi)外迭代法。此算法的創(chuàng)新之處在于，在內(nèi)層迭代上利用隨機(jī)坐標(biāo)下降方法，同時采用Armijo下降準(zhǔn)則對迭代步長進(jìn)行選擇，這樣就可以保證目標(biāo)函數(shù)值隨著迭代步數(shù)的增加而單調(diào)下降。同時在外層迭代上更新積極集和對應(yīng)的非積極集，并采用投影算子，將不在可行域中的數(shù)值解分量投影到可行域邊界上。數(shù)值實驗驗證此算法的高效性。在偏差準(zhǔn)則的收斂條件下，新的積極集內(nèi)外迭代法所用的計算量、迭代步數(shù)和CPU時間都比前人提出的算法更少。在后續(xù)工作中，將考慮稀疏約束條件下的最小二乘問題，并構(gòu)造類似的積極集隨機(jī)迭代算法進(jìn)行求解。

作者貢獻(xiàn)聲明：

鄭 寧：負(fù)責(zé)算法構(gòu)造、數(shù)值實驗和分析、論文撰寫等。

殷俊鋒：負(fù)責(zé)算法設(shè)計和數(shù)值案例分析。