荊燕飛，李彩霞，胡少亮

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

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

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

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

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

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

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

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

1 貪婪雙子空間隨機Kaczmarz方法

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

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

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

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

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

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

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

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

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

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

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

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

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

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

若式（3）的最后兩項的和取到最大值，則x*和x k+1之間的距離可以最小化。基于這個想法，設(shè)計最優(yōu)的貪婪距離選行準則為

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

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

由此可得，對于k=1，2，3…成立

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

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

由指標集的定義，對于，可得

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

算法3貪婪雙子空間隨機Kaczmarz方法。

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

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

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

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

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

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

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

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

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

2 貪婪雙子空間隨機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特殊情況（每次迭代僅選一個工作行）——廣義貪婪隨機Kaczmarz方法，有如下收斂性定理。

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

和

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

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

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

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

通過歸納法，可得

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

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

證明：由算法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）中的誤差與廣義貪婪隨機Kaczmarz方法兩次迭代后獲得的誤差進行比較。設(shè)z是廣義貪婪隨機Kaczmarz方法基于y k的下一次迭代值，則有

由y k的遞推式可得

由μ、ν、ξ的定義可得

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

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

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

對式（19）兩端取期望可得

由定理3可得

再由

因此，可得

再由期望的定義可得

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

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

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

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

基于事實①和②可得

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

浮點運算量方面，原雙子空間隨機Kaczmarz方法每次迭代需要12n+4個浮點運算量。若殘差由遞推公式計算，則貪婪雙子空間隨機Kaczmarz方法每次迭代需要16m+6n+7個浮點運算量。因此，當時，貪婪雙子空間隨機Kaczmarz方法計算量小于原雙子空間隨機Kaczmarz方法。

3 數(shù)值實驗

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

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

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

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

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

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

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

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

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

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

圖1 以A2和A4為系數(shù)矩陣的線性系統(tǒng)的近似解的相對誤差隨著迭代步數(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

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

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

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

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

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

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

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

4 結(jié)論

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

作者貢獻聲明：

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

李彩霞：論文撰寫。

胡少亮：論文修改。