资源｜如何只用NumPy码一个神经网络

选自Towards Data Science，作者：Piotr Skalski，机器之心编译，参与：高璇、张倩。

Keras、TensorFlow、PyTorch 等高级框架可以帮助我们快速构建复杂模型。深入研究并理解其中的理念很有价值。不久前，本文作者发表了一篇文章（参见《资源 | 来自独秀同学的深度网络数学笔记，还不快收藏？》），简明扼要地解释了神经网络的工作原理，但那篇文章偏向于数学理论知识。所以作者打算以一种更实际的方式来跟进这一话题。他们尝试只使用 NumPy 构建一个全运算的神经网络，通过解决简单的分类问题来测试模型，并将其与 Keras 构建的神经网络进行性能比较。

注：本文将包含大量用 Python 编写的代码片段。希望读起来不会太无聊。：)所有源代码都可以在作者的 GitHub 上找到。链接：https://github.com/SkalskiP/ILearnDeepLearning.py

图 1 ：密集神经网络架构

磨刀不误砍柴工

在开始编程之前，需要先整理一个基本的路线图。我们的目标是创建一个程序，该程序能创建一个拥有特定架构（层的数量和大小以及激活函数都是确定的）的密集连接神经网络。图 1 给出了网络的示例。最重要的是，网络必须可训练且能进行预测。

图 2 ：神经网络框图

上图显示了在训练神经网络时需要执行的操作。它还显示了在单次迭代的不同阶段，需要更新和读取多少参数。构建正确的数据结构并熟练地管理其状态是任务中最困难的部分之一。

图 3 ：l 层的权值矩阵 W 和偏置向量 b 的维数。

神经网络层初始化

首先初始化每一层的权值矩阵 W 和偏置向量 b。在图 3 中。先准备一个为系数分配适当维数的清单。上标 [l] 表示当前层的索引 (从 1 数起)，值 n 表示给定层中的单位数。假设描述 NN 架构的信息将以类似 Snippet 1 的列表形式传递到程序中，列表的每一项是一个描述单个网络层基本参数的字典：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"},
]

Snippet 1：包含描述特定神经网络参数的列表。该列表对应图 1 所示的 NN。

如果你对这个话题很熟悉，你可能已经在脑海中听到一个焦虑的声音：「嘿，嘿！这里有问题！有些领域是不必要的……」是的，这次你内心的声音是对的。前一层输出的向量是下一层的输入，所以实际上只知道一个向量的大小就足够了。但我特意使用以下符号来保持所有层之间目标的一致性，使那些刚接触这一课题的人更容易理解代码。

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

Snippet 2：初始化权值矩阵和偏置向量值的函数。

最后是这一部分最主要的任务——层参数初始化。看过 Snippet 2 上的代码并对 NumPy 有一定经验的人会发现，矩阵 W 和向量 b 被小的随机数填充。这种做法并非偶然。权值不能用相同的数字初始化，不然会出现「对称问题」。如果所有权值一样，不管输入 X 是多少，隐藏层中的所有单位都相同。在某种程度上，我们在初始阶段就会陷入死循环，无论训练模型时间多长、网络多深都无法逃脱。线性代数是不会被抵消的。

在第一次迭代中，使用较小的数值可以提高算法效率。通过图 4 所示的 sigmoid 函数图可以看到，对于较大数值，它几乎是平的，这十分影响 NN 的学习速度。总之，使用小随机数进行参数初始化是一种简单的方法，能保证我们的算法有足够好的起点。准备好的参数值存储在带有唯一标定其父层的 python 字典中。字典在函数末尾返回，因此算法的下一步是访问它的内容。

图 4：算法中使用的激活函数。

激活函数

我们将使用的函数中，有几个函数非常简单但功能强大。激活函数可以写在一行代码中，但却能使神经网络表现出自身所需的非线性性能和可表达性。「没有它们，我们的神经网络就会变成由多个线性函数组合而成的线性函数。」可选激活函数很多，但在这个项目中，我决定使用这两种——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;

Snippet 3：ReLU 和 Sigmoid 激活函数及其导数。

前向传播

设计好的神经网络有一个简单的架构。信息以 X 矩阵的形式沿一个方向传递，穿过隐藏的单元，从而得到预测向量 Y_hat。为了便于阅读，我将前向传播分解为两个单独的函数——对单个层进行前向传播和对整个 NN 进行前向传播。

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

Snippet 4：单层前向传播步骤

这部分代码可能是最容易理解的。给定上一层的输入信号，我们计算仿射变换 Z，然后应用选定的激活函数。通过使用 NumPy，我们可以利用向量化——一次性对整个层和整批示例执行矩阵运算。这减少了迭代次数，大大加快了计算速度。除了计算矩阵 A，我们的函数还返回一个中间值 Z。作用是什么呢？答案如图 2 所示。我们需要在反向传播中用到 Z。

图 5 ：在前向传播中使用的单个矩阵的维数。

使用预设好的一层前向函数后，就可以轻松地构建整个前向传播。这个函数稍显复杂，它的作用不仅是预测，还要管理中间值的集合。它返回 Python 字典，其中包含为特定层计算的 A 和 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

Snippnet 5：完整前向传播步骤

损失函数

为了观察进度，保证正确方向，我们通常需要计算损失函数的值。「一般来说，损失函数用来表征我们与『理想』解决方案的距离。」我们根据要解决的问题来选择损失函数，像 Keras 这样的框架会有多种选择。因为我计划测试我们的 NN 在两类点上的分类，所以选择二进制交叉熵，它定义如下。为了获得更多学习过程的信息，我决定引入一个计算准确率的函数。

Snippnet 6：损失函数和准确率计算

反向传播

许多缺乏经验的深度学习爱好者认为反向传播是一种难以理解的算法。微积分和线性代数的结合常常使缺乏数学基础的人望而却步。所以如果你无法马上理解，也不要担心。相信我，我们都经历过这个过程。

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

Snippnet 7：单层反向传播步骤

人们常常混淆反向传播与梯度下降，但实际上这是两个独立的问题。前者的目的是有效地计算梯度，而后者是利用计算得到的梯度进行优化。在 NN 中，我们计算关于参数的代价函数梯度（之前讨论过），但是反向传播可以用来计算任何函数的导数。这个算法的本质是在已知各个函数的导数后，利用微分学中的链式法则计算出结合成的函数的导数。对于一层网络，这个过程可用下面的公式描述。本文主要关注的是实际实现，故省略推导过程。通过公式可以看出，预先记住中间层的 A 矩阵和 Z 矩阵的值是十分必要的。

图 6：一层中的前向和反向传播。

就像前向传播一样，我决定将计算分为两个独立的函数。第一个函数（Snippnet7）侧重一个单独的层，可以归结为用 NumPy 重写上面的公式。第二个表示完全反向传播，主要在三个字典中读取和更新值。然后计算预测向量（前向传播结果）的代价函数导数。这很简单，它只是重述了下面的公式。然后从末端开始遍历网络层，并根据图 6 所示的图计算所有参数的导数。最后，函数返回 python 字典，其中就有我们想求的梯度。

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

Snippnet 8：全反向传播步骤

更新参数值

该方法的目标是利用梯度优化来更新网络参数，以使目标函数更接近最小值。为了实现这项任务，我们使用两个字典作为函数参数：params_values 存储参数的当前值；grads_values 存储根据参数计算出的代价函数导数。虽然该优化算法非常简单，只需对每一层应用下面的方程即可，但它可以作为更高级优化器的一个良好起点，所以我决定使用它，这也可能是我下一篇文章的主题。

def update(params_values, grads_values, nn_architecture, learning_rate):
 for layer_idx, layer in enumerate(nn_architecture):
 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;

Snippnet 9：利用梯度下降更新参数值

组合成型

任务中最困难的部分已经过去了，我们已经准备好了所有必要的函数，现在只需把它们按正确的顺序组合即可。为了更好地理解操作顺序，需要对照图 2 的表。该函数经过训练和期间的权值变化返回了最优权重。只需要使用接收到的权重矩阵和一组测试数据即可运行完整的前向传播，从而进行预测。

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

Snippnet 10：训练模型

David vs Goliath

现在可以检验我们的模型在简单的分类问题上的表现了。我生成了一个由两类点组成的数据集，如图 7 所示。然后让模型学习对两类点分类。为了便于比较，我还在高级框架中编写了 Keras 模型。两种模型具有相同的架构和学习速率。尽管如此，这样对比还是稍有不公，因为我们准备的测试太过于简单。最终，NumPy 模型和 Keras 模型在测试集上的准确率都达到了 95%，但是我们的模型需要多花几十倍的时间才能达到这样的准确率。在我看来，这种状态主要是由于缺乏适当的优化。

图 7：测试数据集

图 8：两种模型实现的分类边界可视化