機器學習演算法實踐-嶺回歸和LASSO

02-03

前言

繼續線性回歸的總結, 本文主要介紹兩種線性回歸的縮減(shrinkage)方法的基礎知識: 嶺回歸(Ridge Regression)和LASSO(Least Absolute Shrinkage and Selection Operator)並對其進行了Python實現。同時也對一種更為簡單的向前逐步回歸計算回歸係數的方法進行了相應的實現。

正文

通過上一篇《機器學習演算法實踐-標準與局部加權線性回歸》中標準線性回歸的公式 $w = (X^{T}X)^{-1}X^{T}y$ 中可以看出在計算回歸係數的時候我們需要計算矩陣 $X^TX$ 的逆，但是如果該矩陣是個奇異矩陣，則無法對其進行求解。那麼什麼情況下該矩陣會有奇異性呢?

X本身存在線性相關關係(多重共線性), 即非滿秩矩陣。如果數據的特徵中存在兩個相關的變數，即使並不是完全線性相關，但是也會造成矩陣求逆的時候造成求解不穩定。
當數據特徵比數據量還要多的時候, 即 $m < n$ , 這時候矩陣 $X$ 是一個矮胖型的矩陣，非滿秩。

對於上面的兩種情況，我們需要對最初的標準線性回歸做一定的變化使原先無法求逆的矩陣變得非奇異，使得問題可以穩定求解。我們可以通過縮減的方式來處理這些問題例如嶺回歸和LASSO.

中心化和標準化

這裡先介紹下數據的中心化和標準化，在回歸問題和一些機器學習演算法中通常要對原始數據進行中心化和標準化處理，也就是需要將數據的均值調整到0，標準差調整為1, 計算過程很簡單就是將所有數據減去平均值後再除以標準差:

$x_i^{『} = frac{x_i - mu}{sigma}$

這樣調整後的均值:

$mu^{『} = frac{(sum_{i=1}^{n}x_i)/n - mu}{sigma} = 0$

調整後的標準差:

$sigma^{『} = frac{(x_i - mu)^2/n}{sigma^2} = frac{sigma^2}{sigma^2} = 1$

之所以需要進行中心化其實就是個平移過程，將所有數據的中心平移到原點。而標準化則是使得所有數據的不同特徵都有相同的尺度Scale, 這樣在使用梯度下降法以及其他方法優化的時候不同特徵參數的影響程度就會一致了。

如下圖所示，可以看出得到的標準化數據在每個維度上的尺度是一致的(圖片來自網路，侵刪)

嶺回歸(Ridge Regression)

標準最小二乘法優化問題:

$f(w) = sum_{i=1}^{m} (y_i - x_{i}^{T}w)^2$

也可以通過矩陣表示:

$f(w) = (y - Xw)^{T}(y - Xw)$

得到的回歸係數為:

$hat{w} = (X^{T}X)^{-1}X^{T}y$

這個問題解存在且唯一的條件就是XX列滿秩: $rank(X) = dim(X)$ .

即使 $X$ 列滿秩，但是當數據特徵中存在共線性，即相關性比較大的時候，會使得標準最小二乘求解不穩定, $X^TX$ 的行列式接近零，計算 $X^TX$ 的時候誤差會很大。這個時候我們需要在cost function上添加一個懲罰項 $lambdasum_{i=1}^{n}w_{i}^2$ ，稱為L2正則化。

這個時候的cost function的形式就為:

$f(w) = sum_{i=1}^{m} (y_i - x_{i}^{T}w)^2 + lambdasum_{i=1}^{n}w_{i}^{2}$

通過加入此懲罰項進行優化後，限制了回歸係數wiwi的絕對值，數學上可以證明上式的等價形式如下:

$f(w) = sum_{i=1}^{m} (y_i - x_{i}^{T}w)^2$

$s.t. sum_{i=1}^{n}w_{j}^2 le t$

其中t為某個閾值。

將嶺回歸係數用矩陣的形式表示:

$hat{w} = (X^{T}X + lambda I)^{-1}X^{T}y$

可以看到，就是通過將 $X^TX$ 加上一個單位矩陣是的矩陣變成非奇異矩陣並可以進行求逆運算。

嶺回歸的幾何意義

以兩個變數為例, 殘差平方和可以表示為 $w_1, w_2$ 的一個二次函數，是一個在三維空間中的拋物面，可以用等值線來表示。而限制條件 $w_1^2 + w_2^2 < t$ ，相當於在二維平面的一個圓。這個時候等值線與圓相切的點便是在約束條件下的最優點，如下圖所示，

嶺回歸的一些性質

當嶺參數 $lambda = 0$ 時，得到的解是最小二乘解
當嶺參數 $lambda$ 趨向更大時，嶺回歸係數 $w_i$ 趨向於0，約束項 $t$ 很小

嶺回歸的Python實現

通過矩陣的形式計算 $hat{w}$ , 可以很簡單的實現

def ridge_regression(X, y, lambd=0.2):n 獲取嶺回歸係數n n XTX = X.T*Xn m, _ = XTX.shapen I = np.matrix(np.eye(m))n w = (XTX + lambd*I).I*X.T*yn return wn

嶺跡圖

可以知道求得的嶺係數 $w_i$ 是嶺參數 $lambda$ 的函數，不同的 $lambda$ 得到不同的嶺參數 $w_i$ , 因此我們可以增大 $lambda$ 的值來得到嶺回歸係數的變化，以及嶺參數的變化軌跡圖(嶺跡圖), 不存在奇異性時，嶺跡圖應穩定的逐漸趨向於0。

通過嶺跡圖我們可以:

觀察較佳的 $lambda$ 取值
觀察變數是否有多重共線性

繪製嶺跡圖

上面我們通過函數ridge_regression實現了計算嶺回歸係數的計算，我們使用《機器學習實戰》中的鮑魚年齡的數據來進行計算並繪製不同 $lambda$ 的嶺參數變化的軌跡圖。數據以及完整代碼詳見 https://github.com/PytLab/MLBox/tree/master/linear_regression

選取30組不同的λλ來獲取嶺係數矩陣包含30個不同的嶺係數。

def ridge_traj(X, y, ntest=30):n 獲取嶺軌跡矩陣n n _, n = X.shapen ws = np.zeros((ntest, n))n for i in range(ntest):n w = ridge_regression(X, y, lambd=exp(i-10))n ws[i, :] = w.Tn return wsn

繪製嶺軌跡圖

if __main__ == __name__:n ntest = 30n # 載入數據n X, y = load_data(abalone.txt)n # 中心化 & 標準化n X, y = standarize(X), standarize(y)n # 繪製嶺軌跡n ws = ridge_traj(X, y, ntest)n fig = plt.figure()n ax = fig.add_subplot(111)n lambdas = [i-10 for i in range(ntest)]n ax.plot(lambdas, ws)n plt.show()n

上圖繪製了回歸係數 $w_i$ 與 $log(lambda)$ 的關係，在最左邊 $lambda$ 係數最小時，可以得到所有係數的原始值(與標準線性回歸相同); 而在右邊，係數全部縮減為0, 從不穩定趨於穩定；為了定量的找到最佳參數值，還需要進行交叉驗證。要判斷哪些變數對結果的預測最具影響力，可以觀察他們的係數大小即可。

LASSO

嶺回歸限定了所有回歸係數的平方和不大於 $t$ , 在使用普通最小二乘法回歸的時候當兩個變數具有相關性的時候，可能會使得其中一個係數是個很大正數，另一個係數是很大的負數。通過嶺回歸的 $sum_{i=1}^{n} w_i le t$ 的限制，可以避免這個問題。

LASSO(The Least Absolute Shrinkage and Selection Operator)是另一種縮減方法，將回歸係數收縮在一定的區域內。LASSO的主要思想是構造一個一階懲罰函數獲得一個精鍊的模型, 通過最終確定一些變數的係數為0進行特徵篩選。

LASSO的懲罰項為:

$sum_{i=1}^{n} vert w_i vert le t$

與嶺回歸的不同在於，此約束條件使用了絕對值的一階懲罰函數代替了平方和的二階函數。雖然只是形式稍有不同，但是得到的結果卻又很大差別。在LASSO中，當λλ很小的時候，一些係數會隨著變為0而嶺回歸卻很難使得某個係數恰好縮減為0. 我們可以通過幾何解釋看到LASSO與嶺回歸之間的不同。

LASSO的幾何解釋

同樣以兩個變數為例，標準線性回歸的cost function還是可以用二維平面的等值線表示，而約束條件則與嶺回歸的圓不同，LASSO的約束條件可以用方形表示，如下圖:

相比圓，方形的頂點更容易與拋物面相交，頂點就意味著對應的很多係數為0，而嶺回歸中的圓上的任意一點都很容易與拋物面相交很難得到正好等於0的係數。這也就意味著，lasso起到了很好的篩選變數的作用。

LASSO回歸係數的計算

雖然懲罰函數只是做了細微的變化，但是相比嶺回歸可以直接通過矩陣運算得到回歸係數相比，LASSO的計算變得相對複雜。由於懲罰項中含有絕對值，此函數的導數是連續不光滑的，所以無法進行求導並使用梯度下降優化。本部分使用坐標下降發對LASSO回歸係數進行計算。

坐標下降法是每次選擇一個維度的參數進行一維優化，然後不斷的迭代對多個維度進行更新直到函數收斂。SVM對偶問題的優化演算法SMO也是類似的原理，這部分的詳細介紹我在之前的一篇博客中進行了整理，參考《機器學習演算法實踐-SVM中的SMO演算法》。

下面我們分別對LASSO的cost function的兩部分求解:

1) RSS部分

$RSS(w) = sum_{i=1}^{m}(y_i - sum_{j=1}^{n}x_{ij}w_j)^2$

求導:

$frac{partial RSS(w)}{partial w_k} = -2sum_{i=1}^{m}x_{ik}(y_i - sum_{j=1}^{n}x_{ij}w_j)$

$= -2sum_{i=1}^{m}(x_{ik}y_i - x_{ik}sum_{j=1, j ne k}^{n}x_{ij}w_j - x_{ik}^2w_k)$

$= -2sum_{i=1}^{m}x_{ik}(y_i - sum_{j=1, j ne k}^{n}x_{ij}w_{j}) + 2w_ksum_{i=1}^{m}x_{ik}^2$

令 $p_k = sum_{i=1}^{m}x_{ik}(y_i - sum_{j=1, j ne k}^{n}x_{ij}w_{j})$ ， $z_k = sum_{i=1}{m}x_{ik}^2$ 得到:

$frac{partial RSS(w)}{partial w_j} = -2p_k + 2z_kw_k$ $frac{partial RSS(w)}{partial w_j} = -2p_k + 2z_kw_k$

2）正則項

關於懲罰項的求導我們需要使用subgradient，可以參考LASSO（least absolute shrinkage and selection operator）回歸中如何用梯度下降法求解？

$lambda frac{partial sum_{i=1}^{n}vert w_j vert}{partial w_k} = begin{cases} -lambda & w_k < 0 [-lambda, lambda] & w_k = 0 lambda & w_k > 0 end{cases}$

這樣整體的偏導數:

$frac{partial f(w)}{partial w_k} = 2z_kw_k - 2p_k + begin{cases} -lambda & w_k < 0 [-lambda, lambda] & w_k = 0 lambda & w_k > 0 end{cases}$

$= begin{cases} 2z_kw_k - 2p_k - lambda & w_k < 0 [-2p_k - lambda, -2p_k + lambda] & w_j = 0 2z_kw_k - 2p_k + lambda & w_k > 0 end{cases}$

令 $frac{partial f(w)}{partial w_k} = 0$ 得到

$hat{w_k} = begin{cases} (p_k + lambda/2)/z_k & p_k < -lambda/2 0 & -lambda/2 le p_k le lambda/2 (p_k - lambda/2)/z_k & p_k > lambda/2 end{cases}$

通過上面的公式我們便可以每次選取一維進行優化並不斷跌打得到最優回歸係數。

LASSO的Python實現

根據上面代碼我們實現梯度下降法並使用其獲取LASSO回歸係數。

def lasso_regression(X, y, lambd=0.2, threshold=0.1):n 通過坐標下降(coordinate descent)法獲取LASSO回歸係數n n # 計算殘差平方和n rss = lambda X, y, w: (y - X*w).T*(y - X*w)n # 初始化回歸係數w.n m, n = X.shapen w = np.matrix(np.zeros((n, 1)))n r = rss(X, y, w)n # 使用坐標下降法優化回歸係數wn niter = itertools.count(1)n for it in niter:n for k in range(n):n # 計算常量值z_k和p_kn z_k = (X[:, k].T*X[:, k])[0, 0]n p_k = 0n for i in range(m):n p_k += X[i, k]*(y[i, 0] - sum([X[i, j]*w[j, 0] for j in range(n) if j != k]))n if p_k < -lambd/2:n w_k = (p_k + lambd/2)/z_kn elif p_k > lambd/2:n w_k = (p_k - lambd/2)/z_kn else:n w_k = 0n w[k, 0] = w_kn r_prime = rss(X, y, w)n delta = abs(r_prime - r)[0, 0]n r = r_primen print(Iteration: {}, delta = {}.format(it, delta))n if delta < threshold:n breakn return wn

我們選取 $lambda = 10$ , 收斂閾值為0.1來獲取回歸係數

if __main__ == __name__:n X, y = load_data(abalone.txt)n X, y = standarize(X), standarize(y)n w = lasso_regression(X, y, lambd=10)n y_prime = X*wn # 計算相關係數n corrcoef = get_corrcoef(np.array(y.reshape(1, -1)),n np.array(y_prime.reshape(1, -1)))n print(Correlation coefficient: {}.format(corrcoef))n

迭代了150步收斂到0.1，計算相對比較耗時:

Iteration: 146, delta = 0.1081124857935265nIteration: 147, delta = 0.10565615985365184nIteration: 148, delta = 0.10326058648411163nIteration: 149, delta = 0.10092418256476776nIteration: 150, delta = 0.09864540659987142nCorrelation coefficient: 0.7255254877587117n

LASSO回歸係數軌跡

類似嶺軌跡，我們也可以改變λλ的值得到不同的回歸係數，通過作圖可以看到回歸係數的軌跡

ntest = 30n# 繪製軌跡nws = lasso_traj(X, y, ntest)nfig = plt.figure()nax = fig.add_subplot(111)nlambdas = [i-10 for i in range(ntest)]nax.plot(lambdas, ws)nplt.show()n

得到的軌跡圖如下:

通過與嶺軌跡圖進行對比發現，隨著λλ的增大，係數逐漸趨近於0，但是嶺回歸沒有係數真正為0，而lasso中不斷有係數變為0. 迭代過程中輸出如下圖:

逐步向前回歸

LASSO計算複雜度相對較高，本部分稍微介紹一下逐步向前回歸，他屬於一種貪心演算法，給定初始係數向量，然後不斷迭代遍歷每個係數，增加或減小一個很小的數，看看代價函數是否變小，如果變小就保留，如果變大就捨棄，然後不斷迭代直到回歸係數達到穩定。

下面給出實現

def stagewise_regression(X, y, eps=0.01, niter=100):n 通過向前逐步回歸獲取回歸係數n n m, n = X.shapen w = np.matrix(np.zeros((n, 1)))n min_error = float(inf)n all_ws = np.matrix(np.zeros((niter, n)))n # 計算殘差平方和n rss = lambda X, y, w: (y - X*w).T*(y - X*w)n for i in range(niter):n print({}: w = {}.format(i, w.T[0, :]))n for j in range(n):n for sign in [-1, 1]:n w_test = w.copy()n w_test[j, 0] += eps*signn test_error = rss(X, y, w_test)n if test_error < min_error:n min_error = test_errorn w = w_testn all_ws[i, :] = w.Tn return all_wsn

我們去變化量為0.005，迭代步數為1000次，得到回歸係數隨著迭代次數的變化曲線:

逐步回歸演算法的主要有點在於他可以幫助人們理解現有的模型並作出改進。當構建了一個模型後，可以運行逐步回歸演算法找出重要的特徵，即使停止那些不重要特徵的收集。

總結

本文介紹了兩種回歸中的縮減方法，嶺回歸和LASSO。兩種回歸均是在標準線性回歸的基礎上加上正則項來減小模型的方差。這裡其實便涉及到了權衡偏差(Bias)和方差(Variance)的問題。方差針對的是模型之間的差異，即不同的訓練數據得到模型的區別越大說明模型的方差越大。而偏差指的是模型預測值與樣本數據之間的差異。所以為了在過擬合和欠擬合之前進行權衡，我們需要確定適當的模型複雜度來使得總誤差最小。

參考

《Machine Learning in Action》
機器學習中的Bias(偏差)，Error(誤差)，和Variance(方差)有什麼區別和聯繫？
Lasso回歸的坐標下降法推導
數據什麼時候需要做中心化和標準化處理？