Lasso回歸算法： 坐標軸下降法與最小角回歸法小結

　　　　前面的文章對線性回歸做了一個小結，文章在這： 線性回歸原理小結。里面對線程回歸的正則化也做了一個初步的介紹。提到了線程回歸的L2正則化-Ridge回歸，以及線程回歸的L1正則化-Lasso回歸。但是對於Lasso回歸的解法沒有提及，本文是對該文的補充和擴展。以下都用矩陣法表示，如果對於矩陣分析不熟悉，推薦學習張賢達的《矩陣分析與應用》。

1. 回顧線性回歸　

　　　　首先我們簡要回歸下線性回歸的一般形式： 

　　　　\(h_\mathbf{\theta}(\mathbf{X}) = \mathbf{X\theta}\)

　　　　需要極小化的損失函數是： 

　　　　\(J(\mathbf\theta) = \frac{1}{2}(\mathbf{X\theta} - \mathbf{Y})^T(\mathbf{X\theta} - \mathbf{Y})\)

　　　　如果用梯度下降法求解，則每一輪\(\theta\)迭代的表達式是： 

　　　　\(\mathbf\theta= \mathbf\theta - \alpha\mathbf{X}^T(\mathbf{X\theta} - \mathbf{Y})\)

　　　　其中\(\alpha\)為步長。

　　　　如果用最小二乘法，則\(\theta\)的結果是：

　　　　\( \mathbf{\theta} = (\mathbf{X^{T}X})^{-1}\mathbf{X^{T}Y} \)

2. 回顧Ridge回歸

　　　　由於直接套用線性回歸可能產生過擬合，我們需要加入正則化項，如果加入的是L2正則化項，就是Ridge回歸，有時也翻譯為脊回歸。它和一般線性回歸的區別是在損失函數上增加了一個L2正則化的項，和一個調節線性回歸項和正則化項權重的系數\(\alpha\)。損失函數表達式如下：

　　　　\(J(\mathbf\theta) = \frac{1}{2}(\mathbf{X\theta} - \mathbf{Y})^T(\mathbf{X\theta} - \mathbf{Y}) + \frac{1}{2}\alpha||\theta||_2^2\)

　　　　其中\(\alpha\)為常數系數，需要進行調優。\(||\theta||_2\)為L2范數。

　　　　Ridge回歸的解法和一般線性回歸大同小異。如果采用梯度下降法，則每一輪\(\theta\)迭代的表達式是：

　　　　\(\mathbf\theta= \mathbf\theta -\beta (\mathbf{X}^T(\mathbf{X\theta} - \mathbf{Y}) + \alpha\theta)\)

　　　　其中\(\beta\)為步長。

　　　　如果用最小二乘法，則\(\theta\)的結果是：

　　　　\(\mathbf{\theta = (X^TX + \alpha E)^{-1}X^TY}\) 

　　　　其中E為單位矩陣。

　　　　Ridge回歸在不拋棄任何一個變量的情況下，縮小了回歸系數，使得模型相對而言比較的穩定，但這會使得模型的變量特別多，模型解釋性差。有沒有折中一點的辦法呢？即又可以防止過擬合，同時克服Ridge回歸

模型變量多的缺點呢？有，這就是下面說的Lasso回歸。

3. 初識Lasso回歸　

　　　　Lasso回歸有時也叫做線性回歸的L1正則化，和Ridge回歸的主要區別就是在正則化項，Ridge回歸用的是L2正則化，而Lasso回歸用的是L1正則化。Lasso回歸的損失函數表達式如下：　

　　　　\(J(\mathbf\theta) = \frac{1}{2n}(\mathbf{X\theta} - \mathbf{Y})^T(\mathbf{X\theta} - \mathbf{Y}) + \alpha||\theta||_1\)

　　　　其中n為樣本個數，\(\alpha\)為常數系數，需要進行調優。\(||\theta||_1\)為L1范數。　　　

　　　　Lasso回歸使得一些系數變小，甚至還是一些絕對值較小的系數直接變為0，因此特別適用於參數數目縮減與參數的選擇，因而用來估計稀疏參數的線性模型。

 

　　　　但是Lasso回歸有一個很大的問題，導致我們需要把它單獨拎出來講，就是它的損失函數不是連續可導的，由於L1范數用的是絕對值之和，導致損失函數有不可導的點。也就是說，我們的最小二乘法，梯度下降法，牛頓法與擬牛頓法對它統統失效了。那我們怎么才能求有這個L1范數的損失函數極小值呢？

 

　　　　OK，本章主角，兩種全新的求極值解法坐標軸下降法（coordinate descent）和最小角回歸法（ Least Angle Regression， LARS）該隆重出場了。　　　　　　　　　　

4. 用坐標軸下降法求解Lasso回歸

　　　　坐標軸下降法顧名思義，是沿着坐標軸的方向去下降，這和梯度下降不同。梯度下降是沿着梯度的負方向下降。不過梯度下降和坐標軸下降的共性就都是迭代法，通過啟發式的方式一步步迭代求解函數的最小值。

　　　　坐標軸下降法的數學依據主要是這個結論（此處不做證明）：一個可微的凸函數\(J(\theta)\), 其中\(\theta\)是nx1的向量，即有n個維度。如果在某一點\(\overline\theta\)，使得\(J(\theta)\)在每一個坐標軸\(\overline\theta_i\)(i = 1,2,...n)上都是最小值，那么\(J(\overline\theta_i)\)就是一個全局的最小值。

　　　　於是我們的優化目標就是在\(\theta\)的n個坐標軸上(或者說向量的方向上)對損失函數做迭代的下降，當所有的坐標軸上的\(\theta_i\)(i = 1,2,...n)都達到收斂時，我們的損失函數最小，此時的\(\theta\)即為我們要求的結果。

　　　　下面我們看看具體的算法過程：

　　　　1. 首先，我們把\(\theta\)向量隨機取一個初值。記為\(\theta^{(0)}\) ，上面的括號里面的數字代表我們迭代的輪數，當前初始輪數為0.

　　　　2. 對於第k輪的迭代。我們從\(\theta_1^{(k)}\)開始，到\(\theta_n^{(k)}\)為止，依次求\(\theta_i^{(k)}\)。\(\theta_i^{(k)}\)的表達式如下：

　　　　\( \theta_i^{(k)}  \in \underbrace{argmin}_{\theta_i} J(\theta_1^{(k)}, \theta_2^{(k)}, ... \theta_{i-1}^{(k)}, \theta_i, \theta_{i+1}^{(k-1)}, ..., \theta_n^{(k-1)})\) 

　　　　也就是說\( \theta_i^{(k)} \)是使\(J(\theta_1^{(k)}, \theta_2^{(k)}, ... \theta_{i-1}^{(k)}, \theta_i, \theta_{i+1}^{(k-1)}, ..., \theta_n^{(k-1)})\)最小化時候的\(\theta_i\)的值。此時\(J(\theta)\)只有\( \theta_i^{(k)} \)是變量，其余均為常量，因此最小值容易通過求導或者一維搜索求得。

　　　　如果上面這個式子不好理解，我們具體一點，在第k輪，\(\theta\)向量的n個維度的迭代式如下：

　　　　\( \theta_1^{(k)}  \in \underbrace{argmin}_{\theta_1} J(\theta_1, \theta_2^{(k-1)}, ... , \theta_n^{(k-1)})\) 

　　　　\( \theta_2^{(k)}  \in \underbrace{argmin}_{\theta_2} J(\theta_1^{(k)}, \theta_2, \theta_3^{(k-1)}... , \theta_n^{(k-1)})\) 

　　　　...

　　　　\( \theta_n^{(k)}  \in \underbrace{argmin}_{\theta_n} J(\theta_1^{(k)}, \theta_2^{(k)}, ... , \theta_{n-1}^{(k)}, \theta_n)\) 

　　　　3. 檢查\(\theta^{(k)}\)向量和\(\theta^{(k-1)}\)向量在各個維度上的變化情況，如果在所有維度上變化都足夠小，那么\(\theta^{(k)}\)即為最終結果，否則轉入2，繼續第k+1輪的迭代。

　　　　以上就是坐標軸下降法的求極值過程，可以和梯度下降做一個比較：

　　　　a) 坐標軸下降法在每次迭代中在當前點處沿一個坐標方向進行一維搜索 ，固定其他的坐標方向，找到一個函數的局部極小值。而梯度下降總是沿着梯度的負方向求函數的局部最小值。

　　　　b) 坐標軸下降優化方法是一種非梯度優化算法。在整個過程中依次循環使用不同的坐標方向進行迭代，一個周期的一維搜索迭代過程相當於一個梯度下降的迭代。

　　　　c) 梯度下降是利用目標函數的導數來確定搜索方向的，該梯度方向可能不與任何坐標軸平行。而坐標軸下降法法是利用當前坐標方向進行搜索，不需要求目標函數的導數，只按照某一坐標方向進行搜索最小值。

　　　　d) 兩者都是迭代方法，且每一輪迭代，都需要O(mn)的計算量(m為樣本數，n為系數向量的維度)

5. 用最小角回歸法求解Lasso回歸

　　　　第四節介紹了坐標軸下降法求解Lasso回歸的方法，此處再介紹另一種常用方法， 最小角回歸法(Least Angle Regression， LARS)。

　　　　在介紹最小角回歸前，我們先看看兩個預備算法，好吧，這個算法真沒有那么好講。

5.1 前向選擇（Forward Selection）算法

　　　　第一個預備算法是前向選擇（Forward Selection）算法。

　　　　前向選擇算法的原理是是一種典型的貪心算法。要解決的問題是對於:

　　　　\(\mathbf{Y = X\theta}\)這樣的線性關系，如何求解系數向量\(\mathbf{\theta}\)的問題。其中\(\mathbf{Y}\)為 mx1的向量，\(\mathbf{X}\)為mxn的矩陣，\(\mathbf{\theta}\)為nx1的向量。m為樣本數量，n為特征維度。

　　　　把 矩陣\(\mathbf{X}\)看做n個mx1的向量\(\mathbf{X_i}\)(i=1,2,...n)，在\(\mathbf{Y}\)的\(\mathbf{X}\)變量\(\mathbf{X_i}\)(i =1,2,...n)中，選擇和目標\(\mathbf{Y}\)最為接近(余弦距離最大)的一個變量\(\mathbf{X_k}\)，用\(\mathbf{X_k}\)來逼近\(\mathbf{Y}\),得到下式：

　　　　\(\overline{\mathbf{Y}} = \mathbf{X_k\theta_k}\)

　　　　其中： \(\mathbf{\theta_k}= \mathbf{\frac{<X_k, Y>}{||X_k||_2}}\)

　　　　即：\(\overline{\mathbf{Y}}\) 是 \(\mathbf{Y}\)在 \(\mathbf{X_k}\)上的投影。那么，可以定義殘差(residual):   \(\mathbf{Y_{yes}} = \mathbf{Y - \overline{Y}}\)。由於是投影，所以很容易知道 \(\mathbf{Y_{yes}} 和\mathbf{X_k}\)是正交的。再以\(\mathbf{Y_{yes}}\)為新的因變量，去掉\(\mathbf{X_k}\)后，剩下的自變量的集合\(\mathbf{X_i}\),i=1,2,3...k−1,k+1,...n}為新的自變量集合，重復剛才投影和殘差的操作，直到殘差為0，或者所有的自變量都用完了，才停止算法。

　　　　當\(\mathbf{X}\)只有2維時，例子如上圖，和\(\mathbf{Y}\)最接近的是\(\mathbf{X_1}\)，首先在\(\mathbf{X_1}\)上面投影，殘差如上圖長虛線。此時\(X_1\theta_1\)模擬了\(\mathbf{Y}\)，\(\theta_1\)模擬了\(\mathbf{ \theta}\)(僅僅模擬了一個維度)。接着發現最接近的是\(\mathbf{X_2}\)，此時用殘差接着在\(\mathbf{X_2}\)投影，殘差如圖中短虛線。由於沒有其他自變量了，此時\(X_1\theta_1+X_2\theta_2\)模擬了\(\mathbf{Y}\),對應的模擬了兩個維度的\(\theta\)即為最終結果，此處\(\theta\)計算設計較多矩陣運算，這里不討論。

　　　　此算法對每個變量只需要執行一次操作，效率高，速度快。但也容易看出，當自變量不是正交的時候，由於每次都是在做投影，所有算法只能給出一個局部近似解。因此，這個簡單的算法太粗糙，還不能直接用於我們的Lasso回歸。

 

5.2 前向梯度（Forward Stagewise）算法

　　　　第二個預備算法是前向梯度（Forward Stagewise）算法。

　　　　前向梯度算法和前向選擇算法有類似的地方，也是在\(\mathbf{Y}\)的\(\mathbf{X}\)變量\(\mathbf{X_i}\)(i =1,2,...n)中，選擇和目標\(\mathbf{Y}\)最為接近(余弦距離最大)的一個變量\(\mathbf{X_k}\)，用\(\mathbf{X_k}\)來逼近\(\mathbf{Y}\)，但是前向梯度算法不是粗暴的用投影，而是每次在最為接近的自變量\(\mathbf{X_t}\)的方向移動一小步，然后再看殘差\(\mathbf{Y_{yes}}\)和哪個\(\mathbf{X_i}\)(i =1,2,...n)最為接近。此時我們也不會把\(\mathbf{X_t}\) 去除，因為我們只是前進了一小步，有可能下面最接近的自變量還是\(\mathbf{X_t}\)。如此進行下去，直到殘差\(\mathbf{Y_{yes}} \)減小到足夠小，算法停止。

　　　　

　　　　當\(\mathbf{X}\)只有2維時，例子如上圖，和\(\mathbf{Y}\)最接近的是\(\mathbf{X_1}\)，首先在\(\mathbf{X_1}\)上面走一小段距離，此處\(\varepsilon\)為一個較小的常量，發現此時的殘差還是和\\(\mathbf{X_1}\)最接近。那么接着沿\(\mathbf{X_1}\)走，一直走到發現殘差不是和\(\mathbf{X_1}\)最接近，而是和\(\mathbf{X_2}\)最接近，此時殘差如上圖長虛線。接着沿着\(\mathbf{X_2}\)走一小步，發現殘差此時又和\(\mathbf{X_1}\)最接近，那么開始沿着\(\mathbf{X_1}\)走，走完一步后發現殘差為0，那么算法停止。此時\(\mathbf{Y}\)由剛才所有的所有步相加而模擬，對應的算出的系數\(\theta\)即為最終結果。此處\(\theta\)計算設計較多矩陣運算，這里不討論。

　　　　當算法在\(\varepsilon\) 很小的時候，可以很精確的給出最優解，當然，其計算的迭代次數也是大大的增加。和前向選擇算法相比，前向梯度算法更加精確，但是更加復雜。

　　　　有沒有折中的辦法可以綜合前向梯度算法和前向選擇算法的優點，做一個折中呢？有！這就是終於要出場的最小角回歸法。

5.3  最小角回歸(Least Angle Regression， LARS)算法

　　　　好吧，最小角回歸(Least Angle Regression， LARS)算法終於出場了。最小角回歸法對前向梯度算法和前向選擇算法做了折中，保留了前向梯度算法一定程度的精確性，同時簡化了前向梯度算法一步步迭代的過程。具體算法是這樣的：　

　　　　首先，還是找到與因變量\(\mathbf{Y}\)最接近或者相關度最高的自變量\(\mathbf{X_k}\)，使用類似於前向梯度算法中的殘差計算方法，得到新的目標\(\mathbf{Y_{yes}}\)，此時不用和前向梯度算法一樣小步小步的走。而是直接向前走直到出現一個\(\mathbf{X_t}\)，使得\(\mathbf{X_t}\)和\(\mathbf{Y_{yes}}\)的相關度和\(\mathbf{X_k}\)與\(\mathbf{Y_{yes}}\)的相關度是一樣的，此時殘差\(\mathbf{Y_{yes}}\)就在\(\mathbf{X_t}\)和\(\mathbf{X_k}\)的角平分線方向上，此時我們開始沿着這個殘差角平分線走，直到出現第三個特征\(\mathbf{X_p}\)和\(\mathbf{Y_{yes}}\)的相關度足夠大的時候，即\(\mathbf{X_p}\)到當前殘差\(\mathbf{Y_{yes}}\)的相關度和\(\theta_t\)，\(\theta_k\)與\(\mathbf{Y_{yes}}\)的一樣。將其也叫入到\(\mathbf{Y}\)的逼近特征集合中，並用\(\mathbf{Y}\)的逼近特征集合的共同角分線，作為新的逼近方向。以此循環，直到\(\mathbf{Y_{yes}}\)足夠的小，或者說所有的變量都已經取完了，算法停止。此時對應的系數\(\theta\)即為最終結果。

 

　　　　當\(\theta\)只有2維時，例子如上圖，和\(\mathbf{Y}\)最接近的是\(\mathbf{X_1}\)，首先在\(\mathbf{X_1}\)上面走一段距離，一直到殘差在\(\mathbf{X_1}\)和\(\mathbf{X_2}\)的角平分線上，此時沿着角平分線走，直到殘差最夠小時停止，此時對應的系數\(\beta\)即為最終結果。此處\(\theta\)計算設計較多矩陣運算，這里不討論。

　　　　最小角回歸法是一個適用於高維數據的回歸算法，其主要的優點有：

　　　　1）特別適合於特征維度n 遠高於樣本數m的情況。

　　　　2）算法的最壞計算復雜度和最小二乘法類似，但是其計算速度幾乎和前向選擇算法一樣

　　　　3）可以產生分段線性結果的完整路徑，這在模型的交叉驗證中極為有用

　　　　主要的缺點是：

　　　　由於LARS的迭代方向是根據目標的殘差而定，所以該算法對樣本的噪聲極為敏感。

6. 總結

　　　　Lasso回歸是在ridge回歸的基礎上發展起來的，如果模型的特征非常多，需要壓縮，那么Lasso回歸是很好的選擇。一般的情況下，普通的線性回歸模型就夠了。

　　　　另外，本文對最小角回歸法怎么求具體的\(\theta\)參數值沒有提及，僅僅涉及了原理，如果對具體的算計推導有興趣，可以參考Bradley Efron的論文《Least Angle Regression》，網上很容易找到。 　　　　

 

（歡迎轉載，轉載請注明出處。歡迎溝通交流： [e-mail]）