大型稀疏線性系統(tǒng)的一類含參數(shù)的貪心隨機(jī)Kaczmarz算法

劉 永,,崔蓉蓉

(1.常州信息職業(yè)技術(shù)學(xué)院,江蘇常州213164;2.上海大學(xué)理學(xué)院,上海200444;3.鹽城師范學(xué)院數(shù)學(xué)與統(tǒng)計學(xué)院,江蘇鹽城224002)

對于大規(guī)模線性系統(tǒng)

式中:A∈Cm×n為復(fù)數(shù)域上的m×n矩陣;b∈Cm為復(fù)數(shù)域上的m維已知向量;x∈Cn為復(fù)數(shù)域上的n維未知向量.Kaczmarz(K)算法[1]是計算其逼近解的一種非常有效的迭代算法.K算法由于其簡單且收斂速度快,故被廣泛應(yīng)用于圖像重構(gòu)[2-4]、信號處理[5]和分布式計算[6]等領(lǐng)域.K算法的迭代公式為

式中:(·)?為矩陣或向量的共軛轉(zhuǎn)置;‖·‖2為向量的Euclidean范數(shù);A(ik)為矩陣A的第ik行;b(ik)為向量b的第ik個元素,ik=(k mod m)+1.從ik的選擇方式可以看出,K算法的收斂速率依賴于矩陣A的行順序,有時一個特定的行順序可能會導(dǎo)致算法收斂非常慢[7-8].

為了提高K算法的收斂速率,Strohmer等[9]在2009年提出了隨機(jī)Kaczmarz(randomized Kaczmarz,RK)算法,該算法依概率P(row=ik)=的大小隨機(jī)選擇系統(tǒng)(1)中系數(shù)矩陣A的某一行,并按迭代式(2)進(jìn)行計算,其中表示矩陣A的Frobenius范數(shù).Strohmer等還首次證明了當(dāng)相容線性系統(tǒng)(1)中系數(shù)矩陣A為列滿秩且m≥n,或行滿秩且m≤n時,RK算法依期望指數(shù)收斂到系統(tǒng)(1)的唯一最小二乘解或唯一最小范數(shù)解,其收斂速率為

式中:x?為系統(tǒng)(1)的解;λmin(A?A)和tr(A?A)分別為矩陣A?A的最小非零特征值和跡.從RK算法的行選擇方式可以看出,其收斂速率并不取決于方程的個數(shù),甚至可以不需要知道整個方程組,只是通過隨機(jī)選取整個方程組的一部分進(jìn)行計算就可以得到系統(tǒng)(1)的解.因此,該算法就特別適用于規(guī)模比較大的系統(tǒng),單從這點來看,RK算法優(yōu)于K算法.但是RK算法中行的隨機(jī)選擇方式也不一定是最優(yōu)的[10].例如,當(dāng)為常數(shù)時,RK算法將等概率地選取A的行,此時算法的收斂速率將大打折扣.2018年,Bai等[11]提出了一種新的概率準(zhǔn)則,即在RK算法的每個迭代步中,盡可能地選取線性系統(tǒng)(1)的殘差向量中較大項所對應(yīng)的行進(jìn)行迭代計算,并在此基礎(chǔ)上構(gòu)造了貪心隨機(jī)Kaczmarz(greedy randomized Kaczmarz,GRK)算法,該算法在理論和實驗上都比RK算法收斂得更快.

本工作在GRK算法的迭代公式中引入松弛因子ω,構(gòu)造了求解系統(tǒng)(1)的含參數(shù)的GRK算法.此外,本工作還給出了該算法的收斂性證明及收斂速率的一個上界.最后,通過數(shù)值實驗證明了在適當(dāng)選擇松弛參數(shù)的情況下,本工作所提出的算法的收斂上界可小于文獻(xiàn)[11]中算法的收斂上界.

1 貪心隨機(jī)Kaczmarz算法

算法1貪心隨機(jī)Kaczmarz(GRK)算法[11].

輸入:A,b,?和x0

輸出:x?

1:for k=0,1,···,??1 do

2:計算

3:構(gòu)造集合

5:依概率Pr(row=ik)=選擇ik∈Uk,

6:計算xk+1=.

7:end for

算法1之所以比RK算法能更高效地計算線性系統(tǒng)(1)的解,主要有如下2個方面的原因.

(1)在算法1的每一步迭代過程中,第3步中的集合Uk都是非空的,這是算法1能成功的保證.為了理解這一點,不妨設(shè)

那么,可以直接得到

式中:?k如算法1的第2步所定義.于是對任意的k∈{0,1,···}都有

換句話說,集合Uk都是非空的.這樣,算法1便能持續(xù)迭代下去.

(2)算法1比RK算法收斂速度更快,主要是因為該算法引入了一個實用的、合適的概率準(zhǔn)則(算法1的第4～5步),能在算法的每個迭代步中盡可能地選取線性系統(tǒng)(1)的殘差向量中較大項所對應(yīng)的行進(jìn)行迭代計算,這與線性方程組的求解思想是一致的,即先迭代計算方程組殘差向量中最大項所在的行,這樣才能以最快的速度達(dá)到方程組求解的精度要求.Bai等在理論上也證明了GRK算法的收斂速率快于RK算法,具體可參見文獻(xiàn)[11]中的注解3.2.

2 含參數(shù)的貪心隨機(jī)Kaczmarz算法

含參數(shù)的GRK算法(GRK(ω))實際上是在算法1的第6步中引入松弛因子ω,也即在GRK算法中采用如下迭代公式:

進(jìn)行計算,這里ω∈(0,2)為給定的松弛因子.因此,可得如下算法2.

算法2含參數(shù)的貪心隨機(jī)Kaczmarz算法(GRK(ω)).

輸入:A,b,?,ω和x0

輸出:x?

1:for k=0,1,···,??1 do

2:計算

3:構(gòu)造集合

5:依概率Pr(row=ik)=選擇ik∈Uk,

6:計算xk+1=(1?ω)xk+ω.

7:end for

關(guān)于算法2,有如下收斂性定理.

定理1令x?=A?b是相容線性系統(tǒng)(1)的解,這里A?是矩陣A的Moore-Penrose逆,x0∈R(A?),R(A?)表示A?的列空間,ω∈(0,2).那么由算法2迭代生成的序列依期望指數(shù)收斂到x?,且解誤差依期望滿足

和

這里,用Ek表示前k次迭代的期望值,即Ek[·]=Ek[·|i0,i1,···,ik?1],這里il(l=0,1,···,k-1)表示算法在第l次迭代時選擇的是第il行.同時也易知E[Ek[·]]=E[·].

證明 令Pik=,則

因此,

另一方面,由Pik的定義容易驗證,于是可得

對上式兩邊同時取期望得

式(4)的第1個不等式可由不等式(3)直接得到,第2個不等式可由下面的估計式得到:

實際上,由式(4)得到的估計式還只是一個粗糙的結(jié)果,故還可以對?k作進(jìn)一步的估計：

對上式兩邊取全期望可得

最后,可得

3 數(shù)值實驗

用GRK算法和GRK(ω)算法求解相容線性系統(tǒng)Ax=b,A∈Rm×n主要有2個來源:①從文獻(xiàn)[12]中選取6個具有應(yīng)用背景的稀疏矩陣,分別是WorldCities、football及Stranke94(Pajek網(wǎng)絡(luò)問題)和bibd?16?8、relat6以及bibd?17?8(組合問題);②用Matlab函數(shù)randn生成的隨機(jī)矩陣A=randn(4 000,1 500)和A=randn(1 500,4 000).此外,設(shè)定線性系統(tǒng)Ax=b的解x?=randn(n,1),右端向量b=Ax?.本算法的數(shù)值實驗是在內(nèi)存為4G,主頻為2.67 GHz,處理器為Intel(R)Core(TM)i5 CPU的電腦上用Matlab(R2016b)進(jìn)行計算的.考慮到隨機(jī)性,這里所有計算結(jié)果取10次獨立實驗結(jié)果的平均值.所有數(shù)值實驗的初始向量均設(shè)為x0=0,停機(jī)準(zhǔn)則設(shè)為相對解誤差(relative solution error,RSE),滿足

對于第一類矩陣,表1中列出了GRK算法和GRK(ω)算法在求解相容線性系統(tǒng)Ax=b達(dá)到精度要求時所需的迭代步(用IT表示)和計算時間(用CPU表示).同時,表1給出了測試矩陣的一些基本信息,如稠密性

和矩陣條件數(shù)cond(A).表1中,加速比

表示GRK(ω)算法相對于GRK算法在計算時間這一指標(biāo)上的加速程度,其數(shù)值越大表明GRK(ω)算法對同一測試矩陣加速得越快.從表1中可以看出,對于矩陣A,無論m≥n還是m＜n,當(dāng)選取合適的參數(shù)ω時,GRK(ω)算法在迭代次數(shù)和計算時間方面都優(yōu)于GRK算法,并且有顯著的提速,加速比(speed-up)最高可取3.03(見矩陣Stranke94).

表1 在不同系數(shù)矩陣下求解線性系統(tǒng)(1)的GRK算法和GRK(ω)算法的數(shù)值結(jié)果Table 1 Numerical results of GRK and GRK(ω)algorithms for solving linear system(1)under different coefficient matrices A

對于第二類矩陣,分別描繪了GRK算法和GRK(ω)算法對于2種不同的系數(shù)矩陣A,log10(RSE)隨迭代步和計算時間的變化曲線圖(見圖1和2).從圖中可以看出,本算法是可行的，并且在選擇適當(dāng)?shù)膮?shù)后比GRK算法收斂得更快.

圖1 當(dāng)線性系統(tǒng)(1)的系數(shù)矩陣A=randn(4 000,1 500)時,GRK算法和GRK(1.5)算法的log10(RSE)的變化曲線Fig.1 Curves of log10(RSE)for GRK and GRK(1.5)algorithms when the coefficient matrix of linear system(1)is A=randn(4 000,1 500)

圖2 當(dāng)線性系統(tǒng)(1)的系數(shù)矩陣A=randn(1 500,4 000)時,GRK算法和GRK(1.5)算法的log10(RSE)的變化曲線Fig.2 Curves of log10(RSE)for GRK and GRK(1.5)algorithms when the coefficient matrix of linear system(1)is A=randn(1 500,4 000)

4 結(jié)束語

本工作在貪心隨機(jī)Kaczmarz(GRK)算法的基礎(chǔ)上,提出了一類含參數(shù)的GRK(GRK(ω))算法,分析了算法的收斂性.數(shù)值結(jié)果表明,在選擇適當(dāng)參數(shù)ω后,GRK(ω)算法在迭代步和計算時間方面明顯優(yōu)于GRK算法.在數(shù)值實驗中,對于不同的系數(shù)矩陣參數(shù)ω是通過多次實驗使得迭代解在達(dá)到停機(jī)準(zhǔn)則時所需迭代步相對最少來確定的,自適應(yīng)的參數(shù)選擇策略(ωk隨迭代步的變化而變化)可能會獲得更好的數(shù)值結(jié)果.