到目前為止，我們已經研究了梯度下降演算法、人工神經網路以及反向傳播演算法，他們各自肩負重任：

• 梯度下降演算法：機器自學習的演算法框架；
• 人工神經網路：“萬能函式”的形式表達；
• 反向傳播演算法：計算人工神經網路梯度下降的高效方法；

基於它們，我們已經具備了構建具有相當實用性的智慧程式的核心知識。它們來之不易，從上世紀40年代人工神經元問世，到80年代末反向傳播演算法被重新應用，歷經了近半個世紀。然而，實現它們並進行復雜的數字手寫體識別任務，只需要74行Python程式碼（忽略空行和註釋）。要知道如果採用程式設計的方法（非學習的方式）來挑戰這個任務，是相當艱難的。

本篇將分析這份Python程式碼“network.py”，它基於NumPy，在對50000張影象學習後，即能夠識別0~9手寫體數字，正確率達到95%以上。強烈建議暫時忘記TF，用心感受凝結了人類文明結晶的滄桑演算法。程式碼來自Micheal Nielsen的《Neural Networks and Deep Learning》，略有修改（格式或環境匹配）。

MNIST資料集

MNIST

早在1998年，在AT&T貝爾實驗室的Yann LeCun就開始使用人工神經網路挑戰數字手寫體識別，用於解決當時銀行支票以及郵局信件郵編自動識別的需求。而現今，數字手寫體識別，已經成了機器學習的入門實驗案例。演算法實驗使用最廣泛的資料集就是MNIST，由Yann LeCun提供下載。它包含了60000張訓練圖片集，以及10000張測試圖片集，最初來源於NIST（National Institute of Standards and Technology，美國國家標準與技術研究院）資料庫。

樣本影象

如上圖所示，MNIST中的影象是灰度影象，畫素值為0的表示白色，為1的表示黑色，中間值是各種灰色。每張樣本影象的大小是28x28，具有784個畫素。

訓練集與測試集

MNIST中的60000張訓練影象掃描自250個人的手寫樣本，他們一半是美國人口普查局的員工，一半是大學生。10000張測試影象來自另外250個人（儘管也是出自美國人口普查局和高校）。可是為什麼要這麼做呢？答案是為了泛化（Generalization）。

人們希望學習訓練集（training set）後獲得的模型，能夠識別出從未見過的樣本，這種能力就是泛化能力，通俗的說，就是舉一反三。人類大腦就具有相當好的泛化能力，一個兩歲小孩在見過少量的鴨子圖片後，即可辨認出他從未見過的各種形態的鴨子。

基於這種考慮，測試集（test set）不會參於模型的訓練，而是特意被留出以測試模型的泛化效能。周志華的西瓜書中有一個比方：如果讓學生複習的題目，就是考試的考題，那麼即便他們考了100分，也不能保證他們真的學會了。

標籤

Label

上圖中右側的部分，稱為標籤（Label），是和樣本資料中的每張圖片一一對應的，由人工進行標註。標籤是資料集必不可少的一部分。模型的訓練過程，就是不斷的使識別結果趨近於標籤的過程。基於標籤的學習，稱為有監督學習。

驗證集與超引數

來自Micheal Nielsen的程式碼，又把60000張訓練集進行了進一步的劃分，其中50000張作為訓練集，10000張作為驗證集（validation set）。所以程式碼使用MNIST資料集與Yann LeCun的是有些區別的。

模型的引數是由訓練資料自動調整的，其他不被學習演算法覆蓋的引數，比如神經網路中的學習率、隨機梯度下降演算法中的mini batch的大小等，它們都被稱為超引數。驗證集被劃分出來就是用於評估模型的泛化能力，並以此為依據優化超引數的。

這裡容易產生一個疑問：評估模型的泛化能力，不是測試集要做的事情嗎？

測試集的確是用於評估模型的泛化能力的，但是理想情況下是用於最終評測。也就是說，測試集產生的任何結果和反饋，都不應該用於改善模型，以避免模型對測試集產生過擬合。那麼從訓練集劃分出驗證集，就沒有這個限制了，一方面驗證集不參與訓練，可以評估模型的泛化能力，另一方面，可以從評估的結果來進一步改善模型的網路架構、超引數。

驗證資料不是MNIST規範的一部分，但是留出驗證資料已經成了一種預設的做法。

Python必知必會：張量構建

本篇使用與《Neural Networks and Deep Learning》示例程式碼一致的Python版本：

• Python 2.7.x，使用了conda建立了專用虛擬環境，具體方法參考1 Hello, TensorFlow；
• NumPy版本：1.13.1。

作為AI時代頭牌語言的Python，具有非常好的生態環境，其中數值演算法庫NumPy做矩陣操作、集合操作，基本都是“一刀斃命”。為了能順暢的分析接下來的Python程式碼，我挑選了1處程式碼重點看下，略作修改（前兩層神經元數量）可以單獨執行。

第1行：匯入numpy並啟用np作為別名。

第2行：是一個數組定義，其中包含了3個元素。

第3行：

• 先看 sizes[1:]，它表示sizes的一個子陣列，包含元素從原陣列的下標1開始，直到原陣列最後1個元素，它的值可以算出是[15, 10];
• 然後是NumPy的隨機數生成方法random.randn，它生成的隨機數，符合均值為0，標準差為1的標準正態分佈；
• random.randn方法的引數，描述生成張量的形狀，例如random.randn(2,3)會生成秩為2，形狀為shape[2,3]的張量，是一個矩陣：array([[-2.17399771, 0.20546498, -1.2405749 ], [-0.36701965, 0.12564214, 0.10203605]])，關於張量請參考2 TensorFlow核心基礎；
• 第3行整體來看才能感受到Python和NumPy的威力：方法引數的引數化，即呼叫randn方法時可傳入變數：randn(y, 1)，而變數y遍歷集合sizes[1:]，效果等同於[randn(15, 1), randn(10, 1)]；

第4行：

• 先看sizes[:-1]表示其包含的元素從原陣列的第1個開始，直到原陣列的最後1個的前一個（倒數第2個），此時sizes[:-1]是[8, 15]；
• 第4行randn的兩個引數都是變數y和x，此時出現的zip方法，限制了兩個變數是同步自增的，效果等同於[randn(15, 8), randn(10, 15)]。

矩陣與神經網路

分析了前面4行程式碼，我們知道了如何高效的定義矩陣，但是和神經網路的構建有什麼關係呢？下面給出網路結構與矩陣結構的對應關係。

3層感知器

上面的神經網路結構即可描述為：sizes = [8, 15, 10]，第一層輸入層8個神經元，第二層隱藏層15個神經元，第三層輸出層10個神經元。

第一層是輸入層，沒有權重和偏置。

第二層的權重和偏置為：

第2層神經元的權重和偏置

第三層的權重和偏置為：

第3層神經元的權重和偏置

回看第3行程式碼，其等價於[randn(15, 1), randn(10, 1)]，相當於把網路中的兩層偏置矩陣放在一起了：

回看第4行程式碼，其等價於[randn(15, 8), randn(10, 15)]，相當於把網路中的兩層權重矩陣放在一起了：

而這4個矩陣本身，就代表了想要構建的神經網路模型，它們中的元素，構成了神經網路的所有可學習引數（不包括超引數）。當明瞭了神經網路與矩陣群的對映關係，在你的腦中即可想象出資料在網路中的層層流動，直到最後的輸出的形態。

隨機梯度下降演算法框架

整個神經網路程式的骨架，就是梯度下降演算法本身，在network.py中，它被抽象成了一個單獨的函式SDG（Stochastic Gradient Descent）：

函式體的實現，非常清晰，有兩層迴圈組成。外層是資料集的迭代（epoch）；內層是隨機梯度下降演算法中小批量集合的迭代，每個批量（batch）都會計算一次梯度，進行一次全體引數的更新（一次更新就是一個step）：

BP

可以想象self.update_mini_batch(mini_batch, eta)中的主要任務就是獲得每個引數的偏導數，然後進行更新，求取偏導數的程式碼即：

delta_nabla_b，delta_nabla_w = self.backprop(x, y)

反向傳播演算法（BP）的實現也封裝成了函式backprop：

識別率

執行程式碼，在Python命令列輸入以下程式碼：

import mnist_loader
import network
training_data, validation_data, test_data = mnist_loader.load_data_wrapper()
net = network.Network([784, 30, 10])net.SGD(training_data, 30, 10, 3.0, test_data=test_data)

上面程式碼中的mnist_loader負責MNIST資料的讀取，這部分程式碼在這裡下載，為了適配資料集的相對路徑做了微調。

接下來，定義了一個3層的神經網路：

• 輸入層784個神經元（對應28x28的數字手寫體影象）；
• 隱藏層30個神經元；
• 輸出層10個神經元（對應10個手寫體數字）。

最後是梯度下降法的設定：

• epoch：30次；
• batch：10個樣本影象；
• 學習率：3.0。

程式碼開始執行，30次迭代學習後，識別準確率即可達到95%。這個識別率是未去逐個優化超引數，就能輕鬆得到的，可以把它當做一個基線水準，在此基礎上再去慢慢接近NN的極限（99.6%以上）。

執行結果如下：

識別準確率

附完整程式碼

"""
network.py
~~~~~~~~~~
A module to implement the stochastic gradient descent learning
algorithm for a feedforward neural network. Gradients are calculated
using backpropagation. Note that I have focused on making the code
simple, easily readable, and easily modifiable. It is not optimized,
and omits many desirable features.
"""
# Libraries
# Standard library
import random

# Third-party libraries
import numpy as np

class Network(object):
    def __init__(self, sizes):
        """The list ``sizes`` contains the number of neurons in the
        respective layers of the network. For example, if the list
        was [2, 3, 1] then it would be a three-layer network, with the
        first layer containing 2 neurons, the second layer 3 neurons,
        and the third layer 1 neuron. The biases and weights for the
        network are initialized randomly, using a Gaussian
        distribution with mean 0, and variance 1. Note that the first
        layer is assumed to be an input layer, and by convention we
        won't set any biases for those neurons, since biases are only
        ever used in computing the outputs from later layers."""

        self.num_layers = len(sizes)
        self.sizes = sizes
        self.biases = [np.random.randn(y, 1) for y in sizes[1:]]
        self.weights = [np.random.randn(y, x)
                        for x, y in zip(sizes[:-1], sizes[1:])]

    def feedforward(self, a):
        """Return the output of the network if ``a`` is input."""
        for b, w in zip(self.biases, self.weights):
            a = sigmoid(np.dot(w, a) + b)
        return a

    def SGD(self, training_data, epochs, mini_batch_size, eta,
            test_data=None):

        """Train the neural network using mini-batch stochastic
        gradient descent. The ``training_data`` is a list of tuples
        ``(x, y)`` representing the training inputs and the desired
        outputs. The other non-optional parameters are
        self-explanatory. If ``test_data`` is provided then the
        network will be evaluated against the test data after each
        epoch, and partial progress printed out. This is useful for
        tracking progress, but slows things down substantially."""
        if test_data:
            n_test = len(test_data)
        n = len(training_data)
        for j in range(epochs):
            random.shuffle(training_data)
            mini_batches = [
                training_data[k:k + mini_batch_size]
                for k in range(0, n, mini_batch_size)]
            for mini_batch in mini_batches:
                self.update_mini_batch(mini_batch, eta)
            if test_data:
                print("Epoch {0}: {1} / {2}".format(
                    j, self.evaluate(test_data), n_test))
            else:
                print("Epoch {0} complete".format(j))

    def update_mini_batch(self, mini_batch, eta):
        """Update the network's weights and biases by applying
        gradient descent using backpropagation to a single mini batch.
        The ``mini_batch`` is a list of tuples ``(x, y)``, and ``eta``
        is the learning rate."""
        nabla_b = [np.zeros(b.shape) for b in self.biases]
        nabla_w = [np.zeros(w.shape) for w in self.weights]
        for x, y in mini_batch:
            delta_nabla_b, delta_nabla_w = self.backprop(x, y)
         nabla_b = [nb + dnb for nb, dnb in zip(nabla_b, delta_nabla_b)]
         nabla_w = [nw + dnw for nw, dnw in zip(nabla_w, delta_nabla_w)]
       self.weights = [w - (eta / len(mini_batch)) * nw for w, nw in zip(
            self.weights, nabla_w)]
        self.biases = [b - (eta / len(mini_batch)) * nb for b, nb in zip(
            self.biases, nabla_b)]

    def backprop(self, x, y):
        """Return a tuple ``(nabla_b, nabla_w)`` representing the
        gradient for the cost function C_x. ``nabla_b`` and
        ``nabla_w`` are layer-by-layer lists of numpy arrays, similar
        to ``self.biases`` and ``self.weights``."""
        nabla_b = [np.zeros(b.shape) for b in self.biases]
        nabla_w = [np.zeros(w.shape) for w in self.weights]
        # feedforward
        activation = x
    activations = [x]  # list to store all the activations, layer by layer
        zs = []  # list to store all the z vectors, layer by layer
        for b, w in zip(self.biases, self.weights):
            z = np.dot(w, activation) + b
            zs.append(z)
            activation = sigmoid(z)
            activations.append(activation)
        # backward pass
        delta = self.cost_derivative(activations[-1], y) * 
            sigmoid_prime(zs[-1])
        nabla_b[-1] = delta
        nabla_w[-1] = np.dot(delta, activations[-2].transpose())
        # Note that the variable l in the loop below is used a little
        # differently to the notation in Chapter 2 of the book. Here,
        # l = 1 means the last layer of neurons, l = 2 is the
        # second-last layer, and so on. It's a renumbering of the
        # scheme in the book, used here to take advantage of the fact
        # that Python can use negative indices in lists.
        for l in range(2, self.num_layers):
            z = zs[-l]
            sp = sigmoid_prime(z)
           delta = np.dot(self.weights[-l + 1].transpose(), delta) * sp
           nabla_b[-l] = delta
           nabla_w[-l] = np.dot(delta, activations[-l - 1].transpose())
        return (nabla_b, nabla_w)

    def evaluate(self, test_data):
        """Return the number of test inputs for which the neural
        network outputs the correct result. Note that the neural
        network's output is assumed to be the index of whichever
        neuron in the final layer has the highest activation."""

        test_results = [(np.argmax(self.feedforward(x)), y)
                        for (x, y) in test_data]
        return sum(int(x == y) for (x, y) in test_results)

    def cost_derivative(self, output_activations, y):
        """Return the vector of partial derivatives partial C_x /
        partial a for the output activations."""
        return (output_activations - y)

# Miscellaneous functions
def sigmoid(z):
    """The sigmoid function."""
    return 1.0 / (1.0 + np.exp(-z))

def sigmoid_prime(z):
    """Derivative of the sigmoid function."""    return sigmoid(z) * (1 - sigmoid(z))