[譯] Scratch 平臺的神經網路實現（R 語言）

Cup發表於2019-02-27

原文地址：Neural Networks from Scratch (in R)

原文作者：Ilia Karmanov

譯文出自：掘金翻譯計劃

本文永久連結：github.com/xitu/gold-m…

譯者：CACppuccino

校對者：Isvih

Scratch 平臺的神經網路實現（R 語言）

這篇文章是針對那些有統計或者經濟學背景的人們，幫助他們通過 R 語言上的 Scratch 平臺更好地學習和理解機器學習知識。

Andrej Karpathy 在 CS231n 課程中這樣說道：

“我們有意識地在設計課程的時候，於反向傳播演算法的程式設計作業中包含了對最底層的資料的計算要求。學生們需要在原始的 numpy 庫中使資料在各層中正向、反向傳播。一些學生因而難免在課程的留言板上抱怨（這些複雜的計算）”

如果框架已經為你完成了反向傳播演算法（BP 演算法）的計算，你又何苦折磨自己而不去探尋更多有趣的深度學習問題呢？

import keras
model = Sequential()
model.add(Dense(512, activation=’relu’, input_shape=(784,)))
model.add(Dense(10, activation=’softmax’))
model.compile(loss=’categorical_crossentropy’, optimizer=RMSprop())
model.fit()複製程式碼

Karpathy教授,將“智力上的好奇”或者“你可能想要晚些提升核心演算法”的論點抽象出來，認為計算實際上是一種洩漏抽象（譯者注：“抽象洩漏”是軟體開發時，本應隱藏實現細節的抽象化不可避免地暴露出底層細節與侷限性。抽象洩露是棘手的問題，因為抽象化本來目的就是向使用者隱藏不必要公開的細節–維基百科）：

“人們很容易陷入這樣的誤區中-認為你可以簡單地將任意的神經層組合在一起然後反向傳播演算法會‘令它們自己在你的資料上工作起來’。”

因此，我寫這篇文章的目的有兩層：

理解神經網路背後的抽象洩漏（通過在 Scratch 平臺上操作），而這些東西的重要性恰恰是我開始所忽略的。這樣如果我的模型沒有達到預期的學習效果，我可以更好地解決問題，而不是盲目地改變優化方案（甚至更換學習框架）。
一個深度神經網路（DNN），一旦被拆分成塊，對於 AI 領域之外的人們也再也不是一個黑箱了。相反，對於大多數有基本的統計背景的人來說，是一個個非常熟悉的話題的組合。我相信他們只需要學習很少的一些（只是那些如何將這一塊塊知識組合一起）知識就可以在一個全新的領域獲得不錯的洞察力。

從線性迴歸開始，藉著 R-notebook，通過解決一系列的數學和程式設計問題直至瞭解深度神經網路（DNN）。希望能夠藉此展示出來，你所需學習的新知識其實只有很少的一部分。

筆記

github.com/ilkarman/De…
github.com/ilkarman/De…
github.com/ilkarman/De…
github.com/ilkarman/De…

一、線性迴歸（見筆記(github-ipynb)）

在 R 中解決最小二乘法的計算器的閉包解決方案只需如下幾行：

# Matrix of explanatory variables
X <- as.matrix(X)
# Add column of 1s for intercept coefficient
intcpt <- rep(1, length(y))
# Combine predictors with intercept
X <- cbind(intcpt, X)
# OLS (closed-form solution)
beta_hat <- solve(t(X) %*% X) %*% t(X) %*% y複製程式碼

變數 beta_hat 所形成的向量包含的數值，定義了我們的“機器學習模型”。線性迴歸是用來預測一個連續的變數的（例如：這架飛機會延誤多久）。在預測分類的時候（例如：這架飛機會延誤嗎-會/不會），我們希望我們的預測能夠落在0到1之間，這樣我們可以將其轉換為各個種類的事件發生的可能性（根據所給的資料）。

當我們只有兩個互斥的結果時我們將使用一個二項邏輯迴歸。當候選結果（或者分類）多於兩個時，即多項互斥（例如：這架飛機延誤時間可能在5分鐘內、5-10分鐘或多於10分鐘），我們將使用多項邏輯迴歸（或者“Softmax 迴歸”）（譯者注：Softmax 函式是邏輯函式的一種推廣，更多知識見知乎）。在這種情況下許多類別不是互斥的（例如：這篇文章中的“R”，“神經網路”和“統計學”），我們可以採用二項式邏輯迴歸（譯者注：不是二項邏輯迴歸）。

另外，我們也可以用梯度下降（GD）這種迭代法來替代我們上文提到的閉包方法。整個過程如下：

從隨機地猜測權重開始
將所猜測的權重值代入損失函式中
將猜測值移向梯度的相反方向移動一小步（即我們所謂的“學習頻率”）
重複上述步驟 N 次

GD 僅僅使用了 Jacobian 矩陣 (而不是 Hessian 矩陣)，不過我們知道，當我們的損失函式為凸函式時，所有的極小值即（區域性最小值）為（全域性）最小值，因此 GD 總能夠收斂至全域性最小值。

線性迴歸中所用的損失函式是均方誤差函式：

要使用 GD 方法我們只需要找出 beta_hat 的偏導數（即 `delta`/梯度）

在 R 中實現方法如下：

# Start with a random guess
beta_hat <- matrix(0.1, nrow=ncol(X_mat))
  # Repeat below for N-iterations
  for (j in 1:N)
  {
    # Calculate the cost/error (y_guess - y_truth)
    residual <- (X_mat %*% beta_hat) - y
    # Calculate the gradient at that point
    delta <- (t(X_mat) %*% residual) * (1/nrow(X_mat))
    # Move guess in opposite direction of gradient
    beta_hat <- beta_hat - (lr*delta)
  }複製程式碼

200次的迭代之後我們會得到和閉包方法一樣的梯度與引數。除了這代表著我們的進步意外（我們使用了 GD），這個迭代方法在當閉包方法因矩陣過大，而無法計算矩陣的逆的時候，也非常有用（因為有記憶體的限制）。

第二步 – 邏輯迴歸 (見筆記(github-ipynb))

邏輯迴歸即一種用來解決二項分類的線性迴歸方法。它與標準的線性迴歸主要的兩種不同在於：

我們使用一種稱為 logistic-sigmoid 的 ‘啟用’/連結函式來將輸出壓縮至 0 到 1 的範圍內
不是最小化損失的方差而是最小化伯努利分佈的負對數似然

其它的都保持不變。

我們可以像這樣計算我們的啟用函式：

sigmoid <- function(z){1.0/(1.0+exp(-z))}複製程式碼

我們可以在 R 中這樣建立對數似然函式：

log_likelihood <- function(X_mat, y, beta_hat)
{
  scores <- X_mat %*% beta_hat
  ll <- (y * scores) - log(1+exp(scores))
  sum(ll)
}複製程式碼

這個損失函式（邏輯損失或對數損失函式）也叫做交叉熵損失。交叉熵損失根本上來講是對“意外”的一種測量，並且會成為所有接下來的模型的基礎，所以值得多花一些時間。

如果我們還像以前一樣建立最小平方損失函式，由於我們目前擁有的是一個非線性啟用函式（sigmoid），那麼損失函式將因不再是凸函式而使優化變得困難。

我們可以為兩個分類設立自己的損失函式。當 y=1 時，我們希望我們的損失函式值在預測值接近0的時候變得非常高，在接近1的時候變得非常低。當 y=0 時，我們所期望的與之前恰恰相反。這導致了我們有了如下的損失函式：

這裡的損失函式中的 delta 與我們之前的線性迴歸中的 delta 非常相似。唯一的不同在於我們在這裡將 sigmoid 函式也應用在了預測之中。這意味著邏輯迴歸中的梯度下降函式也會看起來很相似：

logistic_reg <- function(X, y, epochs, lr)
{
  X_mat <- cbind(1, X)
  beta_hat <- matrix(1, nrow=ncol(X_mat))
  for (j in 1:epochs)
  {
    # For a linear regression this was:
    # 1*(X_mat %*% beta_hat) - y
    residual <- sigmoid(X_mat %*% beta_hat) - y
    # Update weights with gradient descent
    delta <- t(X_mat) %*% as.matrix(residual, ncol=nrow(X_mat)) *  (1/nrow(X_mat))
    beta_hat <- beta_hat - (lr*delta)
  }
  # Print log-likliehood
  print(log_likelihood(X_mat, y, beta_hat))
  # Return
  beta_hat
}複製程式碼

三、Softmax 迴歸函式（無筆記）

邏輯迴歸的推廣即為多項邏輯迴歸（也稱為 ‘softmax 函式’），是對兩項以上的分類進行預測的。我尚未在 R 中建立這個例子，因為下一步的神經網路中也有一些東西簡化之後與之相似，然而為了完整起見，如果你仍然想要建立它的話，我還是要強調一下這裡主要的不同。

首先，我們不再用 sigmoid 函式來講我們所得的值壓縮在 0 至 1 之間：

我們用 softmax 函式來將 n 個值的和壓縮至 1：

這樣意味著每個類別所得的值，可以根據所給的條件，被轉化為該類的概率。同時也意味著當我們希望提高某一分類的權重來提高它所獲得的概率的時候，其它分類的出現概率會有所下降。也就是說，我們的各個類別是互斥的。

其次，我們使用一個更加通用的交叉熵損失函式：

要想知道為什麼-記住對於二項分類（如之前的例子）我們有兩個類別：j = 2，在每個類別是互斥的，a1 + a2 = 1 且 y 是一位有效編碼（one-hot）所以 y1+y2=1，我們可以將通用公式重寫為：
（譯者注：one-hot是將分類的特徵轉化為更加適合分類和迴歸演算法的資料格式（Quora-Håkon Hapnes Strand），中文資料可見此）

這與我們剛開始的等式是相同的。然而，我們現在將 j=2 的條件放寬。這裡的交叉熵損失函式可以被看出來有著與二項分類的邏輯輸出的交叉熵有著相同的梯度。

然而，即使梯度有著相同的公式，也會因為啟用函式代入了不同的值而不一樣（用了 softmax 而不是邏輯中的 sigmoid）。

在大多數的深度學習框架中，你可以選擇‘二項交叉熵（binary_crossentropy）’或者‘分類交叉熵（categorical_crossentropy）’損失函式。這取決於你的最後一層神經包含的是 sigmoid 還是 softmax 啟用函式，相對應著，你可以選擇‘二項交叉熵（binary_crossentropy）’或者‘分類交叉熵（categorical_crossentropy）’。而由於梯度相同，神經網路的訓練並不會被影響，然而所得到的損失（或評測值）會由於搞混它們而錯誤。

之所以要涉及到 softmax 是因為大多數的神經網路，會在各個類別互斥的時候，用 softmax 層作為最後一層（讀出層），用多項交叉熵（也叫分類交叉熵）損失函式，而不是用 sigmoid 函式搭配二項交叉熵損失函式。儘管多項 sigmoid 也可以用於多類別分類（並且會被用於下個例子中），但這總體上僅用於多項不互斥的時候。有了 softmax 作為輸出，由於輸出的和被限制為 1，我們可以直接將輸出轉化為概率。

四、神經網路（見筆記(github-ipynb))）

一個神經網路可以被看作為一系列的邏輯迴歸堆疊在一起。這意味著我們可以說，一個邏輯迴歸實際上是一個（帶有 sigmoid 啟用函式）無隱藏層的神經網路。

隱藏層，使神經網路具有非線性且導致了用於通用近似定理所描述的特性。該定理宣告，一個神經網路和一個隱藏層可以逼近任何線性或非線性的函式。而隱藏層的數量可以擴充套件至上百層。

如果將神經網路看作兩個東西的結合會很有用：1）很多的邏輯迴歸堆疊在一起形成‘特徵生成器’ 2）一個 softmax 迴歸函式構成的單個讀出層。近來深度學習的成功可歸功於‘特徵生成器’。例如：在以前的計算機視覺領域，我們需要痛苦地宣告我們需要找到各種長方形，圓形，顏色和結合方式（與經濟學家們如何決定哪些相互作用需要用於線性迴歸中相似）。現在，隱藏層是對決定哪個特徵（哪個‘相互作用’）需要提取的優化器。很多的深度學習實際上是通過用一個訓練好的模型，去掉讀出層，然後用那些特徵作為輸入（或者是促進決策樹（boosted decision-trees））來生成的。

隱藏層同時也意味著我們的損失函式在引數中不是一個凸函式，我們不能夠通過一個平滑的山坡來到達底部。我們會用隨機梯度下降（SGD）而不是梯度下降（GD），不像我們之前在邏輯迴歸中做的一樣，這樣基本上在每一次小批量（mini-batch）（比觀察總數小很多）被在神經網路中傳播後都會重編觀察（隨機）並更新梯度。這裡有很多 SGD 的替代方法，Sebastian Ruder 為我們做了很多工作。我認為這確實是個迷人的話題，不過卻超出這篇博文所討論的範圍了，很遺憾。簡要來講，大多數優化方法是一階的（包括 SGD，Adam，RMSprop和 Adagrad）因為計算二階函式的計算難度過高。然而，一些一階方法有一個固定的學習頻率（SGD）而有一些擁有適應性學習頻率（Adam），這意味著我們通過成為損失函式所更新權重的‘數量’-將會在開始有巨大的變化而隨著我們接近目標而逐漸變小。

需要弄清楚的一點是，最小化訓練資料上的損失並非我們的主要目標-理論上我們希望最小化‘不可見的’（測試）資料的損失；因此所有的優化方法都代表著已經一種假設之下，即訓練資料的的低損失會以同樣的（損失）分佈推廣至‘新’的資料。這意味著我們可能更青睞於一個有著更高的訓練資料損失的神經網路；因為它在驗證資料上的損失很低（即那些未曾被用於訓練的資料）-我們則會說該神經網路在這種情況下‘過度擬合’了。這裡有一些近期的論文聲稱，他們發現了很多很尖的最小值點，所以適應性優化方法並不像 SGD 一樣能夠很好的推廣。（譯者注：即演算法在一些驗證資料中表現地出奇的差）

之前我們需要將梯度反向傳播一層，現在一樣，我們也需要將其反向傳播過所有的隱藏層。關於反向傳播演算法的解釋，已經超出了本文的範圍，然而理解這個演算法卻是十分必要的。這裡有一些不錯的資源可能對各位有所幫助。

我們現在可以在 Scratch 平臺上用 R 通過四個函式建立一個神經網路了。

我們首先初始化權重：

neuralnetwork <- function(sizes, training_data, epochs, mini_batch_size, lr, C, verbose=FALSE, validation_data=training_data)

由於我們將引數進行了複雜的結合，我們不能簡單地像以前一樣將它們初始化為 1 或 0，神經網路會因此而在計算過程中卡住。為了防止這種情況，我們採用高斯分佈（不過就像那些優化方法一樣，這也有許多其他的方法）：

biases <- lapply(seq_along(listb), function(idx){
    r <- listb[[idx]]
    matrix(rnorm(n=r), nrow=r, ncol=1)
    })

weights <- lapply(seq_along(listb), function(idx){
    c <- listw[[idx]]
    r <- listb[[idx]]
    matrix(rnorm(n=r*c), nrow=r, ncol=c)
    })複製程式碼

我們使用隨機梯度下降（SGD）作為我們的優化方法：

    SGD <- function(training_data, epochs, mini_batch_size, lr, C, sizes, num_layers, biases, weights,verbose=FALSE, validation_data)
    {
      # Every epoch
      for (j in 1:epochs){
    # Stochastic mini-batch (shuffle data)
    training_data <- sample(training_data)
    # Partition set into mini-batches
    mini_batches <- split(training_data,
      ceiling(seq_along(training_data)/mini_batch_size))
    # Feed forward (and back) all mini-batches
    for (k in 1:length(mini_batches)) {
      # Update biases and weights
      res <- update_mini_batch(mini_batches[[k]], lr, C, sizes, num_layers, biases, weights)
      biases <- res[[1]]
      weights <- res[[-1]]
    }
      }
      # Return trained biases and weights
      list(biases, weights)
    }複製程式碼

作為 SGD 方法的一部分，我們更新了

 update_mini_batch <- function(mini_batch, lr, C, sizes, num_layers, biases, weights)
 {
   nmb <- length(mini_batch)
   listw <- sizes[1:length(sizes)-1]
   listb <-  sizes[-1]

 # Initialise updates with zero vectors (for EACH mini-batch)
   nabla_b <- lapply(seq_along(listb), function(idx){
 r <- listb[[idx]]
 matrix(0, nrow=r, ncol=1)
   })
   nabla_w <- lapply(seq_along(listb), function(idx){
 c <- listw[[idx]]
 r <- listb[[idx]]
 matrix(0, nrow=r, ncol=c)
   })

 # Go through mini_batch
   for (i in 1:nmb){
 x <- mini_batch[[i]][[1]]
 y <- mini_batch[[i]][[-1]]
 # Back propagation will return delta
 # Backprop for each observation in mini-batch
 delta_nablas <- backprop(x, y, C, sizes, num_layers, biases, weights)
 delta_nabla_b <- delta_nablas[[1]]
 delta_nabla_w <- delta_nablas[[-1]]
 # Add on deltas to nabla
 nabla_b <- lapply(seq_along(biases),function(j)
   unlist(nabla_b[[j]])+unlist(delta_nabla_b[[j]]))
 nabla_w <- lapply(seq_along(weights),function(j)
   unlist(nabla_w[[j]])+unlist(delta_nabla_w[[j]]))
   }
   # After mini-batch has finished update biases and weights:
   # i.e. weights = weights - (learning-rate/numbr in batch)*nabla_weights
   # Opposite direction of gradient
   weights <- lapply(seq_along(weights), function(j)
 unlist(weights[[j]])-(lr/nmb)*unlist(nabla_w[[j]]))
   biases <- lapply(seq_along(biases), function(j)
 unlist(biases[[j]])-(lr/nmb)*unlist(nabla_b[[j]]))
   # Return
   list(biases, weights)
 }複製程式碼

我們用來計算 delta 的演算法是反向傳播演算法。

在這個例子中我們使用交叉熵損失函式，產生了以下的梯度：

cost_delta <- function(method, z, a, y) {if (method==`ce`){return (a-y)}}複製程式碼

同時，為了與我們的邏輯迴歸例子保持連續，我們在隱藏層和讀出層上使用 sigmoid 啟用函式：

# Calculate activation function
    sigmoid <- function(z){1.0/(1.0+exp(-z))}
    # Partial derivative of activation function
    sigmoid_prime <- function(z){sigmoid(z)*(1-sigmoid(z))}複製程式碼

如之前所說，一般來講 softmax 啟用函式適用於讀出層。對於隱藏層，線性整流函式（ReLU）更加地普遍，這裡就是最大值函式（負數被看作為0）。隱藏層使用的啟用函式可以被想象為一場扛著火焰同時保持它（梯度）不滅的比賽。sigmoid 函式在0和1處平坦化，成為一個平坦的梯度，相當於火焰的熄滅（我們失去了訊號）。而線性整流函式（ReLU）幫助儲存了這個梯度。

反向傳播函式被定義為：

backprop <- function(x, y, C, sizes, num_layers, biases, weights)複製程式碼

請在筆記中檢視完整的程式碼-然而原則還是一樣的：我們有一個正向傳播，使得我們在網路中將權重傳導過所有神經層，併產生預測值。然後將預測值代入損失梯度函式中並將所有神經層中的權重更新。

這總結了神經網路的建成（搭配上你所需要的儘可能多的隱藏層）。將隱藏層的啟用函式換為 ReLU
函式，讀出層換為 softmax 函式，並且加上 L1 和 L2 的歸一化，是一個不錯的練習。把它在筆記中的 iris 資料集跑一遍，只用一個隱藏層，包含40個神經元，我們就可以在大概30多回合訓練後得到一個96%精確度的神經網路。

筆記中還提供了一個100個神經元的手寫識別系統的例子，來根據28*28畫素的影像預測數字。

五、卷積神經網路（[見筆記(github.com/ilkarman/De…)]）

在這裡，我們只會簡單地測試卷積神經網路（CNN）中的正向傳播。CNN 首次受到關注是因為1998年的LeCun的精品論文。自此之後，CNN 被證實是在影像、聲音、視訊甚至文字中最好的演算法。

影像識別開始時是一個手動的過程，研究者們需要明確影像的哪些位元（特徵）對於識別有用。例如，如果我們希望將一張圖片歸類進‘貓’或‘籃球’，我們可以寫一些程式碼提取出顏色（如籃球是棕色）和形狀（貓有著三角形耳朵）。這樣我們或許就可以在這些特徵上跑一個線性迴歸，來得到三角形個數和影像是貓還是樹的關係。這個方法很受圖片的大小、角度、質量和光線的影響，有很多問題。規模不變的特徵變換(SIFT) 在此基礎上做了大幅提升並曾被用來對一個物體提供‘特徵描述’，這樣可以被用來訓練線性迴歸（或其他的關係型學習器）。然而，這個方法有個一成不變的規則使其不能被為特定的領域而優化。

CNN 卷積神經網路用一種很有趣的方式看待影像（提取特徵）。開始時，他們只觀察影像的很小一部分（每次），比如說一個大小為 5*5 畫素的框（一個過濾器）。2D 用於影像的卷積，是將這個框掃遍整個影像。這個階段會專門用於提取顏色和線段。然而，下一個神經層會轉而關注之前過濾器的結合，因而‘放大來觀察’。在一定數量的層數之後，神經網路會放的足夠大而能識別出形狀和更大的結構。

這些過濾器最終會成為神經網路需要去學習、識別的‘特徵’。接著，它就可以通過統計各個特徵的數量來識別其與影像標籤（如‘籃球’或‘貓’）的關係。這個方法看起來對圖片來講很自然-因為它們可以被拆成小塊來描述（它們的顏色，紋理等）。CNN 看起來在影像分形特徵分析方面會蓬勃發展。這也意味著它們不一定適合其他形式的資料，如 excel 工作單中就沒有固有的樣式：我們可以改變任意幾列的順序而資料還是一樣的——不過在影像中交換畫素點的位置就會導致影像的改變。

在之前的例子中我們觀察的是一個標準的神經網路對手寫字型的歸類。在神經網路中的 i 層的每個神經元，與 j 層的每個神經元相連-我們所框中的是整個影像（譯者注：與 CNN 之前的 5*5 畫素的框不同）。這意味著如果我們學習了數字 2 的樣子，我們可能無法在它被錯誤地顛倒的時候識別出來，因為我們只見過它正的樣子。CNN 在觀察數字 2 的小的位元時並且在比較樣式的時候有很大的優勢。這意味著很多被提取出的特徵對各種旋轉，歪斜等是免疫的（譯者注：即適用於所有變形）。對於更多的細節，Brandon 在這裡解釋了什麼是真正的 CNN。

我們在 R 中如此定義 2D 卷積函式：

convolution <- function(input_img, filter, show=TRUE, out=FALSE)
{
  conv_out <- outer(
    1:(nrow(input_img)-kernel_size[[1]]+1),
    1:(ncol(input_img)-kernel_size[[2]]+1),
    Vectorize(function(r,c) sum(input_img[r:(r+kernel_size[[1]]-1),
                                          c:(c+kernel_size[[2]]-1)]*filter))
  )
}複製程式碼

並用它對一個圖片應用了一個 3*3 的過濾器：

conv_emboss <- matrix(c(2,0,0,0,-1,0,0,0,-1), nrow = 3)
convolution(input_img = r_img, filter = conv_emboss)複製程式碼

你可以檢視筆記來看結果，然而這看起來是從圖片中提取線段。否則，卷積可以‘銳化’一張圖片，就像一個3*3的過濾器：

conv_sharpen <- matrix(c(0,-1,0,-1,5,-1,0,-1,0), nrow = 3)
convolution(input_img = r_img, filter = conv_sharpen)複製程式碼

很顯然我們可以隨機地隨機地初始化一些個數的過濾器（如：64個）：

filter_map <- lapply(X=c(1:64), FUN=function(x){
    # Random matrix of 0, 1, -1
    conv_rand <- matrix(sample.int(3, size=9, replace = TRUE), ncol=3)-2
    convolution(input_img = r_img, filter = conv_rand, show=FALSE, out=TRUE)
})複製程式碼

我們可以用以下的函式視覺化這個 map：

square_stack_lst_of_matricies <- function(lst)
{
    sqr_size <- sqrt(length(lst))
    # Stack vertically
    cols <- do.call(cbind, lst)
    # Split to another dim
    dim(cols) <- c(dim(filter_map[[1]])[[1]],
                   dim(filter_map[[1]])[[1]]*sqr_size,
                   sqr_size)
    # Stack horizontally
    do.call(rbind, lapply(1:dim(cols)[3], function(i) cols[, , i]))
}複製程式碼