機器學習及演算法

林嶔 (Lin, Chin)

Lesson 11 機器學習概論3(正則化與彈性網路)

第一節:正則化(1)

– 由於數學上的相似性,我們已經教過的4種統計模型可以利用相同的技術進行擴展,這其實是相當方便的。

– 在統計問題中,有一系列的統計方法用到了正則化技術,而主要所希望解決的問題是「變數多個案少」的情況,我們將統計領域中常用的正則化技術分為兩類:

  1. L1 Regularization - 使用這個方法的迴歸模型稱為套索迴歸(LASSO regression)

  2. L2 Regularization - 使用這個方法的迴歸模型稱為脊迴歸(Ridge regression)

– 請至這裡下載範例資料

dat <- read.csv("ECG_train.csv", header = TRUE, fileEncoding = 'CP950', stringsAsFactors = FALSE, na.strings = "")

第一節:正則化(2)

\[ \begin{align} \hat{y_i} & = b_{0} + b_{1}x_i \\ loss & = \frac{{1}}{2n}\sum \limits_{i=1}^{n} \left(y_{i} - \hat{y_i}\right)^{2} \end{align} \]

– 而L1 Regularization的意思,就是對\(loss\)的部分做一些修正: (這裡\(\lambda\)是一個懲罰權重)

\[ \begin{align} loss & = \frac{{1}}{2n}\sum \limits_{i=1}^{n} \left(y_{i} - \hat{y_i}\right)^{2} + \lambda ( |b_0| + |b_1|) \end{align} \]

– L2 Regularization與L1 Regularization非常相似,只是對\(loss\)的修正方式不同:

\[ \begin{align} loss & = \frac{{1}}{2n}\sum \limits_{i=1}^{n} \left(y_{i} - \hat{y_i}\right)^{2} + \frac{\lambda}{2} ( b_0^2 + b_1^2) \end{align} \]

– 而我們再考慮一下他的物理意義,權重為0代表該項輸入並不重要,因此也代表該變項可有可無,所以使用正則化技術會讓許多變項的權重為0,這對於變項很多需要選變項時非常好用。

第一節:正則化(3)

– 整合預測函數與損失函數

\[loss = \frac{{1}}{2n}\sum \limits_{i=1}^{n} \left(y_{i} - \left(b_{0} + b_{1}x_{i}\right)\right)^{2} + \lambda ( |b_0| + |b_1|)\]

\(b_0\)以及\(b_1\)的偏導函數(過程略)

\[ \begin{align} \frac{\partial}{\partial b_0}loss & = \frac{{1}}{n} \sum \limits_{i=1}^{n}\left(y_{i} - \left(b_{0} + b_{1}x_{i}\right)\right) + \lambda \frac{b_0}{|b_0|} \\ \frac{\partial}{\partial b_1}loss & = \frac{{1}}{n} \sum \limits_{i=1}^{n}\left(y_{i} - \left(b_{0} + b_{1}x_{i}\right)\right)\cdot x_{i} + \lambda \frac{b_1}{|b_1|} \end{align} \]

– 另外,他似乎只提供了方向,而不提供大小。這是因為對於絕對值函數而言,導函數不是\(1\)就是\(-1\)

– 因此在訓練的時候,假設某個參數\(b\)的重要性實在太低,那由原始損失函數部分(左邊)所提供的梯度就會太小,而該參數\(b\)的數值就會在\(\lambda\)\(-\lambda\)之間震盪。

– 由於我們非常清楚這個介於在\(\lambda\)\(-\lambda\)之間是完全沒有意義的,因而能將所有介於\(\lambda\)\(-\lambda\)之間的參數\(b\)指定為\(0\),這樣暗示著這個參數所對應到的特徵\(x\)是完全沒有意義的。

– 從以上觀點來看,L1 Regularization具有很強的「稀疏化」特性,可以拿來做變項挑選。

第一節:正則化(4)

– 這個問題則可以被L2 Regularization解決,我們再看看他對於\(b_0\)以及\(b_1\)的偏導函數為何:

\[loss = \frac{{1}}{2n}\sum \limits_{i=1}^{n} \left(y_{i} - \left(b_{0} + b_{1}x_{i}\right)\right)^{2} + \frac{\lambda}{2} ( b_0^2 + b_1^2)\]

\(b_0\)以及\(b_1\)的偏導函數(過程略)

\[ \begin{align} \frac{\partial}{\partial b_0}loss & = \frac{{1}}{n} \sum \limits_{i=1}^{n}\left(y_{i} - \left(b_{0} + b_{1}x_{i}\right)\right) + \lambda b_0\\ \frac{\partial}{\partial b_1}loss & = \frac{{1}}{n} \sum \limits_{i=1}^{n}\left(y_{i} - \left(b_{0} + b_{1}x_{i}\right)\right)\cdot x_{i} + \lambda b_1 \end{align} \]

– 很明顯的,當\(\lambda\)很大的時候,那所有權重都將會很難離開0,但與L1 Regularization正則化不同,單純使用L2 Regularization經常會造成很多不重要的參數\(b\)有個很小的數值,卻沒有辦法捨去。

練習1:在Cox proportional hazard model中套用正則化的技術(1)

– 我們只使用前50筆資料加快他的速度:

subdat <- dat[!(dat[,"time"] %in% c(NA, 0)) & !(dat[,"death"] %in% NA) & !(dat[,"AGE"] %in% NA) & !(dat[,"Rate"] %in% NA),]
subdat <- subdat[1:50,]
survival_encode <- function (event, time) {
  
  sub_dat <- data.frame(id = 1:length(event), event = event, time = time, stringsAsFactors = FALSE)
  sub_dat <- sub_dat[order(sub_dat[,'event'], decreasing = TRUE),]
  sub_dat <- sub_dat[order(sub_dat[,'time']),]
  
  label <- array(0L, dim = rep(nrow(sub_dat), 2))
  mask <- array(0L, dim = rep(nrow(sub_dat), 2))
  
  diag(label) <- 1L
  
  for (m in 1:nrow(sub_dat)) {
    
    mask[m,m:nrow(sub_dat)] <- 1L
    
    if (sub_dat[m,'event'] == 0) {
      
      mask[m,] <- 0L
      
    } else if (m > 1) {
      
      zero_pos <- which(sub_dat[1:(m-1),'time'] == sub_dat[m,'time'])
      
      if (length(zero_pos) > 0) {
        
        mask[zero_pos,m] <- 0L
        
      }
      
    }
    
  }
  
  original_pos <- order(sub_dat[,'id'])
  
  return(list(label = label[,original_pos,drop = FALSE], mask = mask[,original_pos,drop = FALSE]))
  
}

練習1:在Cox proportional hazard model中套用正則化的技術(2)

X_mat <- matrix(as.matrix(subdat[,c("AGE", "Rate")]), ncol = 2)
Y_list <- survival_encode(event = subdat[,'death'], time = subdat[,'time'])
Y_mat <- Y_list$label
M_mat <- Y_list$mask

num.iteration <- 2000 
lr <- 0.01
B_mat <- matrix(0, nrow = ncol(X_mat), ncol = 1)

for (i in 1:num.iteration) {
  
  H_mat <- matrix(X_mat %*% B_mat, nrow = nrow(Y_mat), ncol = ncol(Y_mat), byrow = TRUE)

  P11_mat <- M_mat * exp(H_mat)
  P12_mat <- P11_mat %*% X_mat
  P13_mat <- P11_mat %*% matrix(1, nrow = nrow(X_mat), ncol = ncol(X_mat))
  P14_mat <- apply(P12_mat / P13_mat, 2, sum, na.rm = TRUE)
  
  P21_mat <- M_mat * Y_mat * exp(H_mat)
  P22_mat <- P21_mat %*% X_mat
  P23_mat <- P21_mat %*% matrix(1, nrow = nrow(X_mat), ncol = ncol(X_mat))
  P24_mat <- apply(P22_mat / P23_mat, 2, sum, na.rm = TRUE)

  grad_B <- P14_mat - P24_mat
  
  B_mat <- B_mat - lr * grad_B / nrow(X_mat)
  
}

print(B_mat)
##             [,1]
## [1,] 0.037072488
## [2,] 0.008077607

練習1答案(1)

– 這是使用L1 Regularization的結果:

lambda <- 0.005
num.iteration <- 2000 
lr <- 0.01
B_mat <- matrix(0, nrow = ncol(X_mat), ncol = 1)

for (i in 1:num.iteration) {
  
  H_mat <- matrix(X_mat %*% B_mat, nrow = nrow(Y_mat), ncol = ncol(Y_mat), byrow = TRUE)

  P11_mat <- M_mat * exp(H_mat)
  P12_mat <- P11_mat %*% X_mat
  P13_mat <- P11_mat %*% matrix(1, nrow = nrow(X_mat), ncol = ncol(X_mat))
  P14_mat <- apply(P12_mat / P13_mat, 2, sum, na.rm = TRUE)
  
  P21_mat <- M_mat * Y_mat * exp(H_mat)
  P22_mat <- P21_mat %*% X_mat
  P23_mat <- P21_mat %*% matrix(1, nrow = nrow(X_mat), ncol = ncol(X_mat))
  P24_mat <- apply(P22_mat / P23_mat, 2, sum, na.rm = TRUE)

  grad_B <- P14_mat - P24_mat
  
  B_mat <- B_mat - lr * grad_B / nrow(X_mat) - lambda * sign(B_mat)
  
}

print(B_mat)
##             [,1]
## [1,] 0.023285025
## [2,] 0.003745863

– 使用L1 Regularization的結果告訴我們其實Rate這個變項不重要!

練習1答案(2)

lambda <- 0.05
num.iteration <- 2000 
lr <- 0.01
B_mat <- matrix(0, nrow = ncol(X_mat), ncol = 1)

for (i in 1:num.iteration) {
  
  H_mat <- matrix(X_mat %*% B_mat, nrow = nrow(Y_mat), ncol = ncol(Y_mat), byrow = TRUE)

  P11_mat <- M_mat * exp(H_mat)
  P12_mat <- P11_mat %*% X_mat
  P13_mat <- P11_mat %*% matrix(1, nrow = nrow(X_mat), ncol = ncol(X_mat))
  P14_mat <- apply(P12_mat / P13_mat, 2, sum, na.rm = TRUE)
  
  P21_mat <- M_mat * Y_mat * exp(H_mat)
  P22_mat <- P21_mat %*% X_mat
  P23_mat <- P21_mat %*% matrix(1, nrow = nrow(X_mat), ncol = ncol(X_mat))
  P24_mat <- apply(P22_mat / P23_mat, 2, sum, na.rm = TRUE)

  grad_B <- P14_mat - P24_mat
  
  B_mat <- B_mat - lr * grad_B / nrow(X_mat) - lambda * B_mat
  
}

print(B_mat)
##             [,1]
## [1,] 0.031659045
## [2,] 0.008512594

– 你可以試試看使用不同的\(\lambda\)

第二節:彈性網路(1)

– 答案當然是可以的,先讓我們以視覺化的方式觀察他們分別的意義。

– 以剛剛的資料為例,讓我們看看「主損失函數」、「L1 Regularization」、「L2 Regularization」他們分別的損失梯度圖,其中白色的點是在無正則化下的最優解。

– 讓我們將它合併,左邊是全部使用L1 Regularization、中間是全部使用L2 Regularization、右邊是各用一半: