基于隨機投影和OMP的大規模稀疏線性回歸算法*

祁德濤，楊 健，鄒仕祥，郝 琳

（1.陸軍工程大學，江蘇 南京 210007；2.96037部隊55分隊，江蘇 南京 210007）

0 引言

隨著大數據時代的蓬勃發展，互聯網技術突飛猛進，網絡帶寬、數據存儲等技術等的提升讓各行各業的數據都呈指數式增長，大數據分析與利用逐漸成為各領域角逐的焦點。在當前市場環境下，各行各業對其掌握的數據進行預測具有極高實用價值，大規模稀疏線性回歸分析作為數據預測最基本手段也是當前學者們研究的熱點。目前稀疏線性回歸算法大多是基于梯度下降的算法，其在大規模問題求解上具有局限性。本文研究大規模稀疏線性問題高效求解方案。假設給定數據矩陣及其對應的觀測值,其中元素xi∈?d表示數據矩陣第i個樣本的特征，yi是其對應的觀測值。稀疏線性回歸要解決的問題是學習稀疏線性函數f(x)=XTβ使得每一個yi都能夠被f(x)很好的近似，其問題的核心是求解由模型系數稀疏向量。本文研究數據特征維度大于樣本個數的情形，即d＞n。

對于稀疏學習領域，解決稀疏向量β求解問題可建模為經典的LASSO問題：

其中，λ為稀疏模型的正則化系數，||●||是向量l2范數。

近年來，很多求解式（1）優化問題方法嘗試在非光滑凸優化算法上做出改進，提出了很多高效的優化算法。文獻[1]中，作者提出適用于求解非光滑凸優化問題的修正鄰近梯度法，算法的特點是采用自適應步長參數,并且該算法的線性收斂性不需要目標函數的強凸性作為前提。文獻[2]中，為了解決更為復雜的可分塊的凸優化lasso問題，采用分快思想，將目標函數轉化為多個凸函數之和，借助鄰近算子，利用交替鄰近梯度法有效的解決了這一類問題。文獻[3]提出了基于同倫策略的近端梯度下降算法，利用迭代策略，將問題轉化為求解一系列形如式（1）的子問題。以上算法再解決過程中雖然都能達到一定的收斂速度，但由于其迭代都是以當前梯度為基準，每次迭代過程中對梯度的求解需要計算形如XXTβt的矩陣相乘，其計算量為O(nd)，當數據樣本與特征數都較大時，迭代運算復雜度過高。

本文以降維為出發點，借助了隨機投影算法在大數據應用中的高效降維優勢。由于隨機矩陣算法在大數據應用中的高效性，現在的隨機矩陣算法很多均被應用到了最小二乘[4]和低秩矩陣近似[5]等問題的快速求解上。最近有研究結果表明，相比傳統的求解算法（比如，Cholesky分解，QR分解），基于隨機矩陣的最小二乘問題近似求解以及矩陣的低秩近似算法的算法復雜度取得明顯的降低。

具體的，近幾年學者們提出了很多隨機矩陣算法[6]用以加速大規模數據數據分析中的相關問題。他們通常通過構建數據矩陣X的低秩近似矩陣縮小原問題的規模，然后通過求解小規模問題獲得原問題近似解。目前的隨機矩陣算法根據構建隨機矩陣的方式，主要分為兩大類：基于隨機采樣的方法和基于隨機投影的方法?；陔S機采樣的方法[7]中，構建的隨機矩陣的每一列都是根據特定的概率從原始數據矩陣抽樣得來。而基于隨機投影的方法[8]中，隨機矩陣的列元素都是原始數據矩陣列元素的線性組合。由于隨機矩陣算法相對是一個較新的研究內容，在大規模稀疏學習模型訓練算法中，暫時還沒較多研究成果，最近的研究成果中文獻[9]則利用隨機投影的辦法，結合基于同倫策略的近端梯度下降算法設計了一個復合的稀疏線性回歸模型訓練算法。利用隨機投影首先降低了數據矩陣的規模，在一定程度上提高了算法整體效率。

綜上所述，本文設計了一種基于隨機投影的加速稀疏線性回歸問題求解方法。算法分為兩步：首先利用隨機投影構建數據矩陣X的低秩近似矩陣其中從而使得梯度計算涉及的矩陣乘法XTXβt近似為WTQTQWβt。通過此方式，計算復雜度由O(nd)簡化為O(dk+nk)，當k＜＜min(n,d)時，效率提升非常明顯；第二步，通過正交匹配追蹤算法，利用最小二乘估計方法迭代地逼近已知支撐集下的稀疏向量。對于K-稀疏向量，OMP算法只需要K次迭代就能找到稀疏向量所有稀疏位置，最優地恢復稀疏向量。

通過上述的分析，本文的貢獻主要以下兩點：（1）提出利用經典的OMP算法求解稀疏線性回歸問題，避免了基于梯度下降之類算法迭代計算梯度高復雜度的問題；（2）結合隨機投影與OMP算法，提出面向稀疏線性回歸問題的基于隨機投影的OMP算法，大幅度提高了大規模稀疏線性回歸問題求解效率。

1 基于隨機投影的OMP算法

基于當前大規模數據稀疏線性回歸算法存在的弊端，本文結合隨機矩陣算法在大規模數據分析的優勢以及OMP算法的優勢，提出基于隨機投影的OMP算法。首先利用隨機矩陣算法對大規模數據進行k-低秩近似降維，達到減小規模的目的?？紤]到本文待解決問題是欠定的，如果采用隨機采樣的方式，會導致問題更加欠定，因此本文利用隨機投影的辦法，算法第二步將k-低秩近似降維的數據矩陣和對應的觀測值作為輸入，利用OMP算法求解稀疏向量β。據目前所知，這是第一次將隨機投影算法與OMP算法結合。算法流程如下：

基于隨機投影的正交匹配追蹤算法

1：輸入：數據矩陣x、響應向量y、參數λ0和k

3：采樣d×k隨機矩陣Z,其中Zij服從正態分布N(0,1)

4:計算矩陣XTZ的QR分解，即XTZ=QR

5：計算x的低秩近似矩陣=(XQ)QT

6:計算稀疏向量β：

7：初始化β0=0

8：fort=0 tok-1 do

9:初始化殘差r0=y

11：end for

12:輸出β

2 算法理論分析

2.1 矩陣低秩近似的隨機投影算法

考慮到本文與最小二乘問題和低秩矩陣近似的隨機投影算法關聯較大，本節首先介紹隨機投影算法在經典最小二乘和低秩矩陣近似問題的原理。

引理2.1 對任意ε，δ∈(0,1)和正整數n，令d為正整數且滿足：

定義A為上述隨機投影，對于Γ上含有I個點的任意集合I，至少以1-δ的概率，對所有的，都有下述不等式成立，

這就是Johnson-Lindenstrauss引理，引理表明，高維數據能夠通過滿足一定性質的Lipschitz映射，將其映射到低維空間中，并且能較高概率近似地保持原始數據的幾何形狀。

針對最小二乘問題中隨機投影的應用：

2.2 正交匹配追蹤算法

本文第一部分利用隨機投影預處理達到了降維的目的，然后利用OMP算法解決稀疏線性回歸問題。之前的學者已經證明了OMP算法在壓縮感知重構中良好的表現，OMP算法能夠快速有效的在有限的步數內最優的恢復感知矩陣，那么基于壓縮感知問題和稀疏線性回歸的相似性，算法是否也能快速有效的在有限的步數內恢復出稀疏向量。本文的理論推理和實驗結果最終證明其是可行的。

OMP算法是一種貪婪算法，貪婪算法是基于l0優化的啟發式恢復算法。相比于凸優化算法，它具有簡單、快速、高效的特點，其核心思想是貪婪的選取當前匹配測試中最匹配的原子，并在追蹤下一次測試前去掉已選擇原子帶來的影響，直到算法收斂。它主要包括四個關鍵步驟：匹配測試、貪婪選取、追蹤準則和停止準則。匹配追蹤（MP）是最簡單的貪婪算法，它每次迭代選擇最匹配的原子，其追蹤準則只是簡單地利用相干投影系數不斷地逼近稀疏信號，因而收斂速度慢、重構精度差。OMP算法則在追蹤時用最小二乘估計的方法逼近已知支撐集下的稀疏信號，而最小二乘是當前支撐集下的正交的最優化估計，從而避免了在下次匹配測試中重復地選擇已選擇過的原子。對于K-稀疏信號，OMP算法只需要K次迭代就能找到所有的稀疏位置，并最優地恢復出信號。

假設對于K-稀疏向量β，其對應的稀疏支撐集為S={S1,…,SK}∈{1,…,N}和對應的數據矩陣X，由稀疏恢復公式得：

匹配測試是算法的核心，它根據殘差與數據矩陣列原子的匹配程度來貪婪的選擇稀疏原子，對于第t次選取，匹配測試為：

匹配測試能夠準確的選擇正確的稀疏位置的保證為|Tt∈S|＞＞|Tt?S|。從上式可知，列原子的匹配程度與信號的稀疏度、噪聲成分和數據矩陣列原子間的非相干程度有關。如果數據矩陣是N×N正交矩陣，那么原子間相干系數為：

因此在一定的信噪比下，|Tt∈S|＞＞|Tt?S|必然成立。其中：

對于稀疏信號的數據矩陣來說。M×N的數據矩陣是冗余的，其列原子并不是正交的，而且具有較小相干系數。為保證匹配測試能夠選擇正確的稀疏位置，數據矩陣的列原子必須是近似正交的，在感知矩陣良好的非相干性保證下，匹配測試的成功概率隨著信號稀疏度和信號噪聲能量的增加而不斷衰減。因而，對于同一數據矩陣，貪婪算法的恢復性能隨著稀疏度的增加和信號信噪比的減小而減少。

本小結最后，對算法進行簡要的總結：（1）算法更新過程中嚴格控制追蹤準則避免出現當前最匹配原子地重復選取，對于K-稀疏信號，算法只需要迭代K次就能找到完整的稀疏位置，不再需要步長參數。（2）與PGH算法以及目前流行地基于梯度下降一類地算法相比，本算法更新過程中避免了迭代過程中的梯度計算，同時充分利用了隨機投影在大規模數據矩陣的高效性，算法整體效率非常高。

3 實 驗

上一節論述了算法的理論可行性。接下來通過實驗驗證算法的準確性。本文的實驗結果將與文獻[11]中PGH算法做對比。

PGH算法：當前最優，具有指數收斂速度

實驗主要關注如下三點：（1）算法的收斂速度；（2）待解決模型稀疏度與算法迭代次數的的關系；（3）一定噪聲水平下算法的收斂效果。

3.1 合成數據集上的實驗

數據集：本實驗采用合成數據集，具體產生方式如下。數據矩陣X：首先構造一個服從[-1，1]均勻分布的隨機矩陣，然后對其進行QR分解。即U=QR，其中，。隨后選取Q的前k項得到新矩陣通過該計算得到低秩矩陣最后在低秩矩陣上加上高斯噪聲即得數據矩陣X=X0+eX。

稀疏向量β1：根據數據矩陣的特征數d，設置稀疏向量并隨即從其中選取s個元素設為1，其余均為零。

響應值y：響應值y通過如下公式計算獲得y=XTβ1+ey。其中噪聲ey服從高斯隨機分布

參數設置：實驗中，具體的參數設置如下。對于訓練數據集，d=10000，N=5000，k=500，s=25取噪聲ex,ey的標準差為1。在PGH算法中，設置λopt=0.001，Lmin=10-9，δ=0.2。在OMP算法中，Lmin=0.001，γ=10-5。

實驗目標：通過數據矩陣X和響應值y中恢復稀疏向量β。

算法評估標準：實驗需要重點評估學習出來的稀疏向量β的準確性，所以采用相對于最優解的誤差||β1-β||2，殘差以及運行時間這三個衡量指標。每個實驗重復100次，并將均值作為最終結果。

實驗結果：如圖1、2，展示了噪聲水平設置為1時，隨著迭代次數的增加，稀疏向量誤差、響應值殘差的變化情況。圖1、2可以看出OMP算法與PGH算法相比，β以越來越快的速度提升質量，殘差也相應的以越來越快的速度減小并在較小迭代次數內收斂。圖3、4展示了運行時間與稀疏向量誤差、觀測值殘差的關系。明顯地發現OMP算法在整體時間效率上有絕對的優勢，以更短的時間收斂，同時結合圖1、2也說明OMP算法單次迭代效率更高，計算更加簡單。

圖1 迭代過程中稀疏向量相對誤差變化

圖2 迭代過程中殘差變化

圖3 算法運行時間與稀疏向量相對誤差關系

圖4 算法運行時間與殘差關系

為進一步說明對于K-稀疏向量OMP算法只需要進行K次迭代，實驗設置了不同稀疏度的稀疏向量。如下圖5展示，證明了OMP算法的迭代次數與稀疏度緊密相關，因此OMP算法求解小稀疏度模型會有更大的速度優勢。

最后，為了觀察不同噪聲水平對兩種算法運行結果的影響，實驗設計了不同水平的噪聲進行比較，相關結果以列表形式展示如表1所示。

實驗結果表明對于不同噪聲水平，OMP算法始終保持它的效率優勢，且其殘差和誤差均比PGH低，說明其具備一定抗噪能力，收斂效果有保證。

圖5 不同稀疏度迭代次數情況

表1 當d=10000,N=5000,k=500時的結果

3.2 實際數據集上的實驗

數據集：本文采用Minst數據集進行實驗。MNIST數據集是由手寫體數字的掃描圖片組成的，數據集包含訓練樣本60 000個和測試樣本10 000個。每個樣本是一張28×28的圖片，其維度是784。本實驗，隨機地從訓練數據集中采集10 000個樣本構成數據，然后從測試集中隨機地抽取一個樣本作為響應值

實驗目標：學習出一個稀疏向量β使得y能夠被XTβ較好的近似。

參數設置：實驗中，k=100，其余參數設置如下。在PGH算法中，λopt=0.001，Lmin=10-9，δ=0.2，Tmax=500。在OMP算法中，Lmin=0.001，γ=10-5。每次實驗均重復100次，取均值作為實驗結果。

算法評估指標：對于此實驗，由于β是未知的，所以采用響應值殘差和時間效率作為評估標準，并將實驗結果以表格形式展示如表2所示。

表2 當T=500時的結果

實驗結論：通過對實驗結果的分析，發現OMP算法在時間效率上有絕對的優勢。同時，求得的殘差與PGH算法相近，說明OMP算法在恢復效果上也是可靠的。

4 結語

本文提出了面向稀疏線性回歸的基于隨機投影的正交匹配追蹤算法。算法首先利用高斯隨機投影獲得k-低秩近似矩陣，達到降維的目的，降低計算量。然后，提出將正交匹配追蹤算法用于求解稀疏線性回歸問題。通過理論分析和實驗結果對比證明了該算法恢復準確，同時與基于梯度下降一類的算法相比有非常好的時間優勢，特別是針對小稀疏度模型問題，算法在恢復效果和恢復效率上都有顯著優勢。