基於numpy的前饋神經網路(feedforwardneuralnetwork)

***
程式碼部分可以直接通過Jupyter Notebook來檢視

---

這幾天在上Andrew Ng教授開的Coursera系列課程Deep Learning，總覺得光是看視訊和做作業還不夠，還是得自己動手寫寫程式碼，親自實現課程裡提到的演算法內容，於是便有了這篇部落格，作為自己入門深度學習的里程碑吧。

前饋神經網路

機器學習有兩個基本問題，一是迴歸，二是分類，神經網路大多用於解決分類問題，前饋神經網路（feedforward neural network）是整個神經網路家族中較為常見和較為基礎的一種，如下圖右上角的DFF所示。圖片來源是Cheat Sheets for AI, Neural Networks, Machine Learning, Deep Learning & Big Data。

神經網路中的基本元素是神經元，每層都有一定數量的神經元，神經元組合的多樣性決定了神經網路的豐富性。下面是一個簡單的前饋神經網路，總共有三層，從左到右分別是輸入層、隱層和輸出層，輸入層的x1和x2表示這個樣本只有兩個特徵（自變數），因為輸入層通常不計入內，所以這是一個兩層的神經網路，第一層有4個神經元，第二層只有1個。注意，隱層可以不止一層，隱層設定得越多，整個神經網路越龐大。

這個神經網路的工作原理是，給定一個樣本的資料，將資料傳輸到第一層，進行線性變換和啟用變換，得到加工過後的資料，這份新資料傳到第二層，作為第二層的輸入，接著進行線性變換和啟用變換，又得到一份新的資料，因為第二層是最後一層了（如果不止兩層就一直進行這樣的操作直到抵達最後一層為止），所以最終的輸出作為我們對該樣本的預測值y_hat。

每一個神經元如同工廠的流水車間的機器，它重複做著四件事情：【接受上一層資料作為輸入>>線性變換>>啟用變換>>輸出資料到下一層】，每個神經元中有三個組成部分：權重(weight)矩陣W，偏置(bias)向量b，以及啟用函式(activation function) g，用公式表達為下圖，其中上標(i)表示這是第i個樣本資料，上標[1]和[2]分別表示這是神經網路的第一層與第二層：

$z^{[1](i)} = W^{[1]} x^{(i)} + b^{[1]} $
$a^{[1](i)} = anh(z^{[1](i)})$
$z^{[2](i)} = W^{[2]} a^{[1](i)} + b^{[2]}$
$hat{y}^{(i)} = a^{[2](i)} = sigma(z^{ [2](i)})$
$y^{(i)}_{prediction} = egin{cases} 1 & mbox{if } a^{[2](i)} > 0.5 \ 0 & mbox{otherwise } end{cases}$

公式(1)中，通過簡單的線性變換得到了z[1]，z稱為prev_activation，接著通過啟用變換（這裡用的啟用函式是tanh函式，下面會講到），得到了a[1]，稱為activation，公式(3)(4)表達的是第二層的線性變換和啟用變換，和第一層大同小異，只不過第二層的啟用函式用的不是tanh函式而是sigmoid函式。

啟用函式有下面三種，都是執行了非線性變換，實現的效果都是將prev_activation轉化為activation。

$ sigmoid(z) = frac{1}{1+e^{-z}} $
$ tanh(z) = frac{e^{z} – e^{-z}}{e^{z} + e^{-z}} = 2 sigmoid(2z) – 1 $
$
relu(z) =
max(0, z) =
begin{cases}

z, & z>0 \

0, & z leq 0
end{cases}
$

每一個神經元都可以從上面三種啟用函式中選取一個作為自己的啟用函式g，經驗表明，使用tanh函式的效果總是碾壓使用sigmoid函式，所以人們大多使用tanh作為啟用函式，近年來人們發現了relu函式，發現它的效能比tanh更好，relu成為了廣受歡迎的啟用函式。既然sigmoid效能最差，為什麼還要介紹它？在一開始的時候說到，神經網路通常用於分類，比方說，給定一張圖片，去識別預測它是不是一隻貓：

我們的返回值應該是範圍在0~1之間的概率值，sigmoid的函式範圍是(0, 1), tanh的範圍是(-1, 1), relu的是[0, +∞)，使用sigmoid顯然更合適些。所以通常一個神經網路的配置是，中間的隱層的所有神經元使用tanh或者relu作為啟用函式，輸出層的神經元使用sigmoid。

之所以要線上性變換之後進行非線性變化，是因為，如果沒有非線性變換，純粹使用線性變化的話，不管使用了多少層的線性變換，最終的結果通過合併同類項之後仍然是線性變換，100層的神經網路和1層的沒有任何差別。神經網路從本質上來說就是一系列的非線性變換。

講完了基本概念，我們現在的問題是，如何訓練這個神經網路。

和其他的監督學習方法一樣，神經網路同樣是需要定義一個損失函式(cost function)，通過對它進行最小化，得到最優的引數W和b，即每個神經元的權重矩陣和偏置向量。神經網路用的損失函式是交叉熵：

$ J = -frac{1}{m} sum_{i = 1}^{m}{ [ y^{(i)}log(hat{y}^{(i)}) + (1-y^{(i)}) log(1-hat{y}^{(i)}) ] } $

其中m為樣本總數，y_hat為最終預測值，是通過每一層的神經元進行層層加工處理得到的，L表示整個神經網路的層數，即[L]表示最後一層：

$ Z^{[1]} = W^{[1]}X + b^{[1]} $
$ A^{[1]} = g^{[1]}(Z^{[1]}) $
$ Z^{[2]} = W^{[2]}X + b^{[2]} $
$ A^{[2]} = g^{[2]}(Z^{[2]}) $
$ … $
$ A^{[L]} = g^{[L]}(Z^{[L]}) = hat{Y}$

優化損失函式的方法是梯度下降法，梯度下降的原理這裡就不說了，就貼一下引數的更新公式：

$ W^{[l]} = W^{[l]} – alpha ext{ } dW^{[l]}$
$ b^{[l]} = b^{[l]} – alpha ext{ } db^{[l]}$
$ alpha:learning ext{ }rate $
$ dW ext{ | } db: gradients$

以前在看書的時候，看到神經網路部分，總是在說“神經網路訓練使用的是反向傳播(back propagation)演算法”，一直搞不明白這四個字究竟是什麼意思，現在總算是有些理解了。

既然有反向傳播，就有正向傳播，所謂正向，就是沿著神經網路從輸入層到輸出層，就是沿著每個神經元執行【接受上一層資料作為輸入>>線性變換>>啟用變換>>輸出資料到下一層】的方向，就是下圖中從左到右的方向，而反向傳播，指的是從輸出層到輸入層，即從右到左。

所謂的傳播，其實就是把整個神經網路運算的過程抽象成了流程圖的形式。正向傳播時，資料通過第一層神經元，被處理成了新的資料，傳輸到下一層，再加工成新資料，再傳到下一層，這便是資料的傳播，而在反向傳播時，傳播的是梯度，從最後一層（輸出層）出發，首先根據正向傳播得到的預測值，計算得出最後一個神經元處的A(activation)的梯度，進而得到Z(prev_activation)的梯度，根據公式再算出W, b以及倒數第二層的A的梯度，就像多米諾骨牌，不斷進行下去，直到將梯度傳播到第一層的神經網路處。

在正向傳播時，每個神經元執行了兩步操作，線性變化和啟用變換，線性變換是由W[l], b[l], A[l-1]得到了z[l]，啟用變化是由z[l]得到了A[l]；而在反向傳播時，在每個神經元中，根據相應的啟用函式由dA[l]得到dz[l]，再由dz[l]得到dW[l], db[l], dA[l-1]，這可以看做是線性變化和啟用變換的反向操作。

由此得到了每一層的神經元的W和b的梯度，引數W和b也由此得以更新，反向傳播時傳播的是梯度，也可以理解為傳播的是誤差，是當前模型的不足，引數W和b藉助這個資訊來進行修正，從而優化模型。於是整個神經網路就按照下圖的流程進行迭代，直到達到理想的模型效果：

所以整個神經網路的操作流程是這樣的：

• 確定神經網路的層數和每層的神經元數
• 初始化引數W和b

• 進行迴圈迭代

  • 前向傳播，得到預測值
  • 計算損失函式
  • 反向傳播，得到引數W和b的梯度
  • 對引數W和b進行更新
• 對新資料做預測

---

程式碼實現

課程裡有給出前饋神經網路的python程式碼，但是我覺得寫得過於繁瑣（簡單的幾步操作居然寫了300多行？），太注重形式，太刻意地去遵循程式碼規範，反倒增加了閱讀負擔，不太合我心意。

於是我按照自己對演算法的理解寫了一份程式碼，進行了精簡，只保留最重要的部分，並且加入了L2正則化的部分。為了測試程式碼的正確性，這裡使用的神經網路的層數和神經元數和Coursera上的是一樣的，共有四層，每層的神經元數量分別是20/ 7/ 5/ 1，用的資料也是課程的資料。

下面進入程式碼部分。

首先要呼叫numpy。為什麼要使用numpy？因為不想使用一層又一層的for迴圈，使用向量化的計算能夠大幅度地降低運算時間。

準備好啟用函式，隱層統一使用relu，輸出層使用sigmoid：

def sigmoid(z):
    ```
    z為prev_activation, size為 nl * m
    ```
    return 1 / (1 + np.exp(-z))

def relu(z):
    ```
    z為prev_activation, size為 nl * m
    ```
    return np.maximum(0, z)

設定好神經網路的結構：

n0, m = X.shape
n1 = 20
n2 = 7
n3 = 5
n4 = 1
layers_dims = [n0, n1, n2, n3, n4] # [12288, 20, 7, 5, 1] 
L = len(layers_dims) - 1  # 4層神經網路，不計輸入層

構建神經網路模型：

##### neural network model
def neural_network(X, Y, learning_rate=0.01, num_iterations=2000, lambd=0):
    m = X.shape[1]
    ### initialize forward propagation
    param_w = [i for i in range(L+1)]
    param_b = [i for i in range(L+1)]
    np.random.seed(10)
    for l in range(1, L+1):
        if l < L:
            param_w[l] = np.random.randn(layers_dims[l], layers_dims[l - 1]) * np.sqrt(2 / layers_dims[l - 1])
        if l == L:
            param_w[l] = np.random.randn(layers_dims[l], layers_dims[l - 1]) * 0.01
        param_b[l] = np.zeros((layers_dims[l], 1))

    activations = [X, ] + [i for i in range(L)]
    prev_activations = [i for i in range(L+1)]

    dA = [i for i in range(L+1)]
    dz = [i for i in range(L+1)]
    dw = [i for i in range(L+1)]
    db = [i for i in range(L+1)]

    for i in range(num_iterations):
        ### forward propagation
        for l in range(1, L+1):
            prev_activations[l] = np.dot(param_w[l], activations[l-1]) + param_b[l]
            if l < L:
                activations[l] = relu(prev_activations[l])
            else:
                activations[l] = sigmoid(prev_activations[l])

        cross_entropy_cost = -1/m * (np.dot(np.log(activations[L]), Y.T) 
                                     + np.dot(np.log(1-activations[L]), 1-Y.T))
        regularization_cost = 0
        for l in range(1, L+1):
            regularization_cost += np.sum(np.square(param_w[l])) * lambd/(2*m)
        cost = cross_entropy_cost + regularization_cost

        ### initialize backward propagation
        dA[L] =  np.divide(1-Y, 1-activations[L]) - np.divide(Y, activations[L])
        assert dA[L].shape == (1, m)

        ### backward propagation
        for l in reversed(range(1, L+1)):
            if l == L:
                dz[l] = dA[l] * activations[l] * (1-activations[l])
            else:
                dz[l] = dA[l].copy()
                dz[l][prev_activations[l] <= 0] = 0

            dw[l] = 1/m * np.dot(dz[l], activations[l-1].T) + param_w[l] * lambd/m
            db[l] = 1/m * np.sum(dz[l], axis=1, keepdims=True)
            dA[l-1] = np.dot(param_w[l].T, dz[l])

            assert dz[l].shape == prev_activations[l].shape
            assert dw[l].shape == param_w[l].shape
            assert db[l].shape == param_b[l].shape
            assert dA[l-1].shape == activations[l-1].shape

            param_w[l] = param_w[l] - learning_rate * dw[l]
            param_b[l] = param_b[l] - learning_rate * db[l]

        if i % 100 == 0:
            print("cost after iteration {}: {}".format(i, cost))

Andrew Ng教授是用一個dict來儲存每層神經元的引數的，比如說，在呼叫第三層的引數W3和b3時，他的寫法是:parameters[`W` + str(3)]和parameters[`b` + str(3)]，這樣寫沒有錯誤，雖然直觀，但是很麻煩，我的做法是，分別使用list來儲存W和b，根據位置讀取，對應上面的就是param_w[3]和param_b[3]。

注意到，python（以及其他程式語言）是從0開始計數的，而不是從1開始，這意味著，當神經網路共有4層的時候，我的param_w和param_b的長度是5，我對activation、prev_activation以及各個梯度dA/dz/dw/db都是用長度為5的list來儲存的，而Ng教授記錄這些變數的長度都是4。

我認為我的寫法是更容易理解的寫法，因為整個神經網路包括輸入層在內一共有5層，但是因為習慣上我們不計輸入層，所以這是個4層網路，如果我們把輸入層稱為第0層，輸出層稱為第4層，沒有什麼問題，並且恰好符合了程式語言從0開始計數的習慣。所以在我的寫法下，activations[3]就表示第三層的輸出，param_b[2]表示第二層的偏置向量，不會有什麼誤解，很直觀。

而對於Ng教授的做法，獲取每個引數時是基於字串的，比如說grads["db" + str(4)]表示第四層的b的梯度，而在用for迴圈遍歷每一層神經元的時候，又是基於位置的，這麼一來，你就會在+1 和 -1之後迷失自我，即便你確保了自己沒有出錯，你也已經花費了不少精力來判斷這個地方究竟是該+1還是-1還是保持原樣。

下面是我在做Ng教授佈置的程式設計作業時，需要填寫的反向傳播的部分，要求我在裡面填入5行程式碼來完成，這導致了我在 l 和 l+1和 l-1三者之間糾結猶豫了很久才終於填寫正確。我個人認為不是一段user-friendly的程式碼。

for l in reversed(range(L - 1)):
    # lth layer: (RELU -> LINEAR) gradients.
    # Inputs: "grads["dA" + str(l + 1)], current_cache". 
    # Outputs: "grads["dA" + str(l)] , grads["dW" + str(l + 1)] , grads["db" + str(l + 1)]

    ### START CODE HERE ### (approx. 5 lines)
    current_cache = caches[l]
    dA_prev_temp, dW_temp, db_temp = linear_activation_backward(grads[`dA` + str(l + 1)], current_cache, activation="relu")
    grads["dA" + str(l)] = dA_prev_temp
    grads["dW" + str(l + 1)] = dW_temp
    grads["db" + str(l + 1)] = db_temp
    ### END CODE HERE ###

Ng教授的程式碼中，對於不同的子函式，層數L時而等於5，時而等於4，很是混亂，而我的層數L始終等於4，符合python從0開始計數的習慣，使得整個演算法寫起來很快，也很好理解，一看就明白當前迭代到了哪一層。

另外，在寫演算法的時候，要重點關注每個引數的size，注意看它是幾乘幾的矩陣，很多時候出bug都來自於引數的size弄錯了，下面是各引數的size表，裡面的12288和209表示我們使用的訓練集的資料有209個樣本（圖片），每個樣本有12288個特徵（每張圖片是64畫素*64畫素的，所以有64*64*3=12288）。

關於size的規律是，梯度和引數的size相同，同一層的W和dW的size相同，b和db的size相同，其實也很好理解，看引數更新的公式就知道了；另外，同一層的A/z/dA/dz的size都是相同的，這個要記住。可以這麼理解，原資料的size是p×m，即m個p維的樣本（注意在神經網路中，每一列代表一個樣本，和其他機器學習方法中每一行代表一個樣本是反過來的），經過每一層的神經元處理之後，得到了新的資料，可以理解為是降維的過程，由一開始的12288×209的資料，通過第一層後變為了20×209的資料，即把高維的12288個特徵濃縮為20個新特徵，再濃縮為7個特徵，最後得到了1×209的預測值。

最後是預測函式，將訓練好的神經網路用於一套新的資料上：

def predict(X_new, parameters, threshold=0.5):
    param_w = parameters["param_w"]
    param_b = parameters["param_b"]

    activations = [X_new, ] + [i for i in range(L)]
    prev_activations = [i for i in range(L + 1)]
    m = X_new.shape[1]

    for l in range(1, L + 1):
        prev_activations[l] = np.dot(param_w[l], activations[l - 1]) + param_b[l]
        if l < L:
            activations[l] = relu(prev_activations[l])
        else:
            activations[l] = sigmoid(prev_activations[l])
    prediction = (activations[L] > threshold).astype("int")
    return prediction

---

經驗與收穫

• 這是我第一次用python寫一個具體的演算法，雖然一直在python做資料分析，但是一直沒有寫過一個邏輯完整的、能應用於實際場合的演算法（自己動手寫完整演算法的經歷，在今天之前，是用R語言寫了逐步迴歸和隨機森林），今天算是填補了我在python上的這塊空白。
• 使用assert語句來確保每個引數的shape正確，能夠減少出現bug的機率
• numpy中有一些操作和函式是element-wise的，要留意

• 整個演算法沒一會兒就寫完了，但是訓練的時候cost一直沒有收斂，檢查了老半天，才發現是由兩個問題導致的：

  • relu的導函式寫錯了，返回值是0與1沒錯，但是我用來判斷的引數竟然是A而不是z，真是寫得莫名其妙；
  • 初始化不正確（我後來才發現初始化是神經網路的關鍵），對於權重W我一直使用的寫法是param_w[l] = np.random.randn(layers_dims[l], layers_dims[l - 1]) * 0.01，導致cost沒能實現收斂，應該是由vanishing gradients導致的，後來我使用了He初始化和Xavier初始化，即把上面的0.01改成np.sqrt(2 / layers_dims[l - 1])或者np.sqrt(1 / layers_dims[l - 1])，效果明顯。
• 怎麼知道自己的最終演算法寫對了？我和課程使用的是同一套資料，如果使用我的演算法，在相同的神經網路結構下，和Ng教授使用的演算法得到的精確度差不多，就說明沒問題了。
• 比起寫程式碼，理解程式碼背後的演算法內容更重要。