只用NumPy實現神經網路
作者：Piotr Skalski
編譯：weakish
【編者按】和Piotr Skalski一起，基於NumPy手寫神經網路，通過親自動手實踐，加深對神經網路內部機制的理解。
Keras、TensorFlow、PyTorch等高層框架讓我們可以快速搭建複雜模型。然而，花一點時間瞭解下底層概念是值得的。前不久我發過一篇文章，以簡單的方式解釋了神經網路是如何工作的。但這是一篇高度理論性的文章，主要以數學為主（數學是神經網路超能力的來源）。我打算以實踐的方式繼續這一主題，在這篇文章中，我們將嘗試僅僅基於NumPy建立一個可以工作的神經網路。最後，我們將測試一下建立的模型——用它來解決簡單的分類問題，並和基於Keras搭建的神經網路比較一下。
顯然，今天的文章會包含很多Python程式碼，我希望你不會因此覺得枯燥。我同時也把文中所有的程式碼發在GitHub上： ofollow,noindex">SkalskiP/ILe arnDeepLearning.py
概覽
在開始程式設計之前，先讓我們準備一份基本的路線圖。我們的目標是建立一個特定架構（層數、層大小、啟用函式）的密集連線神經網路。然後訓練這一神經網路並做出預測。
上面的示意圖展示了訓練網路時進行的操作，以及單次迭代不同階段需要更新和讀取的引數。
初始化神經網路層
讓我們從初始化每一層的權重矩陣 W 和偏置向量 b 開始。下圖展示了網路層l的權重矩陣和偏置向量，其中，上標 [l] 表示當前層的索引，n表示給定層中的神經元數量。
我們的程式也將以類似的列表形式描述神經網路架構。列表的每一項是一個字典，描述單個網路層的基本引數： input_dim 是網路層輸入的訊號向量的大小， output_dim 是網路層輸出的啟用向量的大小， activation 是網路層所用的啟用函式。
nn_architecture = [
{"input_dim": 2, "output_dim": 4, "activation": "relu"},
{"input_dim": 4, "output_dim": 6, "activation": "relu"},
{"input_dim": 6, "output_dim": 6, "activation": "relu"},
{"input_dim": 6, "output_dim": 4, "activation": "relu"},
{"input_dim": 4, "output_dim": 1, "activation": "sigmoid"},
]
如果你熟悉這一主題，你的腦海中大概已經迴盪起焦急的聲音：“喂喂！搞錯了！這裡有些是不必要的……”這一次，你內心的聲音是對的，一個網路層的輸出向量同時也是下一層的輸入，所以其實只用列出兩者之一就夠了。不過，我故意使用這樣的表示法，使每一層的格式保持一致，也讓初次接觸這一主題的人更容易理解程式碼。
def init_layers(nn_architecture, seed = 99):
np.random.seed(seed)
number_of_layers = len(nn_architecture)
params_values = {}
for idx, layer in enumerate(nn_architecture):
layer_idx = idx + 1
layer_input_size = layer["input_dim"]
layer_output_size = layer["output_dim"]
params_values['W' + str(layer_idx)] = np.random.randn(
layer_output_size, layer_input_size) * 0.1
params_values['b' + str(layer_idx)] = np.random.randn(
layer_output_size, 1) * 0.1
return params_values
上面的程式碼初始化了網路層的引數。注意我們用隨機的小數字填充矩陣 W 和向量 b 。這並不是偶然的。權重值無法使用相同的數字初始化，否則會造成 破壞性的對稱問題 。 基本上，如果權重都一樣，不管輸入X是什麼，隱藏層的所有單元也都一樣 。這樣，我們就會陷入初始狀態，不管訓練多久，網路多深，都無望擺脫。線性代數不會原諒我們。
小數值增加了演算法的效率。我們可以看看下面的sigmoid函式影象，大數值處的函式影象幾乎是扁平的，這會對神經網路的學習速度造成顯著影響。所有引數使用小隨機數是一個簡單的方法，但它保證了演算法有一個 足夠好 的開始。
啟用函式
啟用函式只需一行程式碼就可以定義，但它們給神經網路帶來了非線性和所需的表達力。 “沒有它們，神經網路將變成線性函式的組合，也就是單個線性函式。” 啟用函式有很多種，但在這個專案中，我決定使用其中兩種——sigmoid和ReLU。為了同時支援前向傳播和反向傳播，我們還需要準備好它們的導數。
def sigmoid(Z): return 1/(1+np.exp(-Z)) def relu(Z): return np.maximum(0,Z) def sigmoid_backward(dA, Z): sig = sigmoid(Z) return dA * sig * (1 - sig) def relu_backward(dA, Z): dZ = np.array(dA, copy = True) dZ[Z <= 0] = 0; return dZ;
前向傳播
我們設計的神經網路有一個簡單的架構。輸入矩陣 X 傳入網路，沿著隱藏單元傳播，最終得到預測向量 Y_hat 。為了讓程式碼更易讀，我將前向傳播拆分成兩個函式——單層前向傳播，和整個神經網路前向傳播。
def single_layer_forward_propagation(A_prev, W_curr, b_curr, activation="relu"):
Z_curr = np.dot(W_curr, A_prev) + b_curr
if activation is "relu":
activation_func = relu
elif activation is "sigmoid":
activation_func = sigmoid
else:
raise Exception('Non-supported activation function')
return activation_func(Z_curr), Z_curr
這部分程式碼大概是最直接，最容易理解的。給定來自上一層的輸入訊號，我們計算仿射變換 Z ，接著應用選中的啟用函式。基於NumPy，我們可以對整個網路層和整批樣本一下子進行矩陣操作，無需迭代，這大大加速了計算。除了計算結果外，函式還返回了一個反向傳播時需要用到的中間值 Z 。
基於單層前向傳播函式，編寫整個前向傳播步驟很容易。這是一個略微複雜一點的函式，它的角色不僅是進行預測，還包括組織中間值。
def full_forward_propagation(X, params_values, nn_architecture):
memory = {}
A_curr = X
for idx, layer in enumerate(nn_architecture):
layer_idx = idx + 1
A_prev = A_curr
activ_function_curr = layer["activation"]
W_curr = params_values["W" + str(layer_idx)]
b_curr = params_values["b" + str(layer_idx)]
A_curr, Z_curr = single_layer_forward_propagation(A_prev, W_curr, b_curr, activ_function_curr)
memory["A" + str(idx)] = A_prev
memory["Z" + str(layer_idx)] = Z_curr
return A_curr, memory
損失函式
損失函式可以監測進展，確保我們向著正確的方向移動。 “一般來說，損失函式是為了顯示我們離‘理想’解答還有多遠。” 損失函式根據我們計劃解決的問題而選用，Keras之類的框架提供了很多選項。因為我計劃將神經網路用於二元分類問題，我決定使用交叉熵：
為了取得更多關於學習過程的資訊，我決定另外實現一個計算精確度的函式。
def get_cost_value(Y_hat, Y): m = Y_hat.shape[1] cost = -1 / m * (np.dot(Y, np.log(Y_hat).T) + np.dot(1 - Y, np.log(1 - Y_hat).T)) return np.squeeze(cost) def get_accuracy_value(Y_hat, Y): Y_hat_ = convert_prob_into_class(Y_hat) return (Y_hat_ == Y).all(axis=0).mean()
反向傳播
不幸的是，很多缺乏經驗的深度學習愛好者都覺得反向傳播很嚇人，難以理解。微積分和線性代數的組合經常會嚇退那些沒有經過紮實的數學訓練的人。所以不要過於擔心你現在還不能理解這一切。相信我，我們都經歷過這個過程。
def single_layer_backward_propagation(dA_curr, W_curr, b_curr, Z_curr, A_prev, activation="relu"):
m = A_prev.shape[1]
if activation is "relu":
backward_activation_func = relu_backward
elif activation is "sigmoid":
backward_activation_func = sigmoid_backward
else:
raise Exception('Non-supported activation function')
dZ_curr = backward_activation_func(dA_curr, Z_curr)
dW_curr = np.dot(dZ_curr, A_prev.T) / m
db_curr = np.sum(dZ_curr, axis=1, keepdims=True) / m
dA_prev = np.dot(W_curr.T, dZ_curr)
return dA_prev, dW_curr, db_curr
人們經常搞混反向傳播和梯度下降，但事實上它們不一樣。前者是為了高效地計算梯度，後者則是為了基於計算出的梯度進行優化。在神經網路中，我們計算損失函式在引數上的梯度，但反向傳播可以用來計算任何函式的導數。 反向傳播演算法的精髓在於遞迴地使用求導的鏈式法則，通過組合導數已知的函式，計算函式的導數 。下面的公式描述了單個網路層上的反向傳播過程。由於本文的重點在實際實現，所以我將省略求導過程。從公式上我們可以很明顯地看到，為什麼我們需要在前向傳播時記住中間層的 A 、 Z 矩陣的值。
和前向傳播一樣，我決定將計算拆分成兩個函式。之前給出的是單個網路層的反向傳播函式，基本上就是以NumPy方式重寫上面的數學公式。而定義完整反向傳播過程的函式，主要是讀取、更新三個字典中的值。
def full_backward_propagation(Y_hat, Y, memory, params_values, nn_architecture):
grads_values = {}
m = Y.shape[1]
Y = Y.reshape(Y_hat.shape)
dA_prev = - (np.divide(Y, Y_hat) - np.divide(1 - Y, 1 - Y_hat));
for layer_idx_prev, layer in reversed(list(enumerate(nn_architecture))):
layer_idx_curr = layer_idx_prev + 1
activ_function_curr = layer["activation"]
dA_curr = dA_prev
A_prev = memory["A" + str(layer_idx_prev)]
Z_curr = memory["Z" + str(layer_idx_curr)]
W_curr = params_values["W" + str(layer_idx_curr)]
b_curr = params_values["b" + str(layer_idx_curr)]
dA_prev, dW_curr, db_curr = single_layer_backward_propagation(
dA_curr, W_curr, b_curr, Z_curr, A_prev, activ_function_curr)
grads_values["dW" + str(layer_idx_curr)] = dW_curr
grads_values["db" + str(layer_idx_curr)] = db_curr
return grads_values
基於單個網路層的反向傳播函式，我們從最後一層開始迭代計算所有引數上的導數，並最終返回包含所需梯度的python字典。
更新引數值
反向傳播是為了計算梯度，以根據梯度進行優化，更新網路的引數值。為了完成這一任務，我們將使用兩個字典作為函式引數： params_values ，其中儲存了當前引數值； grads_values ，其中儲存了用於更新引數值所需的梯度資訊。現在我們只需在每個網路層上應用以下等式即可。這是一個非常簡單的優化演算法，但我決定使用它作為更高階的優化演算法的起點（大概會是我下一篇文章的主題）。
Line"/>
def update(params_values, grads_values, nn_architecture, learning_rate): for idx, layer in enumerate(nn_architecture): layer_idx = idx + 1 params_values["W" + str(layer_idx)] -= learning_rate * grads_values["dW" + str(layer_idx)] params_values["b" + str(layer_idx)] -= learning_rate * grads_values["db" + str(layer_idx)] return params_values;
整合一切
萬事俱備只欠東風。最困難的部分已經完成了——我們已經準備好了所需的函式，現在只需以正確的順序把它們放到一起。
def train(X, Y, nn_architecture, epochs, learning_rate): params_values = init_layers(nn_architecture, 2) cost_history = [] accuracy_history = [] for i in range(epochs): Y_hat, cashe = full_forward_propagation(X, params_values, nn_architecture) cost = get_cost_value(Y_hat, Y) cost_history.append(cost) accuracy = get_accuracy_value(Y_hat, Y) accuracy_history.append(accuracy) grads_values = full_backward_propagation(Y_hat, Y, cashe, params_values, nn_architecture) params_values = update(params_values, grads_values, nn_architecture, learning_rate) return params_values, cost_history, accuracy_history
大衛對戰歌利亞
該是看看我們的模型能不能解決一個簡單的分類問題的時候了。我生成了一個包含兩個分類的資料點的資料集，如下圖所示。
作為對比，我用高層框架Keras搭建了一個模型。兩個模型採用相同的架構和學習率。不過這仍然是一場不公平的較量，因為我們的模型用的都是最簡單的方案。最終，基於NumPy和Keras的模型在測試集上都取得了類似的精確度（95%）。不過，我們的模型收斂的速度要慢很多倍。在我看來，這主要是因為我們的模型缺乏恰當的優化演算法。
https://www.zhihu.com/video/1036214027331997696
再見
希望我的文章拓展了你的視野，增加了你對神經網路內部運作機制的理解——這將是對我投入精力撰寫本文的最好回報。我在編寫程式碼和說明的時候也學到了很多東西，親自動手實踐能讓你學到很多東西。