求解帶擾動的線性方程組的貪婪隨機(jī)Kaczmarz方法

巫文婷

（北京理工大學(xué)數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，北京 100081）

對于系數(shù)矩陣為A∈Cm×n且右端向量為b∈Cm的大規(guī)模相容線性代數(shù)方程組

的求解，Kaczmarz方法［1］是經(jīng)典的行處理迭代方法，在信號與圖像處理領(lǐng)域有著廣泛的應(yīng)用。其每步迭代只需按照給定的循環(huán)順序選取系數(shù)矩陣的某一行，并將當(dāng)前迭代向量正交投影至由該行所形成的超平面上。Strohmer等［2］提出按照與系數(shù)矩陣每一行的歐氏范數(shù)平方成比例的概率準(zhǔn)則隨機(jī)選取系數(shù)矩陣的行，得到了收斂更為快速的隨機(jī)Kaczmarz方法。若用(·)*表示相應(yīng)矩陣或向量的共軛轉(zhuǎn)置，則當(dāng)初始迭代向量在A*的列空間中時(shí)，隨機(jī)Kaczmarz方法期望線性收斂［2-5］到線性代數(shù)方程組（1）的最小歐氏范數(shù)解x?=A?b，其中A?表示系數(shù)矩陣A的Moore-Penrose偽逆。當(dāng)線性代數(shù)方程組（1）的右端向量發(fā)生擾動時(shí)，Needell［6］給出了隨機(jī)Kaczmarz方法的期望解誤差的上界，并說明了隨著迭代步數(shù)的增長，隨機(jī)Kaczmarz方法的期望解誤差會以線性速率下降至一個(gè)誤差閾值。之后，Zouzias等［7］對Needell所提出的期望解誤差的上界進(jìn)行了改進(jìn)。

影響隨機(jī)Kaczmarz方法收斂速率的關(guān)鍵在于其中所蘊(yùn)含的用于選取每步迭代所需調(diào)用的系數(shù)矩陣行的概率準(zhǔn)則。為了提高隨機(jī)Kaczmarz方法的收斂速率，Bai等［8-9］提出了一個(gè)可以獲取每步迭代中殘向量的模較大分量的概率準(zhǔn)則，并基于該概率準(zhǔn)則構(gòu)造出了貪婪隨機(jī)Kaczmarz方法。當(dāng)初始迭代向量在A*的列空間時(shí)，貪婪隨機(jī)Kaczmarz方法所產(chǎn)生的迭代序列收斂到線性代數(shù)方程組（1）的最小范數(shù)解x?=A?b且具有期望線性收斂速率

當(dāng)線性代數(shù)方程組（1）的右端向量b被加上一個(gè)不為零的擾動向量r∈Cm時(shí)，實(shí)際求解的線性代數(shù)方程組問題將變?yōu)?/p>

其中，y=b+r。若對任意矩陣G∈Cm×n和任意向量u∈Cm，用G(i)和u(i)分別表示矩陣G的第i行和向量u的第i個(gè)分量，由于貪婪隨機(jī)Kaczmarz方法的第k步迭代將當(dāng)前迭代向量x k投影至超平面=上，而線性代數(shù)方程組（1）的最小范數(shù)解x?=A?b滿足Ax?=y-r，故其并不能收斂到x?。在這種情況下，本文對貪婪隨機(jī)Kaczmarz方法的期望解誤差進(jìn)行分析，得到了其上界，證明了隨著迭代步數(shù)的增長，由貪婪隨機(jī)Kaczmarz方法所產(chǎn)生的迭代解與原線性代數(shù)方程組（1）的最小范數(shù)解x?之間的誤差以線性速率下降至一個(gè)給定閾值。數(shù)值實(shí)驗(yàn)表明本文所給出的閾值能夠很好地估計(jì)貪婪隨機(jī)Kaczmarz方法的迭代解誤差所能達(dá)到的最小值。

1 貪婪隨機(jī)Kaczmarz方法

不失一般性，在本文的討論中，總是假設(shè)系數(shù)矩陣A不存在零行，即A(i)≠0，i=1，2，…，m。當(dāng)線性代數(shù)方程組（1）的右端向量發(fā)生擾動時(shí)，求解線性代數(shù)方程組（2）的貪婪隨機(jī)Kaczmarz方法［8］如下：

輸入：A，y，?與x0

輸出：x?

1：fork=0，1，…，?-1 do

2：計(jì)算

3： 確定正整數(shù)指標(biāo)集

4： 計(jì)算向量r?k∈Cm的第i個(gè)分量

7：endfor

令R(A)和R(A)⊥分別為系數(shù)矩陣A的像空間和其像空間的正交補(bǔ)子空間，則有r=r R(A)+r R(A)⊥，其中r R(A)和r R(A)⊥分別表示向量r在R(A)和R(A)⊥上的正交投影。若線性代數(shù)方程組（1）右端向量的擾動r在系數(shù)矩陣A的像空間R(A)中，則線性代數(shù)方程組（2）等價(jià)于線性代數(shù)方程組

易知線性代數(shù)方程組（3）是相容的且其最小范數(shù)解為=A?(b+r R(A))。此時(shí)，若初始迭代向量在A*的列空間中，則求解線性代數(shù)方程組（2）的貪婪隨機(jī)Kaczmarz方法所產(chǎn)生的迭代序列期望線性收斂到x??。

2 誤差分析

與求解原線性代數(shù)方程組（1）的貪婪隨機(jī)Kaczmarz方法的分析［8］類似，根據(jù)求解線性代數(shù)方程組（2）的貪婪隨機(jī)Kaczmarz方法中εk的定義可知，對于k=1，2，…，由于

有

而對于k=0，有

因此，可得引理1。

引理1求解系數(shù)矩陣為A∈Cm×n且右端向量為y∈Cm的線性代數(shù)方程組（2）的貪婪隨機(jī)Kaczmarz方法的概率準(zhǔn)則中的量εk，k=0，1，2，…，滿足

和

當(dāng)線性代數(shù)方程組（1）右端向量的擾動r在R(A)中時(shí)，若初始迭代向量在A*的列空間中，則求解線性代數(shù)方程組（2）的貪婪隨機(jī)Kaczmarz方法所產(chǎn)生的迭代向量期望線性收斂到線性代數(shù)方程組（3）的最小范數(shù)解x??。對于一般的擾動情形，求解線性代數(shù)方程組（2）的貪婪隨機(jī)Kaczmarz方法所產(chǎn)生的迭代向量與x??之間的期望誤差的上界可以由引理2給出。

引理2如果初始迭代向量x0∈Cn在A*的列空間中，則貪婪隨機(jī)Kaczmarz方法求解擾動后的線性代數(shù)方程組（2）所產(chǎn)生的迭代序列與線性代數(shù)方程組（3）的最小范數(shù)解之間的期望誤差滿足

其中

證明：由求解線性代數(shù)方程組（2）的貪婪隨機(jī)Kaczmarz方法的定義和Ax??=b+r R(A)可知

其中I表示具有合適階數(shù)的單位矩陣，則有

用Ek表示固定前k步迭代的條件期望，即Ek[·]=E[·|i0，i1，…，ik-1]，其中il(l=0，1，…，k-1)為第l步迭代所選取的行指標(biāo)，則對等式（5）兩邊取條件期望可得

由于

可知

則根據(jù)正整數(shù)指標(biāo)集Uk的定義可得

對于在矩陣A*的像空間R(A*)中的任意向量u，有不等式

成立。利用該不等式，通過引理1、向量b+r R(A)-Ax k=A(-x k)∈R(A)和向量r R(A)⊥的正交性、x k-∈R(A*)以 及可 知，對 于k=1，2，…，有

而對于k=0，有

對不等式兩邊取期望可得對于k=1，2，…，有

而對于k=0，有

因此，由于0＜α＜1，β≥0，對于k=1，2，…，有

由Jensen不等式可知

基于引理2，關(guān)于求解擾動后的線性代數(shù)方程組（2）的貪婪隨機(jī)Kaczmarz方法所產(chǎn)生的迭代近似解與原線性代數(shù)方程組（1）的最小范數(shù)解x?之間的期望誤差的上界，可以給出如下定理。

定理1如果初始迭代向量x0∈Cn在A*的列空間中且為線性代數(shù)方程組（3）的最小范數(shù)解，則貪婪隨機(jī)Kaczmarz方法求解擾動后的線性代數(shù)方程組（2）所產(chǎn)生的迭代序列與原線性代數(shù)方程組（1）的最小范數(shù)解x?=A?b之間的期望誤差滿足

其中，α、α0和β如（4）式中定義。

證明：因?yàn)?x?∈R(A*)，則通過不等式（6）可知

又因?yàn)橛?/p>

成立，則由引理2可得

定理1說明了當(dāng)?shù)綌?shù)k趨于無窮時(shí)，貪婪隨機(jī)Kaczmarz方法求解擾動后的線性代數(shù)方程組（2）所產(chǎn)生的迭代序列的期望相對解誤差以的線性速率下降至某個(gè)閾值，并給出了該閾值的估計(jì)

若r R(A)⊥為零，則擾動后的線性代數(shù)方程組（2）即為相容的線性代數(shù)方程組（3）。因此，求解線性代數(shù)方程組（2）的貪婪隨機(jī)Kaczmarz方法的迭代序列收斂到線性代數(shù)方程組（3）的最小范數(shù)解，且其與原線性代數(shù)方程組（1）的最小范數(shù)解x?的差趨于-x?。此時(shí)，如果初始迭代向量在A*的列空間中，則求解線性代數(shù)方程組（2）的貪婪隨機(jī)Kaczmarz方法所產(chǎn)生的迭代序列滿足

若r R(A)為零，則線性代數(shù)方程組（3）即為線性代數(shù)方程組（1），且線性代數(shù)方程組（3）的最小范數(shù)解即為x?=A?b。此時(shí)，如果初始迭代向量在A*的列空間中，則求解線性代數(shù)方程組（2）的貪婪隨機(jī)Kaczmarz方法所產(chǎn)生的迭代序列滿足

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

通過數(shù)值實(shí)驗(yàn)測試貪婪隨機(jī)Kaczmarz方法求解帶右端項(xiàng)擾動的線性代數(shù)方程組（2）時(shí)的數(shù)值表現(xiàn)，并將所估計(jì)的閾值與貪婪隨機(jī)Kaczmarz方法所產(chǎn)生的真實(shí)迭代解誤差進(jìn)行比較。

所測試的系數(shù)矩陣A分為兩類：一類為利用MATLAB函數(shù)randn隨機(jī)產(chǎn)生的元素服從標(biāo)準(zhǔn)正態(tài)分布的矩陣，矩陣維數(shù)分別為200×100 000和100000×200；另一類為取自稀疏矩陣庫［10］的稀疏矩陣bibd_16_8和ash958，矩陣維數(shù)分別為120×12870和958×292。

線性代數(shù)方程組（1）的右端向量b取為Ax*，其中x*∈Rn是利用MATLAB函數(shù)randn隨機(jī)產(chǎn)生的。線性代數(shù)方程組（1）的右端項(xiàng)的擾動向量r分為三類：一類為隨機(jī)擾動，一類為在系數(shù)矩陣的像空間中的擾動，另一類為在系數(shù)矩陣的像空間的正交補(bǔ)空間中的擾動。這三類擾動均由MATLAB函數(shù)randn生成，并且其歐氏范數(shù)為原始右端向量b的歐氏范數(shù)的0.0005倍。由于200×100000的隨機(jī)矩陣和稀疏矩陣bibd_16_8為行滿秩矩陣，故其所對應(yīng)的右端項(xiàng)擾動r一定在其像空間中。所有計(jì)算均從初始向量x0=0開始。

圖1描繪了當(dāng)線性代數(shù)方程組（1）的系數(shù)矩陣為隨機(jī)產(chǎn)生的矩陣且右端項(xiàng)擾動r為隨機(jī)擾動時(shí)，貪婪隨機(jī)Kaczmarz方法重復(fù)運(yùn)行50次所產(chǎn)生的相對解誤差ERS的中值和閾值τ分別取以10為底的對數(shù)之后相對于迭代步數(shù)的曲線，其中

從圖1可以看出，貪婪隨機(jī)Kaczmarz方法所產(chǎn)生的相對解誤差下降至10-3左右之后不再減小，且閾值τ能夠很好地給出該最小值的估計(jì)。特別地，當(dāng)系數(shù)矩陣為200×100 000的隨機(jī)矩陣時(shí)，τ與真實(shí)相對解誤差所能達(dá)到的最小值非常接近。

圖1 當(dāng)m=200，n=100 000和m=100 000，n=200時(shí)，lg E RS和lgτ相對于迭代步數(shù)的曲線Fig.1 Pictures of lg E RS and lgτversus the iteration step when m=200,n=100000 and m=100000,n=200

圖2描繪了對于100 000×200的隨機(jī)產(chǎn)生的系數(shù)矩陣，當(dāng)線性代數(shù)方程組（1）的右端項(xiàng)擾動r分別在系數(shù)矩陣的像空間和系數(shù)矩陣的像空間的正交補(bǔ)空間中時(shí)，貪婪隨機(jī)Kaczmarz方法重復(fù)運(yùn)行50次所產(chǎn)生的相對解誤差的中值和閾值τ分別取以10為底的對數(shù)之后相對于迭代步數(shù)的曲線。從圖2可以看出，貪婪隨機(jī)Kaczmarz方法所產(chǎn)生的相對解誤差下降至10-3左右之后不再減小，且閾值τ能夠很好地給出該最小值的估計(jì)，特別當(dāng)擾動向量r在系數(shù)矩陣的像空間中時(shí)。

圖2 當(dāng)m=100 000，n=200且r∈R(A)和r∈R(A)⊥時(shí)，lg E RS和lgτ相對于迭代步數(shù)的曲線Fig.2 Pictures of lg E RS and lgτversus the iteration step when m=100 000，n=200,and r∈R(A)and r∈R(A)⊥

當(dāng)線性代數(shù)方程組（1）的系數(shù)矩陣為稀疏矩陣bibd_16_8和ash958時(shí)，圖3描繪了右端項(xiàng)擾動r為隨機(jī)擾動時(shí)，貪婪隨機(jī)Kaczmarz方法重復(fù)運(yùn)行50次所產(chǎn)生的相對解誤差的中值和閾值τ分別取以10為底的對數(shù)之后相對于迭代步數(shù)的曲線；圖4描繪了右端項(xiàng)擾動r分別在系數(shù)矩陣的像空間和系數(shù)矩陣的像空間的正交補(bǔ)空間中時(shí)的相應(yīng)曲線。類似地，從圖3和圖4可以看出，貪婪隨機(jī)Kaczmarz方法所產(chǎn)生的相對解誤差下降至10-3左右之后不再減小，且閾值τ能夠很好地給出該最小值的估計(jì)，特別當(dāng)系數(shù)矩陣為bibd_16_8和當(dāng)系數(shù)矩陣為ash958且擾動向量r在系數(shù)矩陣的像空間中時(shí)。

圖3 當(dāng)矩陣為bibd_16_8和ash958時(shí)，lg E RS和lgτ相對于迭代步數(shù)的曲線Fig.3 Pictures of lg E RS and lgτversus the iteration step when the matrices are bibd_16_8 and ash958

圖4 當(dāng)矩陣為ash958且r∈R(A)和r∈R(A)⊥時(shí)，lg E RS與lgτ相對于迭代步數(shù)的曲線Fig.4 Pictures of lg E RS and lgτversus the iteration step when the matrix is ash958,and r∈R(A)and r∈R(A)⊥

4 結(jié)論

當(dāng)線性代數(shù)方程組的右端向量發(fā)生擾動時(shí)，證明了貪婪隨機(jī)Kaczmarz方法的期望解誤差以線性速率下降至一個(gè)給定閾值。數(shù)值實(shí)驗(yàn)表明該閾值能夠很好地估計(jì)貪婪隨機(jī)Kaczmarz方法的迭代解誤差所能達(dá)到的最小值?？梢园l(fā)現(xiàn)當(dāng)擾動不是很大且對解的精度要求不是很高時(shí)，貪婪隨機(jī)Kaczmarz方法仍然能夠很好地給出原線性代數(shù)方程組最小范數(shù)解的近似。