第一節：邏輯斯回歸的限制(1)

目前我們學的都到的都是線性模型，其中最好的彈性網路其實也只是邏輯斯回歸的加強版而已。

– 下面是一個能夠對任意2維分布做邏輯斯迴歸運算並分類的函數，請大家試試看不同分布的狀況下邏輯斯迴歸的分類能力。

LOGISTIC_PLOT <- function (x1, x2, y, fomula) {
  
  require(scales)
  require(plot3D)
  
  model = glm(fomula, family = 'binomial')
  
  x1_seq = seq(min(x1), max(x1), length.out = 100)
  x2_seq = seq(min(x2), max(x2), length.out = 100)
  
  z_matrix = sapply(x2_seq, function(x) {1/(1+exp(-predict(model, data.frame(x1 = x1_seq, x2 = x))))})
  
  
  image2D(z = z_matrix,
          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)
  
}

set.seed(0)
x1 = rnorm(50, sd = 1) 
x2 = rnorm(50, sd = 1) 
lr1 = 1 + 0.5 * x1 + 2 * x2
y = lr1 > 0

LOGISTIC_PLOT(x1 = x1, x2 = x2, y = y, fomula = y ~ x1 + x2)

第一節：邏輯斯回歸的限制(2)

讓我們試試看不同分布的效果吧：

set.seed(0)
x1 = rnorm(50, sd = 1) 
x2 = rnorm(50, sd = 1) 
lr1 = 1 + 0.5 * x1 + 0.7 * x2 + 3 * x1 * x2
y = lr1 > 0

LOGISTIC_PLOT(x1 = x1, x2 = x2, y = y, fomula = y ~ x1 + x2)

這個情境下，似乎並不是分割得非常好，我們可以透過在迴歸式中加入交互作用來解決：

LOGISTIC_PLOT(x1 = x1, x2 = x2, y = y, fomula = y ~ x1 + x2 + x1:x2)

第一節：邏輯斯回歸的限制(3)

再換個分布吧，我們同樣嘗試具有交互作用項的邏輯斯回歸：

set.seed(0)
x1 = rnorm(50, sd = 1) 
x2 = rnorm(50, sd = 1) 
lr1 = - 0.5 + 0.5 * x1^2 + 0.3 * x2^2 + 0.4 * x1 * x2
y = lr1 > 0

LOGISTIC_PLOT(x1 = x1, x2 = x2, y = y, fomula = y ~ x1 + x2 + x1:x2)

看來加入交互作用項並非萬解，其實這個情境下還要加入平方項：

LOGISTIC_PLOT(x1 = x1, x2 = x2, y = y, fomula = y ~ x1 + x2 + x1:x2 + poly(x1, 2) + poly(x2, 2))

第二節：神經網路介紹(1)

線性迴歸及邏輯斯回歸最大的問題在於其是線性預測方法，他的預測式限制了他的發展潛力。當遇到「線性不可分割」的問題時，他們將需要「特徵工程」的協助才能做出完美預測。

– 線性迴歸

\[\hat{y} = f(x) = b_{0} + b_{1}x_1 + b_{2}x_2\]

– 邏輯斯回歸

\[log(\frac{{p}}{1-p}) = b_{0} + b_{1}x_1 + b_{2}x_2\]

儘管我們剛剛透過「特徵工程」做出了完美預測，但在低維空間中我們還能透過肉眼敏銳的觀察，但到高維空間中這將顯得不現實，我們迫切的需要一種方法能夠幫助我們「自動」做特徵工程。

第二節：神經網路介紹(2)

這個問題從1950年代後吸引眾多學者研究，而解決的方法其實很多，而效果各異。為了更好的與後續深度學習的課程做結合，我們這裡從David Rumelhart、Geoffrey Hinton、Ronald Williams在1986年的研究開始講起（事實上這並非當時的最佳解法）。

F01

以剛剛的問題來說，我們不禁思考為什麼辨認貓狗對於人類來說這麼容易的事情，讓電腦做起來卻如此困難?

– 這時候生物學家開始踏入了這個領域，如果我們能了解大腦是如何運作的，那是否能模仿其結構進行分類?

第二節：神經網路介紹(3)

大腦是神經系統的一部分，而他是由多個神經元互相結合而成的。(下圖為Ramón y Cajal在1905年所畫下的神經細胞樣貌)

F03

– 神經細胞的構造如下，不論是何種神經元皆可分成：接收區、觸發區、傳導區和輸出區。

F04

神經細胞傳導訊息的過程是透過電位的變化，在接收區會接收到一個電位改變的訊號，再交給突觸來處理。

F05

透過樹突(dendrite)能接收上一個神經元的訊息，而有些會在接收訊息後產生抑制性作用，有些會產生興奮性作用，然後這些訊號再透過神經元整合，之後再透過軸突(axon)將訊號傳導出去。
我們根據這樣的生物學知識，開始來用電腦模擬一個簡單的神經元。

第二節：神經網路介紹(4)

第一代的人工智慧是1958年由Frank Rosenblatt所發展的感知機(perceptron)，此時還稱不上是神經「網路」，他的結構如下：

F06

讓我們用數學語言來描述他：

\[ \begin{align} \mbox{weighted sum} & = w_{0} + w_{1}x_1 + w_{2}x_2 + \dots \\ \hat{y} & = step(\mbox{weighted sum}) \end{align} \]

仔細一看，他跟邏輯斯回歸實在是很像，差別只是在output並非是一個機率，而是直接決定是TRUE還是FALSE

– 既然是邏輯斯回歸的簡單版，那就不用談了，我們已經很清楚邏輯斯回歸的極限了。

第二節：神經網路介紹(5)

顯然大腦的構成不是僅僅使用一個神經元，既然邏輯斯回歸的構造與神經元很像，所以我們試著把多個邏輯斯回歸結合在一起吧！

– 我們不要拿太複雜的結構，就用下面這個最簡單的結構，這被稱為多層感知機(Multilayer Perceptron, MLP)

F07

讓我們用數學語言來描述他，為了簡化表達此一預測方程，我們將大量的使用函數組來描述：

線性預測函數L：

\[L^k(x_1, x_2) = w_{0}^k + w_{1}^kx_1 + w_{2}^kx_2\]

邏輯斯轉換函數S：

\[ \begin{align} S(x) & = \frac{{1}}{1+e^{-x}} \end{align} \]

多層感知機預測函數：

\[ \begin{align} h_1 & = S(L^1(x_1, x_2)) \\ h_2 & = S(L^2(x_1, x_2)) \\ o & = S(L^3(h_1, h_2)) \end{align} \]

第三節：反向傳播演算法(1)

有了多層感知機的預測函數之後，我們需要再給他一個損失函數再求解。而剛剛那篇1986年的Nature文章就是結合連鎖率與相關微分技術用於多層感知機，並將他命名為反向傳播演算法(backpropagation algorithm)。在這裡我們使用的是交叉熵函數作為損失函數，我們把求解目標式寫下來：

\[ \begin{align} h_1 & = S(L^1(x_1, x_2)) \\ h_2 & = S(L^2(x_1, x_2)) \\ o & = S(L^3(h_1, h_2)) \\ loss & = CE(y, o) = \frac{{1}}{n}\sum \limits_{i=1}^{n} -\left(y_{i} \cdot log(o_{i}) + (1-y_{i}) \cdot log(1-o_{i})\right) \end{align} \]

我們的目標很明確，找出一組\(w_{0}^1\)、\(w_{1}^1\)、\(w_{2}^1\)、\(w_{0}^2\)、\(w_{1}^2\)、\(w_{2}^2\)、\(w_{0}^3\)、\(w_{1}^3\)、\(w_{2}^3\)，使的整體的loss最小化。我們一樣是使用梯度下降法進行求解，而其實就是要求得上述所有係數的偏導函數。

– 其實這種複雜的微分是有規律可循的，我們把幾個重要的微分工具寫出來：

連鎖率

\[\frac{\partial}{\partial x}h(x) = \frac{\partial}{\partial x}f(g(x)) = \frac{\partial}{\partial g(x)}f(g(x)) \cdot\frac{\partial}{\partial x}g(x)\]

邏輯斯轉換函數S微分(第一課有詳細解法)

\[ \begin{align} \frac{\partial}{\partial x}S(x) & = S(x)(1-S(x)) \end{align} \]

單一交叉熵函數的微分(第一課有詳細解法)

\[ \begin{align} \frac{\partial}{\partial p_i}CE(y_i, p_i) & = \frac{p_i - y_i}{p_i(1-p_i)} \end{align} \]

第三節：反向傳播演算法(2)

下列是\(w_{0}^1\)、\(w_{1}^1\)、\(w_{2}^1\)、\(w_{0}^2\)、\(w_{1}^2\)、\(w_{2}^2\)、\(w_{0}^3\)、\(w_{1}^3\)、\(w_{2}^3\)各自的偏微分：

– 先求離結果比較近的部分(比較深層)

\(w_{0}^3\)的偏導函數：

\[ \begin{align} \frac{\partial}{\partial w_{0}^3}CE(y, o) & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left(\frac{\partial}{\partial w_{0}^3}CE(y_i, o_i)\right) \\ & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left(\frac{\partial}{\partial o_i}CE(y_i, o_i) \cdot \frac{\partial}{\partial w_{0}^3}S(L^3(h_{1i}, h_{2i})) \right) \\ & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left(\frac{o_i-y_i}{o_i(1-o_i)} \cdot \frac{\partial}{\partial L^3(h_{1i}, h_{2i})}S(L^3(h_{1i}, h_{2i})) \cdot \frac{\partial}{\partial w_{0}^3}L^3(h_{1i}, h_{2i}) \right) \\ & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left(\frac{o_i-y_i}{o_i(1-o_i)} \cdot S(L^3(h_{1i}, h_{2i})) \cdot (1 - S(L^3(h_{1i}, h_{2i}))) \cdot \frac{\partial}{\partial w_{0}^3} \left( w_{0}^3 + w_{1}^3h_{1i} + w_{2}^3h_{2i} \right) \right) \\ & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left(\frac{o_i-y_i}{o_i(1-o_i)} \cdot o_i(1-o_i) \right) \\ & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \end{align} \]

\(w_{1}^3\)的偏導函數(過程略)：

\[ \begin{align} \frac{\partial}{\partial w_{1}^3}CE(y, o) & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot h_{1i} \end{align} \]

\(w_{2}^3\)的偏導函數(過程略)：

\[ \begin{align} \frac{\partial}{\partial w_{2}^3}CE(y, o) & = \frac{{1}}{n}\sum \limits_{i=1}^{n}\left( o_i-y_i \right) \cdot h_{2i} \end{align} \]

– 再求離結果比較遠的部分(比較淺層)

\(w_{0}^2\)的偏導函數：

\[ \begin{align} \frac{\partial}{\partial w_{0}^2}CE(y, o) & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot \frac{\partial}{\partial w_{0}^2} \left( w_{0}^3 + w_{1}^3h_{1i} + w_{2}^3h_{2i} \right) \\ & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot w_{2}^3\cdot \frac{\partial}{\partial w_{0}^2} \left( h_{2i} \right) \\ & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot w_{2}^3\cdot \frac{\partial}{\partial w_{0}^2} \left( S(L^2(x_{1i}, x_{2i})) \right) \\ & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot w_{2}^3 \cdot \frac{\partial}{\partial L^2(x_{1i}, x_{2i})} S(L^2(x_{1i}, x_{2i})) \cdot \frac{\partial}{\partial w_{0}^2} L^2(x_{1i}, x_{2i}) \\ & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot w_{2}^3 \cdot h_{2i} (1 - h_{2i}) \cdot \frac{\partial}{\partial w_{0}^2} (w_{0}^2 + w_{1}^2x_{1i} + w_{2}^2x_{2i}) \\ & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot w_{2}^3 \cdot h_{2i} (1 - h_{2i}) \end{align} \]

\(w_{1}^2\)的偏導函數(過程略)：

\[ \begin{align} \frac{\partial}{\partial w_{1}^2}CE(y, o) & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot w_{2}^3 \cdot h_{2i} (1 - h_{2i}) \cdot x_{1i} \end{align} \]

\(w_{2}^2\)的偏導函數(過程略)：

\[ \begin{align} \frac{\partial}{\partial w_{2}^2}CE(y, o) & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot w_{2}^3 \cdot h_{2i} (1 - h_{2i}) \cdot x_{2i} \end{align} \]

\(w_{0}^1\)的偏導函數(過程略)：

\[ \begin{align} \frac{\partial}{\partial w_{0}^1}CE(y, o) & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot w_{1}^3 \cdot h_{1i} (1 - h_{1i}) \end{align} \]

\(w_{1}^1\)的偏導函數(過程略)：

\[ \begin{align} \frac{\partial}{\partial w_{1}^1}CE(y, o) & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot w_{1}^3 \cdot h_{1i} (1 - h_{1i}) \cdot x_{1i} \end{align} \]

\(w_{2}^1\)的偏導函數(過程略)：

\[ \begin{align} \frac{\partial}{\partial w_{2}^1}CE(y, o) & = \frac{{1}}{n}\sum \limits_{i=1}^{n} \left( o_i-y_i \right) \cdot w_{1}^3 \cdot h_{1i} (1 - h_{1i}) \cdot x_{2i} \end{align} \]

練習1：利用梯度下降法獲得一個多層感知機

怕你看不懂公式，幫妳寫成公式的樣子

#Sample generation

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

#Forward

S.fun = function (x, eps = 1e-5) {
  S = 1/(1 + exp(-x))
  S[S < eps] = eps
  S[S > 1 - eps] = 1 - eps
  return(S)
}

h1.fun = function (w10, w11, w12, x1 = x1, x2 = x2) {
  L1 = w10 + w11 * x1 + w12 * x2
  return(S.fun(L1))
}

h2.fun = function (w20, w21, w22, x1 = x1, x2 = x2) {
  L2 = w20 + w21 * x1 + w22 * x2
  return(S.fun(L2))
}

o.fun = function (w30, w31, w32, h1, h2) {
  L3 = w30 + w31 * h1 + w32 * h2
  return(S.fun(L3))
}

loss.fun = function (o, y = y) {
  loss = -1/length(y) * sum(y * log(o) + (1 - y) * log(1 - o))
  return(loss)
}

#Backward

differential.fun.w30 = function(o, y = y) {
  return(-1/length(y)*sum(y-o))
}

differential.fun.w31 = function(o, h1, y = y) {
  return(-1/length(y)*sum((y-o)*h1))
}

differential.fun.w32 = function(o, h2, y = y) {
  return(-1/length(y)*sum((y-o)*h2))
}

differential.fun.w20 = function(o, h2, w32, y = y) {
  return(-1/length(y)*sum((y-o)*w32*h2*(1-h2)))
}

differential.fun.w21 = function(o, h2, w32, y = y, x1 = x1) {
  return(-1/length(y)*sum((y-o)*w32*h2*(1-h2)*x1))
}

differential.fun.w22 = function(o, h2, w32, y = y, x2 = x2) {
  return(-1/length(y)*sum((y-o)*w32*h2*(1-h2)*x2))
}

differential.fun.w10 = function(o, h1, w31, y = y) {
  return(-1/length(y)*sum((y-o)*w31*h1*(1-h1)))
}

differential.fun.w11 = function(o, h1, w31, y = y, x1 = x1) {
  return(-1/length(y)*sum((y-o)*w31*h1*(1-h1)*x1))
}

differential.fun.w12 = function(o, h1, w31, y = y, x2 = x2) {
  return(-1/length(y)*sum((y-o)*w31*h1*(1-h1)*x2))
}

練習1答案(1)

num.iteration = 10000
lr = 0.1
W_matrix = matrix(0, nrow = num.iteration + 1, ncol = 9)
loss_seq = rep(0, num.iteration)
colnames(W_matrix) = c('w10', 'w11', 'w12', 'w20', 'w21', 'w22', 'w30', 'w31', 'w32')

#Start random values
W_matrix[1,] = rnorm(9, sd = 1) 

for (i in 2:(num.iteration+1)) {
  
  #Forward
  
  current_H1 = h1.fun(w10 = W_matrix[i-1,1], w11 = W_matrix[i-1,2], w12 = W_matrix[i-1,3],
                      x1 = x1, x2 = x2)
  
  current_H2 = h2.fun(w20 = W_matrix[i-1,4], w21 = W_matrix[i-1,5], w22 = W_matrix[i-1,6],
                      x1 = x1, x2 = x2)
  
  current_O = o.fun(w30 = W_matrix[i-1,7], w31 = W_matrix[i-1,8], w32 = W_matrix[i-1,9],
                    h1 = current_H1, h2 = current_H2)
  
  loss_seq[i-1] = loss.fun(o = current_O, y = y)
  
  #Backward
  
  W_matrix[i,1] = W_matrix[i-1,1] - lr * differential.fun.w10(o = current_O, h1 = current_H1,
                                       w31 = W_matrix[i-1,8], y = y)
  
  W_matrix[i,2] = W_matrix[i-1,2] - lr * differential.fun.w11(o = current_O, h1 = current_H1,
                                       w31 = W_matrix[i-1,8], y = y, x1 = x1)
  
  W_matrix[i,3] = W_matrix[i-1,3] - lr * differential.fun.w12(o = current_O, h1 = current_H1,
                                       w31 = W_matrix[i-1,8], y = y, x2 = x2)
  
  W_matrix[i,4] = W_matrix[i-1,4] - lr * differential.fun.w20(o = current_O, h2 = current_H2,
                                       w32 = W_matrix[i-1,9], y = y)
    
  W_matrix[i,5] = W_matrix[i-1,5] - lr * differential.fun.w21(o = current_O, h2 = current_H2,
                                       w32 = W_matrix[i-1,9], y = y, x1 = x1)
  
  W_matrix[i,6] = W_matrix[i-1,6] - lr * differential.fun.w22(o = current_O, h2 = current_H2,
                                       w32 = W_matrix[i-1,9], y = y, x2 = x2)
    
  W_matrix[i,7] = W_matrix[i-1,7] - lr * differential.fun.w30(o = current_O, y = y)
    
  W_matrix[i,8] = W_matrix[i-1,8] - lr * differential.fun.w31(o = current_O, h1 = current_H1, y = y)
  
  W_matrix[i,9] = W_matrix[i-1,9] - lr * differential.fun.w32(o = current_O, h2 = current_H2, y = y)
  
}

練習1答案(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, W_list = W_matrix[nrow(W_matrix),]) {
  H1 = h1.fun(w10 = W_list[1], w11 = W_list[2], w12 = W_list[3], x1 = x1, x2 = x2)
  H2 = h2.fun(w20 = W_list[4], w21 = W_list[5], w22 = W_list[6], x1 = x1, x2 = x2)
  O = o.fun(w30 = W_list[7], w31 = W_list[8], w32 = W_list[9], h1 = H1, h2 = H2)
  return(O)
}

z_matrix = sapply(x2_seq, function(x) {pre_func(x1 = x1_seq, x2 = x)})
  
  
image2D(z = z_matrix,
        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)

這真是太有趣了，組合3個邏輯斯回歸讓分割邊界產生曲線，讓我們觀察一下\(loss\)下降的軌跡

plot(loss_seq, type = 'l', main = 'loss', xlab = 'iter.', ylab = 'CE loss')

第四節：矩陣化多層感知機(1)

透過多個人工神經元(邏輯斯迴歸)達到的效果令人驚嘆，但每次增加神經元時都將導致微分過程的變化，從而導致其難以擴展。

– 透過矩陣能夠將多層感知機變成較為簡易的公式組合(以下為個人層級的方程式)：

線性預測函數L之矩陣式(d為輸出維度)：

\[ \begin{align} L^k_d(X, W^k_d) & = XW^k_d \\ X & = \begin{pmatrix} x_{1,1} & x_{1,2} & \cdots & x_{1,m} \end{pmatrix} \\ W^k_d & = \begin{pmatrix} w^k_{0,1} & w^k_{0,2} & \cdots & w^k_{0,d} \\ w^k_{1,1} & w^k_{1,2} & \cdots & w^k_{1,d} \\ w^k_{2,1} & w^k_{2,2} & \cdots & w^k_{2,d} \\ \vdots & \vdots & \ddots & \vdots \\ w^k_{m,1} & w^k_{m,2} & \cdots & w^k_{m,d} \end{pmatrix} \\ \frac{\partial}{\partial W^k_d}L^k_d(X) & = \begin{pmatrix} X^T & & X^T & \cdots & X^T \\ \end{pmatrix} \mbox{ [repeat } d \mbox{ times]} \end{align} \]

邏輯斯轉換函數S：

\[ \begin{align} S(x) & = \frac{{1}}{1+e^{-x}} \\ \frac{\partial}{\partial x}S(x) & = S(x)(1-S(x)) \end{align} \]

多層感知機預測函數之矩陣式，矩陣上標E代表該矩陣的第一欄已被填滿1(這是剛剛的單一隱藏層網路，其隱藏層中神經元數目為d，在剛剛的例子中d=2)：

\[ \begin{align} l_1 & = L^1_d(x^E,W^1_d) \\ h_1 & = S(l_1) \\ l_2 & = L^2_1(h_1^E,W^2_1) \\ o & = S(l_2) \end{align} \]

第四節：矩陣化多層感知機(2)

我們同樣的需要解矩陣的微分方程，求解目標式如下：

\[ \begin{align} l_1 & = L^1_d(x^E,W^1_d) \\ h_1 & = S(l_1) \\ l_2 & = L^2_1(h_1^E,W^2_1) \\ o & = S(l_2) \\ loss & = CE(y, o) = \frac{{1}}{n}\sum \limits_{i=1}^{n} -\left(y \cdot log(o) + (1-y) \cdot log(1-o)\right) \end{align} \]

讓我們分別對這裡所有元素進行微分(符號\(\otimes\)為按元素相乘，符號\(\bullet\)為按矩陣乘法)：

\[ \begin{align} grad.o & = \frac{\partial}{\partial o}loss = \frac{o-y}{o(1-o)} \\ grad.l_2 & = \frac{\partial}{\partial l_2}loss = grad.o \otimes \frac{\partial}{\partial l_2}o= o-y \\ grad.W^2_1 & = \frac{\partial}{\partial W^2_1}loss = grad.l_2 \otimes \frac{\partial}{\partial W^2_1}l_2 = \frac{1}{n} \otimes (h_1^E)^T \bullet grad.l_2\\ grad.h_1^E & = \frac{\partial}{\partial h_1^E}loss = grad.l_2 \otimes \frac{\partial}{\partial h_1^E}l_2 = grad.l_2 \bullet (W^2_1)^T \\ grad.l_1 & = \frac{\partial}{\partial l_1}loss = grad.h_1 \otimes \frac{\partial}{\partial l_1}h_1 = grad.h_1 \otimes h_1 \otimes (1-h_1) \\ grad.W^1_d & = \frac{\partial}{\partial W^1_d}loss = grad.l_1 \otimes \frac{\partial}{\partial W^1_d}l_1 = \frac{1}{n} \otimes (x^E)^T \bullet grad.l_1 \end{align} \]

練習2：利用矩陣式求解並試著增加隱藏層的數目

試著實現矩陣化的MLP吧！

– 我們用程式碼描述上述式子：

#Sample generation

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

#Forward

S.fun = function (x, eps = 1e-5) {
  S = 1/(1 + exp(-x))
  S[S < eps] = eps
  S[S > 1 - eps] = 1 - eps
  return(S)
}

L.fun = function (X, W) {
  X.E = cbind(1, X)
  L = X.E %*% W
  return(L)
}

CE.fun = function (o, y, eps = 1e-5) {
  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, h1) {
  return(grad_h1*(h1*(1-h1)))
}

grad_W1.fun = function (grad_l1, x) {
  x.E = cbind(1, x)
  return(t(x.E) %*% grad_l1/nrow(x))
}

練習2答案(1)

這是一個MLP訓練函數，由於已經矩陣化了，我們可以很方便的指定隱藏層的神經元數量(num.hidden)

MLP_Trainer = function (num.iteration = 500, num.hidden = 2, lr = 0.1, x1 = x1, x2 = x2, y = y) {
  
  #Functions
  
  #Forward
  
  S.fun = function (x, eps = 1e-5) {
    S = 1/(1 + exp(-x))
    S[S < eps] = eps
    S[S > 1 - eps] = 1 - eps
    return(S)
  }
  
  L.fun = function (X, W) {
    X.E = cbind(1, X)
    L = X.E %*% W
    return(L)
  }
  
  CE.fun = function (o, y, eps = 1e-5) {
    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, h1) {
    return(grad_h1*(h1*(1-h1)))
  }
  
  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 = S.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-5)
    
    #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, h1 = current_h1)
    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 = S.fun(x = L.fun(X = new_X, W = W1)), W = W2))
    return(O)
  }
  
  z_matrix = sapply(x2_seq, function(x) {pre_func(x1 = x1_seq, x2 = x)})
  
  par(mfrow = c(1, 2))
  
  image2D(z = z_matrix,
          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)
  
  plot(loss_seq, type = 'l', main = 'loss', xlab = 'iter.', ylab = 'CE loss')
  
}

練習2答案(2)

看結果

MLP_Trainer(num.iteration = 10000, num.hidden = 2, lr = 0.1, x1 = x1, x2 = x2, y = y)

現在讓我們試試看多給他幾個神經元，是不是又切的更好了?

MLP_Trainer(num.iteration = 10000, num.hidden = 5, lr = 0.1, x1 = x1, x2 = x2, y = y)

第五節：使用套件進行多層感知機(1)

– 請至這裡下載範例資料

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

套件「nnet」內的函數「nnet()」可以協助我們快速地進行多層感知機的建立。

library(nnet)

第五節：使用套件進行多層感知機(2)

這個函數能支援3種型態的outcome，其範例如下：

– 這是想要做線性預測的方式：

subdat <- dat[!(dat[,'K'] %in% NA) & !(dat[,'GENDER'] %in% NA) & !(dat[,'AGE'] %in% NA), c('K', 'GENDER', 'AGE')]
subdat[,'GENDER'] <- as.factor(subdat[,'GENDER'])

fit_nnet <- nnet(K ~ GENDER + AGE, data = subdat, size = 5, maxit = 100, linout = TRUE, decay = 0.001, trace = FALSE, MaxNWts = 1e5)
pred_val <- predict(fit_nnet, subdat)

– 這是想要做二元類別預測的方式：

subdat <- dat[!(dat[,'LVD'] %in% NA) & !(dat[,'GENDER'] %in% NA) & !(dat[,'AGE'] %in% NA), c('LVD', 'GENDER', 'AGE')]
subdat[,'GENDER'] <- as.factor(subdat[,'GENDER'])

fit_nnet <- nnet(LVD ~ GENDER + AGE, data = subdat, size = 5, maxit = 100, entropy = TRUE, decay = 0.001, trace = FALSE, MaxNWts = 1e5)
pred_val <- predict(fit_nnet, subdat)

– 這是想要做多元類別預測的方式：

subdat <- dat[!(dat[,'AMI'] %in% NA) & !(dat[,'GENDER'] %in% NA) & !(dat[,'AGE'] %in% NA), c('AMI', 'GENDER', 'AGE')]
subdat[,'GENDER'] <- as.factor(subdat[,'GENDER'])
subdat[,'AMI'] <- as.factor(subdat[,'AMI'])

fit_nnet <- nnet(AMI ~ GENDER + AGE, data = subdat, size = 5, maxit = 100, decay = 0.001, trace = FALSE, MaxNWts = 1e5)
pred_val <- predict(fit_nnet, subdat)

練習3：比較彈性網路與多層感知機的結果

讓我們來進行一個資料科學實驗，比較線性模型與非線性模型在鉀離子預測上的準確度：

subdat <- dat[apply(is.na(dat[,c(-1, -3:-5)]), 1, sum) %in% 0, c(-1, -3:-5)]

subdat[,'GENDER'] <- as.factor(subdat[,'GENDER'])

set.seed(0)
all_idx <- 1:nrow(subdat)

train_idx <- sample(all_idx, nrow(subdat) * 0.6)
valid_idx <- sample(all_idx[!all_idx %in% train_idx], nrow(subdat) * 0.2)
test_idx <- all_idx[!all_idx %in% c(train_idx, valid_idx)]

train_dat <- subdat[train_idx,]
valid_dat <- subdat[valid_idx,]
test_dat <- subdat[test_idx,]

比較一下彈性網路與多層感知機的準確度吧！

練習3答案(1)

這是彈性網路的語法，我們使用\(\alpha = 0.5\)來進行實驗。

– 需要注意一下彈性網路的函數需要特定格式的輸入喔！

library(glmnet)

train_X <- model.matrix(~ ., data = train_dat[,-1])
train_X <- train_X[,-1]
train_Y <- train_dat[,1]

valid_X <- model.matrix(~ ., data = valid_dat[,-1])
valid_X <- valid_X[,-1]
valid_Y <- valid_dat[,1]

glmnet_fit <- cv.glmnet(x = train_X, y = train_Y, family = 'gaussian', alpha = 0.5, nfolds = 5)
pred_valid <- predict(glmnet_fit, valid_X, s = "lambda.1se")
cor(pred_valid, valid_Y)

##       [,1]
## 1 0.333089

– 相關係數是0.333，再讓我們試試看MLP：

library(nnet)

fit_nnet <- nnet(K ~ ., data = train_dat, size = 100, maxit = 100, linout = TRUE, decay = 0.001, trace = FALSE, MaxNWts = 1e5)
pred_valid <- predict(fit_nnet, valid_dat)
cor(pred_valid, valid_Y)

##           [,1]
## [1,] 0.3116847

– 相關係數是0.312，好像沒有比較好。

練習3答案(2)

但在訓練樣本上是確實比較好的，可以看一下這樣的結果：

– 這是彈性網路，結果是0.375

pred_train <- predict(glmnet_fit, train_X, s = "lambda.1se")
cor(pred_train, train_Y)

##        [,1]
## 1 0.3745448

– 這是多層感知機，結果是0.515

pred_train <- predict(fit_nnet, train_dat)
cor(pred_train, train_Y)

##           [,1]
## [1,] 0.5152867

這個叫做過度擬合(overfitting)，在使用非線性模型時真的要非常小心，參數沒有調好就會發生這樣的事情！未來我們會有更多課程來教大家避免這個問題。

F08

課程小結

本週介紹的多層感知機(Multilayer Perceptron, MLP)很明顯是基於統計工具的擴展，因此他其實還可以擴展至「存活分析」之上。

– 很可惜的是，並沒有方便的套件能幫助我們利用MLP做存活分析，這裡有待後續課程的學習。

– 現代的深度學習其實就是基於這個技術，但想要擴展至更深層的網路還需要更多的數學技術，相關突破會在下學期的課程繼續。

現在，我們終於會讓模型「自動」進行特徵工程，這解決了我們很多問題，我們現在再也不需要關心\(x\)與\(y\)之間的關聯性了。

– 隨著模型變得更加複雜，「正則化」技術對我們更加重要了，資料科學實驗的完整流程是非常重要的，你必須調整好最佳的參數才能得到好的結果。

– 過度擬合(overfitting)是非線性模型面臨最重要的問題，我們下學期的課程會強調這個部份的技術！

機器學習及演算法