大型稀疏線性系統的一類含參數的貪心隨機Kaczmarz算法

劉 永,,崔蓉蓉

(1.常州信息職業技術學院,江蘇常州213164;2.上海大學理學院,上海200444;3.鹽城師范學院數學與統計學院,江蘇鹽城224002)

對于大規模線性系統

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

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

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

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

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

1 貪心隨機Kaczmarz算法

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

輸入:A,b,?和x0

輸出:x?

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

2:計算

3:構造集合

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

6:計算xk+1=.

7:end for

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

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

那么,可以直接得到

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

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

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

2 含參數的貪心隨機Kaczmarz算法

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

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

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

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

輸出:x?

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

2:計算

3:構造集合

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

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

7:end for

關于算法2,有如下收斂性定理.

定理1令x?=A?b是相容線性系統(1)的解,這里A?是矩陣A的Moore-Penrose逆,x0∈R(A?),R(A?)表示A?的列空間,ω∈(0,2).那么由算法2迭代生成的序列依期望指數收斂到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)得到的估計式還只是一個粗糙的結果,故還可以對?k作進一步的估計：

對上式兩邊取全期望可得

最后,可得

3 數值實驗

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

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

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

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

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

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

圖1 當線性系統(1)的系數矩陣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 當線性系統(1)的系數矩陣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 結束語

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