本来这是之前课程的作业报告，本来打算多加润色一下再发上来的，但是新的PJ已经出现，后来想想就把这道稚嫩的心路历程直接放在这里也未尝不好吧。
不知道大家看自己过去写的博客会不会有尬得抠脚的感觉，反正我每次都不敢点开，即使当初是很用心地写的。真的打开看了反倒觉得，当初是很认真地在感悟的，并没有想象中那么可笑。这或许也是特别的经历吧。

Part1. 回归问题

实现思路

BP神经网络的实现主要分为两部分：

  1. 前向传播
  2. 后向传播

接下来分别说明实现方法及部分细节。

前向传播

数据生成

Untitled

首先是数据生成。这一部分对应函数 produce_data() 。

考虑到我们希望拟合函数 y=sin(x),x∈[−π,π]y=sin(x), x\in [-\pi, \pi]，我们希望创建出两个数组，X_data为输入，y_data为对应的输出。

对于输入，需要在区间内随机采样，所以使用 np.random.rand() ，美中不足是取到的数据右侧为开区间，不过正态分布本来就几乎取不到端点，影响不大。

对于输出，我们希望有小于0.01的平均误差，也就是加入正态分布的噪点。

结构初始化

Untitled

接下来开始搭建基础的网络结构。这一部分对应函数 generate_parameters()。

首先是输入层。输入的向量有几个维度，输入层就有几个神经元，由此可以确定输入层的神经元个数为 X_train.T[0]。

然后是输出层。我们只需给出一个预测值作为输出，故输出数据是一维的，只需要一个神经元。

而主要的设计部分在于隐藏层。考虑每两层之间的关系。记前一层的输出为a，当前层的输出为b，则该层第0个神经元的输出为 b[0]=w0,0a[0]+w0,1a[1]+...+w0,0a[m]+biasb[0]=w_{0,0} a[ 0] +w_{0,1} a[ 1] +...+w_{0,0} a[ m] +bias，其中m为第前一层的维数。整理下来我们将得到一个向量，向量中的每个元素为一权重和偏置构成的二元组，元素个数为当前层的维数。

Untitled

我们用一个大数组layers来记录权重和梯度等信息。layers由layer组成，layer又由neurons组成，每个neuron是一个python对象，记录了输入数组、权重、输出值以及梯度。另外再使用数组biases记录偏置信息，通过biases[layer_idx]获得该层的偏置。

层层推进

Untitled

如图所示，点乘权重和该层的输入（也就是前一层的输出），然后再应用激活函数即可。

注意到输出层需要输出一个值，而激活函数的值域是有限的，所以在输出层不能再应用激活函数。直接采用其线性输出。

损失计算

使用损失函数计算出与实际标签值的差异即可。

后向传播

现在已经能够算出loss了，接下来做的事情就是返工层层背锅。

Untitled

我们的目标是找到一个类似递推式的东西，从而能够将loss传播回去；那么，需要考虑怎样向前面一层走。

考虑前一层的第i个神经元，用rir_i表示该神经元的输出，而 sjs_j 表示下一层的第j个神经元的输出，则有si=σ(w1r1+w2r2+...+wmrm+b)s_{i} =\sigma ( w_{1} r_{1} +w_{2} r_{2} +...+w_{m} r_{m} +b)。记 zi=w1r1+w2r2+...+wmrm+bz_{i} =w_{1} r_{1} +w_{2} r_{2} +...+w_{m} r_{m} +b ，则有 si=σ(zi)s_i=\sigma(z_i)。

于是接着可以得到：

∂loss∂ri=∑j=1n∂loss∂sj×∂sj∂zi×∂zi∂ri=∑j=1n∂loss∂sj×σ′×wi\frac{\partial loss}{\partial r_{i}} =\sum_{j=1}^{n}\limits\frac{\partial loss}{\partial s{j}} \times \frac{\partial s_{j}}{\partial z_{i}} \times \frac{\partial z_{i}}{\partial r_{i}} =\sum_ {j=1}^{n}\limits\frac{\partial loss}{\partial s{j}} \times \sigma '\times w_{i}

∂loss∂b1=∑j=1n∂loss∂sj×∂sj∂zi×∂zi∂b1=∑j=1n∂loss∂sj×σ′\frac{\partial loss}{\partial b_{1}} =\sum_{j=1}^{n}\limits\frac{\partial loss}{\partial s{j}} \times \frac{\partial s_{j}}{\partial z_{i}} \times \frac{\partial z_{i}}{\partial b_{1}} =\sum_{j=1}^{n}\limits\frac{\partial loss}{\partial s{j}} \times \sigma '

其中第i个神经元所带有的权重还未调整，我们再对 rir_iri 继续求偏导即可。若 ri=p1x1+p2x2+…+pmxm+br_i=p_1 x_1+p_2 x_2+ …+p_m x_m+b ，再对其中第k个权重求偏导会得到 ∂ri∂pk=xk\frac{\partial r_i}{\partial p_k} =x_k ，代入上面的第一式，可知 ∂loss∂pk=∑j=1n∂loss∂sj×σ′×wi×xk\frac{\partial loss}{\partial p_{k}}=\sum_ {j=1}^{n}\limits\frac{\partial loss}{\partial s{j}} \times \sigma '\times w_{i}\times x_k。

分别调整 pkp_k 和 bb 为减去学习率和偏导的乘积的结果，这就是更新后的值。

代码结构

从main进入代码，进入到NeuralNetwork函数中，在这个函数的内部会首先初始化网络，然后在给定epoch数量内反复循环。

在一个epoch中，流程如下：

  1. 打散数据并从中分离出x和y
  2. 对每个x中的标签值进行前向传播
  3. mini-batch：如果这个数据被抽样抽到，就以这个数据进行反向传播。
  4. 计算误差并绘制图表。每50个epoch打印一次loss结果，每1000个epoch使用matplotlib打印一次该epoch的预测结果。

参数调整/实验过程

由于训练结果时好时坏，参数的调整主要依靠玄学。这里仅记录我的做法。

参数个数

由于sin函数是一个很简单的函数，因此我采用了较少的参数个数进行拟合。我是在训练过程中观察图像变化的，图像好像一头降下去一头就会翘起来，或者变得非常平滑类似线性函数，总之形态很怪异。

在增加了参数个数之后，拟合能力得到了显著的提高。

踩坑：激活函数求导

在反向传播求导的过程中需要对激活函数也求导，例如sigmoid()。但是我在理解上犯了一些错误，如下：

1
2
3
4
5
def sigmoid(x):
    return 1 / (1 + np.exp(-x))

def sigmoid_derivative(x):
    return sigmoid(x)*(1-sigmoid(x)) # !! wrong

从公式表达上来看sigmoid_derivative确实是这样，但是这样这个x从哪来呢？从我们保存的值着手，我们只有激活后的值还在。之前我一直填入了错误的参数，导致训练似是而非。

这也引起我的思考，由于有这样的特性，可用的激活函数并不多。

学习率、batch

实际操作过程中发现学习率提高后表现很糟，所以取了0.0001左右的值。batch大小设置为200。

操作发现batch大小取得较大时，拟合的曲线更平滑，同时变化缓慢（毕竟取样点变少了）。batch取得越小，学习率相应也要越小，原因很显然，batch小的时候学习很频繁。

Untitled

【这是一次实验中batch=20时候的情况，batch调小更容易出现尖的图像】

在训练过程中使用了如下优化：

1
2
3
4
5
6
7
learning_rate = 0.0005
if loss < 0.05:
    learning_rate = 0.00005
elif loss < 0.030:
    learning_rate = 0.00001
elif loss < 0.015:
    learning_rate = 0.000002

梯度爆炸

最开始损失函数的结果经常出现nan，发现是梯度爆炸训练不下去了。为了防止这种情况的发生，使用了一个简单朴素的函数 clipped_gradient 来处理，超过给定的阈限值的梯度会被裁剪。

1
2
3
4
5
def clipped_gradient(gradient, threshold=10):
    if abs(gradient) <= abs(threshold):
        return gradient
    else:
        return threshold*np.sign(gradient)

处理之后会发现，小的扰动被遏止了，但是该来的严重梯度爆炸还是会爆炸。梯度爆炸的情况一般总伴随着较高的学习率，在拟合的过程中可以时常见到图像产生巨大的变化，并且图像往往是不太平滑的。观察那些有危险的图像，我认为也有拟合能力不足的原因——梯度被裁剪之后，即使网络不爆炸，往往也不能成功拟合到可接受的值；所以梯度裁剪至多算是一种权宜之计，真正重要的还是调整网络自身的参数使得它们处于一个合理的量级。因此我最终将clipped_gradient删去了。

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
import numpy as np
import matplotlib.pyplot as plt


def tanh(x):
    return (np.exp(x)-np.exp(-x))/(np.exp(x)+np.exp(-x))


def tanh_derivative(output):
    return 1-output**2


def produce_data(num_samples):
    x = np.random.rand(num_samples)*2*np.pi-np.pi
    y = np.sin(x)
    return x, y


def mse_loss(prediction, target):
    return 1/2*np.mean((prediction - target) ** 2)  # mse


class Neuron:
    def __init__(self, weight_size, input=None):
        self.weight = np.random.random(weight_size)*1/weight_size
        self.input = input
        self.output = 0
        self.delta = 0


def initialize_network(layer_dims):
    layers = []
    input_layer = []
    # input layer: with no weight
    for num in range(layer_dims[0]):
        input_layer.append(Neuron(weight_size=0))
    layers.append(input_layer)

    # other layers: with weight size equal to the neurons in the layer before
    for i in range(1, len(layer_dims)):
        layer = []
        for num in range(layer_dims[i]):
            layer.append(
                Neuron(weight_size=layer_dims[i-1], input=layers[i-1]))
        layers.append(layer)

    # biases: negative random numbers.
    biases = -0.1*np.random.random(len(layer_dims))
    return layers, biases


def prop_forward(input, layers, biases):
    for i in range(len(input)):
        layers[0][i].output = input[i]

    for layer_idx in range(1, len(layers)-1):
        bias = biases[layer_idx]
        # get input values from the previous layer
        values = []
        for neuron_idx in range(len(layers[layer_idx-1])):
            values.append(layers[layer_idx-1][neuron_idx].output)
        # calculate
        for neuron_idx in range(len(layers[layer_idx])):
            layers[layer_idx][neuron_idx].output = tanh(np.dot(
                values, layers[layer_idx][neuron_idx].weight)+bias)

    # output layer
    hidden_values = []
    prediction = []
    for neuron_idx in range(len(layers[-2])):
        hidden_values.append(layers[-2][neuron_idx].output)
    for neuron_idx in range(len(layers[-1])):
        layers[-1][neuron_idx].output = np.dot(hidden_values,
                                               layers[-1][neuron_idx].weight)+biases[-1]
        # collect prediction from the output layer
        prediction.append(layers[-1][neuron_idx].output)

    return prediction


def prop_backward(target, layers, biases, learning_rate=0.01):
    # output layer
    layers[-1][0].delta = layers[-1][0].output-target
    biases[-1] -= layers[-1][0].delta*learning_rate
    for neuron_idx in range(len(layers[-2])):
        layers[-2][neuron_idx].delta = layers[-1][0].delta * \
            layers[-1][0].weight[neuron_idx]
        layers[-1][0].weight[neuron_idx] -= layers[-1][0].delta * \
            layers[-1][0].input[neuron_idx].output*learning_rate

    # hidden layers
    for layer_idx in reversed(range(0, len(layers)-2)):
        for neuron_idx in range(len(layers[layer_idx])):
            layers[layer_idx][neuron_idx].delta =\
                np.sum([n.delta*n.weight[neuron_idx]*tanh_derivative(n.output)
                       for n in layers[layer_idx+1]])
        for neuron_idx in range(len(layers[layer_idx+1])):
            neuron = layers[layer_idx+1][neuron_idx]
            neuron.weight -= np.multiply([n.output for n in layers[layer_idx]],
                                         neuron.delta * tanh_derivative(neuron.output)*learning_rate)
            biases[layer_idx] -= neuron.delta * \
                tanh_derivative(neuron.output)*learning_rate


def NeuralNetwork(X_train, y_train, layer_dims, epochs=10000, batch_size=20):
    # 1. initialize the network
    layers, biases = initialize_network(layer_dims)
    # 2. train
    loss = 1
    for epoch in range(epochs):
        learning_rate = 0.0005
        if loss < 0.05:
            learning_rate = 0.00005
        elif loss < 0.030:
            learning_rate = 0.00001
        elif loss < 0.015:
            learning_rate = 0.000002
        elif loss < 0.008:
            learning_rate = 0.0000005
        # concatenate and shuffle the data
        train_data = np.column_stack((X_train, y_train))
        np.random.shuffle(train_data)
        # retrieve x and y from concatenated data
        x = train_data.T[0]
        y = train_data.T[1]
        losses = []
        predictions = []
        for i in range(x.__len__()):
            # prop forward
            prediction = prop_forward([x[i]], layers, biases)
            predictions.append(prediction)
            losses.append(mse_loss(prediction, y[i]))
            if i % batch_size == 0:
                prop_backward(y[i], layers, biases, learning_rate)
        # compute loss
        loss = np.mean(losses)
        if (epoch+1) % 1000 == 0:
            plt.clf()
            plt.plot(x, predictions, "ro")
            plt.xlabel("x")
            plt.ylabel("predict")
            plt.show()
        print(f"Epoch {epoch+1}, loss: {loss}")


if __name__ == "__main__":
    X_train, y_train = produce_data(200)

    X_train = X_train.reshape(X_train.shape[0], -1)
    y_train = y_train.T.reshape(X_train.shape[0], -1)

    nn = NeuralNetwork(X_train, y_train,
                       layer_dims=[X_train.T.shape[0], 200, 1])

我必须要说，这个代码收敛不怎么样，如果不幸被某学弟学妹发现了，请万万不要抄我的。自己的作业给我自己写啊混蛋！

实验结果

Untitled

参考

一步步手写神经网络_手写 神经网络-CSDN博客

深度学习——反向传播（Backpropagation）_南方惆怅客的博客-CSDN博客

深度学习中常用的激活函数_回归的激活函数_DeepDriving的博客-CSDN博客

15、优化算法之Mini-batch 梯度下降法 - Hzzhbest - 博客园

Part2. 分类问题

反向传播算法 - 汉字分类

汉字分类任务大多数代码都可以复用之前的网络，相同部分这里不再赘述。主要区别在于：第一，我们读进来的数据是图；第二，不同于之前的回归任务，这是一个分类任务。

以下是我在学习处理汉字分类任务时的过程的笔记，也包括我的一些思考。

one-hot编码

对于分类问题，需要对label进行独热编码。方法是，创建一个Identity矩阵，然后选取相应的行。

1
2
3
onehot_mat = np.eye(12)
# ...
data.append([img_data, onehot_mat[label]])

softmax

由于输出不再是单值，需要设计新的输出层以及损失函数。

考虑这样一个输出作为线性加权的结果：

1
2
[ 0.00104033  0.0021963   0.00333928  0.00163153  0.00411633  0.00242474
  0.00228806  0.00163667  0.00508671  0.00216196 -0.00026968  0.00251382]

对于这样的输出，我们希望加一层激活函数，将这些项映射为概率，且所有项的和为1。为此，需要使用softmax函数。

softmax是相对于（hard）argmax而言的，argmax非黑即白，只会输出0或1，而softmax会输出各个类别的置信度。

softmax的公式为：

softmax(zi)=ezi∑c=1Cezcsoftmax( z_{i}) =\frac{e^{z_{i}}}{\sum _{c=1}^{C}\limits e^{z_{c}}}

其中 ziz_i 是第i个节点的输出值，C是输出节点的个数。

损失函数

在分类问题中，损失函数不宜再使用mse，一般使用交叉熵损失。

交叉熵损失一般具有如下形式： loss=−∑xp(x)log⁡q(x)loss=-\sum_{x}\limits p( x)\log q( x) ，其中p是真实分布，q是拟合分布。

踩坑：反向传播

输出侧的部分有以下结构：

Untitled

——本应这样计算，但是实在太复杂了，因为softmax激活函数其实类似于一个全连接层，和relu或sigmoid这样的激活函数并不相同。

所以我想象将softmax另起一层，从而获得这样的结构：

Untitled

但依旧出现了很多问题。后来我在网上找到了这样一张图：

Untitled

考虑一下，softmax依旧是激活函数，它只不过是在原本前一层线性加权所计算出来的output之上再附加了一层激活函数。所以重新绘制示意图如下：

Untitled

考虑线性输出层的第i个神经元 ziz_i ，有 grad(zi)=∑j=112∂aj∂zi×∂loss∂ajgrad( z_{i}) =\sum _{j=1}^{12}\limits\frac{\partial a{j}}{\partial z_{i}} \times \frac{\partial loss}{\partial a_{j}}，这里先求 ∂loss∂aj\frac{\partial loss}{\partial a_j}的值：

∂loss∂ai=(−∑xtilog⁡(ai))ai′=−tiai\displaystyle \frac{\partial loss}{\partial a_{i}} =\left( -\sum_{x}\limits t{i}\log( a_{i})\right)_{a_i}^{'} =-\frac{t_{i}}{a_{i}}

接下来还需要求的就只有 ∂aj∂zi\frac{\partial a{j}}{\partial z_{i}}，也就是softmax的求导。

考虑两种情况：

若i=j，则有 ——

∂∂zi(ezj∑c=1Cezc)=ezi∗(∑c=1Cezc)−(∑c=1Cezc)′∗ezi(∑c=1Cezc)2=ezi∗(∑c=1Cezc)−(ezi)2(∑c=1Cezc)2=softmax(zi)∗(1−softmax(zi))\frac{\partial }{\partial z_{i}}\left(\frac{e^{z_{j}}}{\sum _{c=1}^{C} e^{z_{c}}}\right) =\frac{e^{z_{i}} *\left(\sum _{c=1}^{C} e^{z_{c}}\right) -\left(\sum _{c=1}^{C} e^{z_{c}}\right)^{'} *e^{z_{i}}}{\left(\sum _{c=1}^{C} e^{z_{c}}\right)^{2}} \\=\frac{e^{z_{i}} *\left(\sum _{c=1}^{C} e^{z_{c}}\right) -\left( e^{z_{i}}\right)^{2}}{\left(\sum _{c=1}^{C} e^{z_{c}}\right)^{2}}\\ =softmax( z_{i}) *( 1-softmax( z_{i}))

若i≠j，则有——

∂∂zi(ezj∑c=1Cezc)=−(∑c=1Cezc)′∗ezi(∑c=1Cezc)2=−ezj∗ezi(∑c=1Cezc)2=−softmax(zi)∗softmax(zj)\frac{\partial }{\partial z_{i}}\left(\frac{e^{z_{j}}}{\sum _{c=1}^{C} e^{z_{c}}}\right) =\frac{-\left(\sum _{c=1}^{C} e^{z_{c}}\right)^{'} *e^{z_{i}}}{\left(\sum _{c=1}^{C} e^{z_{c}}\right)^{2}} =\frac{-e^{z_{j}} *e^{z_{i}}}{\left(\sum _{c=1}^{C} e^{z_{c}}\right)^{2}} \\=-softmax( z_{i}) *softmax( z_{j})

将这些式子代回 grad(zi)=∑j=112∂aj∂zi×∂loss∂ajgrad( z_{i}) =\sum _{j=1}^{12}\limits\frac{\partial a{j}}{\partial z_{i}} \times \frac{\partial loss}{\partial a_{j}}，得到：

grad(zi)=∑j=112∂aj∂zi×∂loss∂aj=∂ai∂zi×∂loss∂ai+∑j≠i12∂aj∂zi×∂loss∂aj=−tiai∗ai(1−ai)+∑j≠i12tjaj∗ai∗aj=−ti+ai∑j=112tj     (sum up to 1)=−ti+ai\begin{array}{l}grad( z_{i}) \\=\sum _{j=1}^{12}\limits\frac{\partial a_{j}}{\partial z_{i}} \times \frac{\partial loss}{\partial a_{j}}\\=\frac{\partial a_{i}}{\partial z_{i}} \times \frac{\partial loss}{\partial a_{i}} +\sum \limits_{j\neq i}^{12}\frac{\partial a_{j}}{\partial z_{i}} \times \frac{\partial loss}{\partial a_{j}}\\=-\frac{t_{i}}{a_{i}} *a_{i}( 1-a_{i}) +\sum\limits _{j\neq i}^{12}\frac{t_{j}}{a_{j}} *a_{i} *a_{j}\\=-t_{i} +a_{i}\sum\limits _{j=1}^{12} t_{j} \ \ \ \ \ ( sum\ up\ to\ 1)\\=-t_{i} +a_{i}\end{array}

故loss对 ziz_i 的偏导为prediction-target。

Untitled

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
import os
import numpy as np
import imageio


def load_images(data_dir):
    data = []
    for i in range(1, 13):
        folder_path = os.path.join(data_dir, str(i))
        for pic in os.listdir(folder_path):
            data.append([os.path.join(folder_path, pic), i])
    return data


def relu(x):
    return np.maximum(0, x)


def relu_derivative(x):
    if x > 0:
        return 1
    else:
        return 0


class Neuron:
    def __init__(self, weight_size, input=None):
        self.weight = np.random.random(weight_size)*1/weight_size
        self.input = input
        self.output = 0
        self.delta = 0


def initialize_network(layer_dims):
    layers = []
    input_layer = []
    # input layer: with no weight
    for num in range(layer_dims[0]):
        input_layer.append(Neuron(weight_size=0))
    layers.append(input_layer)

    # other layers: with weight size equal to the neurons in the layer before
    for i in range(1, len(layer_dims)):
        layer = []
        for num in range(layer_dims[i]):
            layer.append(
                Neuron(weight_size=layer_dims[i-1], input=layers[i-1]))
        layers.append(layer)

    # biases: negative random numbers.
    biases = -0.1*np.random.random(len(layer_dims))
    biases[-1] = 0
    biases[0] = 0
    return layers, biases


def softmax(values):
    exp_array = [np.exp(i) for i in values]
    return [i/np.sum(exp_array) for i in exp_array]


def prop_forward(input, layers, biases):
    for i in range(len(input)):
        layers[0][i].output = input[i]

    for layer_idx in range(1, len(layers)-1):
        bias = biases[layer_idx]
        # get input values from the previous layer
        values = []
        for neuron_idx in range(len(layers[layer_idx-1])):
            values.append(layers[layer_idx-1][neuron_idx].output)
        # calculate
        for neuron_idx in range(len(layers[layer_idx])):
            layers[layer_idx][neuron_idx].output = relu(np.dot(
                values, layers[layer_idx][neuron_idx].weight)+bias)

    # output layer
    hidden_values = []
    prediction = []
    for neuron_idx in range(len(layers[-2])):
        hidden_values.append(layers[-2][neuron_idx].output)
    for neuron_idx in range(len(layers[-1])):
        layers[-1][neuron_idx].output = np.dot(hidden_values,
                                               layers[-1][neuron_idx].weight)
        prediction.append(layers[-1][neuron_idx].output)
    # softmax
    prediction = softmax(prediction)
    for i in range(len(prediction)):
        layers[-1][i].output = prediction[i]
    return prediction


def ce_loss(prediction, target):  # cross entropy
    return -1*np.dot(target, [np.log(i) for i in prediction])


def clipped_gradient(gradient, threshold=10):
    # return gradient
    if abs(gradient) <= abs(threshold):
        return gradient
    else:
        return threshold*np.sign(gradient)


def prop_backward(target, layers, biases, learning_rate=0.01):
    # output layer
    for neuron_idx in range(len(layers[-1])):
        layers[-1][neuron_idx].delta = (layers[-1]
                                        [neuron_idx].output-target)[neuron_idx]
    for neuron_idx in range(len(layers[-2])):
        layers[-2][neuron_idx].delta = np.sum(
            [n.delta*n.weight[neuron_idx]for n in layers[-1]])
    for neuron_idx in range(len(layers[-1])):
        neuron = layers[-1][neuron_idx]
        neuron.weight -= np.multiply([n.output for n in layers[-2]],
                                     neuron.delta * learning_rate)

    # hidden layers
    for layer_idx in reversed(range(0, len(layers)-2)):
        for neuron_idx in range(len(layers[layer_idx])):
            layers[layer_idx][neuron_idx].delta =\
                np.sum([n.delta*n.weight[neuron_idx]*relu_derivative(n.output)
                       for n in layers[layer_idx+1]])
        for neuron_idx in range(len(layers[layer_idx+1])):
            neuron = layers[layer_idx+1][neuron_idx]
            neuron.weight -= np.multiply([n.output for n in layers[layer_idx]],
                                         neuron.delta * relu_derivative(neuron.output)*learning_rate)
            biases[layer_idx] -= neuron.delta * \
                relu_derivative(neuron.output)*learning_rate
    biases[0] = 0


def NeuralNetwork(X_train, y_train, X_test, y_test, layer_dims, epochs=300, batch_size=500):
    # 1. initialize the network
    layers, biases = initialize_network(layer_dims)

    # 2. train
    loss = 100
    for epoch in range(epochs):
        learning_rate = 0.00005
        if loss < 1:
            learning_rate = 0.00001
        elif loss < 0.8:
            learning_rate = 0.000002
        # concatenate and shuffle the data
        train_data = np.column_stack((X_train, y_train))
        np.random.shuffle(train_data)
        # retrieve x and y from concatenated data
        x = train_data.T[0]
        y = train_data.T[1]
        losses = []
        predictions = []
        hit = 0
        for i in range(x.__len__()):
            prediction = prop_forward(x[i], layers, biases)
            predictions.append(prediction)
            if np.argmax(y[i]) == np.argmax(prediction):
                hit += 1
            losses.append(ce_loss(prediction, y[i]))
            if i % batch_size == 0:
                prop_backward(y[i], layers, biases, learning_rate)

        # compute loss
        loss = np.mean(losses)
        print(f"Epoch {epoch+1}, loss: {loss}, accuracy:{hit/len(y)}")
    # 3. test
    test_data = np.column_stack((X_test, y_test))
    test_x = test_data.T[0]
    test_y = test_data.T[1]
    test_losses = []
    test_predictions = []
    test_hit = 0
    for i in range(test_x.__len__()):
        prediction = prop_forward(test_x[i], layers, biases)
        test_predictions.append(prediction)
        if np.argmax(test_y[i]) == np.argmax(prediction):
            test_hit += 1
        test_losses.append(ce_loss(prediction, test_y[i]))
        if i % batch_size == 0:
            prop_backward(test_y[i], layers, biases, learning_rate)

    # compute loss
    test_loss = np.mean(test_losses)
    print(f"Test loss: {test_loss}, accuracy:{test_hit/len(test_y)}")


def image_to_array(data_dir):
    img_list = load_images(data_dir)
    data = []
    onehot_mat = np.eye(12)
    for i in range(len(img_list)):
        img_name = img_list[i][0]
        img_arr = imageio.v2.imread(img_name)
        img_data = 255.0-img_arr.reshape(784)
        label = img_list[i][1]-1
        data.append([img_data, onehot_mat[label]])
    data = np.array(data, dtype=object)
    np.random.shuffle(data)
    x = data.T[0]
    y = data.T[1]
    return x, y


if __name__ == "__main__":
    X_train, y_train = image_to_array("train_data\\train")
    X_test, y_test = image_to_array("train_copy\\train")
    nn = NeuralNetwork(X_train.T, y_train.T, X_test,
                       y_test, layer_dims=[784, 40, 12])

我必须要说，这个代码收敛更加不怎么样，跑起来非常非常非常慢，如果不幸被某学弟学妹发现了，请万万不要抄我的。自己的作业给我自己写啊混蛋！

References

一文详解Softmax函数

分类问题中的损失函数

交叉熵（cross-entropy）损失函数求导过程推导_交叉熵损失函数求导-CSDN博客