手写感知器的反向传播算法 Part2：基于理论推导的代码实践：如何实现自己的第一个神经网络完成手写体识别？（上）

前言：

本文是手写实现BP神经网络的第二篇文章，也是本专栏中最核心的一篇文章，因为本文作者主要介绍的是作者的程序是如何实现的。但是整个网络最大的难度其实是在于计算梯度，搞清楚怎么做反向传播才能写出正确的代码。关于如何计算梯度部分，请看作者的上一篇文章，链接在下面：

手写感知器的反向传播算法 Part1：梯度传播的理论推导：如何计算梯度？-CSDN博客

在上一篇文章之中作者给出了应该如何将梯度运算运用在程序中的大致思路，这篇文章中作者就将一一加以实现。那么闲言少叙，让我们进入正题。记住作者的一句话：掌握了梯度传播的计算方式，那么整个BP神经网络就已经完成了一半了。由于一整个程序的讲解有些过于庞大和冗长了，所以作者这里决定分成上下两个部分。上篇讲的是除了前向推理，反向传播操作的具体计算的代码的所有部分，读者会在这个模块了解应该怎么写出一个神经网络，框架应该是如何的。而具体怎么计算的代码会在下篇进行呈现。

Part1：整体结构简介:

完整的程序由两个部分组成：

第一个部分是对于图片的处理过程，在BP神经网络这个语境下至少包含一个作用：将28*28的图片展开成为784*1的列向量。其余的作用作者之后会详细讲到。

第二个部分就是主程序。主程序的构成主要是三个Part：

Part1. main函数中完成对于网络结构，学习率的确定。在本程序中，注意学习率，网络结构，和每一层的参数都是可以灵活调整的。这就是为什么第一个部分和第二个部分不能合二为一的原因：如果我们固定了网络结构就只需要第二个和第三个Part了

Part2.“BP_Model_for_classification”类。这个类主要是为了完成深度学习的整个流程设计的：在这个类中最主要的两个函数是create_model和train函数，前者是使用第三个类创建一个模型，后者是使用这个创造出来的model，调用其成员函数完成训练时候的前向推理和反向传播等。

Part3:"BP_foward_for_classification"类。这个类的目的是承载所有的计算过程。train函数调用的所有的计算操作都是这个类中的成员函数。这个部分作者会在（下）中详细介绍。

那么下面作者将按照顺序依次给出代码和解析。

Part2：准备部分：图像处理：

思路解析：

1.图像是什么样的？

手写体识别的图片一般都是一张白底黑字的图片，图上实际上只有黑色与白色这两种颜色。

 这个汉字在python中读入成为数字矩阵之后，这个字的数字矩阵只有两个数字，分别是255和0。其中0元素对应的点位是空白，而255元素对应的点位是黑色。

2.我们想要将图片处理成什么样子？

（1）从减小计算量的角度思考：图片的两种色块分别表示为255和0，那么从减小计算量的角度思考那么设定成1和0或者是1和-1这样的操作是和原本的图片效果一致。但是具体设定为哪种下面作者会有一个明确的介绍

（2）从标准化的角度思考：在CV的领域之中，实际使用的时候，torchvision提供了一个transforms类可以对于图片进行操作。其能达成的效果比如裁剪，翻转等可以处理原图片产生一些“新”的数据，其最大的益处之一就是可以扩充数据集。但是在这个程序之中最重要的是一步Normalize操作。这一步操作的目的是对于图片进行正则化，处理成均值为0，方差为1的矩阵。这一步经作者自己尝试之后得到的结论是有助于收敛。关于这个模块具体的实验论证可以看作者后续的文章。这里只给出结论：torchvision中的Normalize将图片处理成了1和-1的数字矩阵。

 于是基于上面这两个部分的分析，我们已经确定了总方针：将图片变成每一列由1和-1组成的列向量，并且在最后一行加上所属于的类别序号。

数据集形式：

数据集形如下：

其中每一个文件夹中对应着一种文字对应的图片。

具体代码：

import numpy as np
from PIL import Image

from main_pic import copy

# import main
np.random.seed(0)
data_all = np.zeros((64*64+1,5620),dtype=float)
count = 0
root = "train_2022/"
for i in range(10):
    for j in range(500):
        img = Image.open(root + f"{i}/{j+1}.bmp").convert("L")
        matrix = np.array(img)
        matrix_reshape = matrix.reshape((-1,1))
        matrix_column = np.zeros((4096,1))
        for p in range(64*64):
            if matrix_reshape[p,0] ==255:
                matrix_column[p,0] = 1
            else:
                matrix_column[p, 0] = -1
        #         既然整个灰度图只有0,255两个类型，那不妨把255变成1，便于运算

        copy(data_all[:,count],matrix_column)
        count += 1
for i in range(5000):
    data_all[-1,i] = i//500 +1

np.savetxt("train_data.csv",data_all,delimiter=",")

基于以上的代码，我们可以得到一个csv文件。在这个文件中前64*64行是展开的图片，最后一行是这个图片中的汉字属于的种类。

Part3：代码正式部分：main函数部分

main函数部分比较简单，其目的就是给我们一个可以设定参数的区域。在作者的这个代码之中，网络结构，learning rate等等参数都可以直接进行调整，于是这里如下设定下面的参数

if __name__ == '__main__':
    train_data = np.loadtxt("train_data.csv", delimiter=",", dtype=float)
    test_data = np.loadtxt("val_data.csv", delimiter=",", dtype=float)
    # read the datasets from local
    layer_num = 2
    neuron_sit = [4096, 1000, 250, 10]
    lr = 0.001
    model_path = "final_test_for_resume.pkl"
    model = BP_Model_for_classification(layer_num, neuron_sit, lr,model_path)
    #     上面完成了整个参数的初始化和所有的准备工作，下面将完成对于整个模型的前向推理操作
    train_loss,train_acc = model.train(train_data, test_data)
    plot(train_acc,train_loss)

注意：这里的neuron_sit是每一层的神经元的情况。这里的4096就是64*64的结果，而这里的10是分类的结果，即将结果分成10个类别。这里的10直接接上softmax得到最后的每一个类别的概率。

上面的BP_Model_for_classification是等会儿要介绍的Part4的内容。通过这个类别我们创建出了整个模型。返回的事实模型本身。通过调用这个模型的成员函数我们完成推理。

Part4：代码正式部分：BP_Model_for_classification类：

1.构造函数：

def __init__(self, layer_num, neuron_sit, lr, last_model = None):
    self.layer_num = layer_num
    self.neuron_sit = neuron_sit
    self.learning_rate = lr
    self.epochs = 20
    self.train_loss = 0
    self.test_loss = 0
    self.mini_batch_num = 20
    if last_model == None:
        self.create_model()
    else:
        with open(model_path, 'rb') as f:
            # with open('my_model_for_classification_one_layer.pkl', 'rb') as f:
            last_model = pickle.load(f)
            self.create_model(last_bias = last_model.model.bias_para,last_weight = last_model.model.weight_para)

解析：

1.从main函数处读出一些程序的参数：learning_rate等

2.定下到底训练多少epoch。其实这一步对于大家来说没有那么有必要。因为如果我们训练的目的就是得到一个好的模型从而能跑出来更高的分数，那么理论上来说跑多少个epoch都是合适的。毕竟如果模型的精度达到了需求，直接打断就可以了。至今为止的最好的模型会背保存在本地。但是作者这里加上这一功能的原因是作者有一个绘图函数，可以反映整个训练过程中的loss和精度的变化。因此少一些的epoch就能确保作者可以得到图像。

3.mini-batch：这个必须的项目其实是大多数初学者根本想不起来做或者不理解应该怎么做的一个内容。batch的作用效果是平均。试想如果每一个样本都进行更新和推断操作，那么不仅会导致运行训练速率极慢，更会导致某些噪声也被训练进入了模型。加上batch之后，每一次反向传播的梯度都变为了20次反向传播累计的梯度的平均值，那么噪声会被平均掉很大一部分，这可以有效提升对于噪声的处理能力，同时也会加快训练速度。

4.last_model：这个部分的功能是断点续训。这个功能在后来的作者来看是非常实用的。在根本不清楚多少epoch会达到过拟合的情况下，epoch的数量不好确定：如果设定的太小了那么就会导致正确率还在上升但是训练已经终止了。这里添加了这个功能就是可以在上一次的训练基础上继续进行训练。正如代码中所写：如果没有从main函数传入last_model的地址，那么就用上次保存的model创建新的model；如果没有则create_model函数创造一个新的。

2.切换训练/测试状态以及创建一个新的模型：

这个部分的代码如下：

    def train_state(self):
        self.model.train_state = True

    def test_state(self):
        self.model.train_state = False

    def create_model(self,last_bias = None, last_weight = None):
        if last_weight is None and last_bias is None:
            self.model = BP_foward_for_classification(self.layer_num, self.neuron_sit, self.learning_rate)
        else:
            self.model = BP_foward_for_classification(self.layer_num, self.neuron_sit, self.learning_rate,last_weight,last_bias)

（1）前两个函数模仿的是成型的神经网络中的xxx.eval()，即切换模式。这里为什么需要这个机制呢？这是因为训练模式和评估模式需要的结果是不同的：训练模式需要每一个状态的计算结果以便于反向传播，而测试模式则只需要最后一个状态的计算结果。因此这里分成两个部分保存。

（2）最后的一个函数就是创建模型的部分。这个部分中的last_weight和last_bias是当我们加载上次训练的模型的时候，从保存模型的pkl文件中读出来的参数。用这两个参数创建出的模型确实可以认定是和之前的模型效果相同的。关于这个模型具体是怎么创建的请读者移步Part5的构造函数。

3.train函数部分:

这个部分是整个程序中最公式化的部分，这一系列的思路希望读者能记住并熟练运用。

Step1:训练准备工作：

x_train, y_train_digit, y_train_matrix, x_test, y_test_digit, y_test_matrix = prepare_dataset(train_data,test_data)
back_num = x_train.shape[1]
num_of_batch = back_num / self.mini_batch_num

这里首先做的是使用prepare_dataset函数准备好数据集。后续的两行代码的目的是计算出数据集中应该分出多少个batch，便于后续取图片和类别信息。

关于prepare_dataset函数效果如下：

prepare_dataset函数代码&解析

def prepare_dataset(train_data, test_data):
    train_x = train_data[:4096, :]
    train_y_digit = train_data[-1, :]
    #     这里考虑到后面的运算中正确答案出现的语境都是one-hot向量，所以将y转化成one-hot向量的形式保存为一个矩阵用于反向传播
    train_y_matrix = np.zeros((10, np.size(train_y_digit)))
    for i in range(np.size(train_y_digit)):
        x_i = train_y_digit[i] - 1
        train_y_matrix[int(x_i), i] = 1
    #     上面得到的矩阵的每一列是每一个one-hot向量
    test_x = test_data[:4096, :]
    test_y_digit = test_data[-1, :]
    #     这里考虑到后面的运算中正确答案出现的语境都是one-hot向量，所以将y转化成one-hot向量的形式保存为一个矩阵用于反向传播
    test_y_matrix = np.zeros((10, np.size(test_y_digit)))
    for i in range(np.size(test_y_digit)):
        x_i = test_y_digit[i] - 1
        test_y_matrix[int(x_i), i] = 1
    return train_x, train_y_digit, train_y_matrix, test_x, test_y_digit, test_y_matrix

关于这几个矩阵的内容，作者在注释里都做了标记。这里作者做一个简要的说明：

train_x和test_x比较明显，就是图片对应的数字矩阵展开之后得到的列向量形成的矩阵

train_y_digit和test_y_digit得到目标就是给出图片对应的类别的序号

而train_y_matrix和test_y_matrix的目标是图片对应的类别的one-hot向量编码。正如作者上一篇文章中所推导的式子一样，在反向传播时候是需要one-hot向量和softmax后的结果作差的。因此这里作者直接在这个函数之中给出这些矩阵以便后续使用。

Step2：在每一个epoch开始之前，对于数据集完成打乱

这个部分的代码如下：

            loss_all = 0
            self.model.create_result_set(x_train, x_test)
            # 每一个epoch都shuffle一次测试集和训练集
            if epoch != 0:
                x_train, y_train_digit, y_train_matrix, x_test, y_test_digit, y_test_matrix = renew_dataset(x_train,
                                                                                                            y_train_digit,
                                                                                                            x_test,
                                                                                                            y_test_digit)

create_result_set函数这里就不细说了，实际上就是创造两个成员变量，分别用来存放训练过程中的每一层的计算结果和测试时的计算结果。下面主要介绍renew_dataset函数。这个函数和上面提到的prepare_dataset效果比较相似。那这里为什么要打乱呢？

为什么要打乱数据集？

这是因为这可以避免一种极端情况。我们假设运气不好，有一个batch中的噪声全都是朝着一个方向的。那么这样的话如果每个epoch都对于这一个batch的结果进行平均后反向传播，那么这一组“坏样本”就可能对于训练造成破坏。每一次训练开始之前对于数据集进行打乱有利于避免极端情况的发生，有助于减少噪声造成的影响。

 下面给出打乱的代码：

renew_dataset函数：

def renew_dataset(x_train, y_train_digit, x_test, y_test_digit):
    train_data = np.vstack((x_train, y_train_digit))
    test_data = np.vstack((x_test, y_test_digit))
    mid_train_data = train_data.T
    mid_test_data = test_data.T
    np.random.shuffle(mid_train_data)
    np.random.shuffle(mid_test_data)
    train_data = mid_train_data.T
    test_data = mid_test_data.T
    return prepare_dataset(train_data, test_data)

为什么要加上一个.T，也就是一个转置呢？这是因为shuffle只能确保按照行进行打乱。所以为了确保每一张图片和其对应的类别可以一一对应，需要先转置后再进行打乱。 

Step3：正向推理得到一个batch中每一个图片对应的类别的矩阵：

for i in range(int(num_of_batch)):
    start = i * self.model.minibatch_num
    end = (i + 1) * self.model.minibatch_num
    # 在train_x中的开始的和结束的序号。
    self.train_state()
    # 2.利用模型做正向推理
    self.model.forward(x_train, neuron_sit, i)
    # 3.计算损失函数
    calc_result = self.model.train_result[:10, -1, :]
    #             损失函数使用cross-entropy-loss
    final_result = cross_entropy_loss(calc_result, y_train_matrix)
    #           这就是损失函数的具体值,loss之前都是对的
    print(f"step {i + 1},train loss {final_result}")
    loss_all += final_result

关于forward函数具体的操作过程都是第三个类中实现。这里读者可以先将这个部分当做一个黑盒，从思路上理解这一段干了些什么。实际上，这一段所做的就是将预测的结果和之前得到的one-hot向量形式表示的ground-truth放入cross-entropy-loss函数中计算总的loss值，并且将这些loss累加起来作为这一epoch训练结果的一个评判标准之一。

Step4：反向传播更新参数：

这一个部分的代码如下：

# 4.反向传播
for j in range(self.model.minibatch_num):
    self.train_state()
    self.model.calc_loss(j, y_train_matrix[:, start:end])
    # 这一步中的loss对应的是偏loss比上偏M，M是未经过sigmoid的z值。和上一个实验相比，主要的不同点在于第一个的信息。
    self.model.calc_delta_weight(j)
    # 用上面的Loss矩阵计算出对应的weight矩阵的减少项
    self.model.calc_delta_bias()
#     由于实际上减去的bias就是loss值本身，所以这里不用加上i这个参数，直接减去loss矩阵就行
self.model.backward()
self.model.clear_grad()
#     这里调整成每一个batch一个反向传播，然后立刻得到新的梯度信息指导下一次梯度下降

这里我们结合上一篇文章中作者推导的笔记来看：

上面的代码中的中间三行分别就对应着：计算loss矩阵，使用loss矩阵计算delta_weight和delta_bias矩阵。这几个矩阵计算出梯度之后下面就是用原先的参数减去这几个梯度矩阵，即对应着上面的backward步骤。最后将所有的梯度全部清零防止对于下一次累计梯度造成影响。至此，所有的训练流程已经完成了一遍。 

Step5：测试集验证：

代码如下：

self.test_state()
self.model.forward(x_test, neuron_sit)
# 一个epoch的训练过后，推理测试集观察正确率
calc_result_test = self.model.test_result[:10, -1, :]
final_result = cross_entropy_loss(calc_result_test, y_test_matrix)
print(f"the test result is {final_result}")
correct_num = 0
predict_result = get_max(calc_result_test)
mid_matrix = predict_result - y_test_digit
for item in mid_matrix:
    for item in item:
        if item == 0:
            correct_num += 1
result_acc = correct_num / 620
print(f"the correctness is {result_acc}")
if result_acc > best_acc:
    # 当正确率高于测试集的时候，保存模型
    print(f"the best acc is altered to{result_acc}")
    best_acc = result_acc
    with open('test_5.pkl', 'wb') as f:
        pickle.dump(model, f)
elif epoch == 19:
    with open('final_test_for_resume.pkl', 'wb') as f:
        pickle.dump(model, f)
return_acc[epoch] = result_acc

这段代码的部分的内容和之前相比多了一个get_max函数。通过这个函数，我们就可以通过softmax函数得到预测具体属于哪一个类别。在这里test_y_digit矩阵就有意义了，两者作差得到的结果中，0的数量就是预测正确的数量。统计0的数量除上总样本量结果就是正确率。

Part5：绘图代码：

这个部分的目的是得到一个正确率，训练误差等信息的直观图示。代码如下：

def plot(acc, loss, mode='train', best_acc_=None):
    plt.figure(figsize=(10, 4))
    plt.suptitle('%s_curve' % mode)
    plt.subplots_adjust(wspace=0.2, hspace=0.2)
    epoch = len(acc)

    plt.subplot(1, 2, 1)
    plt.plot(np.arange(epoch), loss, label='loss')
    plt.xlabel('epoch')
    plt.ylabel('loss')
    plt.legend(loc='upper left')

    plt.subplot(1, 2, 2)
    plt.plot(np.arange(epoch), acc, label='acc')
    if best_acc_ is not None:
        plt.scatter(best_acc_[0], best_acc_[1], c='r')
    plt.xlabel('epoch')
    plt.ylabel('acc')
    plt.legend(loc='upper left')
    plt.savefig('my_own_model_%s.jpg' % mode, bbox_inches='tight')
    plt.show()

做出来的图片如下：

 可以看出相当直观。

总结：

本文章中给出的部分虽然没有涉及到实际上的正向推理或者反向传播的具体计算方式，但是本文作者认为在重要性上不弱于下篇，也就是关于细节上如何实现正向传播和反向传播的操作介绍。因为对于初学者而言，这个范式，也就是怎么把一个程序的框架搭出来，其内在信息是相当丰富的，也是后续在各种实验中都会有用的一种技能。所以建议所有的读者仔细阅读搞懂背后的思想，内化于自身。同时敬请期待我的下篇文章。将上下两篇文章同时阅读理解之后，完成手写BP神经网络完成手写体识别将不再话下。

都看到这里了，不给可怜的大学牲点个赞嘛awa