前言
繼續線性迴歸的總結, 本文主要介紹兩種線性迴歸的縮減(shrinkage)方法的基礎知識: 嶺迴歸(Ridge Regression)和LASSO(Least Absolute Shrinkage and Selection Operator)並對其進行了Python實現。同時也對一種更為簡單的向前逐步迴歸計算迴歸係數的方法進行了相應的實現。
正文
通過上一篇《機器學習算法實踐-標準與局部加權線性迴歸》中標準線性迴歸的公式w=(X^T*X)^(-1)X^T*y中可以看出在計算迴歸係數的時候我們需要計算矩陣X^TX的逆，但是如果該矩陣是個奇異矩陣，則無法對其進行求解。那麼什麼情況下該矩陣會有奇異性呢?
X本身存在線性相關關係(多重共線性), 即非滿秩矩陣。如果數據的特徵中存在兩個相關的變量，即使並不是完全線性相關，但是也會造成矩陣求逆的時候造成求解不穩定。
當數據特徵比數據量還要多的時候, 即m<n, 這時候矩陣XX是一個矮胖型的矩陣，非滿秩。
對於上面的兩種情況，我們需要對最初的標準線性迴歸做一定的變化使原先無法求逆的矩陣變得非奇異，使得問題可以穩定求解。我們可以通過縮減的方式來處理這些問題例如嶺迴歸和LASSO.
中心化和標準化
這裡先介紹下數據的中心化和標準化，在迴歸問題和一些機器學習算法中通常要對原始數據進行中心化和標準化處理，也就是需要將數據的均值調整到0，標準差調整為1, 計算過程很簡單就是將所有數據減去平均值後再除以標準差:
這樣調整後的均值:
調整後的標準差:
之所以需要進行中心化其實就是個平移過程，將所有數據的中心平移到原點。而標準化則是使得所有數據的不同特徵都有相同的尺度Scale, 這樣在使用梯度下降法以及其他方法優化的時候不同特徵參數的影響程度就會一致了。
如下圖所示，可以看出得到的標準化數據在每個維度上的尺度是一致的(圖片來自網絡，侵刪)
嶺迴歸(Ridge Regression)
標準最小二乘法優化問題:
也可以通過矩陣表示:
得到的迴歸係數為:
這個問題解存在且唯一的條件就是XX列滿秩:rank(X)=dim(X).
即使X列滿秩，但是當數據特徵中存在共線性，即相關性比較大的時候，會使得標準最小二乘求解不穩定, XTX的行列式接近零，計算XTX的時候誤差會很大。這個時候我們需要在cost function上添加一個懲罰項，稱為L2正則化。
這個時候的cost function的形式就為:
通過加入此懲罰項進行優化後，限制了迴歸係數wiwi的絕對值，數學上可以證明上式的等價形式如下:
其中t為某個閾值。
將嶺迴歸係數用矩陣的形式表示:
可以看到，就是通過將XTX加上一個單位矩陣是的矩陣變成非奇異矩陣並可以進行求你運算。
嶺迴歸的幾何意義
以兩個變量為例, 殘差平方和可以表示為w1,w2的一個二次函數，是一個在三維空間中的拋物面，可以用等值線來表示。而限制條件w21+w22<t， 相當於在二維平面的一個圓。這個時候等值線與圓相切的點便是在約束條件下的最優點，如下圖所示，
嶺迴歸的一些性質
當嶺參數λ=0時，得到的解是最小二乘解
當嶺參數λ趨向更大時，嶺迴歸係數wi趨向於0，約束項t很小
嶺迴歸的Python實現
通過矩陣的形式計算ŵ , 可以很簡單的實現
defridge_regression(X,y,lambd=0.2):
''' 獲取嶺迴歸係數
'''
XTX = X.T*X
m,_ = XTX.shape
I = np.matrix(np.eye(m))
w = (XTX + lambd*I).I*X.T*y
returnw
嶺跡圖
可以知道求得的嶺係數wi是嶺參數λ的函數，不同的λλ得到不同的嶺參數wi, 因此我們可以增大λλ的值來得到嶺迴歸係數的變化，以及嶺參數的變化軌跡圖(嶺跡圖), 不存在奇異性時，嶺跡圖應穩定的逐漸趨向於0。
通過嶺跡圖我們可以:
觀察較佳的λ取值
觀察變量是否有多重共線性
繪製嶺跡圖
上面我們通過函數ridge_regression實現了計算嶺迴歸係數的計算，我們使用《機器學習實戰》中的鮑魚年齡的數據來進行計算並繪製不同λλ的嶺參數變化的軌跡圖。數據以及完整代碼詳見 https://github.com/PytLab/MLBox/tree/master/linear_regression
選取30組不同的λλ來獲取嶺係數矩陣包含30個不同的嶺係數。
defridge_traj(X,y,ntest=30):
''' 獲取嶺軌跡矩陣
'''
_,n = X.shape
ws = np.zeros((ntest,n))
foriinrange(ntest):
w = ridge_regression(X,y,lambd=exp(i-10))
ws[i, :] = w.T
returnws
繪製嶺軌跡圖
if'__main__' == __name__:
ntest = 30
# 加載數據
X,y = load_data('abalone.txt')
# 中心化 & 標準化
X,y = standarize(X),standarize(y)
# 繪製嶺軌跡
ws = ridge_traj(X,y,ntest)
fig = plt.figure()
ax = fig.add_subplot(111)
lambdas = [i-10foriinrange(ntest)]
ax.plot(lambdas,ws)
plt.show()
上圖繪製了迴歸係數wi與log(λ)的關係，在最左邊λλ係數最小時，可以得到所有係數的原始值(與標準線性迴歸相同); 而在右邊，係數全部縮減為0, 從不穩定趨於穩定；為了定量的找到最佳參數值，還需要進行交叉驗證。要判斷哪些變量對結果的預測最具影響力，可以觀察他們的係數大小即可。
LASSO
嶺迴歸限定了所有迴歸係數的平方和不大於tt, 在使用普通最小二乘法迴歸的時候當兩個變量具有相關性的時候，可能會使得其中一個係數是個很大正數，另一個係數是很大的負數。通過嶺迴歸正則項的限制，可以避免這個問題。
LASSO(The Least Absolute Shrinkage and Selection Operator)是另一種縮減方法，將回歸係數收縮在一定的區域內。LASSO的主要思想是構造一個一階懲罰函數獲得一個精煉的模型, 通過最終確定一些變量的係數為0進行特徵篩選。
LASSO的懲罰項為:
與嶺迴歸的不同在於，此約束條件使用了絕對值的一階懲罰函數代替了平方和的二階函數。雖然只是形式稍有不同，但是得到的結果卻又很大差別。在LASSO中，當λλ很小的時候，一些係數會隨著變為0而嶺迴歸卻很難使得某個係數恰好縮減為0. 我們可以通過幾何解釋看到LASSO與嶺迴歸之間的不同。
LASSO的幾何解釋
同樣以兩個變量為例，標準線性迴歸的cost function還是可以用二維平面的等值線表示，而約束條件則與嶺迴歸的圓不同，LASSO的約束條件可以用方形表示，如下圖:
相比圓，方形的頂點更容易與拋物面相交，頂點就意味著對應的很多係數為0，而嶺迴歸中的圓上的任意一點都很容易與拋物面相交很難得到正好等於0的係數。這也就意味著，lasso起到了很好的篩選變量的作用。
LASSO迴歸係數的計算
雖然懲罰函數只是做了細微的變化，但是相比嶺迴歸可以直接通過矩陣運算得到迴歸係數相比，LASSO的計算變得相對複雜。由於懲罰項中含有絕對值，此函數的導數是連續不光滑的，所以無法進行求導並使用梯度下降優化。本部分使用座標下降發對LASSO迴歸係數進行計算。
座標下降法是每次選擇一個維度的參數進行一維優化，然後不斷的迭代對多個維度進行更新直到函數收斂。SVM對偶問題的優化算法SMO也是類似的原理，這部分的詳細介紹我在之前的一篇博客中進行了整理，參考《機器學習算法實踐-SVM中的SMO算法》。
下面我們分別對LASSO的cost function的兩部分求解：
1）RSS部分
求導:
1）正則項
關於懲罰項的求導我們需要使用subgradient，可以參考LASSO（least absolute shrinkage and selection operator） 迴歸中 如何用梯度下降法求解？
這樣整體的偏導數:
通過上面的公式我們便可以每次選取一維進行優化並不斷跌打得到最優迴歸係數。
LASSO的Python實現
根據上面代碼我們實現梯度下降法並使用其獲取LASSO迴歸係數。
deflasso_regression(X,y,lambd=0.2,threshold=0.1):
''' 通過座標下降(coordinate descent)法獲取LASSO迴歸係數
'''
# 計算殘差平方和
rss = lambdaX,y,w: (y - X*w).T*(y - X*w)
# 初始化迴歸係數w.
m,n = X.shape
w = np.matrix(np.zeros((n,1)))
r = rss(X,y,w)
# 使用座標下降法優化迴歸係數w
niter = itertools.count(1)
forit inniter:
forkinrange(n):
# 計算常量值z_k和p_k
z_k = (X[:,k].T*X[:,k])[0,0]
p_k = 0
foriinrange(m):
p_k += X[i,k]*(y[i,0] - sum([X[i,j]*w[j,0]forjinrange(n)ifj != k]))
ifp_k < -lambd/2:
w_k = (p_k + lambd/2)/z_k
elifp_k > lambd/2:
w_k = (p_k - lambd/2)/z_k
else:
w_k = 0
w[k,0] = w_k
r_prime = rss(X,y,w)
delta = abs(r_prime - r)[0,0]
r = r_prime
print('Iteration: {}, delta = {}'.format(it,delta))
ifdelta < threshold:
break
returnw
我們選取λ=10, 收斂閾值為0.1來獲取迴歸係數
if'__main__' == __name__:
X,y = load_data('abalone.txt')
X,y = standarize(X),standarize(y)
w = lasso_regression(X,y,lambd=10)
y_prime = X*w
# 計算相關係數
corrcoef = get_corrcoef(np.array(y.reshape(1, -1)),
np.array(y_prime.reshape(1, -1)))
print('Correlation coefficient: {}'.format(corrcoef))
迭代了150步收斂到0.1，計算相對比較耗時:
Iteration: 146,delta = 0.1081124857935265
Iteration: 147,delta = 0.10565615985365184
Iteration: 148,delta = 0.10326058648411163
Iteration: 149,delta = 0.10092418256476776
Iteration: 150,delta = 0.09864540659987142
Correlation coefficient: 0.7255254877587117
LASSO迴歸係數軌跡
類似嶺軌跡，我們也可以改變λλ的值得到不同的迴歸係數，通過作圖可以看到迴歸係數的軌跡
ntest = 30
# 繪製軌跡
ws = lasso_traj(X,y,ntest)
fig = plt.figure()
ax = fig.add_subplot(111)
lambdas = [i-10foriinrange(ntest)]
ax.plot(lambdas,ws)
plt.show()
得到的軌跡圖如下:
通過與嶺軌跡圖進行對比發現，隨著λλ的增大，係數逐漸趨近於0，但是嶺迴歸沒有係數真正為0，而lasso中不斷有係數變為0.
逐步向前回歸
LASSO計算複雜度相對較高，本部分稍微介紹一下逐步向前回歸，他屬於一種貪心算法，給定初始係數向量，然後不斷迭代遍歷每個係數，增加或減小一個很小的數，看看代價函數是否變小，如果變小就保留，如果變大就捨棄，然後不斷迭代直到迴歸係數達到穩定。
下面給出實現
defstagewise_regression(X,y,eps=0.01,niter=100):
''' 通過向前逐步迴歸獲取迴歸係數
'''
m,n = X.shape
w = np.matrix(np.zeros((n,1)))
min_error = float('inf')
all_ws = np.matrix(np.zeros((niter,n)))
# 計算殘差平方和
rss = lambdaX,y,w: (y - X*w).T*(y - X*w)
foriinrange(niter):
print('{}: w = {}'.format(i,w.T[0, :]))
forjinrange(n):
forsign in[-1,1]:
w_test = w.copy()
w_test[j,0] += eps*sign
test_error = rss(X,y,w_test)
iftest_error < min_error:
min_error = test_error
w = w_test
all_ws[i, :] = w.T
returnall_ws
我們去變化量為0.005，迭代步數為1000次，得到迴歸係數隨著迭代次數的變化曲線:
逐步迴歸算法的主要有點在於他可以幫助人們理解現有的模型並作出改進。當構建了一個模型後，可以運行逐步迴歸算法找出重要的特徵，即使停止那些不重要特徵的收集。
總結
本文介紹了兩種迴歸中的縮減方法，嶺迴歸和LASSO。兩種迴歸均是在標準線性迴歸的基礎上加上正則項來減小模型的方差。這裡其實便涉及到了權衡偏差(Bias)和方差(Variance)的問題。方差針對的是模型之間的差異，即不同的訓練數據得到模型的區別越大說明模型的方差越大。而偏差指的是模型預測值與樣本數據之間的差異。所以為了在過擬合和欠擬合之前進行權衡，我們需要確定適當的模型複雜度來使得總誤差最小。
End.
運行人員：中國統計網小編（微信號：itongjilove）
微博ID：中國統計網
中國統計網，是國內最早的大數據學習網站，公眾號：中國統計網
//www.itongji.cn