注：正則化是用來防止過擬合的方法。在最開始學習機器學習的課程時，只是覺得這個方法就像某種魔法一樣非常神奇的改變了模型的參數。但是一直也無法對其基本原理有一個透徹、直觀的理解。直到最近再次接觸到這個概念，經過一番苦思冥想后終於有了我自己的理解。
0. 正則化（Regularization ）
前面使用多項式回歸，如果多項式最高次項比較大，模型就容易出現過擬合。正則化是一種常見的防止過擬合的方法，一般原理是在代價函數后面加上一個對參數的約束項，這個約束項被叫做正則化項（regularizer）。在線性回歸模型中，通常有兩種不同的正則化項：
- 加上所有參數（不包括$\theta_0$）的絕對值之和，即$l1$范數，此時叫做Lasso回歸；
- 加上所有參數（不包括$\theta_0$）的平方和，即$l2$范數的平方，此時叫做嶺回歸.
看過不少關於正則化原理的解釋，但是都沒有獲得一個比較直觀的理解。下面用代價函數的圖像以及正則化項的圖像來幫助解釋正則化之所以起作用的原因。
0.1 代價函數的圖像
為了可視化，選擇直線方程進行優化。假設一個直線方程以及代價函數如下：
$\hat{h}_{\theta} = \theta_0 + \theta_1 x$，該方程只有一個特征$x$，兩個參數$\theta_0$和$\theta_1$
$J(\theta) = \frac{1}{m} \sum_{i=1}^{m}{(\theta_0 + \theta_1 x^{(i)} - y^{(i)})^2}$，該代價函數為均方誤差函數（MSE），其中$m$表示樣本量.
為了保持簡單，只取一個樣本點$(1, 1)$代入上面的代價函數方程中，可得$J(\theta) = (\theta_0 + \theta_1 - 1)^2$. 該式是一個二元一次方程，可以在3維空間中作圖（下面利用網站GeoGebra畫出該方程的圖像）：
圖0-1，代入樣本點$(1, 1)$后的代價函數MSE的圖像
由於多個樣本點的代價函數是所有樣本點代價函數之和，且不同的樣本點只是相當於改變了代價函數中兩個變量的參數（此時$\theta_0$和$\theta_1$是變量，樣本點的取值是參數）。因此多樣本的代價函數MSE的圖像只會在圖0-1上發生縮放和平移，而不會發生過大的形變。
對於坐標軸，表示如下：
- 使用$J$軸表示藍色軸線，上方為正向；
- 使用$\theta_1$表示紅色軸線，左邊為正向；
- 使用$\theta_0$表示綠色軸線，指向屏幕外的方向為正向.
此時的函數圖像相當於一條拋物線沿着平面$J = 0$上直線$\theta_0 = - \theta_1$平移后形成的圖像。
0.2 正則化項的圖像
這里使用$L1$范數作為正則化項，加上正則化項之后MSE代價函數變成：
$J(\theta) = \frac{1}{m} \sum_{i=1}^{m}{(\theta_0 + \theta_1 x^{(i)} - y^{(i)})^2} + \lambda ||\theta_1||_1$,
上式中$\lambda$是正則化項的參數，為了簡化取$\lambda = 1$。由於正則化項中始終不包含截距項$\theta_0$，此時的$L1$范數相當於參數$\theta_1$的絕對值，函數圖像如下：
圖0-2，$L1$正則化項的圖像
此時的函數圖像相當於一張對折后，半張開的紙。紙的折痕與平面$J = 0$上$\theta_0$軸重疊。
0.3 代價函數與正則化項圖像的疊加
直接將這兩個圖像放在一起的樣子：
圖0-3，同時顯示代價函數與正則化項的圖像
將兩個方程相加之后，即$J(\theta) = (\theta_0 + \theta_1 - 1)^2 + |\theta_1|$，做圖可以得到下面的圖像：
圖0-4，加入正則化項之后代價函數的圖像
此時的圖像，就像是一個圓錐體被捏扁了之后，立在坐標原點上。觀察添加正則化項前后的圖像，我們會發現：
- 加上正則化項之后，此時損失函數就分成了兩部分：第1項為原來的MSE函數，第2項為正則化項，最終的結果是這兩部分的線性組合;
- 在第1項的值非常小但在第2項的值非常大的區域，這些值會受到正則化項的巨大影響，從而使得這些區域的值變的與正則化項近似：例如原來的損失函數沿$\theta_0 = -\theta_1$，$J$軸方向上的值始終為0，但是加入正則化項$J = |\theta_1|$后，該直線上原來為0的點，都變成了$\theta_1$的絕對值。這就像加權平均值一樣，哪一項的權重越大，對最終結果產生的影響也越大;
- 如果想象一種非常極端的情況：在參數的整個定義域上，第2項的取值都遠遠大於第一項的取值，那么最終的損失函數幾乎100%都會由第2項決定，也就是整個代價函數的圖像會非常類似於$J=|\theta_1|$（圖0-2）而不是原來的MSE函數的圖像（圖0-1）。這時候就相當於$\lambda$的取值過大的情況，最終的全局最優解將會是坐標原點，這就是為什么在這種情況下最終得到的解全都為0.
1. 嶺回歸
嶺回歸與多項式回歸唯一的不同在於代價函數上的差別。嶺回歸的代價函數如下：
$$J(\theta) = \frac{1}{m} \sum_{i=1}^{m}{(y^{(i)} - (w x^{(i)} + b))^2} + \lambda ||w||_2^2 = MSE(\theta) + \lambda \sum_{i = 1}^{n}{\theta_i^2} \ \quad \cdots \ (1 - 1)$$
為了方便計算導數，通常也寫成下面的形式：
$$J(\theta) = \frac{1}{2m} \sum_{i=1}^{m}{(y^{(i)} - (w x^{(i)} + b))^2} + \frac{\lambda}{2} ||w||_2^2 = \frac{1}{2}MSE(\theta) + \frac{\lambda}{2} \sum_{i = 1}^{n}{\theta_i^2} \ \quad \cdots \ (1 - 2)$$
上式中的$w$是長度為$n$的向量，不包括截距項的系數$\theta_0$；$\theta$是長度為$n + 1$的向量，包括截距項的系數$\theta_0$；$m$為樣本數；$n$為特征數.
嶺回歸的代價函數仍然是一個凸函數，因此可以利用梯度等於0的方式求得全局最優解（正規方程）：
$$\theta = (X^T X + \lambda I)^{-1}(X^T y)$$
上述正規方程與一般線性回歸的正規方程相比，多了一項$\lambda I$，其中$I$表示單位矩陣。假如$X^T X$是一個奇異矩陣（不滿秩），添加這一項后可以保證該項可逆。由於單位矩陣的形狀是對角線上為1其他地方都為0，看起來像一條山嶺，因此而得名。
除了上述正規方程之外，還可以使用梯度下降的方式求解（求梯度的過程可以參考一般線性回歸，3.2.2節）。這里采用式子$1 - 2$來求導：
$$\nabla_{\theta} J(\theta) = \frac{1}{m} X^T \cdot (X \cdot \theta - y) + \lambda w \ \quad \cdots \ (1 - 3) $$
因為式子$1- 2$中和式第二項不包含$\theta_0$，因此求梯度后，上式第二項中的$w$本來也不包含$\theta_0$。為了計算方便，添加$\theta_0 = 0$到$w$.
因此在梯度下降的過程中，參數的更新可以表示成下面的公式：
$$\theta = \theta - (\frac{\alpha}{m} X^T \cdot (X \cdot \theta - y) + \lambda w) \ \quad \cdots \ (1 - 4) $$
其中$\alpha$為學習率，$\lambda$為正則化項的參數
1.1 數據以及相關函數
1 import numpy as np 2 import matplotlib.pyplot as plt 3 from sklearn.preprocessing import PolynomialFeatures 4 from sklearn.metrics import mean_squared_error 5 6 data = np.array([[ -2.95507616, 10.94533252], 7 [ -0.44226119, 2.96705822], 8 [ -2.13294087, 6.57336839], 9 [ 1.84990823, 5.44244467], 10 [ 0.35139795, 2.83533936], 11 [ -1.77443098, 5.6800407 ], 12 [ -1.8657203 , 6.34470814], 13 [ 1.61526823, 4.77833358], 14 [ -2.38043687, 8.51887713], 15 [ -1.40513866, 4.18262786]]) 16 m = data.shape[0] # 樣本大小 17 X = data[:, 0].reshape(-1, 1) # 將array轉換成矩陣 18 y = data[:, 1].reshape(-1, 1)
繼續使用多項式回歸中的數據。
1.2 嶺回歸的手動實現
有了上面的理論基礎，就可以自己實現嶺回歸了，下面是Python代碼：
1 # 代價函數 2 def L_theta(theta, X_x0, y, lamb): 3 """ 4 lamb: lambda, the parameter of regularization 5 theta: (n+1)·1 matrix, contains the parameter of x0=1 6 X_x0: m·(n+1) matrix, plus x0 7 """ 8 h = np.dot(X_x0, theta) # np.dot 表示矩陣乘法 9 theta_without_t0 = theta[1:] 10 L_theta = 0.5 * mean_squared_error(h, y) + 0.5 * lamb * np.sum(np.square(theta_without_t0)) 11 return L_theta 12 13 # 梯度下降 14 def GD(lamb, X_x0, theta, y, alpha): 15 """ 16 lamb: lambda, the parameter of regularization 17 alpha: learning rate 18 X_x0: m·(n+1), plus x0 19 theta: (n+1)·1 matrix, contains the parameter of x0=1 20 """ 21 for i in range(T): 22 h = np.dot(X_x0, theta) 23 theta_with_t0_0 = np.r_[np.zeros([1, 1]), theta[1:]] # set theta[0] = 0 24 theta -= (alpha * 1/m * np.dot(X_x0.T, h - y) + lamb*(theta_with_t0_0)) # add the gradient of regularization term 25 if i%50000==0: 26 print(L_theta(theta, X_x0, y, lamb)) 27 return theta 28 29 T = 1200000 # 迭代次數 30 degree = 11 31 theta = np.ones((degree + 1, 1)) # 參數的初始化，degree = 11，一個12個參數 32 alpha = 0.0000000006 # 學習率 33 # alpha = 0.003 # 學習率 34 lamb = 0.0001 35 # lamb = 0 36 poly_features_d = PolynomialFeatures(degree=degree, include_bias=False) 37 X_poly_d = poly_features_d.fit_transform(X) 38 X_x0 = np.c_[np.ones((m, 1)), X_poly_d] # ADD X0 = 1 to each instance 39 theta = GD(lamb=lamb, X_x0=X_x0, theta=theta, y=y, alpha=alpha)
上面第10行對應公式$1-2$，第24行對應公式$1-3$。由於自由度比較大，此時利用梯度下降的方法訓練模型比較困難，學習率稍微大一點就會出現出現損失函數的值越過最低點不斷增長的情況。下面是訓練結束后的參數以及代價函數值：
[[ 1.00078848e+00] [ -1.03862735e-05] [ 3.85144400e-05] [ -3.77233288e-05] [ 1.28959318e-04] [ -1.42449160e-04] [ 4.42760996e-04] [ -5.11518471e-04] [ 1.42533716e-03] [ -1.40265037e-03] [ 3.13638870e-03] [ 1.21862016e-03]] 3.59934190413
從上面的結果看，截距項的參數最大，高階項的參數都比較小。下面是比較原始數據和訓練出來的模型之間的關系：
1 X_plot = np.linspace(-2.99, 1.9, 1000).reshape(-1, 1) 2 poly_features_d_with_bias = PolynomialFeatures(degree=degree, include_bias=True) 3 X_plot_poly = poly_features_d_with_bias.fit_transform(X_plot) 4 y_plot = np.dot(X_plot_poly, theta) 5 plt.plot(X_plot, y_plot, 'r-') 6 plt.plot(X, y, 'b.') 7 plt.xlabel('x') 8 plt.ylabel('y') 9 plt.show()
圖1-1，手動實現嶺回歸的效果
圖中模型與原始數據的匹配度不是太好，但是過擬合的情況極大的改善了，模型變的更簡單了。
1.2 正規方程
下面使用正規方程求解：
其中$\lambda = 10$
1 theta2 = np.linalg.inv(np.dot(X_x0.T, X_x0) + 10*np.identity(X_x0.shape[1])).dot(X_x0.T).dot(y) 2 print(theta2) 3 print(L_theta(theta2, X_x0, y, lamb)) 4 5 X_plot = np.linspace(-3, 2, 1000).reshape(-1, 1) 6 poly_features_d_with_bias = PolynomialFeatures(degree=degree, include_bias=True) 7 X_plot_poly = poly_features_d_with_bias.fit_transform(X_plot) 8 y_plot = np.dot(X_plot_poly, theta2) 9 plt.plot(X_plot, y_plot, 'r-') 10 plt.plot(X, y, 'b.') 11 plt.xlabel('x') 12 plt.ylabel('y') 13 plt.show()
參數即代價函數的值：
[[ 0.56502653] [-0.12459546] [ 0.26772443] [-0.15642405] [ 0.29249514] [-0.10084392] [ 0.22791769] [ 0.1648667 ] [-0.05686718] [-0.03906615] [-0.00111673] [ 0.00101724]] 0.604428719639
從參數來看，截距項的系數減小了，1-7階都有比較大的參數都比較大，后面更高階項的參數越來越小，下面是函數圖像：
圖1-2，使用正規方程求解
從圖中可以看到，雖然模型的自由度沒變，還是11，但是過擬合的程度得到了改善。
1.3 使用scikit-learn
scikit-learn中有專門計算嶺回歸的函數，而且效果要比上面的方法好。使用scikit-learn中的嶺回歸，只需要輸入以下參數：
- alpha: 上面公式中的$\lambda$，正則化項的系數；
- solver: 求解方法；
- X: 訓練樣本；
- y: 訓練樣本的標簽.
1 from sklearn.linear_model import Ridge 2 3 # 代價函數 4 def L_theta_new(intercept, coef, X, y, lamb): 5 """ 6 lamb: lambda, the parameter of regularization 7 theta: (n+1)·1 matrix, contains the parameter of x0=1 8 X_x0: m·(n+1) matrix, plus x0 9 """ 10 h = np.dot(X, coef) + intercept # np.dot 表示矩陣乘法 11 L_theta = 0.5 * mean_squared_error(h, y) + 0.5 * lamb * np.sum(np.square(coef)) 12 return L_theta 13 14 lamb = 10 15 ridge_reg = Ridge(alpha=lamb, solver="cholesky") 16 ridge_reg.fit(X_poly_d, y) 17 print(ridge_reg.intercept_, ridge_reg.coef_) 18 print(L_theta_new(intercept=ridge_reg.intercept_, coef=ridge_reg.coef_.T, X=X_poly_d, y=y, lamb=lamb)) 19 20 X_plot = np.linspace(-3, 2, 1000).reshape(-1, 1) 21 X_plot_poly = poly_features_d.fit_transform(X_plot) 22 h = np.dot(X_plot_poly, ridge_reg.coef_.T) + ridge_reg.intercept_ 23 plt.plot(X_plot, h, 'r-') 24 plt.plot(X, y, 'b.') 25 plt.show()
訓練結束后得到的參數為（分別表示截距，特征的系數；代價函數的值）：
[ 3.03698398] [[ -2.95619849e-02 6.09137803e-02 -4.93919290e-02 1.10593684e-01 -4.65660197e-02 1.06387336e-01 5.14340826e-02 -2.29460359e-02 -1.12705709e-02 -1.73925386e-05 2.79198986e-04]] 0.213877232488
圖1-3，使用scikit-learn訓練嶺回歸
經過與前面兩種方法得到的結果比較，這里得到的曲線更加平滑，不僅降低了過擬合的風險，代價函數的值也非常低。
2. Lasso回歸
Lasso回歸於嶺回歸非常相似，它們的差別在於使用了不同的正則化項。最終都實現了約束參數從而防止過擬合的效果。但是Lasso之所以重要，還有另一個原因是：Lasso能夠將一些作用比較小的特征的參數訓練為0，從而獲得稀疏解。也就是說用這種方法，在訓練模型的過程中實現了降維(特征篩選)的目的。
Lasso回歸的代價函數為：
$$J(\theta) = \frac{1}{2m} \sum_{i=1}^{m}{(y^{(i)} - (w x^{(i)} + b))^2} + \lambda ||w||_1 = \frac{1}{2}MSE(\theta) + \lambda \sum_{i = 1}^{n}{|\theta_i|} \ \quad \cdots \ (2 - 1)$$
上式中的$w$是長度為$n$的向量，不包括截距項的系數$θ_0$, $θ$是長度為$n+1$的向量，包括截距項的系數$θ_0$，$m$為樣本數，$n$為特征數.
$||w||_1$表示參數$w$的$l1$范數，也是一種表示距離的函數。加入$w$表示3維空間中的一個點$(x, y, z)$，那么$||w||_1 = |x| + |y| + |z|$，即各個方向上的絕對值（長度）之和。
式子$2-1$的梯度為：
$$\nabla_{\theta}MSE(\theta) + \lambda \begin{pmatrix} sign(\theta_1) \\ sign(\theta_2) \\ \vdots \\ sign(\theta_n) \end{pmatrix} \quad \cdots \ (2-2)$$
其中$sign(\theta_i)$由$\theta_i$的符號決定: $\theta_i > 0, sign(\theta_i) = 1; \ \theta_i = 0, sign(\theta_i) = 0; \ \theta_i < 0, sign(\theta_i) = -1$.
2.1 Lasso的實現
直接使用scikit-learn中的函數：
可以參考官方文檔，http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Lasso.html
下面模型中的參數alpha就是公式(2-1)中的參數$\lambda$，是正則化項的系數，可以取大於0的任意值。alpha的值越大，對模型中參數的懲罰力度越大，因此會有更多的參數被訓練為0（只對線性相關的參數起作用），模型也就變得更加簡單了。
1 from sklearn.linear_model import Lasso 2 3 lamb = 0.025 4 lasso_reg = Lasso(alpha=lamb) 5 lasso_reg.fit(X_poly_d, y) 6 print(lasso_reg.intercept_, lasso_reg.coef_) 7 print(L_theta_new(intercept=lasso_reg.intercept_, coef=lasso_reg.coef_.T, X=X_poly_d, y=y, lamb=lamb)) 8 9 X_plot = np.linspace(-3, 2, 1000).reshape(-1, 1) 10 X_plot_poly = poly_features_d.fit_transform(X_plot) 11 h = np.dot(X_plot_poly, lasso_reg.coef_.T) + lasso_reg.intercept_ 12 plt.plot(X_plot, h, 'r-') 13 plt.plot(X, y, 'b.') 14 plt.show()
最終獲得的參數以及代價函數的值為：
其中計算代價函數值的函數"L_theta_new"需要修改其中的"L_theta"為"L_theta = 0.5 * mean_squared_error(h, y) + lamb * np.sum(np.abs(coef))"
[ 2.86435179] [ -0.00000000e+00 5.29099723e-01 -3.61182017e-02 9.75614738e-02 1.61971116e-03 -3.42711766e-03 2.78782527e-04 -1.63421713e-04 -5.64291215e-06 -1.38933655e-05 1.02036898e-06]
0.0451291096773
從結果可以看到，截距項的值最大，一次項的系數為0，二次項的系數是剩下的所有項中值最大的，也比較符合數據的真實來源。這里也可以看出來，更高階的項雖然系數都非常小但不為0，這是因為這些項之間的關系是非線性的，無法用線性組合互相表示。
圖2-1，Lasso回歸得到的圖像
圖2-1是目前在$degree=11$的情況下，得到的最好模型。
3. 彈性網絡（ Elastic Net）
彈性網絡是結合了嶺回歸和Lasso回歸，由兩者加權平均所得。據介紹這種方法在特征數大於訓練集樣本數或有些特征之間高度相關時比Lasso更加穩定。
其代價函數為：
$$J(\theta) = \frac{1}{2}MSE(\theta) + r\lambda \sum_{i = 1}^{n}{|\theta_i|} + \frac{1-r}{2} \lambda \sum_{i=1}^{n} {\theta_i^2} \ \quad \cdots \ (3 - 1)$$
其中$r$表示$l1$所占的比例。
使用scikit-learn的實現：
1 from sklearn.linear_model import ElasticNet 2 3 # 代價函數 4 def L_theta_ee(intercept, coef, X, y, lamb, r): 5 """ 6 lamb: lambda, the parameter of regularization 7 theta: (n+1)·1 matrix, contains the parameter of x0=1 8 X_x0: m·(n+1) matrix, plus x0 9 """ 10 h = np.dot(X, coef) + intercept # np.dot 表示矩陣乘法 11 L_theta = 0.5 * mean_squared_error(h, y) + r * lamb * np.sum(np.abs(coef)) + 0.5 * (1-r) * lamb * np.sum(np.square(coef)) 12 return L_theta 13 14 elastic_net = ElasticNet(alpha=0.5, l1_ratio=0.8) 15 elastic_net.fit(X_poly_d, y) 16 print(elastic_net.intercept_, elastic_net.coef_) 17 print(L_theta_ee(intercept=elastic_net.intercept_, coef=elastic_net.coef_.T, X=X_poly_d, y=y, lamb=0.1, r=0.8)) 18 19 X_plot = np.linspace(-3, 2, 1000).reshape(-1, 1) 20 X_plot_poly = poly_features_d.fit_transform(X_plot) 21 h = np.dot(X_plot_poly, elastic_net.coef_.T) + elastic_net.intercept_ 22 plt.plot(X_plot, h, 'r-') 23 plt.plot(X, y, 'b.') 24 plt.show()
得到的結果為：
[ 3.31466833] [ -0.00000000e+00 0.00000000e+00 -0.00000000e+00 1.99874040e-01 -1.21830209e-02 2.58040545e-04 3.01117857e-03 -8.54952421e-04 4.35227606e-05 -2.84995639e-06 -8.36248799e-06] 0.0807738447192
該方法中得到了，更多的0，當然這也跟參數的設置有關。
圖3-1，使用elastic-net得到的結果
4. 正則化項的使用以及l1與l2的比較
根據吳恩達老師的機器學習公開課，建議使用下面的步驟來確定$\lambda$的值：
- 創建一個$\lambda$值的列表，例如$\lambda \in {0, 0.01, 0.02, 0.04, 0.08, 0.16, 0.32, 0.64, 1.28, 2.56, 5.12, 10.24}$;
- 創建不同degree的模型（或改變其他變量）;
- 遍歷不同的模型和不同的$\lambda$值;
- 使用學習到的參數$\theta$（包含正則化項）計算驗證集上的誤差（計算誤差時不包含正則化項），$J_{CV}(\theta)$;
- 選擇在驗證集上誤差最小的參數組合（degree和$\lambda$）;
- 使用選出來的參數和$\lambda$在測試集上測試，計算$J_{test}(\theta)$.
下面通過一張圖像來比較一下嶺回歸和Lasso回歸：
圖4-1，Lasso與嶺回歸的比較（俯瞰圖）
上圖中，左上方表示$l1$（圖中菱形圖案）和代價函數（圖中深色橢圓環）；左下方表示$l2$（橢圓形線圈）和代價函數（圖中深色橢圓環）。同一條線上（或同一個環上），表示對應的函數值相同；圖案中心分別表示$l1, l2$范數以及代價函數的最小值位置。
右邊表示代價函數加上對應的正則化項之后的圖像。添加正則化項之后，會影響原來的代價函數的最小值的位置，以及梯度下降時的路線（如果參數調整合適的話，最小值應該在距離原來代價函數最小值附近且與正則化項的圖像相交，因為此時這兩項在相互約束的情況下都取到最小值，它們的和也最小）。右上圖，顯示了Lasso回歸中參數的變化情況，最終停留在了$\theta_2 = 0$這條線上；右下方的取值由於受到了$l2$范數的約束，也產生了位移。
當正則化項的權重非常大的時候，會產生左側黃色點標識的路線，最終所有參數都為0，但是趨近原點的方式不同。這是因為對於范數來說，原點是它們的最小值點。
修改記錄：
7/14, 2020: 經評論區@小雞快跑learning指出，嶺回歸中加在代價函數后面的各項平方和是l2范數的平方，已更正
Reference
http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Ridge.html
Géron A. Hands-on machine learning with Scikit-Learn and TensorFlow: concepts, tools, and techniques to build intelligent systems[M]. " O'Reilly Media, Inc.", 2017. github
https://www.coursera.org/learn/machine-learning
edx: UCSanDiegoX - DSE220x Machine Learning Fundamentals