求解大型稀疏線性方程組的貪婪距離隨機Kaczmarz方法

杜亦疏，殷俊鋒，張 科

（1.同濟大學(xué)數(shù)學(xué)科學(xué)學(xué)院，上海200092；2.上海海事大學(xué)文理學(xué)院，上海201306）

考慮求解具有如下形式的大規(guī)模稀疏線性方程組

其中系數(shù)矩陣A∈Rm×n可為滿秩或秩虧矩陣，且b∈Rm.當(dāng)線性方程組（1）不相容時，通?？紤]求其最小二乘解

當(dāng)線性方程組（1）相容時，通?？紤]求其最小歐氏范數(shù)解

問題（2）等價于求解等價的正規(guī)方程組

為了求解大型稀疏線性方程組（1），許多經(jīng)典的迭代方法陸續(xù)被提出來并加以深入研究，例如，Kaczmarz迭代方法［1］、超松弛迭代方法（SOR）［2-3］、共軛梯度法（CG）［4］、廣義極小殘量法（GMRES）［5］和穩(wěn)定化的雙共軛梯度法（BiCGSTAB）［6］等.關(guān)于這些迭代法理論及其應(yīng)用的綜述，參見文獻(xiàn)［7-13］.

自從20世紀(jì)30年代波蘭數(shù)學(xué)家Stefan Kaczmarz提出Kaczmarz方法后，該方法迅速得到數(shù)值計算專家和學(xué)者的廣泛關(guān)注，其迭代方法的各種變體和收斂理論獲得了大量的研究和發(fā)展［14-15］.由于該方法是一種具有代表性的行處理迭代方法，其行投影方法很容易在計算機上實現(xiàn)和并行計算，因此被廣泛應(yīng)用在分布式計算、電子計算機斷層掃描、數(shù)字信號處理、人工智能和圖像恢復(fù)等科學(xué)與工程計算領(lǐng)域［16-20］.

Strohmer和Vershynin［21］在2009年首次提出采用隨機的思想來選取工作行，并證明了這種隨機Kaczmarz方法具有線性收斂速率，這再度使得Kaczmarz方法成為數(shù)值代數(shù)領(lǐng)域的一個研究熱點.Bai等［22］結(jié)合貪婪和隨機的數(shù)學(xué)思想，提出了一種新的且更加有效的概率準(zhǔn)則來選取工作行，即貪婪隨機Kaczmarz方法.此后，松弛貪婪隨機Kaczmarz方法［23］和部分隨機拓展的Kaczmarz方法［24］也被陸續(xù)提出并加以深入研究.Popa［25］總結(jié)了各類Kaczmarz方法的收斂速率.更多的隨機Kaczmarz方法及其理論可以參見文獻(xiàn)［26］.另一方面，Bai等［27］利用同樣的數(shù)學(xué)思想，提出了一類貪婪隨機坐標(biāo)下降方法，用以求解列滿秩超定線性代數(shù)方程組的最小二乘問題，可以證明新方法較原有的隨機坐標(biāo)下降方法［28］具有更快的收斂速率.

在研究中發(fā)現(xiàn)，提高隨機Kaczmarz方法收斂速率的重點在于能否提出一類從系數(shù)矩陣中選取工作行的有效概率準(zhǔn)則.在原始的Kaczamrz方法中，工作行被循環(huán)選取，然后將當(dāng)前點投影到工作行所形成的超平面上.在文獻(xiàn)［29］中，貪婪Kaczmarz方法的工作行選取辦法被分為兩類，第一類是選取上一次迭代后各殘差分量中模長最大的行作為本次迭代的工作行，第二類是選取當(dāng)前迭代解離系數(shù)矩陣各行所形成的超平面的距離最大的行作為本次迭代的工作行.受該思想啟發(fā)，筆者提出一類基于距離的貪婪隨機Kaczmarz方法（GDRK）.

1 貪婪距離隨機Kaczmarz方法

在原始的Kaczmarz方法中，工作行被循環(huán)選取，然后將當(dāng)前點投影到工作行所形成的超平面上，更具體地說，如果定義A()i代表系數(shù)矩陣A的第i行，b()i代表向量b的第i個分量，初始向量為x0，則Kaczmarz方法如算法1所示.

算法1［1］Kaczmarz方法.①置k：=0.計算i k=(kmodm)+1.②計算x k+1=x k+③置k=k+1，轉(zhuǎn)步驟①.其中，(?)T代表矩陣或向量的轉(zhuǎn)置.在Kaczmarz方法中，通過按順序選取工作行，然后將當(dāng)前點x k正交投影到超平面A(i k)x k=b(i k)上.

在文獻(xiàn)［14］中，Ansorge首次采用貪婪的思想來選取工作行，即選取上一次迭代后各殘差分量中模長最大的行作為本次迭代的工作行，提出貪婪殘差Kaczmarz方法，并將Cimmino方法與Kaczmarz方法的聯(lián)系做了詳盡的分析，具體過程見算法2.

算法2［14］貪婪殘差Kaczmarz方法.①置k：=0.計 算i k=arg maix②計算x k+1=x k+③置k=k+1，轉(zhuǎn)步驟①.

Nutini等［29］采用另一類貪婪選取工作行的辦法，即選取當(dāng)前迭代解離系數(shù)矩陣各行所形成的超平面的距離最大的行作為本次迭代的工作行，提出貪婪距離Kaczmarz方法.論文分析收斂性并嘗試?yán)谜粓D提出快速的隨機Kaczmarz方法，具體過程見算法3.

算法3［29］貪婪距離Kaczmarz方法.①置k：=0.計 算i k=arg maix②計算x k+1=x k+③置k=k+1，轉(zhuǎn)步驟①.

在文獻(xiàn)［21］中，Strohmer和Vershynin將各行2-范數(shù)與‖ ‖AF比值的平方作為概率，然后根據(jù)此概率隨機選取本次迭代的工作行，提出隨機Kaczmarz方法，具體過程見算法4.

算法4［21］隨機Kaczmarz方法.①置k：=0.根據(jù)概率Pr(row=i k)=選取指標(biāo)i k∈{1，2，…，m}.②計算x k+1=x k+③置k=k+1，轉(zhuǎn)步驟①.

在算法4中，Pr(?)代表選取系數(shù)矩陣第i k行作為本次迭代工作行的概率.

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

定理1［21-22］線性方程組（1）相容，其中系數(shù)矩陣A∈Rm×n且右端項b∈Rm.初始向量x0∈Rn在AT的列空間中，則通過隨機Kaczmarz方法生成的迭代序列在期望上收斂到其最小范數(shù)解x?=A?b，且迭代解序列的誤差的數(shù)學(xué)期望滿足

Bai等［22］提出了一種新的概率準(zhǔn)則，即先選取當(dāng)前點x k離系數(shù)矩陣各行所形成的超平面的距離較大的行所對應(yīng)的指標(biāo)i k，將其放入集合U中，然后定義集合U中每行被選取的概率為各殘差分量模長的平方除以殘差歐式范數(shù)的平方，具體過程見算法5.

算法5［22］貪婪隨機Kaczmarz方法.①置k：=②定義正整數(shù)指標(biāo)集U k=A(i k)x k|2≥εk‖ ‖b-A x k③計算的第i個分 量=④根據(jù)概率Pr(row=i k)=選取指標(biāo)i k∈U k.⑤計算x k+1=x k+⑥置k=k+1，轉(zhuǎn)步驟①.

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

定理2［22］若線性方程組（1）相容，其中系數(shù)矩陣A∈Rm×n且右端項b∈Rm.初始向量x0∈Rn在AT的列空間中，則通過貪婪隨機Kaczmarz方法生成的迭代序列在期望上收斂到其最小范數(shù)解x?=A?b，且迭代解序列的誤差的數(shù)學(xué)期望滿足

和

其中γ=更一般地有

在算法5中，希望指標(biāo)集U k中對應(yīng)較大距離的指標(biāo)i k可以以較大概率被取到.然而在算法5的第④步，選取指標(biāo)i k的概率定義為正比于各殘差分量模長的平方.假設(shè)集合中指標(biāo)i k所對應(yīng)的行具有較大殘量，則它有較大概率被取為工作行；若同時該指標(biāo)對應(yīng)的行范數(shù)很大，則當(dāng)前點投影到此工作行所形成的超平面的距離較小，這與算法5第②步中指標(biāo)集U k的定義相矛盾，此時需要對這2個條件進行平衡.在該想法的啟發(fā)下，提出基于距離的貪婪隨機Kaczmarz方法，具體過程見算法6.

算法6貪婪距離隨機Kaczmarz方法.①置k：=0.計算

②定義正整數(shù)指標(biāo)集

④根據(jù)概率Pr(row=i k)=選取指標(biāo)i k∈U k.⑤計算x k+1=x k+置k=k+1，轉(zhuǎn)步驟①.

由于當(dāng)前點x k到所有系數(shù)矩陣行所形成的超平面的最大距離平方滿足下式：

可知算法6中步驟②指標(biāo)集U k非空，具體參見文獻(xiàn)［22］.

由于計算殘量具有如下遞推式：

其中B=AAT且B(i k)表示矩陣B的第i k列.也就是說，如果在迭代之前已知AAT，可以進一步提高貪婪距離隨機Kaczmarz方法的運行效率，參見文獻(xiàn)［22］.關(guān)于采用隨機策略近似計算矩陣AAT，參見文獻(xiàn)［30］.

2 貪婪距離隨機Kaczmarz方法的收斂性分析

首先說明若貪婪距離隨機Kaczmarz方法收斂，則其一定可以收斂到最小范數(shù)解.假設(shè)x?為滿足相容 系 統(tǒng)（1）的 最 小 范 數(shù) 解，由（2）可 知，x??A?b∈R(AT)，其中，R(AT)表示AT列空間.另一方面，若初始向量x0在AT列空間中，由算法6中步驟⑤可知，算法6生成的迭代序列一定在AT列空間中，故若算法收斂，則一定收斂到相容系統(tǒng)（1）的最小范數(shù)解.

關(guān)于貪婪距離隨機Kaczmarz方法，有如下的收斂性定理.

定理3若線性方程組（1）相容，其中系數(shù)矩陣A∈Rm×n且右端項b∈Rm.初始向量x0∈Rn在AT的列空間中，則通過貪婪隨機Kaczmarz方法生成的迭代序列在期望上收斂到其最小范數(shù)解x?=A?b，且迭代解序列的誤差的數(shù)學(xué)期望滿足

證明：由算法6步驟⑤，可得

即向量x k+1-x k平行于向量(A(i k))T.又由

可得向量x k+1-x?正交于向量(A(i k))T，即正交于向量x k+1-x k，故有下式成立

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

其中，式（5）最后一個不等式需要用到不等式（6），即對于在AT列空間中任意u∈Rn，有

由算法6中εk定義，對于k=1，2，…，有

將不等式（7）和（8）代入不等式（5）中，通過歸納法，有

記定理1中隨機Kaczmarz方法的收斂因子為

定理2中貪婪隨機Kaczmarz方法的收斂因子為

定理3中貪婪距離隨機Kaczmarz方法的收斂因子為

一般而言，收斂因子越小，收斂越快.這里需要特別指出的是，貪婪距離隨機Kaczmarz方法和貪婪隨機Kaczmarz方法在理論上具有相同的收斂因子，但是對于實際問題，由于貪婪距離隨機Kaczmarz方法選取工作行的概率更加符合指標(biāo)集U k定義，故實際計算時新方法有可能收斂更快.

在浮點運算量方面，隨機Kaczmarz方法每次迭代需要4n+1個浮點運算量.若殘差由遞推公式計算，貪婪隨機Kaczmarz方法以及貪婪距離隨機Kaczmarz方法每次迭代需要7m+2n+2個浮點運算量［22］.因此，當(dāng)m＞時，貪婪距離隨機Kaczmarz方法計算量小于隨機Kaczmarz方法.

3 數(shù)值實驗

通過數(shù)值實驗比較隨機Kaczmarz方法（RK）、貪婪隨機Kaczmarz方法（GRK）、貪婪距離隨機Kaczmarz方法（GDRK）的計算效率.這里，迭代步數(shù)（IT）和計算時間（CPU）取50次計算結(jié)果的中位數(shù).

實驗中，通過MATLAB函數(shù)randn隨機生成解x∈Rn，右端項b∈Rm由A給出，通過MATLAB函數(shù)pinv計算最小范數(shù)解x?=A?b.此外，需要指出GRK方法與GDRK方法按照算法5和6中定義的過程精確執(zhí)行，并沒有顯式地計算矩陣B=AAT來計算殘量.在本節(jié)的所有數(shù)值實驗中，均取初始向量為x0=0，停機準(zhǔn)則為近似解的相對誤差

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

相關(guān)信息見表1～4.表1和表3給出了測試矩陣的階數(shù)、稠密度、秩和條件數(shù)，其中稠密度定義為矩陣非零元個數(shù)除以矩陣元素個數(shù).關(guān)于GDRK方法對于RK方法的加速比，其定義為RK方法所需的計算時間除以GDRK方法所需的計算時間.所有測試矩陣都來自于佛羅里達(dá)稀疏矩陣集，見參考文獻(xiàn)［31］.

例1.表1列出了測試的方陣的相關(guān)信息.包含滿 秩 矩 陣Trefethen_150、cage5、Trefethen_300、Stranke94和秩虧矩陣Sandi_authors.其中，矩陣Trefethen_150、Trefethen_300、Stranke94、Sandi_authors為對稱矩陣，矩陣cage5為非對稱矩陣.

表1 方陣信息Tab.1 Information of square matrices for Example 1

對于測試的方陣，表2給出了RK、GRK、GDRK三種方法的迭代步數(shù)（IT）與計算時間（CPU）.

表2 方陣數(shù)值結(jié)果Tab.2 Numerical results of square matrices for Example 1

從表2可以看出，GDRK與GRK方法總能計算出符合精度的解，而RK方法有時在迭代步數(shù)達(dá)到20萬步后仍不能計算出符合精度的解.而且即便RK方法收斂，GDRK與GRK方法在迭代步數(shù)與計算時間上也均少于RK方法.此外，GDRK方法進一步優(yōu)化了GRK方法，其關(guān)于RK方法迭代時間的加速比最大可以達(dá)到36.80（矩陣cage5），最小可以達(dá)到2.20（矩陣Stranke94）.因此，無論對于滿秩還是秩虧的方陣，GDRK方法均收斂，且其迭代步數(shù)與計算時間均優(yōu)于傳統(tǒng)的RK方法與GRK方法.

圖1和圖2是近似解的相對誤差與迭代步數(shù)的關(guān)系.

圖1 例1近似解的相對誤差隨著迭代步數(shù)的變化曲線Fig.1 Relative solution error with the change of IT for Example 1

圖1描繪了矩陣cage5（左）和Stranke94（右）的近似解的相對誤差隨著迭代步數(shù)變化的曲線，進一步驗證了GDRK方法比經(jīng)典的RK方法收斂更快.

圖2 例2近似解的相對誤差隨著迭代步數(shù)變化的曲線Fig.2 Relative solution error with the change of IT for Example 2

例2.表3列出了測試的長方形矩陣的相關(guān)信息.包含行滿秩矩陣refine、bibd_13_6、Trec8、列滿秩矩陣WorldCities和秩虧矩陣GL7d26.

表3 長方形矩陣信息Tab.3 Information of rectangular matrices for Ex?ample 2

對于測試的長方形矩陣，表4給出了RK、GRK、GDRK三種方法的迭代步數(shù)（IT）與計算時間（CPU）.

從表4可以得出和例1類似結(jié)論，即GDRK與GRK方法總能計算出符合精度的解，而RK方法有時在迭代步數(shù)達(dá)到20萬步后仍不能計算出符合精度的解.而且即使RK方法收斂，GDRK與GRK方法在迭代步數(shù)與計算時間上也均少于RK方法.此外，GDRK方法進一步優(yōu)化了GRK方法，其關(guān)于RK方法迭代時間的加速比最大可以達(dá)到57.23（矩陣refine），最小可以達(dá)到4.12（矩陣bibd_13_6）.因此，無論對于滿秩或秩虧的長方形矩陣，GDRK方法均收斂，且其迭代步數(shù)與計算時間均優(yōu)于傳統(tǒng)的RK方法與GRK方法.

表4 長方形矩陣數(shù)值結(jié)果Tab.4 Numerical results of rectangular matrices for Example 2

圖2 描繪了矩陣WorldCities和Trec8的近似解的相對誤差隨著迭代步數(shù)變化的曲線，進一步驗證了GDRK方法比經(jīng)典的RK方法收斂更快.

4 結(jié)論

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