求解大型稀疏線性系統(tǒng)的貪婪雙子空間隨機(jī)Kaczmarz方法

荊燕飛，李彩霞，胡少亮

（1.電子科技大學(xué) 數(shù)學(xué)科學(xué)學(xué)院，四川成都 611731；2.中國工程物理研究院 高性能數(shù)值模擬軟件中心，北京 100088）

考慮求解具有如下形式的相容線性系統(tǒng)：

其中，系數(shù)矩陣A∈Rm×n(m＞n)可為滿秩或秩虧矩陣，且b∈Rm。通?？紤]求式（1）的最小范數(shù)解

當(dāng)系數(shù)矩陣A為列滿秩時(shí)，x*為式（1）的惟一解；當(dāng)線性系統(tǒng)（1）有無窮多組解時(shí)，x*為式（1）的最小范數(shù)解（這里的范數(shù)指歐式范數(shù)）。

為求解線性系統(tǒng)（1），許多迭代方法［1-5］已被開發(fā)和進(jìn)一步研究。其中，Kaczmarz［5］于1937年首次提出的Kaczmarz方法，因其成本低廉、易于操作而迅速得到數(shù)值計(jì)算領(lǐng)域的專家學(xué)者的認(rèn)可和廣泛關(guān)注，進(jìn)而使其獲得巨大的理論發(fā)展［6-9］。因其是一種具有代表性的行處理迭代方法且在計(jì)算機(jī)上易于實(shí)現(xiàn)和并行化，因此被廣泛應(yīng)用于計(jì)算機(jī)斷層掃描［10-13］，圖像重建［14-16］，分布式計(jì)算［17-18］和信號(hào)處理［19-21］等領(lǐng)域。

若以隨機(jī)順序而不是以給定順序選用系數(shù)矩陣的行可以極大地提高Kaczmarz方法的收斂速度［13，16，21］。盡管這些隨機(jī)選行的Kaczmarz方法在應(yīng)用中頗具吸引力，但尚不能保證其收斂速率。Strohmer等［22］首次采用選取工作行的概率與其所對(duì)向量的歐式范數(shù)的平方成正比這樣的選行準(zhǔn)則，提出能保證收斂速率的隨機(jī)Kaczmarz方法，并證明其誤差的期望具有指數(shù)收斂速率。Dai等［23］通過求解最小化收斂速度的上限這個(gè)凸優(yōu)化問題來獲得選擇行的最佳概率分布，提出最優(yōu)的隨機(jī)Kaczmarz方法；Bai等［24］結(jié)合貪婪和隨機(jī)的數(shù)學(xué)思想引入一種有效的行選擇準(zhǔn)則提出貪婪隨機(jī)Kaczmarz方法。此后，松弛貪婪隨機(jī)Kaczmarz方法［25］和貪婪距離隨機(jī)Kaczmarz方法［26］也被提出并深入研究。隨機(jī)Kaczmarz方法的各種變種［27-39］和收斂理論［30-35］獲得了大量研究和發(fā)展。另一方面，Leventhal等［36］采用同樣的隨機(jī)思想，提出了求解最小二乘問題的隨機(jī)坐標(biāo)下降方法。因隨機(jī)Kaczmarz方法和隨機(jī)坐標(biāo)下降方法存在一定的共性，后期有學(xué)者將兩者同時(shí)比對(duì)研究［37-38］。

盡管隨機(jī)Kaczmarz方法適用于任何相容線性系統(tǒng)的求解，但當(dāng)系數(shù)矩陣有許多相關(guān)行時(shí)（常出現(xiàn)于地球物理學(xué)中），收斂可能停滯。為克服這一缺點(diǎn)，Needell和Ward采用等可能地隨機(jī)選擇兩個(gè)不同工作行的選擇準(zhǔn)則，提出了隨機(jī)Kaczmarz方法的雙子空間拓展-雙子空間隨機(jī)Kaczmarz方法［27］，并從理論和數(shù)值上說明了該方法至少與隨機(jī)Kaczmarz方法有相同的指數(shù)收斂速率。此外，當(dāng)系數(shù)矩陣高度相干時(shí)，它能極大程度地提高隨機(jī)Kaczmarz方法的收斂速率。此后同時(shí)選用多個(gè)工作行的隨機(jī)塊方法［39-40］也被陸續(xù)提出并加以深入研究。

Nutini等［41］研究表明，采用非等可能概率準(zhǔn)則選取工作行的Kaczmnarz方法收斂速度至少與采用等可能概率準(zhǔn)則選取工作行的Kaczmarz方法收斂速度一樣快。在原始雙子空間隨機(jī)Kaczamrz方法［27］中，等可能地隨機(jī)選擇兩個(gè)不同的工作行，然后將當(dāng)前解投影到由這兩行所確定的解的超平面上獲得下一個(gè)近似解。基于上一次迭代所生成的殘差或基于當(dāng)前解離系數(shù)矩陣各行所形成的超平面的距離來選取工作行都能較大程度地提高Kaczmarz方法的收斂速率［41］。受該思想啟發(fā)，筆者通過引入控制參數(shù)θ建立了一種基于殘差的準(zhǔn)則來選擇工作行，導(dǎo)出一類貪婪雙子空間隨機(jī)Kaczmarz方法。理論證明新方法收斂到相容線性方程組的最小范數(shù)解，而且新方法的理論收斂因子小于原始雙子空間隨機(jī)Kaczmarz方法的收斂因子。數(shù)值實(shí)驗(yàn)驗(yàn)證了新方法的有效性，其在迭代步數(shù)和計(jì)算時(shí)間上均優(yōu)于原始雙子空間隨機(jī)Kaczmarz方法。

1 貪婪雙子空間隨機(jī)Kaczmarz方法

在經(jīng)典的隨機(jī)Kaczmarz方法（RK）中，Strohmer和Vershynin將各行所對(duì)向量的歐式范數(shù)與‖A‖F(xiàn)比值的平方作為概率，然后根據(jù)此概率分布隨機(jī)選取工作行。具體地說，如果定義A(i)代表系數(shù)矩陣A的第i行，b(i)代表向量b的第i個(gè)分量，初始向量為x0，則隨機(jī)Kaczmarz方法的具體過程見算法1。

算法1［22］隨機(jī)Kaczmarz方法。①置k：=0。②根據(jù)概率選取指標(biāo)ik∈{1，2，…，m}。③計(jì)算。④置k=k+1，轉(zhuǎn)步驟②。

算法1中，P(r=ik)=代表選取矩陣A的第ik行作為本次迭代工作行的概率為

關(guān)于隨機(jī)Kaczmarz方法，有如下收斂性定理。

定理1［22，24］若線性系統(tǒng)（1）相容，其中系數(shù)矩陣A∈Rm×n且 右 端 項(xiàng)b∈Rm。初始向量x0∈Ran(AT)，其中Ran(AT)表示AT的列空間，令x k為通過隨機(jī)Kaczmarz方法生成的第k個(gè)迭代值，則

其中σmin(·)表示矩陣的最小非零奇異值。

Needell和Ward發(fā)現(xiàn)當(dāng)線性系統(tǒng)是超定且高度相干時(shí)，隨機(jī)Kaczmarz方法的求解效率低下甚至無效，因此提出同時(shí)啟用兩個(gè)工作行的雙子空間隨機(jī)Kaczmarz方法（2S-RK）去求解這類特殊系統(tǒng)。

為了簡便，Needell和Ward假設(shè)系數(shù)矩陣是標(biāo)準(zhǔn)化的，這意味著其每行都具有單位歐幾里得范數(shù)，在接下來的部分，仍沿用該假設(shè)。在原始雙子空間隨機(jī)Kaczmarz方法中，等可能地隨機(jī)選取兩個(gè)不同工作行，再將當(dāng)前解投影到由兩個(gè)工作行所確定的超平面上獲得下一個(gè)估計(jì)值，具體過程見算法2。

算法2［31］雙子空間隨機(jī)Kaczmarz方法。①置k：=0。②隨機(jī)均勻地選擇行指標(biāo)r k和sk。③計(jì)算μrk，sk=＜A(rk)，A(sk)＞。④計(jì) 算y k=x k+(b(sk)-A(sk)x k)(A(sk))T。⑤計(jì)算νk=。⑥計(jì)算

并以概率P。⑦計(jì)算x k+1=y k+(βkνk y k)(νk)T。置k=k+1，轉(zhuǎn)步②。

關(guān)于雙子空間隨機(jī)Kaczmarz方法，有如下收斂性定理。

定理2［31］若線性系統(tǒng)（1）相容，其中系數(shù)矩陣A∈Rm×n且 右 端 項(xiàng)b∈Rm。初 始 向 量x0∈Ran(AT)，其中Ran(AT)表示AT的列空間，令x k+1為通過雙子空間隨機(jī)Kaczmarz方法生成的第(k+1)個(gè)迭代值，則

這里，e k=x*-x k，μr，s=＜A(r)，A(s)＞。

通過算法2，可知雙子空間隨機(jī)Kaczmarz方法的相應(yīng)誤差向量的2-范數(shù)的平方滿足

若式（3）的最后兩項(xiàng)的和取到最大值，則x*和x k+1之間的距離可以最小化?；谶@個(gè)想法，設(shè)計(jì)最優(yōu)的貪婪距離選行準(zhǔn)則為

若雙子空間隨機(jī)Kaczmarz方法采用式（4）的選行準(zhǔn)則，則其定能以一個(gè)極快的速率收斂到線性系統(tǒng)（1）的解。因此，可采用式（4）的規(guī)則來建立雙子空間隨機(jī)Kaczmarz方法的最佳版本。但在實(shí)踐中計(jì)算滿足式（4）行指標(biāo)r k和sk的成本非常昂貴。為克服其耗時(shí)的缺點(diǎn)，筆者將通過依次構(gòu)造兩個(gè)具有控制參數(shù)θ的行索引集來構(gòu)造次優(yōu)版本的貪婪雙子空間隨機(jī)Kaczmarz方法（2S-GRK（θ）），具體過程見算法3。

首先，構(gòu)造與x k有關(guān)的指標(biāo)集U k選取指標(biāo)sk∈U k，其中向量的定義見算法3的步3。將當(dāng)前解x k投影到解空間{z|A(sk)z=b(sk)}得到中間解向量y k

由此可得，對(duì)于k=1，2，3…成立

其次，構(gòu)造與y k有關(guān)的指標(biāo)集

并以概率P選取指標(biāo)，其中向量的定義見算法3的步8。將當(dāng)前解y k投影到解空間{z|A(rk)z=b(r k)}得到滿足A(sk)x k+1=b(sk)的解向量x k+1

由指標(biāo)集的定義，對(duì)于，可得

對(duì)于k=1，2，3…，同樣可得和，即 是。同理，由指標(biāo)集U k的定義，對(duì)于?s∈U k，可得

算法3貪婪雙子空間隨機(jī)Kaczmarz方法。

步1置k：=0。計(jì)算A(i)x k|2}。定義正整數(shù)指標(biāo)集U k={s||b(s)-A(s)x k|2≥εk}。 計(jì) 算的 第s個(gè) 分 量

步 3 計(jì)算y k=x k+(b(sk)-A(sk)x k)(A(sk))T，

步4 定義正整數(shù)指標(biāo)集={r||b(r)-A(r)y k|2≥?}。計(jì) 算的 第r個(gè) 分 量

步7 計(jì)算x k+1=y k+(βk-νk y k)(νk)T。

步8置k=k+1，轉(zhuǎn)步1。

此外，由y k和x k+1的遞推式可以推得相對(duì)應(yīng)的殘量的遞推公式如下：

其中B=AAT且B(sk)、B(rk)表示矩陣B的第sk、r k列。若迭代前已知AAT，則可進(jìn)一步提高貪婪雙子空間隨機(jī)Kaczmarz方法的求解性能，參見文獻(xiàn)［24］。關(guān)于采用隨機(jī)策略近似計(jì)算矩陣AAT，參見文獻(xiàn)［42］。

若每次迭代僅選擇一個(gè)工作行，則算法3成為一類廣義貪婪隨機(jī)Kaczmarz（GRK（θ））方法。通過改變參數(shù)θ，可得到一系列常用算法作為特殊情況。

時(shí)，GRK（θ）方法分別簡化為貪婪Kaczmarz方法［41］，貪婪距離隨機(jī)Kaczmarz方法［25］，貪婪隨機(jī)Kaczmarz方法［24］和松弛貪婪隨機(jī)Kaczmarz方法［33］。

2 貪婪雙子空間隨機(jī)Kaczmarz方法的收斂性分析

相容線性系統(tǒng)（1）的最小范數(shù)解x*形如

其中Ran(AT)表示AT的列空間。若初始向量x0∈Ran(AT)，由算法3的步5和步13可知，算法3生成的迭代序列一定在Ran(AT)中，故若算法3收斂，則一定收斂到相容線性系統(tǒng)（1）的最小范數(shù)解。

關(guān)于算法3特殊情況（每次迭代僅選一個(gè)工作行）——廣義貪婪隨機(jī)Kaczmarz方法，有如下收斂性定理。

定理3若線性系統(tǒng)（1）相容，其中系數(shù)矩陣A∈Rm×n且 右 端 項(xiàng)b∈Rm。初 始 向 量x0∈Ran(AT)，令x k+1為通過廣義貪婪隨機(jī)Kaczmarz方法生成的第(k+1)個(gè)迭代值，則有

和

證明：由廣義貪婪隨機(jī)Kaczmarz方法（算法3的前5步）可得

對(duì)等式（9）兩側(cè)取期望，有

其中式（11）最后一個(gè)不等式需要用到不等式（12），即對(duì)于任意u∈Ran(AT)，成立

聯(lián)立等式（11）和不等式（12）可得

通過歸納法，可得

關(guān)于貪婪雙子空間隨機(jī)Kaczmarz方法，有如下的收斂性定理。

定理4若線性系統(tǒng)（1）相容，其中系數(shù)矩陣A∈Rm×n且 右 端 項(xiàng)b∈Rm。初 始 向 量x0∈Ran(AT)，令x k+1為通過貪婪雙子空間隨機(jī)Kaczmarz方法生成的第(k+1)個(gè)迭代值，則

證明：由算法3的步5和步13可得

記ξk=＜A(r k)，νk＞。為了便利，在后續(xù)的證明中，A(r k)、A(sk)、μk、νk、ξk分別簡記為A(r)、A(s)、μ、ν、ξ。由于＜A(s)，ν＞=0，即向量ν正交于向量A(s)，展開可得

在算法3中，一次迭代需執(zhí)行兩次投影。因此，將式（14）中的誤差與廣義貪婪隨機(jī)Kaczmarz方法兩次迭代后獲得的誤差進(jìn)行比較。設(shè)z是廣義貪婪隨機(jī)Kaczmarz方法基于y k的下一次迭代值，則有

由y k的遞推式可得

由μ、ν、ξ的定義可得

聯(lián)立式（15）～（17）可得

利用向量ν與向量A(s)的正交性，展開可得

聯(lián)立式（14）和式（18）可得

對(duì)式（19）兩端取期望可得

由定理3可得

再由

因此，可得

再由期望的定義可得

把式（22）和式（23）帶入式（20）可得式（13）。

把式（16）和式（17）帶入式（23），同時(shí)再結(jié)合式（6）可得

收斂因子越小，方法收斂得越快。通過定理4收斂界可知，參數(shù)θ越小，收斂因子M k(θ)越小，即貪婪雙子空間隨機(jī)Kaczmarz方法收斂速度越快。為便于比較，用θ取值特殊的貪婪雙子空間隨機(jī)Kaczmarz方法與原雙子空間隨機(jī)Kaczmarz方法作比較，但對(duì)實(shí)際問題，因貪婪雙子空間隨機(jī)Kaczmarz方法選行概率更具優(yōu)勢性，對(duì)于?θ∈[0，1]，貪婪雙子空間隨機(jī)Kaczmarz方法都比原雙子空間隨機(jī)Kaczmarz方法收斂更快。

事實(shí) ①令t1，t2，…，tl是在概率分別為p1，p2，…，pl的某個(gè)概率空間上定義的一組隨機(jī)變量。若越大的|ti|所對(duì)應(yīng)的概率pi越大，則值越大；②對(duì)于任意數(shù)組Γ={t1，t2，…，tl}，ti∈R，ave(Γ)表示數(shù)組Γ的平均值。令?={ti≥ave(Γ)，ti∈Γ}，有ave(?)≥ave(Γ)。

基于事實(shí)①和②可得

因此，貪婪雙子空間隨機(jī)Kaczmarz方法收斂因子小于原雙子空間隨機(jī)Kaczmarz方法的收斂因子。

浮點(diǎn)運(yùn)算量方面，原雙子空間隨機(jī)Kaczmarz方法每次迭代需要12n+4個(gè)浮點(diǎn)運(yùn)算量。若殘差由遞推公式計(jì)算，則貪婪雙子空間隨機(jī)Kaczmarz方法每次迭代需要16m+6n+7個(gè)浮點(diǎn)運(yùn)算量。因此，當(dāng)時(shí)，貪婪雙子空間隨機(jī)Kaczmarz方法計(jì)算量小于原雙子空間隨機(jī)Kaczmarz方法。

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

在本節(jié)中，分別通過數(shù)值實(shí)驗(yàn)比較隨機(jī)Kaczmarz方法（RK）、雙子空間隨機(jī)Kaczmarz方法（2S-RK）和貪婪雙子空間隨機(jī)Kaczmarz方法（2SGRK（θ）），并顯示后者無論在迭代步數(shù)（IT）還是計(jì)算時(shí)間（CPU）上均優(yōu)于前兩者。這里，迭代步數(shù)（IT）和計(jì)算時(shí)間（CPU）取30次計(jì)算結(jié)果的平均值。

在實(shí)驗(yàn)中，通過Matlab函數(shù)randn隨機(jī)生成解x*∈Rn，右端項(xiàng)b∈Rm由Ax*給出。此外，需要指出2S-GRK（θ）方法按照算法3中定義的過程精確執(zhí)行，并沒有顯式地計(jì)算矩陣B=AAT來計(jì)算殘量。本節(jié)的所有數(shù)值實(shí)驗(yàn)，均取初始向量為x0=0，停機(jī)準(zhǔn)則為近似解的相對(duì)誤差RSE滿足

或者迭代步數(shù)超過30萬。在實(shí)驗(yàn)表格中，若在30萬步內(nèi)未達(dá)到指定精度，即用“--”表示。

為了測試這些方法，選取了兩類矩陣進(jìn)行實(shí)驗(yàn)。一類是通過Matlab函數(shù)unifrnd生成的500×100的相干矩陣。矩陣元素是在區(qū)間[d，1]上獨(dú)立且均勻分布的隨機(jī)變量。通過改變d的值，可以得到具有不同相干系數(shù)對(duì)的矩陣，相干系數(shù)對(duì)(σ，Δ)定義如下

表1列出了測試的相干矩陣的相關(guān)信息。包含不同相干程度的500×100矩陣。

表1 相干矩陣信息Tab.1 Information of coherent matrices

對(duì)于測試的矩陣，表2給出了RK，2S-RK，2SGRK（θ）3種方法的迭代步數(shù)（IT）與計(jì)算時(shí)間（CPU），以及2S-GRK（θ）方法對(duì)于2S-RK方法的加速比。

表2 相干矩陣數(shù)值結(jié)果Tab.2 Numerical results of coherent matrices

由表2可知，2S-RK與2S-GRK（θ）方法總能計(jì)算出符合精度的解，而RK方法有時(shí)在迭代步數(shù)達(dá)到30萬步后仍不能計(jì)算出符合精度的解。而且即便RK方法收斂，2S-RK與2S-GRK（θ）方法在迭代步數(shù)與計(jì)算時(shí)間上也均少于RK方法。此外，2SGRK（θ）方法進(jìn)一步優(yōu)化了2S-RK方法，其關(guān)于2S-RK方法迭代時(shí)間的加速比最大可以達(dá)到3.64，最小可以達(dá)到2.48。值得注意的是，無論矩陣的相干程度如何，2S-GRK（θ）方法均收斂，且其迭代步數(shù)與計(jì)算時(shí)間均優(yōu)于原始的2S-RK方法。

圖1描繪了相干矩陣A2（圖1a）和A4（圖1b）的近似解的相對(duì)誤差RSE以10為底的對(duì)數(shù)隨著迭代步數(shù)變化的曲線，進(jìn)一步驗(yàn)證了2S-GRK（θ）方法比經(jīng)典的2S-RK方法收斂更快。

圖1 以A2和A4為系數(shù)矩陣的線性系統(tǒng)的近似解的相對(duì)誤差隨著迭代步數(shù)變化的曲線Fig.1 lg(R SE)versus IT for linear systems with coefficient matrices:A2 and A4

表3列出了測試的稀疏矩陣的相關(guān)信息。

表3 稀疏矩陣信息Tab.3 Information of sparsematrices

對(duì)于測試的稀疏矩陣，表4給出了RK、2S-RK、2S-GRK（θ）3種方法的迭代步數(shù)（IT）與計(jì)算時(shí)間（CPU），以及2S-GRK（θ）方法對(duì)于2S-RK方法的加速比。

由表4可得與表2類似結(jié)論，即2S-GRK（θ）與2S-RK方法在迭代步數(shù)與計(jì)算時(shí)間上也均少于RK方法。此外，2S-GRK（θ）方法進(jìn)一步優(yōu)化了2S-RK方法，其關(guān)于2S-RK方法迭代時(shí)間的加速比最大可以達(dá)到6.17，最小可以達(dá)到1.75。因此，無論對(duì)于相干矩陣還是稀疏（不相干）矩陣，2S-GRK（θ）方法均收斂，且其迭代步數(shù)與計(jì)算時(shí)間均優(yōu)于傳統(tǒng)的2S-RK方法。

表4 稀疏矩陣數(shù)值結(jié)果Tab.4 Numerical results of sparse matrices

圖2描繪了矩陣Ash219（圖2a）和Ash958（圖2b）的近似解的相對(duì)誤差以10為底的對(duì)數(shù)隨著迭代步數(shù)變化的曲線，進(jìn)一步驗(yàn)證了2S-GRK（θ）方法方法比經(jīng)典的2S-RK方法收斂更快。

圖2 以Ash219和Ash958為系數(shù)矩陣的線性系統(tǒng)的近似解的相對(duì)誤差隨著迭代步數(shù)變化的曲線Fig.2 lg(R SE)versus IT for linear systems with coefficient matrices:Ash219 and Ash958

對(duì)于4個(gè)測試矩陣A2、A4、Ash219和Ash958，圖3描繪了2S-GRK（θ）方法的計(jì)算時(shí)間（圖3a）和迭代步數(shù)（圖3b）隨控制參數(shù)θ變化的曲線。由圖3可知，只要在計(jì)算時(shí)采用合適的控制參數(shù)θ，2SGRK（θ）方法的性能就會(huì)得到很大提升。結(jié)合表2和表4可知，無論θ∈[0，1]取何值，2S-GRK（θ）方法方法總比經(jīng)典的2S-RK方法收斂更快。

圖3 2S-GRK（θ）方法的計(jì)算時(shí)間和迭代步數(shù)隨控制參數(shù)θ變化的曲線Fig.3 CPU and IT versusθof 2S-GRK(θ)

4 結(jié)論

提出一類貪婪雙子空間隨機(jī)Kaczmarz方法，理論分析證明新方法的收斂性，還表明收斂因子小于傳統(tǒng)的雙子空間隨機(jī)Kaczmarz方法的收斂因子。數(shù)值實(shí)驗(yàn)結(jié)果也表明所提出的新方法在迭代步數(shù)和計(jì)算時(shí)間上均優(yōu)于傳統(tǒng)的雙子空間隨機(jī)Kaczmarz方法。貪婪雙子空間隨機(jī)Kaczmarz方法能夠快速求解大規(guī)模稀疏相容線性方程組，主要原因在于新方法定義的隨機(jī)選取工作行的概率更具優(yōu)勢。如何隨機(jī)選取工作行以加速Kaczmarz方法及其應(yīng)用［43］仍然是一個(gè)值得研究的問題。

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

荊燕飛：核心思想提煉、論文修改。

李彩霞：論文撰寫。

胡少亮：論文修改。