局部極值問題(1)

梯度下降法存在非常多問題，除了最適學習率不好找之外，他並非在所有狀況下都能找到loss最小的位置。

我們一樣先從簡單的問題開始，我們希望求\(y = 2x^2 +3x^3 +x^4\)的極值，但我們先看看函數圖形：

F3_1

我們看的出來，這個函數有2個區域極小值，其中一個是全局極值，我們試著用梯度下降法來求解吧：

differential.fun = function(x) {
  return(4*x + 9*x^2 + 4*x^3)
}

start.value = 0.5 
num.iteration = 30
learning.rate = 0.05

result.x = rep(NA, num.iteration)

for (i in 1:num.iteration) {
  if (i == 1) {
    result.x[1] = start.value
  } else {
    result.x[i] = result.x[i-1] - learning.rate * differential.fun(result.x[i-1])
  }
}

tail(result.x, 1)

## [1] 0.000275395

我們可以發現，當我們的起始值設定在0.5，則會找到右邊的極值出現點：0

start.value = -3 
num.iteration = 30
learning.rate = 0.05

result.x = rep(NA, num.iteration)

for (i in 1:num.iteration) {
  if (i == 1) {
    result.x[1] = start.value
  } else {
    result.x[i] = result.x[i-1] - learning.rate * differential.fun(result.x[i-1])
  }
}

tail(result.x, 1)

## [1] -1.640338

如果我們的起始值設定在-3，則會找到左邊的極值出現點：-1.64

局部極值問題(2)

看來起始值設定在哪直接影響到了結果，你是否害怕線性迴歸、邏輯斯迴歸的結果呢?

– 你可以親自在線性迴歸/邏輯斯迴歸中設定不同的起始值，你會發現線性迴歸、邏輯斯迴歸都沒有這個問題，這是因為他們的函數圖形都比較像\(y = x^2\)這種凸函數

– 而\(y = 2x^2 +3x^3 +x^4\)這種叫做凹凸函數

所以我們了解到梯度下降法這種近似解求法在凹凸函數上，通常僅能求出局部最小值，而局部最小值與求解起始值有相當大的關係。
其實從預測未知結果的角度，局部最佳解所求得的方程式儘管不是預測能力最好的方程式，但仍然有一定的預測能力，所以仍然有其實用性。

– 當然，我們仍然希望盡可能的避免局部極值，統計/電腦學家在這方面下了非常多的努力，而我們就即將介紹一些梯度下降法的變種試圖避免我們的優化器陷入局部極值。

慣性動量(1)

在真實世界中，我們若沿著斜率往下走，會受到慣性定律影響讓我們保留部分當前的速度

F3_2

這個想法激發了我們重新改寫梯度下降法的公式，我們可以為他加入一個慣性動量：

\[x_{\left(epoch:t\right)} = x_{\left(epoch:t - 1\right)} - lr \cdot \frac{\partial}{\partial x}f(x_{\left(epoch:t - 1\right)}) + \rho \cdot \left[x_{\left(epoch:t - 1\right)}-x_{\left(epoch:t - 2\right)}\right]\]

這個式子是甚麼意思呢?那就是我們每次下降的值除了受到當前梯度影響之外，還會受到上一次的慣性力作用。而\(\rho\)是慣性動量的保留係數，原則上是一個介於0到1之間的數值，而最常用的是0.9。

慣性動量(2)

讓我們來試試看如何用這個新公式來找尋極值，我們一樣來找尋\(y = 2x^2 +3x^3 +x^4\)的極值

differential.fun = function(x) {
  return(4*x + 9*x^2 + 4*x^3)
}

start.value = 0.5 
num.iteration = 30
learning.rate = 0.1
momentum = 0.9

result.x = rep(NA, num.iteration)

for (i in 1:num.iteration) {
  if (i == 1) {
    result.x[1] = start.value
  } else if (i == 2) {
    result.x[i] = result.x[i-1] - learning.rate * differential.fun(result.x[i-1])
  } else {
    result.x[i] = result.x[i-1] - learning.rate * differential.fun(result.x[i-1]) + momentum * (result.x[i-1] - result.x[i-2])
  }
}

tail(result.x, 1)

## [1] -1.602265

這下子厲害了，藉由慣性動量的協助下，即使我們起始值設定在0.5，仍然衝過了局部最小值，從而到達局部最大值。

F3_3

慣性動量(3)

慣性動量的好處還不只這些，在第一節課時我們曾經示範當學習律不恰當的時候，會導至結果無法收斂，而慣性動量也能在這種情況下發揮妙用。

– 這個例子是求\(y=x^2\)的極值，但學習率錯誤訂成了1.02，在正常狀況下是無法收斂的，但具有慣性動量將可以收斂。

differential.fun = function(x) {
  return(2*x)
}

start.value = 0.5 
num.iteration = 30
learning.rate = 1.02
momentum = 0.9

result.x = rep(NA, num.iteration)

for (i in 1:num.iteration) {
  if (i == 1) {
    result.x[1] = start.value
  } else if (i == 2) {
    result.x[i] = result.x[i-1] - learning.rate * differential.fun(result.x[i-1])
  } else {
    result.x[i] = result.x[i-1] - learning.rate * differential.fun(result.x[i-1]) + momentum * (result.x[i-1] - result.x[i-2])
  }
}

tail(result.x, 1)

## [1] -0.03117517

F3_4

練習1：變種邏輯斯回歸的優化

假定我們非常清楚x1與x2與y之間的關係如下：

\[log(\frac{{p}}{1-p}) = sin(b_{1}x_1x_2+b_2x_2)\]

下面這張圖是\(b_1 = 0.7\)且\(b_2 = 0.7\)的情形，紅色的部分代表\(y\)確定為0，而藍色的部分代表\(y\)確定為1：

現在讓我們試著透過下面這個樣本解出\(b_1\)且\(b_2\)的數值：

set.seed(0)
x1 = runif(1000, min = -3, max = 3) 
x2 = runif(1000, min = -3, max = 3) 
lr = sin(0.7*x1*x2+0.7*x2) + rnorm(1000, sd = 0.3)
y = lr > 0 + 0L

假定我們使用交叉熵作為這個邏輯斯回歸的損失函數，而整體的函數關係式如下：

\[ \begin{align} L(x_1,x_2) & = sin(b_{1}x_1x_2+b_2x_2) \\ S(x) & = \frac{{1}}{1+e^{-x}} \\ CE(y, p) & = -\left(y_{i} \cdot log(p_{i}) + (1-y_{i}) \cdot log(1-p_{i})\right) \\ loss & = CE(y, S(L(x_1,x_2))) \end{align} \]

– \(b_1\)的偏導函數(過程略)：

\[ \begin{align} \frac{\partial}{\partial b_1}loss_i & = (p_i - y_i) \cdot cos(b_{1}x_1x_2+b_2x_2) \cdot x_{1i}x_{2i} \end{align} \]

– \(b_2\)的偏導函數(過程略)：

\[ \begin{align} \frac{\partial}{\partial b_2}loss_i & = (p_i - y_i) \cdot cos(b_{1}x_1x_2+b_2x_2) \cdot x_{2i} \end{align} \]

程式碼如下，請你試著在優化時增加慣性量以及不增加慣性量的優化結果：

pred.fun = function(b1, b2, x1 = x1, x2 = x2) {
  lr = sin(b1*x1*x2+b2*x2)
  p = 1 / (1 + exp(-lr))
  return(p)
}

loss.fun = function(p, y = y, eps = 1e-8) {
  loss = -mean(y * log(p + eps) + (1 - y) * log(1 - p + eps))
  return(loss)
}

differential.fun.b1 = function(p, b1, b2, x1 = x1, x2 = x2, y = y) {
  return(mean((p-y)*cos(b1*x1*x2+b2*x2)*x1*x2))
}

differential.fun.b2 = function(p, b1, b2, x1 = x1, x2 = x2, y = y) {
  return(mean((p-y)*cos(b1*x1*x2+b2*x2)*x2))
}

練習1答案(程式碼)

下面是這個題目的視覺化函數，這能協助我們了解其學習軌跡：

LOGISTIC_VIS = function (x1, x2, y, num.iteration = 100, lr = 0.1, momentum = 0.9, start_b1 = 0, start_b2 = 0) {

  ans_b1 = rep(start_b1, num.iteration)
  ans_b2 = rep(start_b2, num.iteration)
  ans_loss = rep(0, num.iteration)
  
  pred.fun = function(b1, b2, x1 = x1, x2 = x2) {
    lr = sin(b1*x1*x2+b2*x2)
    p = 1 / (1 + exp(-lr))
    return(p)
  }
  
  loss.fun = function(p, y = y, eps = 1e-8) {
    loss = -mean(y * log(p + eps) + (1 - y) * log(1 - p + eps))
    return(loss)
  }
  
  differential.fun.b1 = function(p, b1, b2, x1 = x1, x2 = x2, y = y) {
    return(mean((p-y)*cos(b1*x1*x2+b2*x2)*x1*x2))
  }
  
  differential.fun.b2 = function(p, b1, b2, x1 = x1, x2 = x2, y = y) {
    return(mean((p-y)*cos(b1*x1*x2+b2*x2)*x2))
  }
  
  for (i in 1:num.iteration) {
    if (i >= 2) {
      grad_b1 = differential.fun.b1(p = current_p, b1 = ans_b1[i-1], b2 = ans_b2[i-1], x1 = x1, x2 = x2, y = y)
      grad_b2 = differential.fun.b2(p = current_p, b1 = ans_b1[i-1], b2 = ans_b2[i-1], x1 = x1, x2 = x2, y = y)
      if (i == 2) {
        momentum_b1 = 0
        momentum_b2 = 0
      } else {
        momentum_b1 = momentum * (ans_b1[i-1] - ans_b1[i-2])
        momentum_b2 = momentum * (ans_b2[i-1] - ans_b2[i-2]) 
      }
      ans_b1[i] = ans_b1[i-1] - lr * grad_b1 + momentum_b1
      ans_b2[i] = ans_b2[i-1] - lr * grad_b2 + momentum_b2
    }
    current_p = pred.fun(b1 = ans_b1[i], b2 = ans_b2[i], x1 = x1, x2 = x2)
    ans_loss[i] = loss.fun(p = current_p, y = y)
  }
  
  par(mfrow = c(1, 2))
  
  require(scales)
  require(plot3D)
  
  x1_seq = seq(min(x1), max(x1), length.out = 300)
  x2_seq = seq(min(x2), max(x2), length.out = 300)
  
  z_matrix = sapply(x2_seq, function(x) {pred.fun(b1 = ans_b1[num.iteration], b2 = ans_b2[num.iteration], x1 = x1_seq, x2 = x)})
  
  image2D(z = z_matrix, main = paste0('b1=', formatC(ans_b1[num.iteration], format = 'f', 1), ';b2=', formatC(ans_b2[num.iteration], format = 'f', 1)),
          x = x1_seq, xlab = 'x1',
          y = x2_seq, ylab = 'x2',
          shade = 0.2, rasterImage = TRUE,
          col = colorRampPalette(c("#FFA0A0", "#FFFFFF", "#A0A0FF"))(100))
  
  points(x1, x2, col = (y + 1)*2, pch = 19, cex = 0.5)
  
  require(grDevices) # for trans3d
  
  b1_seq = seq(-0.5, 2, length.out = 50)
  b2_seq = seq(-0.5, 2, length.out = 50)
  LOSS <- b1_seq %*% t(b2_seq)
  for (i in 1:length(b1_seq)) {
    for (j in 1:length(b2_seq)) {
      LOSS[i,j] = loss.fun(p = pred.fun(b1 = b1_seq[i], b2 = b2_seq[j], x1 = x1, x2 = x2), y = y, eps = 1e-8)
    }
  }
  
  persp(b1_seq, b2_seq, LOSS, theta = 120, phi = 50, expand = 0.8, col = "lightblue",
        ltheta = 120, shade = 0.75, ticktype = "detailed",
        xlab = "b1", ylab = "b2", zlab = "loss", main = paste0('loss=', formatC(ans_loss[num.iteration], 4, format = 'f'))) -> res
  
  lines(trans3d(ans_b1, ans_b2, ans_loss, pmat = res), col = '#AA0000', lwd = 1.5)
  points(trans3d(ans_b1[1], ans_b2[1], ans_loss[1], pmat = res), col = '#AA0000', cex = 1, pch = 16)
  points(trans3d(ans_b1[num.iteration], ans_b2[num.iteration], ans_loss[num.iteration], pmat = res), col = '#FF0000', cex = 1, pch = 16)
  
}

練習1答案(實驗結果)

這是當動量為0的結果：

LOGISTIC_VIS(x1 = x1, x2 = x2, y = y, num.iteration = 100, lr = 0.1, momentum = 0, start_b1 = 0, start_b2 = 0)

這是當動量為0.9的結果：

LOGISTIC_VIS(x1 = x1, x2 = x2, y = y, num.iteration = 100, lr = 0.1, momentum = 0.9, start_b1 = 0, start_b2 = 0)

練習1之引申問題

在剛剛的實驗中，我們發現了慣性動量能協助我們衝破局部極值，但如果函數圖形更加複雜呢?事實上剛剛的例子中\(b_1 = 0.7\)且\(b_2 = 0.7\)的情形是一個精心設計的特例，我們試試看產生數據時把他設定成\(b_1 = 2\)且\(b_2 = 3\)再看看結果
這是當動量為0.9的結果：

set.seed(0)
x1 = runif(1000, min = -3, max = 3) 
x2 = runif(1000, min = -3, max = 3) 
lr = sin(2*x1*x2+3*x2) + rnorm(1000, sd = 0.3)
y = lr > 0 + 0L
LOGISTIC_VIS(x1 = x1, x2 = x2, y = y, num.iteration = 300, lr = 0.1, momentum = 0.9, start_b1 = 0, start_b2 = 0)

由於函數圖形的凹凸程度更加複雜，這導致即使帶著慣性動量都無法衝破局部極值，因此我們一定還需要更多手段來解決局部極值問題！

小批量梯度下降法(1)

在我們之前的例子中，每次進行權重更新時都使用了所有的樣本，但如今真實問題訓練數據都是非常巨大，因此現在的常規做法是在每次運算時只使用小部分的樣本。
從『理論上』，只要我們選擇樣本沒有特別偏好，那麼小的隨機樣本提供梯度的期望值應該與全樣本相同，所以結果應該一樣。

– 然而『實際上』，由於每次抽樣造成的變異，這會使的當前梯度會因為抽樣的關係存在於一定的變異性，這將會導致使用「小批量梯度下降法(Mini-batch Stochastic Gradient Descent, SGD)」與「全樣本梯度下降法」的結果有所不同。但問題在於，你認為「小批量梯度下降法」與「全樣本梯度下降法」哪一種能得到較好的結果?而如果使用「小批量梯度下降法」，那批量「樣本大小」與優化結果的關係又是什麼呢?

– 下面是利用練習1的問題所描繪的損失函數，這6張圖分別是批量大小從小到大其最終損失函數的變異情形，試著觀察並討論此問題：

F3_5

小批量梯度下降法(2)

事實上你應該發現了一件事情，那就是由於每次抽樣所造成損失函數的變異，從而導致原來的局部極值似乎偶爾會消失。而批量樣本越大，局部極值就越不容易消失，所以批量樣本過大似乎對優化不利！(這是否違反你的統計直覺?)

– 當然，過小的批量樣本似乎也不是非常理想，這同時對全局極值也存在著危害，同樣有可能讓優化器錯過全局極值。但整體來說，似乎批量樣本較小較為有利！

讓我們實際做實驗來看看在小批量之下的優化情形吧！

– 這個例子是這是\(b_1 = 0.7\)且\(b_2 = 0.7\)所產生的數據，這在練習1中若設為動量為0，則無法離開局部極值，而我們改用批量大小為30再次試試：

F3_6

– 下面這個例子是這是\(b_1 = 2\)且\(b_2 = 3\)所產生的數據，這在練習1中即使我們使用了動量為0.9的優化過程仍然沒辦法找到答案，但是透過調整批量大小(批量大小 = 30)後，我們將能找到全局極值：

F3_7

小批量梯度下降法(3)

關於批量大小與優化效果的關係也有paper在探討，而著名的Yann LeCun也曾對此發表評論：

F3_8

事實上他是引用Revisiting Small Batch Training for Deep Neural Networks的結果，這篇研究是透過大量的實驗探討「批量大小」與「測試集準確度」之間的關係。

– 最終發現在幾個常見的演算法及資料集中，「批量大小」不超過32時「測試集準確度」較高。

F3_9

你是否為「隨機小批量」優化的威力感到驚訝呢?

練習2：調整批量樣本

在第二課我們發現在多數情形下，隱藏層的神經元給太多有時候會對外推性產生危害，現在讓我們來試試看慣性量與小批量的協助下是否能減少這樣的問題：

set.seed(0)
x1 = rnorm(100, sd = 1) 
x2 = rnorm(100, sd = 1) 
lr1 = - 1.5 + x1^2 + x2^2 + rnorm(100)
y = lr1 > 0 + 0L

test_x1 = rnorm(100, sd = 1) 
test_x2 = rnorm(100, sd = 1) 
lr1 = - 1.5 + test_x1^2 + test_x2^2 + rnorm(100)
test_y = lr1 > 0 + 0L

第二課練習4的程式碼如下，請你把它改成帶有慣性量並且可以指定批量大小，並且實驗看看小的批量大小是否能減少過度擬合的危害：

MLP_Trainer = function (num.iteration = 500, num.hidden = 30, lr = 0.1, x1 = x1, x2 = x2, y = y, test_x1 = NULL, test_x2 = NULL, test_y = NULL) {
  
  #Functions
  
  #Forward
  
  S.fun = function (x, eps = 1e-8) {
    S = 1/(1 + exp(-x))
    S[S < eps] = eps
    S[S > 1 - eps] = 1 - eps
    return(S)
  }
  
  ReLU.fun = function (x) {
    x[x < 0] <- 0
    return(x)
  }
  
  L.fun = function (X, W) {
    X.E = cbind(1, X)
    L = X.E %*% W
    return(L)
  }
  
  CE.fun = function (o, y, eps = 1e-8) {
    loss = -1/length(y) * sum(y * log(o + eps) + (1 - y) * log(1 - o + eps))
    return(loss)
  }
  
  #Backward
  
  grad_o.fun = function (o, y) {
    return((o - y)/(o*(1-o)))
  }
  
  grad_l2.fun = function (grad_o, o) {
    return(grad_o*(o*(1-o)))
  }
  
  grad_W2.fun = function (grad_l2, h1) {
    h1.E = cbind(1, h1)
    return(t(h1.E) %*% grad_l2/nrow(h1))
  }
  
  grad_h1.fun = function (grad_l2, W2) {
    return(grad_l2 %*% t(W2[-1,]))
  }
  
  grad_l1.fun = function (grad_h1, l1) {
    de_l1 = l1
    de_l1[de_l1<0] = 0
    de_l1[de_l1>0] = 1
    return(grad_h1*de_l1)
  }
  
  grad_W1.fun = function (grad_l1, x) {
    x.E = cbind(1, x)
    return(t(x.E) %*% grad_l1/nrow(x))
  }
  
  #Caculating
  
  X_matrix = cbind(x1, x2)
  
  W1_list = list()
  W2_list = list()
  loss_seq = rep(0, num.iteration)
  
  #Start random values
  
  W1_list[[1]] = matrix(rnorm(3*num.hidden, sd = 1), nrow = 3, ncol = num.hidden)
  W2_list[[1]] = matrix(rnorm(num.hidden + 1, sd = 1), nrow = num.hidden + 1, ncol = 1)
  
  for (i in 2:(num.iteration+1)) {
    
    #Forward
    
    current_l1 = L.fun(X = X_matrix, W = W1_list[[i - 1]])
    current_h1 = ReLU.fun(x = current_l1)
    current_l2 = L.fun(X = current_h1, W = W2_list[[i - 1]])
    current_o = S.fun(x = current_l2)
    loss_seq[i-1] = CE.fun(o = current_o, y = y, eps = 1e-8)
    
    #Backward
    
    current_grad_o = grad_o.fun(o = current_o, y = y)
    current_grad_l2 = grad_l2.fun(grad_o = current_grad_o, o = current_o)
    current_grad_W2 = grad_W2.fun(grad_l2 = current_grad_l2, h1 = current_h1)
    current_grad_h1 = grad_h1.fun(grad_l2 = current_grad_l2, W2 = W2_list[[i - 1]])
    current_grad_l1 = grad_l1.fun(grad_h1 = current_grad_h1, l1 = current_l1)
    current_grad_W1 = grad_W1.fun(grad_l1 = current_grad_l1, x = X_matrix)
    W2_list[[i]] = W2_list[[i-1]] - lr * current_grad_W2
    W1_list[[i]] = W1_list[[i-1]] - lr * current_grad_W1
    
  }
  
  require(scales)
  require(plot3D)
  
  x1_seq = seq(min(x1), max(x1), length.out = 100)
  x2_seq = seq(min(x2), max(x2), length.out = 100)
  
  pre_func = function (x1, x2, W1 = W1_list[[length(W1_list)]], W2 = W2_list[[length(W2_list)]]) {
    new_X = cbind(x1, x2)
    O = S.fun(x = L.fun(X = ReLU.fun(x = L.fun(X = new_X, W = W1)), W = W2))
    return(O)
  }
  
  pred_y = pre_func(x1 = x1, x2 = x2)
  MAIN_TXT = paste0('Train-Acc:', formatC(mean((pred_y > 0.5) == y), 2, format = 'f'))
  if (!is.null(test_x1)) {
    pred_test_y = pre_func(x1 = test_x1, x2 = test_x2)
    MAIN_TXT = paste0(MAIN_TXT, '; Test-Acc:', formatC(mean((pred_test_y > 0.5) == test_y), 2, format = 'f'))
  }
  
  z_matrix = sapply(x2_seq, function(x) {pre_func(x1 = x1_seq, x2 = x)})
  
  par(mfrow = c(1, 2))
  
  image2D(z = z_matrix, main = MAIN_TXT,
          x = x1_seq, xlab = 'x1',
          y = x2_seq, ylab = 'x2',
          shade = 0.2, rasterImage = TRUE,
          col = colorRampPalette(c("#FFA0A0", "#FFFFFF", "#A0A0FF"))(100))
  
  points(x1, x2, col = (y + 1)*2, pch = 19, cex = 0.5)
  if (!is.null(test_x1)) {
    points(test_x1, test_x2, col = 'black', bg = c('#C00000', '#0000C0')[(test_y + 1)], pch = 21)
  }
  
  plot(loss_seq, type = 'l', main = 'loss', xlab = 'iter.', ylab = 'CE loss')
  
}

練習2答案(程式碼)

修改後的程式碼如下：

MLP_Trainer = function (num.iteration = 500, num.hidden = 30, batch_size = 30, lr = 0.1, momentum = 0, x1 = x1, x2 = x2, y = y, test_x1 = NULL, test_x2 = NULL, test_y = NULL) {
  
  #Functions
  
  #Forward
  
  S.fun = function (x, eps = 1e-8) {
    S = 1/(1 + exp(-x))
    S[S < eps] = eps
    S[S > 1 - eps] = 1 - eps
    return(S)
  }
  
  ReLU.fun = function (x) {
    x[x < 0] <- 0
    return(x)
  }
  
  L.fun = function (X, W) {
    X.E = cbind(1, X)
    L = X.E %*% W
    return(L)
  }
  
  CE.fun = function (o, y, eps = 1e-8) {
    loss = -1/length(y) * sum(y * log(o + eps) + (1 - y) * log(1 - o + eps))
    return(loss)
  }
  
  #Backward
  
  grad_o.fun = function (o, y) {
    return((o - y)/(o*(1-o)))
  }
  
  grad_l2.fun = function (grad_o, o) {
    return(grad_o*(o*(1-o)))
  }
  
  grad_W2.fun = function (grad_l2, h1) {
    h1.E = cbind(1, h1)
    return(t(h1.E) %*% grad_l2/nrow(h1))
  }
  
  grad_h1.fun = function (grad_l2, W2) {
    return(grad_l2 %*% t(W2[-1,]))
  }
  
  grad_l1.fun = function (grad_h1, l1) {
    de_l1 = l1
    de_l1[de_l1<0] = 0
    de_l1[de_l1>0] = 1
    return(grad_h1*de_l1)
  }
  
  grad_W1.fun = function (grad_l1, x) {
    x.E = cbind(1, x)
    return(t(x.E) %*% grad_l1/nrow(x))
  }
  
  #Caculating
  
  X_matrix = cbind(x1, x2)
  
  W1_list = list()
  W2_list = list()
  loss_seq = rep(0, num.iteration)
  
  #Start random values
  
  W1_list[[1]] = matrix(rnorm(3*num.hidden, sd = 1), nrow = 3, ncol = num.hidden)
  W2_list[[1]] = matrix(rnorm(num.hidden + 1, sd = 1), nrow = num.hidden + 1, ncol = 1)
  
  for (i in 2:(num.iteration+1)) {
    
    idx <- sample(1:length(x1), batch_size)
    
    #Forward
    
    current_l1 = L.fun(X = X_matrix[idx,], W = W1_list[[i - 1]])
    current_h1 = ReLU.fun(x = current_l1)
    current_l2 = L.fun(X = current_h1, W = W2_list[[i - 1]])
    current_o = S.fun(x = current_l2)
    loss_seq[i-1] = CE.fun(o = current_o, y = y[idx], eps = 1e-8)
    
    #Backward
    
    current_grad_o = grad_o.fun(o = current_o, y = y[idx])
    current_grad_l2 = grad_l2.fun(grad_o = current_grad_o, o = current_o)
    current_grad_W2 = grad_W2.fun(grad_l2 = current_grad_l2, h1 = current_h1)
    current_grad_h1 = grad_h1.fun(grad_l2 = current_grad_l2, W2 = W2_list[[i - 1]])
    current_grad_l1 = grad_l1.fun(grad_h1 = current_grad_h1, l1 = current_l1)
    current_grad_W1 = grad_W1.fun(grad_l1 = current_grad_l1, x = X_matrix[idx,])
    
    if (i == 2) {
      W2_list[[i]] = W2_list[[i-1]] - lr * current_grad_W2
      W1_list[[i]] = W1_list[[i-1]] - lr * current_grad_W1
    } else {
      W2_list[[i]] = W2_list[[i-1]] - lr * current_grad_W2 + momentum * (W2_list[[i-1]] - W2_list[[i-2]])
      W1_list[[i]] = W1_list[[i-1]] - lr * current_grad_W1 + momentum * (W1_list[[i-1]] - W1_list[[i-2]])
    }
    

    
  }
  
  require(scales)
  require(plot3D)
  
  x1_seq = seq(min(x1), max(x1), length.out = 100)
  x2_seq = seq(min(x2), max(x2), length.out = 100)
  
  pre_func = function (x1, x2, W1 = W1_list[[length(W1_list)]], W2 = W2_list[[length(W2_list)]]) {
    new_X = cbind(x1, x2)
    O = S.fun(x = L.fun(X = ReLU.fun(x = L.fun(X = new_X, W = W1)), W = W2))
    return(O)
  }
  
  pred_y = pre_func(x1 = x1, x2 = x2)
  MAIN_TXT = paste0('Train-Acc:', formatC(mean((pred_y > 0.5) == y), 2, format = 'f'))
  if (!is.null(test_x1)) {
    pred_test_y = pre_func(x1 = test_x1, x2 = test_x2)
    MAIN_TXT = paste0(MAIN_TXT, '; Test-Acc:', formatC(mean((pred_test_y > 0.5) == test_y), 2, format = 'f'))
  }
  
  z_matrix = sapply(x2_seq, function(x) {pre_func(x1 = x1_seq, x2 = x)})
  
  par(mfrow = c(1, 2))
  
  image2D(z = z_matrix, main = MAIN_TXT,
          x = x1_seq, xlab = 'x1',
          y = x2_seq, ylab = 'x2',
          shade = 0.2, rasterImage = TRUE,
          col = colorRampPalette(c("#FFA0A0", "#FFFFFF", "#A0A0FF"))(100))
  
  points(x1, x2, col = (y + 1)*2, pch = 19, cex = 0.5)
  if (!is.null(test_x1)) {
    points(test_x1, test_x2, col = 'black', bg = c('#C00000', '#0000C0')[(test_y + 1)], pch = 21)
  }
  
  plot(loss_seq, type = 'l', main = 'loss', xlab = 'iter.', ylab = 'CE loss')
  
}

練習2答案(實驗結果)

這是批量大小為全部的狀況

MLP_Trainer(num.iteration = 5000, num.hidden = 30, batch_size = 100, lr = 0.1, momentum = 0.9, x1 = x1, x2 = x2, y = y, test_x1 = test_x1, test_x2 = test_x2, test_y = test_y)

這是批量大小為20的狀況

MLP_Trainer(num.iteration = 5000, num.hidden = 30, batch_size = 20, lr = 0.1, momentum = 0.9, x1 = x1, x2 = x2, y = y, test_x1 = test_x1, test_x2 = test_x2, test_y = test_y)

每次批量用太多樣本確實危害健康吧！思考一下為什麼會這樣。

自適應學習率優化演算法(1)

在剛剛的實驗中我們已經證實了小批量的優化過程較佳，但實驗後期\(loss\)的軌跡讓我們相當害怕，是否有方法讓實驗末期的參數更新趨於穩定?

– 一個合理的方法就是透過調整學習率，在疊代初期給一個較大的學習率，而在後期給一個較小的學習率。

– 另外目前為止，所有參數的更新「共享」了同一個學習率，有沒有可能針對每一個參數設計自己的學習率，從而導致優化過程更好?

這邊我們先介紹一個比較早期的解法，叫做Adagrad，他是由John Duchi et al., 2011所發展，對於\(x\)的更新公式如下：

\[ \begin{align} grad.x_t & = \frac{\partial}{\partial x}f(x_{t - 1}) \\ n_t & = n_{t-1} + grad.x_t^2 \\ x_{t} & = x_{t - 1} - \frac{lr}{\sqrt{n_t+\delta}}\cdot grad.x_t \end{align} \]

在這邊\(n\)所代表的是累積梯度平方和，而\(\delta\)是一個很小的數字避免除以0，而這個方法中的\(n\)為一切的核心，他對於梯度較小的參數會給予較大的學習率，並且會隨著時間讓學習率越來越低，從而讓結果收斂。

– 在R語言中的實現如下，我們以求\(y=x^2\)的極值為例實現下列程式碼：

differential.fun = function(x) {
  return(2*x)
}

start.value = 0.5 
num.iteration = 30
learning.rate = 0.1
eps = 1e-8

result.x = rep(NA, num.iteration)

for (i in 1:num.iteration) {
  if (i == 1) {
    n = 0
    result.x[1] = start.value
  } else {
    grad.x = differential.fun(result.x[i-1])
    n = n + grad.x^2
    result.x[i] = result.x[i-1] - learning.rate/sqrt(n + eps) * grad.x
  }
}

tail(result.x, 1)

## [1] 0.01459057

自適應學習率優化演算法(2)

Adagrad提出的「梯度平方項」概念是相當好的，但仔細觀察Adagrad會發現他的公式有一些問題，隨著累積梯度平方和的增加，那整體學習率將會非常非常靠近0，使訓練提前結束。

– 這裡我們可以對其做一些簡單的擴展，主要的概念是來自於Adadelta，我們可以在累積梯度平方和的部分利用了動量的概念，讓當前累積值會隨著時間做衰減，從而保證訓練不會提前結束

這裡我們給定一個衰減係數\(\beta_2\)介於0至1，我們通常設為0.999，擴展後的公式如下：

\[ \begin{align} grad.x_t & = \frac{\partial}{\partial x}f(x_{t - 1}) \\ n_t & = \beta_2 \cdot n_{t-1} + (1 - \beta_2) \cdot grad.x_t^2 \\ x_{t} & = x_{t - 1} - \frac{lr}{\sqrt{n_t+\delta}}\cdot grad.x_t \end{align} \]

接著，你是否注意到這個優化方式與前面學的SGD相比，似乎忽略了動量因子，因此我們再把式子做一些擴展(\(m_t\)代表當前累積動量，而\(\beta_1\)與之前的動量參數\(\rho\)完全等價，我們通常設為0.9)：

\[ \begin{align} grad.x_t & = \frac{\partial}{\partial x}f(x_{t - 1}) \\ m_t & = \beta_1 \cdot m_{t-1} + (1 - \beta_1) \cdot grad.x_t \\ n_t & = \beta_2 \cdot n_{t-1} + (1 - \beta_2) \cdot grad.x_t^2 \\ x_{t} & = x_{t - 1} - \frac{lr}{\sqrt{n_t+\delta}}\cdot m_t \end{align} \]

– 擴展自此，我們考慮最後一個問題，那就是由於\(m_t\)與\(n_t\)是隨著時間累積的，在訓練剛開始的時候這個數字將會非常小，所以我們將計算\(\hat{m_t}\)與\(\hat{n_t}\)作為更新梯度使用，而他的分母項將會隨著時間越來越接近1：

\[ \begin{align} grad.x_t & = \frac{\partial}{\partial x}f(x_{t - 1}) \\ m_t & = \beta_1 \cdot m_{t-1} + (1 - \beta_1) \cdot grad.x_t \\ n_t & = \beta_2 \cdot n_{t-1} + (1 - \beta_2) \cdot grad.x_t^2 \\ \hat{m_t} & = \frac{m_t}{1 - \beta_1^t} \\ \hat{n_t} & = \frac{n_t}{1 - \beta_2^t} \\ x_{t} & = x_{t - 1} - \frac{lr}{\sqrt{\hat{n_t}+\delta}}\cdot \hat{m_t} \end{align} \]

自適應學習率優化演算法(3)

這個結合了動量因子與學習率遞減的方法被稱作Adam(adaptive moment estimation)，他是由Diederik P. Kingma & Jimmy Lei Ba, 2014所發展，統計至2018年6月已被引用近9000次，可謂當代最流行的優化演算法。

F3_10

下面是我們在R語言中的實現，我們同樣以求\(y=x^2\)的極值為例：

differential.fun = function(x) {
  return(2*x)
}

start.value = 0.5 
num.iteration = 30
learning.rate = 0.1 #注意，一般Adam在優化時通常選擇0.001
beta1 = 0.9
beta2 = 0.999
eps = 1e-8

result.x = rep(NA, num.iteration + 1)

for (i in 0:num.iteration) {
  if (i == 0) {
    m = 0
    n = 0
    result.x[i+1] = start.value
  } else {
    grad.x = differential.fun(result.x[i])
    m = beta1 * m + (1 - beta1) * grad.x
    n = beta2 * n + (1 - beta2) * grad.x^2
    m.hat = m/(1 - beta1^i)
    n.hat = n/(1 - beta2^i)
    result.x[i+1] = result.x[i] - learning.rate/sqrt(n.hat + eps) * m.hat
  }
}

tail(result.x, 1)

## [1] 0.02099521

自適應學習率優化演算法(4)

我們以練習1的例子為例並使用\(b_1 = 2\)且\(b_2 = 3\)產生數據，讓我們比較帶有動量的SGD與Adam在優化時的差異：

– SGD的結果(\(lr = 0.1\), \(\rho = 0.9\))：

F3_7

– Adam的結果(\(lr = 0.1\), \(\beta_1 = 0.9\), \(\beta_2 = 0.999\))：

F3_11

看起來雖然Adam在數學式上看上去挺優美，結果也未必比較好。就目前的實驗數據來看，一般認為Adam的優勢在於對學習率較不敏感，但若指定一個恰當的SGD做優化任務，那SGD的優化結果將會比Adam更好。

– 因此，也有研究在訓練初期使用Adam做快速的優化，而在後期使用SGD做收尾(Nitish Shirish Keskar & Richard Socher, 2017)。事實上「最佳化問題」上一直以來都是數學上的千古難題，幾百年的進展下來並沒有一種方法可以在各種問題上都相較於其他有明顯優勢，我們只能說目前的主流是帶有慣性量的SGD與Adam是最常用的兩個方法，但可以預見未來進一步的研究也許會改寫這個格局。

練習3：使用Adam優化MLP

在練習2的實驗中，我們發現小批量相較於大樣本的優化表現較佳，但損失函數在後期並不穩定。讓我們試著以Adam優化MLP，並看看實驗結果如何。

– 初始程式碼的部分請你參考練習2的答案，用這個程式碼開始改。

– 使用的樣本如下：

set.seed(0)
x1 = rnorm(100, sd = 1) 
x2 = rnorm(100, sd = 1) 
lr1 = - 1.5 + x1^2 + x2^2 + rnorm(100)
y = lr1 > 0 + 0L

test_x1 = rnorm(100, sd = 1) 
test_x2 = rnorm(100, sd = 1) 
lr1 = - 1.5 + test_x1^2 + test_x2^2 + rnorm(100)
test_y = lr1 > 0 + 0L

另外，優化方法百百種，除了我們主要介紹的SGD、Adam以及Adagrad之外，還有許多變形，MxNet開發框架在Amazon的支持下有出一系列的網路教學內容，各位有興趣可以點選超連結參考。

練習3答案(程式碼)

修改後的程式碼如下：

MLP_Trainer = function (num.iteration = 500, num.hidden = 2, batch_size = 30,
                        lr = 0.001, beta1 = 0.9, beta2 = 0.999, optimizer = 'adam',
                        x1 = x1, x2 = x2, y = y,
                        test_x1 = NULL, test_x2 = NULL, test_y = NULL) {
  
  #Functions
  
  #Forward
  
  S.fun = function (x, eps = 1e-8) {
    S = 1/(1 + exp(-x))
    S[S < eps] = eps
    S[S > 1 - eps] = 1 - eps
    return(S)
  }
  
  ReLU.fun = function (x) {
    x[x < 0] <- 0
    return(x)
  }
  
  L.fun = function (X, W) {
    X.E = cbind(1, X)
    L = X.E %*% W
    return(L)
  }
  
  CE.fun = function (o, y, eps = 1e-8) {
    loss = -1/length(y) * sum(y * log(o + eps) + (1 - y) * log(1 - o + eps))
    return(loss)
  }
  
  #Backward
  
  grad_o.fun = function (o, y) {
    return((o - y)/(o*(1-o)))
  }
  
  grad_l2.fun = function (grad_o, o) {
    return(grad_o*(o*(1-o)))
  }
  
  grad_W2.fun = function (grad_l2, h1) {
    h1.E = cbind(1, h1)
    return(t(h1.E) %*% grad_l2/nrow(h1))
  }
  
  grad_h1.fun = function (grad_l2, W2) {
    return(grad_l2 %*% t(W2[-1,]))
  }
  
  grad_l1.fun = function (grad_h1, l1) {
    de_l1 = l1
    de_l1[de_l1<0] = 0
    de_l1[de_l1>0] = 1
    return(grad_h1*de_l1)
  }
  
  grad_W1.fun = function (grad_l1, x) {
    x.E = cbind(1, x)
    return(t(x.E) %*% grad_l1/nrow(x))
  }
  
  #Caculating
  
  X_matrix = cbind(x1, x2)
  
  W1_list = list()
  W2_list = list()
  loss_seq = rep(0, num.iteration)
  
  #Start random values
  
  W1_list[[1]] = matrix(rnorm(3*num.hidden, sd = 1), nrow = 3, ncol = num.hidden)
  W2_list[[1]] = matrix(rnorm(num.hidden + 1, sd = 1), nrow = num.hidden + 1, ncol = 1)
  
  for (i in 2:(num.iteration+1)) {
    
    idx <- sample(1:length(x1), batch_size)
    
    #Forward
    
    current_l1 = L.fun(X = X_matrix[idx,], W = W1_list[[i - 1]])
    current_h1 = ReLU.fun(x = current_l1)
    current_l2 = L.fun(X = current_h1, W = W2_list[[i - 1]])
    current_o = S.fun(x = current_l2)
    loss_seq[i-1] = CE.fun(o = current_o, y = y[idx], eps = 1e-8)
    
    #Backward
    
    current_grad_o = grad_o.fun(o = current_o, y = y[idx])
    current_grad_l2 = grad_l2.fun(grad_o = current_grad_o, o = current_o)
    current_grad_W2 = grad_W2.fun(grad_l2 = current_grad_l2, h1 = current_h1)
    current_grad_h1 = grad_h1.fun(grad_l2 = current_grad_l2, W2 = W2_list[[i - 1]])
    current_grad_l1 = grad_l1.fun(grad_h1 = current_grad_h1, l1 = current_l1)
    current_grad_W1 = grad_W1.fun(grad_l1 = current_grad_l1, x = X_matrix[idx,])
    
    if (optimizer == 'adam') {
      
      if (i == 2) {
        cum_m.W1 = 0
        cum_n.W1 = 0
        cum_m.W2 = 0
        cum_n.W2 = 0
      }
      
      cum_m.W1 = beta1 * cum_m.W1 + (1 - beta1) * current_grad_W1
      cum_n.W1 = beta2 * cum_n.W1 + (1 - beta2) * current_grad_W1^2
      cum_m.W1.hat = cum_m.W1/(1 - beta1^(i-1))
      cum_n.W1.hat = cum_n.W1/(1 - beta2^(i-1))
      
      cum_m.W2 = beta1 * cum_m.W2 + (1 - beta1) * current_grad_W2
      cum_n.W2 = beta2 * cum_n.W2 + (1 - beta2) * current_grad_W2^2
      cum_m.W2.hat = cum_m.W2/(1 - beta1^(i-1))
      cum_n.W2.hat = cum_n.W2/(1 - beta2^(i-1))
      
      W2_list[[i]] = W2_list[[i-1]] - lr * cum_m.W2.hat / sqrt(cum_n.W2.hat + 1e-8)
      W1_list[[i]] = W1_list[[i-1]] - lr * cum_m.W1.hat / sqrt(cum_n.W1.hat + 1e-8)
      
    } else if (optimizer == 'sgd') {
      if (i == 2) {
        W2_list[[i]] = W2_list[[i-1]] - lr * current_grad_W2
        W1_list[[i]] = W1_list[[i-1]] - lr * current_grad_W1
      } else {
        W2_list[[i]] = W2_list[[i-1]] - lr * current_grad_W2 + beta1 * (W2_list[[i-1]] - W2_list[[i-2]])
        W1_list[[i]] = W1_list[[i-1]] - lr * current_grad_W1 + beta1 * (W1_list[[i-1]] - W1_list[[i-2]])
      }
    } else {
      stop('optimizer must be selected from "sgd" or "adam".')
    }
    
    
    
  }
  
  require(scales)
  require(plot3D)
  
  x1_seq = seq(min(x1), max(x1), length.out = 100)
  x2_seq = seq(min(x2), max(x2), length.out = 100)
  
  pre_func = function (x1, x2, W1 = W1_list[[length(W1_list)]], W2 = W2_list[[length(W2_list)]]) {
    new_X = cbind(x1, x2)
    O = S.fun(x = L.fun(X = ReLU.fun(x = L.fun(X = new_X, W = W1)), W = W2))
    return(O)
  }
  
  pred_y = pre_func(x1 = x1, x2 = x2)
  MAIN_TXT = paste0('Train-Acc:', formatC(mean((pred_y > 0.5) == y), 2, format = 'f'))
  if (!is.null(test_x1)) {
    pred_test_y = pre_func(x1 = test_x1, x2 = test_x2)
    MAIN_TXT = paste0(MAIN_TXT, '; Test-Acc:', formatC(mean((pred_test_y > 0.5) == test_y), 2, format = 'f'))
  }
  
  z_matrix = sapply(x2_seq, function(x) {pre_func(x1 = x1_seq, x2 = x)})
  
  par(mfrow = c(1, 2))
  
  image2D(z = z_matrix, main = MAIN_TXT,
          x = x1_seq, xlab = 'x1',
          y = x2_seq, ylab = 'x2',
          shade = 0.2, rasterImage = TRUE,
          col = colorRampPalette(c("#FFA0A0", "#FFFFFF", "#A0A0FF"))(100))
  
  points(x1, x2, col = (y + 1)*2, pch = 19, cex = 0.5)
  if (!is.null(test_x1)) {
    points(test_x1, test_x2, col = 'black', bg = c('#C00000', '#0000C0')[(test_y + 1)], pch = 21)
  }
  
  plot(loss_seq, type = 'l', main = 'loss', xlab = 'iter.', ylab = 'CE loss')
  
}

練習3答案(實驗結果)

這是使用Adam的優化結果：

MLP_Trainer(num.iteration = 2000, num.hidden = 100, batch_size = 20,
            lr = 0.1, beta1 = 0.9, beta2 = 0.999, optimizer = 'adam',
            x1 = x1, x2 = x2, y = y,
            test_x1 = test_x1, test_x2 = test_x2, test_y = test_y)

這是使用SGD的優化結果(之前SGD的momentum等同於beta1)：

MLP_Trainer(num.iteration = 2000, num.hidden = 100, batch_size = 20,
            lr = 0.1, beta1 = 0.9, beta2 = 0.999, optimizer = 'sgd',
            x1 = x1, x2 = x2, y = y,
            test_x1 = test_x1, test_x2 = test_x2, test_y = test_y)

最終結果差別不大，但Adam據說優化效率較高，我們僅觀察100代看看：

MLP_Trainer(num.iteration = 100, num.hidden = 100, batch_size = 20,
            lr = 0.1, beta1 = 0.9, beta2 = 0.999, optimizer = 'adam',
            x1 = x1, x2 = x2, y = y,
            test_x1 = test_x1, test_x2 = test_x2, test_y = test_y)

MLP_Trainer(num.iteration = 100, num.hidden = 100, batch_size = 20,
            lr = 0.1, beta1 = 0.9, beta2 = 0.999, optimizer = 'sgd',
            x1 = x1, x2 = x2, y = y,
            test_x1 = test_x1, test_x2 = test_x2, test_y = test_y)

Adam只花了大概50代\(loss\)就已經下降到底了，而SGD似乎直到100代都還在大幅度震動，這樣看到Adam的優點了嗎！

極深網路的優化求解(1)

課程進行至此，我們尚未實際做過超極深的網路求解過程，事實上這件事情非常的困難。困難的原因在於非線性轉換函數的部分，因此容易在反向傳播的過程中慢慢損失梯度，導至網路最前端的權重無法更新。

– 在第二課我們已經改用ReLU作為非線性轉換函數，此舉雖然已經減少了梯度反向傳播的問題，但面對極深的網路時仍然需要其他幫助。

– 在這裡我們試著比較SGD以及Adam在網路極深時的表現，由於實現的程式碼較為複雜，我們直接提供現成的程式碼進行測試：

DEEP_MLP_Trainer = function (num.iteration = 500, num.hidden = c(10, 10, 10), batch_size = 30,
                             lr = 0.001, beta1 = 0.9, beta2 = 0.999, optimizer = 'adam', eps = 1e-8,
                             x1 = x1, x2 = x2, y = y,
                             test_x1 = NULL, test_x2 = NULL, test_y = NULL) {
  
  #Functions
  
  #Forward
  
  S.fun = function (x, eps = eps) {
    S = 1/(1 + exp(-x))
    S[S < eps] = eps
    S[S > 1 - eps] = 1 - eps
    return(S)
  }
  
  ReLU.fun = function (x) {
    x[x < 0] <- 0
    return(x)
  }
  
  L.fun = function (X, W) {
    X.E = cbind(1, X)
    L = X.E %*% W
    return(L)
  }
  
  CE.fun = function (o, y, eps = eps) {
    loss = -1/length(y) * sum(y * log(o + eps) + (1 - y) * log(1 - o + eps))
    return(loss)
  }
  
  #Backward
  
  grad_o.fun = function (o, y) {
    return((o - y)/(o*(1-o)))
  }
  
  grad_s.fun = function (grad_o, o) {
    return(grad_o*(o*(1-o)))
  }
  
  grad_W.fun = function (grad_l, h) {
    h.E = cbind(1, h)
    return(t(h.E) %*% grad_l/nrow(h))
  }
  
  grad_h.fun = function (grad_l, W) {
    return(grad_l %*% t(W[-1,]))
  }
  
  grad_l.fun = function (grad_h, l) {
    de_l = l
    de_l[de_l<0] = 0
    de_l[de_l>0] = 1
    return(grad_h*de_l)
  }
  
  #initialization
  
  X_matrix = cbind(x1, x2)
  Y_matrix = t(t(y))
  
  W_list = list()
  M_list = list()
  N_list = list()
  
  len_h = length(num.hidden)
  
  for (w_seq in 1:(len_h+1)) {
    if (w_seq == 1) {
      NROW_W = ncol(X_matrix) + 1
      NCOL_W = num.hidden[w_seq]
    } else if (w_seq == len_h+1) {
      NROW_W = num.hidden[w_seq - 1] + 1
      NCOL_W = ncol(Y_matrix)
    } else {
      NROW_W = num.hidden[w_seq - 1] + 1
      NCOL_W = num.hidden[w_seq]
    }
    W_list[[w_seq]] = matrix(rnorm(NROW_W*NCOL_W, sd = 1), nrow = NROW_W, ncol = NCOL_W)
    M_list[[w_seq]] = matrix(0, nrow = NROW_W, ncol = NCOL_W)
    N_list[[w_seq]] = matrix(0, nrow = NROW_W, ncol = NCOL_W)
  }
  
  loss_seq = rep(0, num.iteration)
  
  #Caculating
  
  for (i in 1:num.iteration) {
    
    idx = sample(1:length(x1), batch_size)
    
    #Forward
    
    current_l_list = list()
    current_h_list = list()
    
    for (j in 1:len_h) {
      if (j == 1) {
        current_l_list[[j]] = L.fun(X = X_matrix[idx,], W = W_list[[j]])
      } else {
        current_l_list[[j]] = L.fun(X = current_h_list[[j-1]], W = W_list[[j]])
      }
      current_h_list[[j]] = ReLU.fun(x = current_l_list[[j]])
    }
    current_l_list[[len_h+1]] = L.fun(X = current_h_list[[len_h]], W = W_list[[len_h+1]])
    current_o = S.fun(x = current_l_list[[len_h+1]], eps = eps)
    loss_seq[i] = CE.fun(o = current_o, y = y[idx], eps = eps)
    
    #Backward
    
    current_grad_l_list = list()
    current_grad_W_list = list()
    current_grad_h_list = list()
    
    current_grad_o = grad_o.fun(o = current_o, y = y[idx])
    current_grad_l_list[[len_h+1]] = grad_s.fun(grad_o = current_grad_o, o = current_o)
    current_grad_W_list[[len_h+1]] = grad_W.fun(grad_l = current_grad_l_list[[len_h+1]], h = current_h_list[[len_h]])
    
    for (j in len_h:1) {
      current_grad_h_list[[j]] = grad_h.fun(grad_l = current_grad_l_list[[j+1]], W = W_list[[j+1]])
      current_grad_l_list[[j]] = grad_l.fun(grad_h = current_grad_h_list[[j]], l = current_l_list[[j]])
      if (j != 1) {
        current_grad_W_list[[j]] = grad_W.fun(grad_l = current_grad_l_list[[j]], h = current_h_list[[j - 1]])
      } else {
        current_grad_W_list[[j]] = grad_W.fun(grad_l = current_grad_l_list[[j]], h = X_matrix[idx,])
      }
    }
    
    if (optimizer == 'adam') {
      
      for (j in 1:(len_h+1)) {
        M_list[[j]] = beta1 * M_list[[j]] + (1 - beta1) * current_grad_W_list[[j]]
        N_list[[j]] = beta2 * N_list[[j]] + (1 - beta2) * current_grad_W_list[[j]]^2
        M.hat = M_list[[j]]/(1 - beta1^i)
        N.hat = N_list[[j]]/(1 - beta2^i)
        W_list[[j]] = W_list[[j]] - lr*M.hat/sqrt(N.hat+eps)
      }
      
    } else if (optimizer == 'sgd') {
      
      for (j in 1:(len_h+1)) {
        M_list[[j]] = beta1 * M_list[[j]] + lr * current_grad_W_list[[j]]
        W_list[[j]] = W_list[[j]] - M_list[[j]]
      }
      
    } else {
      stop('optimizer must be selected from "sgd" or "adam".')
    }
    
  }
  
  require(scales)
  require(plot3D)
  
  x1_seq = seq(min(x1), max(x1), length.out = 100)
  x2_seq = seq(min(x2), max(x2), length.out = 100)
  
  pre_func = function (x1, x2) {
    new_X = cbind(x1, x2)
    
    current_l_list = list()
    current_h_list = list()
    
    for (j in 1:len_h) {
      if (j == 1) {
        current_l_list[[j]] = L.fun(X = new_X, W = W_list[[j]])
      } else {
        current_l_list[[j]] = L.fun(X = current_h_list[[j-1]], W = W_list[[j]])
      }
      current_h_list[[j]] = ReLU.fun(x = current_l_list[[j]])
    }
    
    current_l_list[[len_h+1]] = L.fun(X = current_h_list[[len_h]], W = W_list[[len_h+1]])
    current_o = S.fun(x = current_l_list[[len_h+1]], eps = eps)
    
    return(current_o)
  }
  
  pred_y = pre_func(x1 = x1, x2 = x2)
  MAIN_TXT = paste0('Train-Acc:', formatC(mean((pred_y > 0.5) == y), 2, format = 'f'))
  if (!is.null(test_x1)) {
    pred_test_y = pre_func(x1 = test_x1, x2 = test_x2)
    MAIN_TXT = paste0(MAIN_TXT, '; Test-Acc:', formatC(mean((pred_test_y > 0.5) == test_y), 2, format = 'f'))
  }
  
  z_matrix = sapply(x2_seq, function(x) {pre_func(x1 = x1_seq, x2 = x)})
  
  par(mfrow = c(1, 2))
  
  image2D(z = z_matrix, main = MAIN_TXT,
          x = x1_seq, xlab = 'x1',
          y = x2_seq, ylab = 'x2',
          shade = 0.2, rasterImage = TRUE,
          col = colorRampPalette(c("#FFA0A0", "#FFFFFF", "#A0A0FF"))(100))
  
  points(x1, x2, col = (y + 1)*2, pch = 19, cex = 0.5)
  if (!is.null(test_x1)) {
    points(test_x1, test_x2, col = 'black', bg = c('#C00000', '#0000C0')[(test_y + 1)], pch = 21)
  }
  
  plot(loss_seq, type = 'l', main = 'loss', xlab = 'iter.', ylab = 'CE loss')
  
}

極深網路的優化求解(2)

現在，讓我們做個實驗看看5個隱藏層的深度網路優化情形：

set.seed(0)
x1 = rnorm(100, sd = 1) 
x2 = rnorm(100, sd = 1) 
lr1 = - 1.5 + x1^2 + x2^2 + rnorm(100)
y = lr1 > 0 + 0L

test_x1 = rnorm(100, sd = 1) 
test_x2 = rnorm(100, sd = 1) 
lr1 = - 1.5 + test_x1^2 + test_x2^2 + rnorm(100)
test_y = lr1 > 0 + 0L

– 這是Adam的結果(這裡我們都使用一般常規默認參數\(lr=0.001\))

DEEP_MLP_Trainer(num.iteration = 2000, num.hidden = c(10, 10, 10, 10, 10), batch_size = 20,
                 lr = 0.001, beta1 = 0.9, beta2 = 0.999, optimizer = 'adam',
                 x1 = x1, x2 = x2, y = y,
                 test_x1 = test_x1, test_x2 = test_x2, test_y = test_y)

– 這是sgd的結果(一般默認\(lr=0.05\))

DEEP_MLP_Trainer(num.iteration = 2000, num.hidden = c(10, 10, 10, 10, 10), batch_size = 20,
                 lr = 0.05, beta1 = 0.9, beta2 = 0.999, optimizer = 'sgd',
                 x1 = x1, x2 = x2, y = y,
                 test_x1 = test_x1, test_x2 = test_x2, test_y = test_y)

極深網路的優化求解(3)

我們可以發現sgd似乎優化失敗了，讓我們做一系列的實驗來看看他的最適學習率：

– \(lr = 0.01\)

DEEP_MLP_Trainer(num.iteration = 2000, num.hidden = c(10, 10, 10, 10, 10), batch_size = 20,
                 lr = 0.01, beta1 = 0.9, beta2 = 0.999, optimizer = 'sgd',
                 x1 = x1, x2 = x2, y = y,
                 test_x1 = test_x1, test_x2 = test_x2, test_y = test_y)

– \(lr = 0.001\)

DEEP_MLP_Trainer(num.iteration = 2000, num.hidden = c(10, 10, 10, 10, 10), batch_size = 20,
                 lr = 0.001, beta1 = 0.9, beta2 = 0.999, optimizer = 'sgd',
                 x1 = x1, x2 = x2, y = y,
                 test_x1 = test_x1, test_x2 = test_x2, test_y = test_y)

– \(lr = 0.0001\)

DEEP_MLP_Trainer(num.iteration = 2000, num.hidden = c(10, 10, 10, 10, 10), batch_size = 20,
                 lr = 0.0001, beta1 = 0.9, beta2 = 0.999, optimizer = 'sgd',
                 x1 = x1, x2 = x2, y = y,
                 test_x1 = test_x1, test_x2 = test_x2, test_y = test_y)

我們可以觀察到，SGD對於學習率相當的敏感，讓我們再試試看Adam

極深網路的優化求解(4)

這是Adam的優化結果：