王 澤，殷俊鋒

（同濟大學數(shù)學科學學院，上海 200092）

求解線性方程組Ax=b的稀疏解在圖像重構(gòu)、信號處理和機器學習中具有廣泛應用［1-4］。通過引入l1-范數(shù)，求解線性方程組的稀疏解可以轉(zhuǎn)化為求解式（1）正則化最小二乘問題：

求解該問題常見的算法包括基追蹤方法、Bregman方法、共軛梯度方法等［5-7］。

自Strohmer和Vershynin［8］證明了隨機Kaczmarz算法具有線性收斂率以來，許多專家學者對隨機Kaczmarz算法進行了深入的研究［9-11］。Kaczmarz算法求解的是最小二乘解，通常是稠密的，為了求解線性方程組的稀疏解，文獻［12］提出稀疏Kaczmarz算法，該算法在原有的Kaczmarz算法基礎之上引入了軟閾值函數(shù)，解決了最小二乘解稠密性的問題。

在線性方程組相容的情況下，文獻［13］證明了稀疏Kaczmarz算法迭代收斂到問題（1）的唯一解。令右端噪聲項bδ滿足‖b?bδ‖2≤δ時，線性方程組可能不相容，文獻［14-15］證明了隨機Kaczmarz算法的迭代值依期望收斂，并且與無噪聲情況下有相同的收縮因子。文獻［16］在稀疏Kaczmarz算法基礎之上，提出了隨機稀疏Kaczmarz算法，并在有噪聲干擾和無噪聲干擾的情況下給出了隨機稀疏Kaczmarz算法的收斂性分析。當系數(shù)矩陣的行范數(shù)相差較小時，隨機Kaczmarz算法將等概率選取A的行，此時隨機Kaczmarz算法的收斂速率較慢，因此Bai和Wu［17］結(jié)合貪婪和隨機思想，提出一種新的概率準則來選取工作行，形成貪婪隨機Kaczmarz算法。更多貪婪隨機Kaczmarz算法的研究參見文獻［18-21］。

在稀疏隨機Kaczmarz算法的基礎之上［21］，受貪婪算法的啟發(fā)，本文提出稀疏貪婪隨機Kaczmarz算法，給出稀疏貪婪隨機Kaczmarz算法在有噪聲干擾和無噪聲干擾情況下的收斂性分析，并以大量的數(shù)值實驗驗證本文算法的計算效率。

符號標記如下，定義為線性方程組Ax=b的解，其中A∈Rm×n，b∈Rm，a i表示系數(shù)矩陣A的第i行。支持集S=supp()={j∈{1，…，n}≠0}是列向量x中非零元的下標所構(gòu)成的集合，A S是A的子矩陣，由支持集S中的指標選取系數(shù)矩陣A的列構(gòu)成，a iS表示的是A S的第i行，x S是在支持集S限制下的列向量；定義(A)=min{σmin(A J)∣J?{1，…，n}，A J≠0}為子矩陣A J的最小奇異值，其中A J為指標集J中A的列構(gòu)成的子矩 陣，令；⊙為哈達瑪積（Hadamard product）；為對支持集大小的估計值；w j為稀疏隨機Kaczmarz算法第j次迭代的權(quán)重值，w j∈Rm×1。

1 稀疏貪婪隨機Kaczmarz算法

當Ax=b的解x稀疏時，文獻［21］提出了稀疏隨機Kaczmarz算法，加快了算法的收斂速率（算法1）。

算法1稀疏隨機Kaczmarz算法［21］。①輸入A∈Rm×n，b∈Rm，最大迭代數(shù)M和估計的支持集的大小。②輸出x j。③初始化S={1，…，n}，x0=0，j=0。④當j≤M時，置j=j+1。⑤選擇行向量a i，i∈{1，…，n}，每一行對應的概率為。⑥確定估計的支持集⑦生成權(quán)重值⑨轉(zhuǎn)步驟④。

稀疏隨機Kaczmarz算法可以用比隨機Kaczmarz算法更少的迭代步數(shù)找到最小二乘解。由于支持集和稀疏度都是未知的，因此稀疏隨機Kaczmarz算法從支持集中所有元素的稀疏度的初始估計開始，然后在每一次迭代中，稀疏隨機Kaczmarz算法通過去掉向量x中數(shù)量級最小的元素下標來更新估計的支持集。該算法第j次迭代的加權(quán)準則為

其中，j為迭代步數(shù)。當j→∞時，w j⊙a i→a iS，因此原線性方程組退化成

由于κ(A S)≤κ(A)，因此稀疏隨機Kaczmarz算法的收斂因子小于隨機Kaczmarz算法。但是由于最小二乘解總是稠密的，該算法在求解超定線性方程組時，雖然在稀疏解零元對應位置上的元素很小，但仍然不等于零，因此在文獻中提出稀疏Kaczmarz算法，文獻［16］在稀疏Kaczmarz算法的基礎上，提出了隨機稀疏Kaczmarz算法（算法2）。

算法2隨機稀疏Kaczmarz算法［16］。①輸入A∈Rm×n，b∈Rm，最大迭代數(shù)M。②輸出x k。③初始化：x0=x*0=0。④置k=0，當k≤M?1時。⑤選擇行向量a ik，ik∈{1，…，n}，每一行對應的概率為x k+1=Sλ(x*k+1)。⑧轉(zhuǎn)步驟④。

算 法 2中λ＞0，Sλ(x)=max{|x|?λ，0}?sign(x)。結(jié)合稀疏隨機Kaczmarz算法和隨機稀疏Kaczmarz算法的思想，受貪婪算法的啟發(fā)，本文提出稀疏貪婪隨機Kaczmarz算法，在保證算法能夠計算出稀疏解的同時，通過貪婪的概率準則來選擇工作行從而達到加快算法收斂速度的目的。算法3給出稀疏貪婪隨機Kaczmarz算法。

算法3稀疏貪婪隨機Kaczmarz算法。①輸入A∈Rm×n，b∈Rm，最大迭代數(shù)M和估計的支持集的大小k。②輸出x k。③初始化S={1，…，n}，x0=x*0=0。④置k=0時，當k≤M?1時。⑤計算

⑥決定正整數(shù)指標集

⑦計算向量的第i行

?計算

?計算x k+1=Sλ(x*k+1)。?轉(zhuǎn)步驟④。

2 稀疏貪婪隨機Kaczmarz算法的收斂性分析

為了證明求解正則基追蹤問題（1）稀疏貪婪隨機Kaczmarz算法依期望線性收斂，先回顧一些基本的概念和性質(zhì)［22］。

令f：Rn→R是凸函數(shù)，用?f(x)表示f在x∈Rn的次微分

其中，?f(x)是一個非空緊凸集。

定義1凸函數(shù)f：Rn→R，如果存在α＞0，使得對任意的x，y∈Rn且x*∈?f(x)，有

那么f是強凸的，當和α的具體值相關時，則稱f是α-強凸的。

定義2［22］如果f：Rn→R是α-強凸，那么它的共軛函數(shù)f*(x*)：=supx∈Rn x*，x?f(x)可微且有1/α?Lipschitz連續(xù)梯度，例如

有式（4）不等式成立：

定義3令f：Rn→R是強凸函數(shù)，那么x，y∈Rn之間的Bregman距離由f和一個次梯度x*∈?f(x)定義

如果f是可微的，那么有?f(x)={?f(x)}，因為D f(x，y)=f(y)?f(x)?，所以可以簡化來寫Df(x，y)=Dx*f(x，y)。

下面的引理都遵循強凸性的假設［13］，給出了隨機方法收斂性分析所需的Bregman距離的關鍵性質(zhì)。

引理1［13］令f：Rn→R是α強凸函數(shù)，對任意的x，y∈Rn和x*∈?f(x)、y*∈?f(y)有

因此

對于序列x k和xk*∈?f(x k)，(x k，y)的有界性意味著x k和x*k有界。如果f有一個L-Lipschitz連續(xù)梯度，那么也有

定義4令f：Rn→R是強凸函數(shù)，且C?Rn是一個非空閉凸集，則x到C上的Bregman投影是關于f和x*∈?f(x)滿足下式的唯一點對于可微的f，可以簡化為ΠC(x)和distf(x，C)。

引理2［13］令f：Rn→R是強凸函數(shù)，x到C上的Bregman投影點∈C是與f和x*∈?f(x)相關，當且僅當存在一個∈?f()滿足式（5）或式（6）：

可以叫這樣的點是的可允次梯度。

引理3表明，Bregman投影到仿射子空間和半空間可以通過求解可微共軛函數(shù)f*的無約束優(yōu)化問題來計算。

引理3［13］令f：Rn→R是α-強凸函數(shù)，A∈Rm×n，b∈Rm，u∈Rn和β∈R。

（1）x∈Rn到仿射空間L(A，b)：={x∈Rn|Ax=b}≠?的Bregman投影是

其中∈Rm是下式的解：

更重要的是，根據(jù)引理2，：=x*?AT是的可允次梯度。如果A行滿秩，那么對任意的y∈L(A，b)有

（2）x∈Rn到超平面H(u，β)：={x∈ Rn∣u，x=β}滿足u≠0的Bregman投影是

其中∈R是下式的解：

更重要的是：=x*??u是的可允次梯度，對任意的y∈H(u，β)有

如果x不是在半空間H≤(u，β)：={x∈Rn∣u，x≤β}，那么有必要滿足＞0，(x)=且上面的不等式對任意的y∈H≤(u，β)都成立。

引理4［16］對于任意的x∈Rn滿 足?f(x)∩R(AT)≠? 和任意的x*=ATy∈?f(x)∩R(AT)，有

定理1（無噪聲情況）令稀疏貪婪隨機Kaczmarz算法所產(chǎn)生的迭代序列x k依期望線性收斂到唯一解，其收縮因子為

從而

且

證明根據(jù)引理3的（2）中的估計，有式（8）成立：

對式（8）兩邊同時取條件期望有

根據(jù)式（2）和式（3）的定義，有不等式（10）成立：

將式（10）代入到式（9）中，從而有

由式（11）的推導，可以得出

根據(jù)引理4有

根據(jù)引理1，令α=1可得

證畢。

定理2（有噪聲情況）已知右段噪聲項bδ∈Rm滿足‖b?bδ‖2≤δ。如果稀疏貪婪隨機Kaczmarz算法的迭代值x k由bδ替代b來計算，那么有噪聲干擾和無噪聲干擾有相同的收縮因子

證明令

和式（8）相類似有

將式（14）進行簡單的變換等價于式（15）：

而

所以

類似于無噪聲情況的證明方法可以得到

通過歸納得出

由引理1，令α=1，且，有

證畢。

注意到稀疏貪婪隨機Kaczmarz算法的收縮因子為

文獻［16］定理3.2中隨機稀疏Kaczmarz算法的收縮因子為

因

由此可以得出稀疏貪婪隨機Kaczmarz算法的收縮因子小于隨機稀疏Kaczmarz算法的收縮因子。一般而言，收縮因子越小，收斂速度越快，因此理論分析表明稀疏貪婪隨機Kaczmarz算法收斂速度會比原有的隨機稀疏Kaczmarz算法的收斂速度更快。

3 數(shù)值實驗

通過隨機生成的高斯矩陣和矩陣市場中的矩陣2個數(shù)值實驗算例來比較稀疏貪婪隨機Kaczmarz（SGRK）算法、隨機稀疏Kaczmarz（RaSK）算法和稀疏隨機Kaczmarz（SRK）算法的收斂速度。為了能夠使三者收斂到相同解上進行比較，在本次數(shù)值實驗中，在稀疏隨機Kaczmarz算法最后一步加上了軟閾值函數(shù)Sλ(x)參與迭代以保證求得稀疏解。對相同參數(shù)的同一個線性方程組Ax=b做了20次試驗，得到迭代步數(shù)（IT）和CPU計算時間的平均值。解的相對誤差（RSE）定義為

例1利用Matlab軟件中的randn函數(shù)隨機生成系數(shù)矩陣A和x n×1，其中有k×n個非零元，之后利用b=A×x得到相容線性方程組，初始向量x0=0。取k=2×k，通過k來設置解的稀疏度，其中k?是對稀疏度的初始估計值。

當解稀疏度為0.2、0.4、0.6以及0.8時，表1和表2分別給出了系數(shù)矩陣A的維數(shù)為1000×150和4000×300時的SGRK算法、RaSK算法和SRK算法的迭代步數(shù)和計算時間。

當稀疏度為0.2和0.8時，圖1和圖2分別給出系數(shù)矩陣A的維數(shù)為1 000×150和4 000×300時3種算法近似解的相對誤差隨迭代步數(shù)變化的曲線。

通過表1和表2中的實驗數(shù)據(jù)以及圖1和圖2中的變化曲線可以看出，隨著稀疏解非零元個數(shù)的不斷增加，SRK算法的收斂曲線不斷向RaSK算法的收斂曲線靠攏，這是因為SRK算法在對應零元位置上乘以權(quán)重，加速減小零元位置上的元素，因此在零元較多的情況下，SRK算法的收斂速度更快。但是，由于SRK算法的計算過程較為復雜，因此在計算時間上一開始低于RaSK算法，后來在迭代步數(shù)不明顯占優(yōu)的情況下，計算時間要多于RaSK算法。SGRK算法在2種算法的基礎之上，每次以貪婪的概率準則隨機選取行指標進行迭代，因此SGRK算法的計算時間和迭代步數(shù)都要優(yōu)于其他2種算法。

圖1 當矩陣A的維數(shù)為1 000×150時SGRK、RaSK和SRK算法的收斂曲線Fig.1 Convergence curves of SGRK,RaSK and SRK methods for Example 1 at a matrix A of 1 000×150

圖2 當矩陣A的維數(shù)為4 000×300時SGRK、RaSK和SRK算法的收斂曲線Fig.2 Convergence curves of SGRK,RaSK and SRK methods for Example 1 at a matrix A of 4 000×300

表1 m=1 000、n=150時SGRK、RaSK和SRK的迭代步數(shù)和計算時間Tab.1 Iteration steps and computation time of SGRK,RaSK and SRK（m=1 000、n=150）

表2 m=4 000、n=300時SGRK、RaSK和SRK的迭代步數(shù)和計算時間Tab.2 Iteration steps and computation time of SGRK,RaSK and SRK（m=4 000、n=300）

例2選取了矩陣市場中的一些矩陣進行數(shù)值實驗。為了得到超定的系數(shù)矩陣，將所有的欠定矩陣取轉(zhuǎn)置之后再做數(shù)值實驗，取稀疏度估計為k=2k，稀疏度k=0.2n。表3列出的是測試矩陣的相關信息。

表3 矩陣信息Tab.3 Matrix information for Example 2

對于測試的矩陣，表4給出了SGRK算法、RaSK算法和SRK算法的迭代步數(shù)和運行時間。

圖3對應于表4，給出了3種算法近似解的相對誤差隨著迭代步數(shù)的變化曲線。通過對比可以發(fā)現(xiàn)，在無噪聲干擾的情況下，SGRK算法的迭代步數(shù)和CPU運行時間都優(yōu)于原來的RaSK算法和SRK算法。

圖3 SGRK、RaSK和SRK算法在稀疏度k=0.2 n時的收斂曲線Fig.3 Convergence curves of SGRK,RaSK and SRK methods convergence curve for Examples 2 at a sparsity k=0.2 n

表4 SGRK、RaSK和SRK的迭代步數(shù)和計算時間Tab.4 Iteration steps and computation time of SGRK,RaSK and SRK

圖4和圖5描述SGRK算法、SRK算法和RaSK算法在有噪聲干擾的系統(tǒng)中的一些數(shù)值實驗結(jié)果。用Matlab軟件中的randn隨機生成高斯矩陣和范數(shù)為0.02的獨立高斯噪聲進行數(shù)值實驗。

圖4和圖5中水平的實線是SGRK算法的相對誤差閾值，水平的虛線是RaSK算法的相對誤差閾值，圖4是在有噪聲干擾的情況下，解的稀疏度為0.2時3種算法的近似解的相對誤差隨迭代步數(shù)的變化曲線。其中圖4a描述系數(shù)矩陣A的維數(shù)為1 000×150時的收斂情況，圖4b描述的是系數(shù)矩陣A的維數(shù)為4 000×300時的收斂情況。

圖5是在有噪聲干擾的情況下，解的稀疏度為0.6時3種算法的近似解的相對誤差隨迭代步數(shù)的變化曲線。其中圖5a和圖5b分別畫出的是系數(shù)矩陣A的維數(shù)為1 000×150和4 000×300時3種算法的收斂曲線。

由圖4和圖5可以看出，在有噪聲干擾的情況下，SGRK算法優(yōu)于其他2種算法最先達到穩(wěn)定閾值。另外，在圖4和圖5中分別用水平的實線和虛線畫出了SGRK算法的相對誤差閾值和RaSK算法的相對誤差閾值，通過數(shù)值實驗表明，由定理2中SGRK算法推導出來的閾值更接近實際情況。

圖4 SGRK、RaSK和SRK算法在稀疏度k=0.2 n時誤差閾值和實際誤差比較Fig.4 Comparison of error threshold and actual error of SGRK,RaSK and SRK methods at a sparsity k=0.2 n

圖5 SGRK、RaSK和SRK算法在稀疏度k=0.6 n時誤差閾值和實際誤差比較Fig.5 Comparison of error threshold and actual error of SGRK,RaSK and SRK methods at a sparsity k=0.6 n

4 結(jié)論

在隨機稀疏Kaczmarz算法和稀疏隨機Kaczmarz算法的基礎上，受貪婪算法的啟發(fā)，提出稀疏貪婪隨機Kaczmarz算法，給出了新算法的收斂性分析，且在有噪聲干擾和無噪聲干擾的情況下，通過理論證明了新算法的收縮因子小于隨機稀疏Kaczmarz算法的收縮因子。數(shù)值實驗表明所提出的新算法在迭代步數(shù)和計算時間上均優(yōu)于傳統(tǒng)的隨機稀疏Kaczmarz算法。

作者貢獻聲明：

王 澤：算法設計者和算法研究的執(zhí)行人，構(gòu)造新的算法，給出收斂性證明，完成數(shù)值實驗和數(shù)據(jù)分析、論文初稿的寫作。

殷俊鋒：研究的構(gòu)思者及負責人，指導實驗設計，數(shù)據(jù)分析，論文寫作與修改。