注：正則化是用來防止過擬合的方法。在最開始學習機器學習的課程時，只是覺得這個方法就像某種魔法一樣非常神奇的改變了模型的引數。但是一直也無法對其基本原理有一個透徹、直觀的理解。直到最近再次接觸到這個概念，經過一番苦思冥想後終於有了我自己的理解。

0. 正則化（Regularization ）

前面使用多項式迴歸，如果多項式最高次項比較大，模型就容易出現過擬合。正則化是一種常見的防止過擬合的方法，一般原理是在代價函式後面加上一個對引數的約束項，這個約束項被叫做正則化項（regularizer）。線上性迴歸模型中，通常有兩種不同的正則化項：

• 加上所有引數（不包括θ0

）的絕對值之和，即l1

• 範數，此時叫做Lasso迴歸；
• 加上所有引數（不包括θ0
• ）的平方和，即l2
  • 範數，此時叫做嶺迴歸.

  看過不少關於正則化原理的解釋，但是都沒有獲得一個比較直觀的理解。下面用代價函式的影象以及正則化項的影象來幫助解釋正則化之所以起作用的原因。

  0.1 代價函式的影象

  為了視覺化，選擇直線方程進行優化。假設一個直線方程以及代價函式如下：

  h^θ=θ0+θ1x

  ，該方程只有一個特徵x，兩個引數θ0和θ1

  J(θ)=1m∑mi=1(θ0+θ1x(i)−y(i))2

  ，該代價函式為均方誤差函式（MSE），其中m

  表示樣本量.

  為了保持簡單，只取一個樣本點(1,1)

  代入上面的代價函式方程中，可得J(θ)=(θ0+θ1−1)2

  . 該式是一個二元一次方程，可以在3維空間中作圖（下面利用網站

  GeoGebra畫出該方程的影象）：

  

  圖0-1，代入樣本點(1,1)

  後的代價函式MSE的影象

  由於多個樣本點的代價函式是所有樣本點代價函式之和，且不同的樣本點只是相當於改變了代價函式中兩個變數的引數（此時θ0

  和θ1

  是變數，樣本點的取值時引數）。因此多樣本的代價函式MSE的影象只會在圖0-1上發生縮放和平移，而不會發生過大的形變。

  對於座標軸，表示如下：

  • 使用J
• 軸表示藍色軸線，上方為正向；
• 使用θ1
• 表示紅色軸線，左邊為正向；
• 使用θ0

• 表示綠色軸線，指向螢幕外的方向為正向.

此時的函式影象相當於一條拋物線沿著平面J=0

上直線θ0=−θ1

平移後形成的影象。

0.2 正則化項的影象

這裡使用L1

範數作為正則化項，加上正則化項之後MSE代價函式變成：

J(θ)=1m∑mi=1(θ0+θ1x(i)−y(i))2+λ||θ1||1

,

上式中λ

是正則化項的引數，為了簡化取λ=1。由於正則化項中始終不包含截距項θ0，此時的L1範數相當於引數θ1

的絕對值，函式影象如下：

圖0-2，L1

正則化項的影象

此時的函式影象相當於一張對摺後，半張開的紙。紙的摺痕與平面J=0

上θ0

軸重疊。

0.3 代價函式與正則化項影象的疊加

直接將這兩個影象放在一起的樣子：

圖0-3，同時顯示代價函式與正則化項的影象

將兩個方程相加之後，即J(θ)=(θ0+θ1−1)2+|θ1|

，做圖可以得到下面的影象：

圖0-4，加入正則化項之後代價函式的影象

此時的影象，就像是一個圓錐體被捏扁了之後，立在座標原點上。觀察新增正則化項前後的影象，我們會發現：

• 加上正則化項之後，此時損失函式就分成了兩部分：第1項為原來的MSE函式，第2項為正則化項，最終的結果是這兩部分的線性組合;
• 在第1項的值非常小但在第2項的值非常大的區域，這些值會受到正則化項的巨大影響，從而使得這些區域的值變的與正則化項近似：例如原來的損失函式沿θ0=−θ1

，J軸方向上的值始終為0，但是加入正則化項J=|θ1|後，該直線上原來為0的點，都變成了θ1

• 的絕對值。這就像加權平均值一樣，哪一項的權重越大，對最終結果產生的影響也越大;
• 如果想象一種非常極端的情況：在引數的整個定義域上，第2項的取值都遠遠大於第一項的取值，那麼最終的損失函式幾乎100%都會由第2項決定，也就是整個代價函式的影象會非常類似於J=|θ1|
• （圖0-2）而不是原來的MSE函式的影象（圖0-1）。這時候就相當於λ
  • 的取值過大的情況，最終的全域性最優解將會是座標原點，這就是為什麼在這種情況下最終得到的解全都為0.

  1. 嶺迴歸

  嶺迴歸與多項式迴歸唯一的不同在於代價函式上的差別。嶺迴歸的代價函式如下：

  J(θ)=1m∑i=1m(y(i)−(wx(i)+b))2+λ||w||22=MSE(θ)+λ∑i=1nθ2i ⋯ (1−1)

  為了方便計算導數，通常也寫成下面的形式：

  J(θ)=12m∑i=1m(y(i)−(wx(i)+b))2+λ2||w||22=12MSE(θ)+λ2∑i=1nθ2i ⋯ (1−2)

  上式中的w

  是長度為n的向量，不包括截距項的係數θ0；θ是長度為n+1的向量，包括截距項的係數θ0；m為樣本數；n

  為特徵數.

  嶺迴歸的代價函式仍然是一個凸函式，因此可以利用梯度等於0的方式求得全域性最優解（正規方程）：

  θ=(XTX+λI)−1(XTy)

  上述正規方程與一般線性迴歸的正規方程相比，多了一項λI

  ，其中I表示單位矩陣。假如XTX

  是一個奇異矩陣（不滿秩），新增這一項後可以保證該項可逆。由於單位矩陣的形狀是對角線上為1其他地方都為0，看起來像一條山嶺，因此而得名。

  除了上述正規方程之外，還可以使用梯度下降的方式求解（求梯度的過程可以參考一般線性迴歸，3.2.2節）。這裡採用式子1−2

  來求導：

  ∇θJ(θ)=1mXT⋅(X⋅θ−y)+λw ⋯ (1−3)

  因為式子1−2

  中和式第二項不包含θ0，因此求梯度後，上式第二項中的w本來也不包含θ0。為了計算方便，新增θ0=0到w

  .

  因此在梯度下降的過程中，引數的更新可以表示成下面的公式：

  θ=θ−(αmXT⋅(X⋅θ−y)+λw) ⋯ (1−4)

  其中α

  為學習率，λ

  為正則化項的引數

  1.1 資料以及相關函式

• import numpy as np
  import matplotlib.pyplot as plt
  from sklearn.preprocessing import PolynomialFeatures
  from sklearn.metrics import mean_squared_error
  
  data = np.array([[ -2.95507616,  10.94533252],
         [ -0.44226119,   2.96705822],
         [ -2.13294087,   6.57336839],
         [  1.84990823,   5.44244467],
         [  0.35139795,   2.83533936],
         [ -1.77443098,   5.6800407 ],
         [ -1.8657203 ,   6.34470814],
         [  1.61526823,   4.77833358],
         [ -2.38043687,   8.51887713],
         [ -1.40513866,   4.18262786]])
  m = data.shape[0]  # 樣本大小
  X = data[:, 0].reshape(-1, 1)  # 將array轉換成矩陣
  y = data[:, 1].reshape(-1, 1)

  繼續使用多項式迴歸中的資料。

  1.2 嶺迴歸的手動實現

   有了上面的理論基礎，就可以自己實現嶺迴歸了，下面是Python程式碼：

• # 代價函式
  def L_theta(theta, X_x0, y, lamb):
      """
      lamb: lambda, the parameter of regularization
      theta: (n+1)·1 matrix, contains the parameter of x0=1
      X_x0: m·(n+1) matrix, plus x0
      """
      h = np.dot(X_x0, theta)  # np.dot 表示矩陣乘法
      theta_without_t0 = theta[1:]
      L_theta = 0.5 * mean_squared_error(h, y) + 0.5 * lamb * np.sum(np.square(theta_without_t0))
      return L_theta
  
  # 梯度下降
  def GD(lamb, X_x0, theta, y, alpha):
      """
      lamb: lambda, the parameter of regularization
      alpha: learning rate
      X_x0: m·(n+1), plus x0
      theta: (n+1)·1 matrix, contains the parameter of x0=1
      """
      for i in range(T):
          h = np.dot(X_x0, theta) 
          theta_with_t0_0 = np.r_[np.zeros([1, 1]), theta[1:]]  # set theta[0] = 0
          theta -= (alpha * 1/m * np.dot(X_x0.T, h - y) + lamb*(theta_with_t0_0))  # add the gradient of regularization term
          if i%50000==0:
              print(L_theta(theta, X_x0, y, lamb))
      return theta
  
  T = 1200000  # 迭代次數
  degree = 11
  theta = np.ones((degree + 1, 1))  # 引數的初始化，degree = 11，一個12個引數
  alpha = 0.0000000006 # 學習率
  # alpha = 0.003  # 學習率
  lamb = 0.0001
  # lamb = 0
  poly_features_d = PolynomialFeatures(degree=degree, include_bias=False)
  X_poly_d = poly_features_d.fit_transform(X)
  X_x0 = np.c_[np.ones((m, 1)), X_poly_d]  # ADD X0 = 1 to each instance
  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

  從上面的結果看，截距項的引數最大，高階項的引數都比較小。下面是比較原始資料和訓練出來的模型之間的關係：

• X_plot = np.linspace(-2.99, 1.9, 1000).reshape(-1, 1)
  poly_features_d_with_bias = PolynomialFeatures(degree=degree, include_bias=True)
  X_plot_poly = poly_features_d_with_bias.fit_transform(X_plot)
  y_plot = np.dot(X_plot_poly, theta)
  plt.plot(X_plot, y_plot, 'r-')
  plt.plot(X, y, 'b.')
  plt.xlabel('x')
  plt.ylabel('y')
  plt.show()

  

   圖1-1，手動實現嶺迴歸的效果

  圖中模型與原始資料的匹配度不是太好，但是過擬合的情況極大的改善了，模型變的更簡單了。

  1.2 正規方程

  下面使用正規方程求解：

   其中λ=10

• theta2 = np.linalg.inv(np.dot(X_x0.T, X_x0) + 10*np.identity(X_x0.shape[1])).dot(X_x0.T).dot(y)
  print(theta2)
  print(L_theta(theta2, X_x0, y, lamb))
  
  X_plot = np.linspace(-3, 2, 1000).reshape(-1, 1)
  poly_features_d_with_bias = PolynomialFeatures(degree=degree, include_bias=True)
  X_plot_poly = poly_features_d_with_bias.fit_transform(X_plot)
  y_plot = np.dot(X_plot_poly, theta2)
  plt.plot(X_plot, y_plot, 'r-')
  plt.plot(X, y, 'b.')
  plt.xlabel('x')
  plt.ylabel('y')
  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: 上面公式中的λ
• ，正則化項的係數；
• solver: 求解方法；
• X: 訓練樣本；
• y: 訓練樣本的標籤.

• from sklearn.linear_model import Ridge
  
  # 代價函式
  def L_theta_new(intercept, coef, X, y, lamb):
      """
      lamb: lambda, the parameter of regularization
      theta: (n+1)·1 matrix, contains the parameter of x0=1
      X_x0: m·(n+1) matrix, plus x0
      """
      h = np.dot(X, coef) + intercept  # np.dot 表示矩陣乘法
      L_theta = 0.5 * mean_squared_error(h, y) + 0.5 * lamb * np.sum(np.square(coef))
      return L_theta
  
  lamb = 10
  ridge_reg = Ridge(alpha=lamb, solver="cholesky")
  ridge_reg.fit(X_poly_d, y)
  print(ridge_reg.intercept_, ridge_reg.coef_)
  print(L_theta_new(intercept=ridge_reg.intercept_, coef=ridge_reg.coef_.T, X=X_poly_d, y=y, lamb=lamb))
  
  X_plot = np.linspace(-3, 2, 1000).reshape(-1, 1)
  X_plot_poly = poly_features_d.fit_transform(X_plot)
  h = np.dot(X_plot_poly, ridge_reg.coef_.T) + ridge_reg.intercept_
  plt.plot(X_plot, h, 'r-')
  plt.plot(X, y, 'b.')
  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(θ)=12m∑i=1m(y(i)−(wx(i)+b))2+λ||w||1=12MSE(θ)+λ∑i=1n|θi| ⋯ (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

  的梯度為：

  ∇θMSE(θ)+λ⎛⎝⎜⎜⎜⎜sign(θ1)sign(θ2)⋮sign(θn)⎞⎠⎟⎟⎟⎟⋯ (2−2)

   其中sign(θi)

  由θi的符號決定: θi>0,sign(θi)=1; θi=0,sign(θi)=0; θi<0,sign(θi)=−1

  .

  2.1 Lasso的實現

  直接使用scikit-learn中的函式：

  可以參考官方文件，http://scikit-learn.org/stable/modules/generated/sklearn.linear_model.Lasso.html

  下面模型中的引數alpha就是公式(2-1)中的引數λ

  ，是正則化項的係數，可以取大於0的任意值。alpha的值越大，對模型中引數的懲罰力度越大，因此會有更多的引數被訓練為0（只對線性相關的引數起作用），模型也就變得更加簡單了。

• from sklearn.linear_model import Lasso
  
  lamb = 0.025
  lasso_reg = Lasso(alpha=lamb)
  lasso_reg.fit(X_poly_d, y)
  print(lasso_reg.intercept_, lasso_reg.coef_)
  print(L_theta_new(intercept=lasso_reg.intercept_, coef=lasso_reg.coef_.T, X=X_poly_d, y=y, lamb=lamb))
  
  X_plot = np.linspace(-3, 2, 1000).reshape(-1, 1)
  X_plot_poly = poly_features_d.fit_transform(X_plot)
  h = np.dot(X_plot_poly, lasso_reg.coef_.T) + lasso_reg.intercept_
  plt.plot(X_plot, h, 'r-')
  plt.plot(X, y, 'b.')
  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(θ)=12MSE(θ)+rλ∑i=1n|θi|+1−r2λ∑i=1nθ2i ⋯ (3−1)

  其中r

  表示l1

  所佔的比例。

  使用scikit-learn的實現：

• from sklearn.linear_model import ElasticNet
  
  # 代價函式
  def L_theta_ee(intercept, coef, X, y, lamb, r):
      """
      lamb: lambda, the parameter of regularization
      theta: (n+1)·1 matrix, contains the parameter of x0=1
      X_x0: m·(n+1) matrix, plus x0
      """
      h = np.dot(X, coef) + intercept  # np.dot 表示矩陣乘法
      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))
      return L_theta
  
  elastic_net = ElasticNet(alpha=0.5, l1_ratio=0.8)
  elastic_net.fit(X_poly_d, y)
  print(elastic_net.intercept_, elastic_net.coef_)
  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))
  
  X_plot = np.linspace(-3, 2, 1000).reshape(-1, 1)
  X_plot_poly = poly_features_d.fit_transform(X_plot)
  h = np.dot(X_plot_poly, elastic_net.coef_.T) + elastic_net.intercept_
  plt.plot(X_plot, h, 'r-')
  plt.plot(X, y, 'b.')
  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的比較

  根據吳恩達老師的機器學習公開課，建議使用下面的步驟來確定λ

  的值：

• 建立一個λ
• 值的列表，例如λ∈0,0.01,0.02,0.04,0.08,0.16,0.32,0.64,1.28,2.56,5.12,10.24
• ;
• 建立不同degree的模型（或改變其他變數）;
• 遍歷不同的模型和不同的λ
• 值;
• 使用學習到的引數θ
• （包含正則化項）計算驗證集上的誤差（計算誤差時不包含正則化項），JCV(θ)
• ;
• 選擇在驗證集上誤差最小的引數組合（degree和λ
• ）;
• 使用選出來的引數和λ
• 在測試集上測試，計算Jtest(θ)
  1. .

  下面通過一張影象來比較一下嶺迴歸和Lasso迴歸：

  

   圖4-1，Lasso與嶺迴歸的比較（俯瞰圖）

  上圖中，左上方表示l1

  （圖中菱形圖案）和代價函式（圖中深色橢圓環）；左下方表示l2（橢圓形線圈）和代價函式（圖中深色橢圓環）。同一條線上（或同一個環上），表示對應的函式值相同；圖案中心分別表示l1,l2

  範數以及代價函式的最小值位置。

  右邊表示代價函式加上對應的正則化項之後的影象。新增正則化項之後，會影響原來的代價函式的最小值的位置，以及梯度下降時的路線（如果引數調整合適的話，最小值應該在距離原來代價函式最小值附近且與正則化項的影象相交，因為此時這兩項在相互約束的情況下都取到最小值，它們的和也最小）。右上圖，顯示了Lasso迴歸中引數的變化情況，最終停留在了θ2=0

  這條線上；右下方的取值由於受到了l2

  範數的約束，也產生了位移。

  當正則化項的權重非常大的時候，會產生左側黃色點標識的路線，最終所有引數都為0，但是趨近原點的方式不同。這是因為對於範數來說，原點是它們的最小值點。