神经网络基础

最近给本科生当机器学习课程的助教,给他们出的作业题需要看这些图,懒得放本地了,直接放博客里。发现jupyter导出markdown好方便,放到博客里面正好,改都不用改。

原来就想过一个问题,为什么我写出来的神经网络不收敛,loss会像火箭一样直接飞了。后来看了一些教程,发现有人在做梯度下降的时候,把梯度除以了梯度的二范数,我尝试之后发现还真好使了,在实验的时候发现是因为没有对数据集进行归一化,如果所有的数据都是很大的数,那么在反向传播的时候,计算出来的梯度的数量级会很大,这就导致更新得到的参数的数量级也很大,预测出的偏差就更大了,然后循环往复,如果给梯度除以一个梯度的二范数,其实就相当于把梯度的数量级降了,这样就可以训练了。但实际上还是将原始数据归一化比较好,对原始数据归一化还能让梯度下降的方向更多。如果数据都是正数,那下降方向会少很多,下降的时候会出现zig-zag现象。

第二题:神经网络:线性回归

实验内容:

  1. 学会梯度下降的基本思想
  2. 学会使用梯度下降求解线性回归
  3. 了解归一化处理的作用

线性回归

我们来完成最简单的线性回归,上图是一个最简单的神经网络,一个输入层,一个输出层,没有激活函数。
我们记输入为$X \in \mathbb{R}^{n \times m}$,输出为$Z \in \mathbb{R}^{n}$。输入包含了$n$个样本,$m$个特征,输出是对这$n$个样本的预测值。
输入层到输出层的权重和偏置,我们记为$W \in \mathbb{R}^{m}$和$b \in \mathbb{R}$。
输出层没有激活函数,所以上面的神经网络的前向传播过程写为:

我们使用均方误差作为模型的损失函数

我们通过调整参数$W$和$b$来降低均方误差,或者说是以降低均方误差为目标,学习参数$W$和参数$b$。当均方误差下降的时候,我们认为当前的模型的预测值$Z$与真值$y$越来越接近,也就是说模型正在学习如何让自己的预测值变得更准确。

在前面的课程中,我们已经学习了这种线性回归模型可以使用最小二乘法求解,最小二乘法在求解数据量较小的问题的时候很有效,但是最小二乘法的时间复杂度很高,一旦数据量变大,效率很低,实际应用中我们会使用梯度下降等基于梯度的优化算法来求解参数$W$和参数$b$。

梯度下降

梯度下降是一种常用的优化算法,通俗来说就是计算出参数的梯度(损失函数对参数的偏导数的导数值),然后将参数减去参数的梯度乘以一个很小的数(下面的公式),来改变参数,然后重新计算损失函数,再次计算梯度,再次进行调整,通过一定次数的迭代,参数就会收敛到最优点附近。

在我们的这个线性回归问题中,我们的参数是$W$和$b$,使用以下的策略更新参数:

其中,$\alpha$ 是学习率,一般设置为0.1,0.01等。

接下来我们会求解损失函数对参数的偏导数。

损失函数MSE记为:

其中,$Z \in \mathbb{R}^{n}$是我们的预测值,也就是神经网络输出层的输出值。这里我们有$n$个样本,实际上是将$n$个样本的预测值与他们的真值相减,取平方后加和。

我们计算损失函数对参数$W$的偏导数,根据链式法则,可以将偏导数拆成两项,分别求解后相乘:

这里我们以矩阵的形式写出推导过程,感兴趣的同学可以尝试使用单个样本进行推到,然后推广到矩阵形式

同理,求解损失函数对参数$b$的偏导数:

因为参数$b$对每个样本的损失值都有贡献,所以我们需要将所有样本的偏导数都加和。

其中,$\frac{\partial \mathrm{loss}}{\partial W} \in \mathbb{R}^{m}$,$\frac{\partial \mathrm{loss}}{\partial b} \in \mathbb{R}$,求解得到的梯度的维度与参数一致。

完成上式两个梯度的计算后,就可以使用梯度下降法对参数进行更新了。

训练神经网络的基本思路:

  1. 首先对参数进行初始化,对参数进行随机初始化(也就是取随机值)
  2. 将样本输入神经网络,计算神经网络预测值 $Z$
  3. 计算损失值MSE
  4. 通过 $Z$ 和 $y$ ,以及 $X$ ,计算参数的梯度
  5. 使用梯度下降更新参数
  6. 循环1-5步,在反复迭代的过程中可以看到损失值不断减小的现象,如果没有下降说明出了问题

接下来我们来实现这个最简单的神经网络。

1. 导入数据

使用kaggle房价数据,选3列作为特征

1
2
3
4
5
6
7
8
9
10
11
12
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline

# 读取数据
data = pd.read_csv('data/kaggle_house_price_prediction/kaggle_hourse_price_train.csv')

# 使用这3列作为特征
features = ['LotArea', 'BsmtUnfSF', 'GarageArea']
target = 'SalePrice'
data = data[features + [target]]

2. 数据预处理

40%做测试集,60%做训练集

1
2
from sklearn.model_selection import train_test_split
trainX, testX, trainY, testY = train_test_split(data[features], data[target], test_size = 0.4, random_state = 32)

训练集876个样本,3个特征,测试集584个样本,3个特征

1
trainX.shape, trainY.shape, testX.shape, testY.shape

3. 参数初始化

这里,我们要初始化参数$W$和$b$,其中$W \in \mathbb{R}^m$,$b \in \mathbb{R}$,初始化的策略是将$W$初始化成一个随机数矩阵,参数$b$为0。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
def initialize(m):
'''
参数初始化,将W初始化成一个随机向量,b是一个长度为1的向量

Parameters
----------
m: int, 特征数

Returns
----------
W: np.ndarray, shape = (m, ), 参数W

b: np.ndarray, shape = (1, ), 参数b

'''

# 指定随机种子,这样生成的随机数就是固定的了,这样就可以与下面的测试样例进行比对
np.random.seed(32)

W = np.random.normal(size = (m, )) * 0.01

b = np.zeros((1, ))

return W, b

4. 前向传播

这里,我们要完成输入矩阵$X$在神经网络中的计算,也就是完成 $Z = XW + b$ 的计算。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
def forward(X, W, b):
'''
前向传播,计算Z = XW + b

Parameters
----------
X: np.ndarray, shape = (n, m),输入的数据

W: np.ndarray, shape = (m, ),权重

b: np.ndarray, shape = (1, ),偏置

Returns
----------
Z: np.ndarray, shape = (n, ),线性组合后的值

'''

# 完成Z = XW + b的计算
# YOUR CODE HERE
Z = np.dot(X, W) + b

return Z
1
2
3
4
# 测试样例
Wt, bt = initialize(trainX.shape[1])
tmp = forward(trainX, Wt, bt)
print(tmp.mean()) # -28.37377

5. 损失函数

接下来编写损失函数,我们以均方误差(MSE)作为损失函数,需要大家实现MSE的计算:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
def mse(y_true, y_pred):
'''
MSE,均方误差

Parameters
----------
y_true: np.ndarray, shape = (n, ),真值

y_pred: np.ndarray, shape = (n, ),预测值

Returns
----------
loss: float,损失值

'''

# 计算MSE
# YOUR CODE HERE
loss = ((y_true - y_pred) ** 2).sum() / len(y_true)

return loss
1
2
3
4
# 测试样例
Wt, bt = initialize(trainX.shape[1])
tmp = mse(trainY, forward(trainX, Wt, bt))
print(tmp) # 39381033680.5

6. 反向传播

这里我们要完成梯度的计算,也就是计算出损失函数对参数的偏导数的导数值:

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
def compute_gradient(X, Z, y_true):
'''
计算梯度

Parameters
----------
X: np.ndarray, shape = (n, m),输入的数据

Z: np.ndarray, shape = (n, ),线性组合后的值

y_true: np.ndarray, shape = (n, ),真值

Returns
----------
dW, np.ndarray, shape = (m, ), 参数W的梯度

db, np.ndarray, shape = (1, ), 参数b的梯度

'''

n = len(y_true)

# 计算W的梯度
# YOUR CODE HERE
dW = np.dot(X.T, (Z - y_true)) * 2 / n

# 计算b的梯度
# YOUR CODE HERE
db = (Z - y_true).sum() / n

return dW, db
1
2
3
4
5
6
7
# 测试样例
Wt, bt = initialize(trainX.shape[1])
Zt = forward(trainX, Wt, bt)
dWt, dbt = compute_gradient(trainX, Zt, trainY)
print(dWt.shape) # (3,)
print(dWt.mean()) # -1532030241.25
print(dbt.mean()) # -182154.277882

7. 梯度下降

这部分需要实现梯度下降的函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
def update(dW, db, W, b, learning_rate):
'''
梯度下降,参数更新,不需要返回值,W和b实际上是以引用的形式传入到函数内部,
函数内改变W和b会直接影响到它们本身

Parameters
----------
dW, np.ndarray, shape = (m, ), 参数W的梯度

db, np.ndarray, shape = (1, ), 参数b的梯度

W: np.ndarray, shape = (m, ),权重

b: np.ndarray, shape = (1, ),偏置

learning_rate, float,学习率

'''
# 更新W
W -= learning_rate * dW

# 更新b
b -= learning_rate * db
1
2
3
4
5
6
7
8
9
10
11
12
# 测试样例
Wt, bt = initialize(trainX.shape[1])
print(Wt.mean()) # 0.00405243937693
print(bt.mean()) # 0.0

Zt = forward(trainX, Wt, bt)
dWt, dbt = compute_gradient(trainX, Zt, trainY)
update(dWt, dbt, Wt, bt, 0.01)

print(Wt.shape) # (3,)
print(Wt.mean()) # 15320302.4166
print(bt.mean()) # 1821.54277882

完成整个参数更新的过程,先计算梯度,再更新参数,将compute_gradient和update组装在一起。

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
def backward(X, Z, y_true, W, b, learning_rate):
'''
使用compute_gradient和update函数,先计算梯度,再更新参数

Parameters
----------
X: np.ndarray, shape = (n, m),输入的数据

Z: np.ndarray, shape = (n, ),线性组合后的值

y_true: np.ndarray, shape = (n, ),真值

W: np.ndarray, shape = (m, ),权重

b: np.ndarray, shape = (1, ),偏置

learning_rate, float,学习率

'''
# 计算参数的梯度
# YOUR CODE HERE
dW, db = compute_gradient(X, Z, y_true)

# 更新参数
# YOUR CODE HERE
update(dW, db, W, b, learning_rate)
1
2
3
4
5
6
7
8
9
10
11
# 测试样例
Wt, bt = initialize(trainX.shape[1])
print(Wt.mean()) # 0.00405243937693
print(bt.mean()) # 0.0

Zt = forward(trainX, Wt, bt)
backward(trainX, Zt, trainY, Wt, bt, 0.01)

print(Wt.shape) # (3,)
print(Wt.mean()) # 15320302.4166
print(bt.mean()) # 1821.54277882

8. 训练

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
def train(trainX, trainY, testX, testY, W, b, epochs, learning_rate = 0.01, verbose = False):
'''
训练,我们要迭代epochs次,每次迭代的过程中,做一次前向传播和一次反向传播,更新参数
同时记录训练集和测试集上的损失值,后面画图用。然后循环往复,直到达到最大迭代次数epochs

Parameters
----------
trainX: np.ndarray, shape = (n, m), 训练集

trainY: np.ndarray, shape = (n, ), 训练集标记

testX: np.ndarray, shape = (n_test, m),测试集

testY: np.ndarray, shape = (n_test, ),测试集的标记

W: np.ndarray, shape = (m, ),参数W

b: np.ndarray, shape = (1, ),参数b

epochs: int, 要迭代的轮数

learning_rate: float, default 0.01,学习率

verbose: boolean, default False,是否打印损失值

Returns
----------
training_loss_list: list(float),每迭代一次之后,训练集上的损失值

testing_loss_list: list(float),每迭代一次之后,测试集上的损失值

'''
training_loss_list = []
testing_loss_list = []

for epoch in range(epochs):

# 这里我们要将神经网络的输出值保存起来,因为后面反向传播的时候需要这个值
Z = forward(trainX, W, b)

# 计算训练集的损失值
training_loss = mse(trainY, Z)

# 计算测试集的损失值
testing_loss = mse(testY, forward(testX, W, b))

# 将损失值存起来
training_loss_list.append(training_loss)
testing_loss_list.append(testing_loss)

# 打印损失值,debug用
if verbose:
print('epoch %s training loss: %s'%(epoch+1, training_loss))
print('epoch %s testing loss: %s'%(epoch+1, testing_loss))
print()

# 反向传播,参数更新
backward(trainX, Z, trainY, W, b, learning_rate)

return training_loss_list, testing_loss_list
1
2
3
4
5
6
7
8
9
10
11
# 测试样例
Wt, bt = initialize(trainX.shape[1])
print(Wt.mean()) # 0.00405243937693
print(bt.mean()) # 0.0

training_loss_list, testing_loss_list = train(trainX, trainY, testX, testY, Wt, bt, 2, learning_rate = 0.01, verbose = False)

print(training_loss_list) # [39381033680.460075, 3.3902307664083424e+23]
print(testing_loss_list) # [38555252685.093872, 4.1516070070405267e+23]
print(Wt.mean()) # -5.70557904608e+13
print(bt.mean()) # -4412133889.08

9. 检查

编写一个绘制损失值变化曲线的函数

一般我们通过绘制损失函数的变化曲线来判断模型的拟合状态。

一般来说,随着迭代轮数的增加,训练集的loss在下降,而测试集的loss在上升,这说明我们正在不断地让模型在训练集上表现得越来越好,在测试集上表现得越来越糟糕,这就是过拟合的体现。

如果训练集loss和测试集loss共同下降,这就是我们想要的结果,说明模型正在很好的学习。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
def plot_loss_curve(training_loss_list, testing_loss_list):
'''
绘制损失值变化曲线

Parameters
----------
training_loss_list: list(float),每迭代一次之后,训练集上的损失值

testing_loss_list: list(float),每迭代一次之后,测试集上的损失值

'''
plt.figure(figsize = (10, 6))
plt.plot(training_loss_list, label = 'training loss')
plt.plot(testing_loss_list, label = 'testing loss')
plt.xlabel('epoch')
plt.ylabel('loss')
plt.legend()

上面这些函数就是完成整个神经网络需要的函数了

函数名 功能
initialize 参数初始化
forward 给定数据,计算神经网络的输出值
mse 给定真值,计算神经网络的预测值与真值之间的差距
backward 计算参数的梯度,并实现参数的更新
compute_gradient 计算参数的梯度
update 参数的更新
backward 计算参数梯度,并且更新参数
train 训练神经网络
plot_loss_curve 绘制损失函数的变化曲线

我们使用参数初始化函数和训练函数,完成神经网络的训练。

1
2
3
4
5
6
7
8
# 特征数m
m = trainX.shape[1]

# 参数初始化
W, b = initialize(m)

# 训练20轮,学习率为0.01
training_loss_list, testing_loss_list = train(trainX, trainY, testX, testY, W, b, 20, learning_rate = 0.01, verbose = True)

绘制损失值的变化曲线

1
plot_loss_curve(training_loss_list, testing_loss_list)

通过打印损失的信息我们可以看到损失值持续上升,这就说明哪里出了问题。但是如果所有的测试样例都通过了,就说明我们的实现是没有问题的。运行下面的测试样例,观察哪里出了问题。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
# 测试样例
Wt, bt = initialize(trainX.shape[1])
print('epoch 0, W:', Wt) # [-0.00348894 0.00983703 0.00580923]
print('epoch 0, b:', bt) # [ 0.]
print()

Zt = forward(trainX, Wt, bt)
dWt, dbt = compute_gradient(trainX, Zt, trainY)
print('dWt:', dWt) # [ -4.18172940e+09 -2.19880296e+08 -1.94481031e+08]
print('db:', dbt) # -182154.277882
print()

update(dWt, dbt, Wt, bt, 0.01)
print('epoch 1, W:', Wt) # [ 41817293.96016914 2198802.97412493 1944810.31544994]
print('epoch 1, b:', bt) # [ 1821.54277882]

可以看到,我们最开始的参数都是在 $10^{-3}$ 这个数量级上,而第一轮迭代时计算出的梯度的数量级在 $10^8$ 左右,这就导致使用梯度下降更新的时候,让参数变成了 $10^6$ 这个数量级左右(学习率为0.01)。产生这样的问题的主要原因是:我们的原始数据 $X$ 没有经过适当的处理,直接扔到了神经网络中进行训练,导致在计算梯度时,由于 $X$ 的数量级过大,导致梯度的数量级变大,在参数更新时使得参数的数量级不断上升,导致参数无法收敛。

解决的方法也很简单,对参数进行归一化处理,将其标准化,使均值为0,缩放到 $[-1, 1]$附近。

10. 标准化处理

标准化处理和第一题一样

1
2
3
4
from sklearn.preprocessing import StandardScaler
stand = StandardScaler()
trainX_normalized = stand.fit_transform(trainX)
testX_normalized = stand.transform(testX)

重新训练模型,这次我们迭代40轮,学习率设置为0.1

1
2
3
m = trainX.shape[1]
W, b = initialize(m)
training_loss_list, testing_loss_list = train(trainX_normalized, trainY, testX_normalized, testY, W, b, 40, learning_rate = 0.1, verbose = False)

打印损失值变化曲线

1
plot_loss_curve(training_loss_list, testing_loss_list)

计算测试集上的MSE

1
2
prediction = forward(testX_normalized, W, b)
mse(testY, prediction) ** 0.5

第三题:神经网络:对数几率回归

实验内容:

  1. 完成对数几率回归
  2. 使用梯度下降求解模型参数
  3. 绘制模型损失值的变化曲线
  4. 调整学习率和迭代轮数,观察损失值曲线的变化
  5. 按照给定的学习率和迭代轮数,初始化新的参数,绘制其训练集和测试集损失值的变化曲线,完成表格内精度的填写

对数几率回归,二分类问题的分类算法,属于线性模型中的一种,我们可以将其抽象为最简单的神经网络。

只有一个输入层和一个输出层,还有一个激活函数,$\rm sigmoid$,简记为$\sigma$。
我们设输入为$X \in \mathbb{R}^{n \times m}$,输入层到输出层的权重为$W \in \mathbb{R}^{m}$,偏置$b \in \mathbb{R}$。

激活函数

这个激活函数,会将输出层的神经元的输出值转换为一个 $(0, 1)$ 区间内的数。

因为是二分类问题,我们设类别为0和1,我们将输出值大于0.5的样本分为1类,输出值小于0.5的类分为0类。

前向传播

其中,$O \in \mathbb{R}^{n}$为输出层的结果,$\sigma$为$\rm sigmoid$激活函数。

注意:这里我们其实是做了广播,将$b$复制了$n-1$份后拼接成了维数为$n$的向量。

所以对数几率回归就可以写为:

损失函数

使用对数损失函数,因为对数损失函数较其他损失函数有更好的性质,感兴趣的同学可以去查相关的资料。

针对二分类问题的对数损失函数:

在这个对数几率回归中,我们的损失函数对所有样本取个平均值:

注意,这里我们的提到的$\log$均为$\ln$,在numpy中为np.log

因为我们的类别只有0和1,所以在这个对数损失函数中,要么前一项为0,要么后一项为0。

如果当前样本的类别为0,那么前一项就为0,损失函数变为 $- \log{(1 - \hat{y})}$ ,因为我们的预测值 $0 < \hat{y} < 1$ ,所以 $0 < 1 - \hat{y} < 1$ ,$- \log{(1 - \hat{y})} > 0$ ,为了降低损失值,模型需要让预测值 $\hat{y}$不断地趋于0。

同理,如果当前样本的类别为1,那么降低损失值就可以使模型的预测值趋于1。

参数更新

求得损失函数对参数的偏导数后,我们就可以使用梯度下降进行参数更新:

其中,$\alpha$ 是学习率,一般设置为0.1,0.01等。

经过一定次数的迭代后,参数会收敛至最优点。这种基于梯度的优化算法很常用,训练神经网络主要使用这类优化算法。

反向传播

我们使用梯度下降更新参数$W$和$b$。为此需要求得损失函数对参数$W$和$b$的偏导数,根据链式法则有:

这里我们一项一项求,先求第一项:

第二项:

第三项:

综上:

同理,求$\rm loss$对$b$的偏导数:

注意,由于$b$是被广播成$n \times K$的矩阵,因此实际上$b$对每个样本的损失都有贡献,因此对其求偏导时,要把$n$个样本对它的偏导数加和。

这样,我们就得到了损失函数对参数的偏导数,然后就可以使用梯度下降算法更新参数

1. 导入数据集

1
2
3
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline

我们生成半月形数据

1
2
from sklearn.datasets import make_moons
X, y = make_moons(n_samples = 2000, noise = 0.3, random_state=0)

选择40%的数据作为测试集,60%作为训练集

1
2
3
4
from sklearn.model_selection import train_test_split
trainX, testX, trainY, testY = train_test_split(X, y, test_size = 0.4, random_state = 32)
trainY = trainY
testY = testY
1
trainX.shape, trainY.shape, testX.shape, testY.shape

2. 数据预处理

使用和第一题一样的预处理方式

1
2
3
4
from sklearn.preprocessing import StandardScaler
s = StandardScaler()
trainX = s.fit_transform(trainX)
testX = s.transform(testX)

3. 定义神经网络

3.1 参数初始化

我们需要对神经网络的参数进行初始化,这个网络中只有两个参数,一个$W \in \mathbb{R}^{m}$,一个$b \in \mathbb{R}$。初始化的时候,我们将参数W随机初始化,参数b初始化为0。为什么要对神经网络的参数进行随机初始化,感兴趣的同学可以去查相关的资料。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
def initialize(m):
'''
初始化参数W和参数b

Returns
----------
W: np.ndarray, shape = (m, ),参数W

b: np.ndarray, shape = (1, ),参数b

'''
np.random.seed(32)
W = np.random.normal(size = (m, )) * 0.01
b = np.zeros((1, ))
return W, b
1
2
3
4
# 测试样例
Wt, bt = initialize(trainX.shape[1])
print(Wt.shape) # (2,)
print(bt.shape) # (1,)

3.2 前向传播

接下来我们要定义神经网络前向传播的过程。

首先计算$Z = XW + b$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
def linear_combination(X, W, b):
'''
完成Z = XW + b的计算

Parameters
----------
X: np.ndarray, shape = (n, m),输入的数据

W: np.ndarray, shape = (m, ),权重

b: np.ndarray, shape = (1, ),偏置

Returns
----------
Z: np.ndarray, shape = (n, ),线性组合后的值

'''

Z = np.dot(X, W) + b # YOUR CODE HERE

return Z
1
2
3
# 测试样例
Wt, bt = initialize(trainX.shape[1])
linear_combination(trainX, Wt, bt).shape #(1200,)

接下来实现激活函数$\rm sigmoid$

1
2
3
4
5
6
7
8
9
10
11
12
13
def my_sigmoid(x):
'''
simgoid 1 / (1 + exp(-x))

Parameters
----------
X: np.ndarray, 待激活的值

'''
# YOUR CODE HERE
activations = 1 / (1 + np.exp(-x))

return activations
1
2
3
4
# 测试样例
Wt, bt = initialize(trainX.shape[1])
Zt = linear_combination(trainX, Wt, bt)
my_sigmoid(Zt).mean() # 0.49999

在实现$\rm sigmoid$的时候,可能会遇到上溢(overflow)的问题,可以看到$\rm sigmoid$中有一个指数运算

当$x$很大的时候,我们使用numpy.exp(x)会直接溢出

1
np.exp(1e56)
1
my_sigmoid(np.array([-1e56]))

虽说程序没有报错,只是抛出了warning,但还是应该解决一下。

解决这种问题的方法有很多,比如,我们可以将$\rm sigmoid$进行变换:

其中,$\mathrm{tanh}(x) = \frac{\mathrm{sinh}(x)}{\mathrm{cosh}(x)} = \frac{e^x - e^{-x}}{e^x + e^{-x}}$

转换成这种形式后,我们就可以直接利用numpy.tanh完成$\rm sigmoid$的计算,就不会产生上溢的问题了。

除此以外,最好的解决方法是使用scipy中的expit函数,完成$\rm sigmoid$的计算。我们现在做的都是神经网络底层相关的运算,很容易出现数值不稳定性相关的问题,最好的办法就是使用别人已经实现好的函数,这样就能减少我们很多的工作量,同时又快速地完成任务。

1
from scipy.special import expit
1
2
def sigmoid(X):
return expit(X)
1
2
# 测试样例
sigmoid(np.array([-1e56]))

接下来完成整个前向传播的函数,也就是 $Z = XW+b$ 和 $\hat{y} = \mathrm{sigmoid}(Z)$

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
def forward(X, W, b):
'''
完成输入矩阵X到最后激活后的预测值y_pred的计算过程

Parameters
----------
X: np.ndarray, shape = (n, m),数据,一行一个样本,一列一个特征

W: np.ndarray, shape = (m, ),权重

b: np.ndarray, shape = (1, ),偏置

Returns
----------
y_pred: np.ndarray, shape = (n, ),模型对每个样本的预测值

'''
# 求Z
Z = linear_combination(X, W, b) # YOUR CODE HERE

# 求激活后的预测值
y_pred = sigmoid(Z) # YOUR CODE HERE

return y_pred
1
2
3
# 测试样例
Wt, bt = initialize(trainX.shape[1])
forward(trainX, Wt, bt).mean() # 0.4999(没有四舍五入)

接下来完成损失函数的编写,我们使用的是对数损失,这里需要注意的一个问题是:

在这个对数损失中,$\hat{y}$中不能有$0$和$1$,如果有$0$,那么损失函数中的前半部分,$\log{0}$就会出错,如果有$1$,那么后半部分$\log{(1-1)}$就会出错。

所以我们要先将$\hat{y}$中的$0$和$1$改变一下,把$0$变成一个比较小但是大于$0$的数,把$1$变成小于$1$但是足够大的数。使用numpy.clip函数就可以作到这点。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
def logloss(y_true, y_pred):
'''
给定真值y,预测值y_hat,计算对数损失并返回

Parameters
----------
y_true: np.ndarray, shape = (n, ), 真值

y_pred: np.ndarray, shape = (n, ),预测值

Returns
----------
loss: float, 损失值

'''
# 下面这句话会把y_pred里面小于1e-10的数变成1e-10,大于1 - 1e-10的数变成1 - 1e-10
y_hat = np.clip(y_pred, 1e-10, 1 - 1e-10)

# 求解对数损失
loss = - np.sum(y_true * np.log(y_hat) + (1 - y_true) * np.log(1 - y_hat)) / len(y_true) # YOUR CODE HERE

return loss
1
2
3
# 测试样例
Wt, bt = initialize(trainX.shape[1])
logloss(trainY, forward(trainX, Wt, bt)) # 0.69740

3.3 反向传播

我们接下来要完成损失函数对参数的偏导数的计算

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
def compute_gradient(y_true, y_pred, X):
'''
给定预测值y_pred,真值y_true,传入的输入数据X,计算损失函数对参数W的偏导数的导数值dW,以及对b的偏导数的导数值db

Parameters
----------
y_true: np.ndarray, shape = (n, ), 真值

y_pred: np.ndarray, shape = (n, ),预测值

X: np.ndarray, shape = (n, m),数据,一行一个样本,一列一个特征

Returns
----------
dW: np.ndarray, shape = (m, ), 损失函数对参数W的偏导数

db: float, 损失函数对参数b的偏导数

'''
# 求损失函数对参数W的偏导数的导数值
dW = np.dot(X.T, (y_pred - y_true)) / len(y_pred) # YOUR CODE HERE

# 求损失函数对参数b的偏导数的导数值
db = np.sum(y_pred - y_true) / len(y_pred) # YOUR CODE HERE

return dW, db
1
2
3
4
5
6
# 测试样例
Wt, bt = initialize(trainX.shape[1])
dWt, dbt = compute_gradient(trainY, forward(trainX, Wt, bt), trainX)
print(dWt.shape) # (2, )
print(dWt.sum()) # 0.04625
print(dbt) # 0.00999

3.4 参数更新

给定学习率,结合上一步求出的偏导数,完成梯度下降的更新公式

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
def update(W, b, dW, db, learning_rate):
'''
梯度下降,给定参数W,参数b,以及损失函数对他们的偏导数,使用梯度下降更新参数W和参数b

Parameters
----------
W: np.ndarray, shape = (m, ),参数W

b: np.ndarray, shape = (1, ),参数b

dW: np.ndarray, shape = (m, ), 损失函数对参数W的偏导数

db: float, 损失函数对参数b的偏导数

learning_rate, float,学习率

'''
# 对参数W进行更新
W -= learning_rate * dW

# 对参数b进行更新
b -= learning_rate * db # YOUR CODE HERE
1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 测试样例
Wt, bt = initialize(trainX.shape[1])
print(Wt) # [-0.00348894 0.00983703]
print(bt) # [ 0.]
print()

dWt, dbt = compute_gradient(trainY, forward(trainX, Wt, bt), trainX)
print(dWt) # [-0.28650366 0.33276308]
print(dbt) # 0.00999999939463
print()

update(Wt, bt, dWt, dbt, 0.01)
print(Wt) # [-0.00062391 0.0065094 ]
print(bt) # [ -9.99999939e-05]

我们来完成整个反向传播和更新参数的函数

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
def backward(y_true, y_pred, X, W, b, learning_rate):
'''
反向传播,包含了计算损失函数对各个参数的偏导数的过程,以及梯度下降更新参数的过程

Parameters
----------
y_true: np.ndarray, shape = (n, ), 真值

y_pred: np.ndarray, shape = (n, ),预测值

X: np.ndarray, shape = (n, m),数据,一行一个样本,一列一个特征

W: np.ndarray, shape = (m, ),参数W

b: np.ndarray, shape = (1, ),参数b

dW: np.ndarray, shape = (m, ), 损失函数对参数W的偏导数

db: float, 损失函数对参数b的偏导数

learning_rate, float,学习率

'''
# 求参数W和参数b的梯度
dW, db = compute_gradient(y_true, y_pred, X)

# 梯度下降
update(W, b, dW, db, learning_rate)
1
2
3
4
5
6
7
8
9
10
11
# 测试样例
Wt, bt = initialize(trainX.shape[1])
y_predt = forward(trainX, Wt, bt)
loss_1 = logloss(trainY, y_predt)
print(loss_1) # 0.697403529518

backward(trainY, y_predt, trainX, Wt, bt, 0.01)

y_predt = forward(trainX, Wt, bt)
loss_2 = logloss(trainY, y_predt)
print(loss_2) # 0.695477626714

4. 训练函数的编写

我们已经实现了完成训练需要的子函数,接下来就是组装了

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
def train(trainX, trainY, testX, testY, W, b, epochs, learning_rate = 0.01, verbose = False):
'''
训练,我们要迭代epochs次,每次迭代的过程中,做一次前向传播和一次反向传播
同时记录训练集和测试集上的损失值,后面画图用

Parameters
----------
trainX: np.ndarray, shape = (n, m), 训练集

trainY: np.ndarray, shape = (n, ), 训练集标记

testX: np.ndarray, shape = (n_test, m),测试集

testY: np.ndarray, shape = (n_test, ),测试集的标记

W: np.ndarray, shape = (m, ),参数W

b: np.ndarray, shape = (1, ),参数b

epochs: int, 要迭代的轮数

learning_rate: float, default 0.01,学习率

verbose: boolean, default False,是否打印损失值

Returns
----------
training_loss_list: list(float),每迭代一次之后,训练集上的损失值

testing_loss_list: list(float),每迭代一次之后,测试集上的损失值

'''

training_loss_list = []
testing_loss_list = []

for i in range(epochs):

# 计算训练集前向传播得到的预测值
# YOUR CODE HERE
train_y_pred = forward(trainX, W, b)

# 计算当前训练集的损失值
# YOUR CODE HERE
training_loss = logloss(trainY, train_y_pred)

# 计算测试集前向传播得到的预测值
# YOUR CODE HERE
test_y_pred = forward(testX, W, b)

# 计算当前测试集的损失值
testing_loss = logloss(testY, test_y_pred)

if verbose == True:
print('epoch %s, training loss:%s'%(i + 1, training_loss))
print('epoch %s, testing loss:%s'%(i + 1, testing_loss))
print()

# 保存损失值
training_loss_list.append(training_loss)
testing_loss_list.append(testing_loss)

# 反向传播更新参数
# YOUR CODE HERE
backward(trainY, train_y_pred, trainX, W, b, learning_rate)

return training_loss_list, testing_loss_list
1
2
3
4
5
# 测试样例
Wt, bt = initialize(trainX.shape[1])
training_loss_list, testing_loss_list = train(trainX, trainY, testX, testY, Wt, bt, 2, 0.1)
print(training_loss_list) # [0.69740352951773121, 0.67843729060725722]
print(testing_loss_list) # [0.69743661286103986, 0.67880126235588389]

5. 绘制模型损失值变化曲线

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
def plot_loss_curve(training_loss_list, testing_loss_list):
'''
绘制损失值变化曲线

Parameters
----------
training_loss_list: list(float),每迭代一次之后,训练集上的损失值

testing_loss_list: list(float),每迭代一次之后,测试集上的损失值

'''
plt.figure(figsize = (10, 6))
plt.plot(training_loss_list, label = 'training loss')
plt.plot(testing_loss_list, label = 'testing loss')
plt.xlabel('epoch')
plt.ylabel('loss')
plt.legend()

6. 预测

接下来编写一个预测的函数,事实上,$\rm sigmoid$输出的是当前这个样本为正例的概率,也就是说,这个输出值是一个0到1的值,一般我们将大于0.5的值变成1,小于0.5的值变成0,也就是说,如果当前输出的概率值大于0.5,那我们认为这个样本的类别就是1,否则就是0,这样输出的就是类标了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
def predict(X, W, b):
'''
预测,调用forward函数完成神经网络对输入X的计算,然后完成类别的划分,大于0.5的变为1,小于等于0.5的变为0

Parameters
----------
X: np.ndarray, shape = (n, m), 训练集

W: np.ndarray, shape = (m, 1),参数W

b: np.ndarray, shape = (1, ),参数b

Returns
----------
prediction: np.ndarray, shape = (n, 1),预测的标记

'''
prediction = (forward(testX, W, b) > 0.5).astype('uint8') # YOUR CODE HERE
return prediction
1
2
3
4
5
# 测试样例
from sklearn.metrics import accuracy_score
Wt, bt = initialize(trainX.shape[1])
predictiont = predict(testX, Wt, bt)
accuracy_score(testY, predictiont) # 0.16250000000000001

7. 训练一个神经网络

我们的学习率是0.01,迭代200轮

1
2
W, b = initialize(trainX.shape[1])
training_loss_list, testing_loss_list = train(trainX, trainY, testX, testY, W, b, 200, 0.01)

计算测试集精度

1
2
prediction = predict(testX, W, b)
accuracy_score(testY, prediction) # 0.83625000000000005

绘制损失值变化曲线

1
plot_loss_curve(training_loss_list, testing_loss_list)

test:初始化新的参数,学习率和迭代轮数按下表设置,绘制其训练集和测试集损失值的变化曲线,完成表格内精度的填写

双击此处填写
学习率 迭代轮数 测试集精度
0.0001 200 0.3325
0.1 1000 0.84
1
2
3
4
5
6
# YOUR CODE HERE
W2, b2 = initialize(trainX.shape[1])
training_loss_list, testing_loss_list = train(trainX, trainY, testX, testY, W2, b2, 200, 0.0001)
prediction = predict(testX, W2, b2)
print(accuracy_score(prediction, testY))
plot_loss_curve(training_loss_list, testing_loss_list)
1
2
3
4
5
6
# YOUR CODE HERE
W2, b2 = initialize(trainX.shape[1])
training_loss_list, testing_loss_list = train(trainX, trainY, testX, testY, W2, b2, 1000, 0.1)
prediction = predict(testX, W2, b2)
print(accuracy_score(prediction, testY))
plot_loss_curve(training_loss_list, testing_loss_list)

第四题:神经网络:三层感知机

实现内容:

  1. 实现一个三层感知机
  2. 对手写数字数据集进行分类
  3. 绘制损失值变化曲线

在这道题中,我们要实现一个三层感知机

前向传播

我们实现一个最简单的三层感知机,一个输入层,一个隐藏层,一个输出层,隐藏层单元个数为$h$个,输出层有$K$个单元。

  1. 我们将第一层的输入,定义为$X \in \mathbb{R}^{n \times m}$,n个样本,m个特征。
  2. 输入层到隐藏层之间的权重(weight)与偏置(bias),分别为$W_1 \in \mathbb{R}^{m \times h}$,$b_1 \in \mathbb{R}^{1 \times h}$。
  3. 隐藏层到输出层的权重和偏置分为别$W_2 \in \mathbb{R}^{h \times K}$,$b_2 \in \mathbb{R}^{1 \times K}$。

隐藏层的激活函数选用ReLU

我们用$H_1$表示第一个隐藏层的输出值,$O$表示输出层的输出值,这样,前向传播即可定义为

其中,$H_1 \in \mathbb{R}^{n \times h}$,$O \in \mathbb{R}^{n \times K}$。

注意:这里我们其实是做了广播,将$b_1$复制了$n-1$份后拼接成了维数为$n \times h$的矩阵,同理,$b_2$也做了广播,拼成了$n \times K$的矩阵。

最后一层的输出,使用softmax函数激活,得到神经网络计算出的各类的概率值:

其中,$\hat{y_i}$表示第$i$类的概率值,也就是输出层第$i$个神经元经$\mathrm{softmax}$激活后的值。

损失函数

损失函数使用交叉熵损失函数:

这样,$n$个样本的平均损失为:

注意,这里我们的提到的$\log$均为$\ln$,在numpy中为np.log

反向传播

我们使用梯度下降训练模型,求解方式就是求出损失函数对参数的偏导数,即参数的梯度,然后将参数减去梯度乘以学习率,进行参数的更新。

其中,$\alpha$是学习率。

在这道题中,交叉熵损失函数的求导比较麻烦,我们先求神经网络的输出层的偏导数,写成链式法则的形式:

首先求解第一项:

然后求解第二项,因为$\hat{y_k}$的分母是$\sum_k \exp{(O_k)}$,里面包含$O_i$,所以每一个$\hat{y_k}$的分母都包含$O_i$,这就要求反向传播的时候需要考虑这$K$项,将这$K$项的偏导数加在一起。

这$K$项分别为:$\frac{\exp{(O_1)}}{\sum_k \exp{(O_k)}}$,$\frac{\exp{(O_2)}}{\sum_k \exp{(O_k)}}$,…,$\frac{\exp{(O_i)}}{\sum_k \exp{(O_k)}}$,…,$\frac{\exp{(O_k)}}{\sum_k \exp{(O_k)}}$。

显然,这里只有分子带有$O_i$的这项与其他的项不同,因为分子和分母同时包含了$O_i$,而其他的项只有分母包含了$O_i$。

这就需要在求解$\frac{\partial \hat{y}}{\partial O_i}$的时候分两种情况讨论

  1. 分子带$O_i$
  2. 分子不带$O_i$

第一种情况,当分子含有$O_i$时:

第二种情况,当分子不含$O_i$时,我们用$j$表示当前项的下标:

这样,$\mathrm{loss}$对$O_i$的偏导数即为:

由于我们处理的多类分类任务,一个样本只对应一个标记,所以$\sum^K_{k = 1} y_k = 1$,上式在这种问题中,即可化简为:

将其写成矩阵表达式:

也就是说,我们的损失函数对输出层的$K$个神经单元的偏导数为$\mathrm{softmax}$激活值减去真值。

接下来我们需要求损失函数对参数$W_2$和$b_2$的偏导数

其中,$\frac{\partial loss}{\partial W_2} \in \mathbb{R}^{h \times K}$,$\frac{\partial loss}{\partial b_2} \in \mathbb{R}^{1 \times K}$。
注意,由于$b_2$是被广播成$n \times K$的矩阵,因此实际上$b_2$对每个样本的损失都有贡献,因此对其求偏导时,要把$n$个样本对它的偏导数加和。

同理,我们可以求得$\mathrm{loss}$对$W_1$和$b_1$的偏导数:

由于我们使用的是$\mathrm{ReLU}$激活函数,它的偏导数为:

所以上式为:

其中,${W_1}_{ij}$表示矩阵$W_1$第$i$行第$j$列的值,${Z}_{ij}$表示矩阵$Z$第$i$行第$j$列的值。
同理:

其中,$\frac{\partial loss}{\partial W_1} \in \mathbb{R}^{m \times h}$,$\frac{\partial loss}{\partial b_1} \in \mathbb{R}^{1 \times h}$。

参数更新

求得损失函数对四个参数的偏导数后,我们就可以使用梯度下降进行参数更新:

其中,$\alpha$是学习率

以上内容,就是一个三层感知机的前向传播与反向传播过程。

1. 导入数据

使用第一题的手写数字数据集

1
2
import matplotlib.pyplot as plt
%matplotlib inline
1
from time import time
1
2
3
import numpy as np
from sklearn.datasets import load_digits
from sklearn.model_selection import train_test_split

40%做测试集,60%做训练集

1
trainX, testX, trainY, testY = train_test_split(load_digits()['data'], load_digits()['target'], test_size = 0.4, random_state = 32)
1
trainX.shape, trainY.shape, testX.shape, testY.shape

2. 数据预处理

使用和第一题一样的标准化处理方法

1
2
3
4
from sklearn.preprocessing import StandardScaler
s = StandardScaler()
trainX = s.fit_transform(trainX)
testX = s.transform(testX)

接下来还要处理输出。
我们的神经网络是针对每个样本,输出其分别属于$K$类的概率,我们要找最大的那个概率,对应的是哪个类。
我们当前的trainY和testY,每个样本都是一个类标,我们需要将其变成one_hot编码,也就是,假设当前样本的类别是3,我们需要把它变成一个长度为10的向量,其中第4个元素为1,其他元素都为0。得到的矩阵分别记为trainY_mat和testY_mat。
这样,模型训练完成后,会针对每个样本输出十个数,分别代表这个样本属于$0,1,…,9$的概率,那我们只要取最大的那个数的下标,就知道模型认为这个样本是哪类了。

1
2
3
4
5
trainY_mat = np.zeros((len(trainY), 10))
trainY_mat[np.arange(0, len(trainY), 1), trainY] = 1

testY_mat = np.zeros((len(testY), 10))
testY_mat[np.arange(0, len(testY), 1), testY] = 1
1
trainY_mat.shape, testY_mat.shape

3. 参数初始化

这题和上一题的区别是,我们把参数用dict存起来

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
def initialize(h, K):
'''
参数初始化

Parameters
----------
h: int: 隐藏层单元个数

K: int: 输出层单元个数

Returns
----------
parameters: dict,参数,键是"W1", "b1", "W2", "b2"

'''
np.random.seed(32)
W_1 = np.random.normal(size = (trainX.shape[1], h)) * 0.01
b_1 = np.zeros((1, h))

np.random.seed(32)
W_2 = np.random.normal(size = (h, K)) * 0.01
b_2 = np.zeros((1, K))

parameters = {'W1': W_1, 'b1': b_1, 'W2': W_2, 'b2': b_2}

return parameters
1
2
3
4
5
6
# 测试样例
parameterst = initialize(50, 10)
print(parameterst['W1'].shape) # (64, 50)
print(parameterst['b1'].shape) # (1, 50)
print(parameterst['W2'].shape) # (50, 10)
print(parameterst['b2'].shape) # (1, 10)

4. 前向传播

完成Z的计算

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
def linear_combination(X, W, b):
'''
计算Z,Z = XW + b

Parameters
----------
X: np.ndarray, shape = (n, m),输入的数据

W: np.ndarray, shape = (m, h),权重

b: np.ndarray, shape = (1, h),偏置

Returns
----------
Z: np.ndarray, shape = (n, h),线性组合后的值

'''

# Z = XW + b
Z = np.dot(X, W) + b

return Z
1
2
3
4
5
# 测试样例
parameterst = initialize(50, 10)
Zt = linear_combination(trainX, parameterst['W1'], parameterst['b1'])
print(Zt.shape) # (1078, 50)
print(Zt.mean()) # -5.27304442123e-19

$\rm ReLU$激活函数

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
def ReLU(X):
'''
ReLU激活函数

Parameters
----------
X: np.ndarray,待激活的矩阵

Returns
----------
activations: np.ndarray, 激活后的矩阵

'''
activations = X.copy()
activations[activations < 0] = 0
return activations
1
2
3
4
5
6
7
8
9
# 测试样例
parameterst = initialize(50, 10)
Zt = linear_combination(trainX, parameterst['W1'], parameterst['b1'])
Ht = ReLU(Zt)
print(Ht.mean()) # 0.0304

Ot = linear_combination(Ht, parameterst['W2'], parameterst['b2'])
print(Ot.shape) # (1078, 10)
print(Ot.mean()) # 0.0006

$\rm softmax$激活

1
2
3
4
5
6
def my_softmax(O):
'''
softmax激活
'''
denominator = np.exp(O).sum(axis = 1, keepdims = True)
return np.exp(O) / denominator
1
2
3
4
5
# 测试样例
test1 = np.array([[-1e32, -1e32, -1e32]])
test2 = np.array([[1e32, 1e32, 1e32]])
print(my_softmax(test1))
print(my_softmax(test2))

这里,其实是有数值计算上的问题的,假设,我们最后的输出有三个数,每个数都特别小,理论上来说,通过$\rm softmax$激活后,三个值都是$\frac{1}{3}$。但实际上就不是这样了,实际上会导致分母为0,除法就不能做了。如果每个数都特别大,会导致做指数运算的时候上溢。

我们需要用其他的方法来实现$\rm softmax$。

我们将传入$\rm softmax$的向量,每个元素减去他们中的最大值,即

这个式子是成立的,感兴趣的同学可以证明一下上面的式子。

当我们做了这样的变换后,向量$O$中的最大值就变成了0,就不会上溢了,而分母中最少有一项为1,也不会出现下溢导致分母为0的问题了。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
def softmax(O):
'''
softmax激活函数

Parameters
----------
O: np.ndarray,待激活的矩阵

Returns
----------
activations: np.ndarray, 激活后的矩阵

'''

# YOUR CODE HEER
t = O - np.max(O, axis = 1, keepdims = True)
denominator = np.exp(t).sum(axis = 1, keepdims = True)

activations = np.exp(t) / denominator
return activations
1
2
3
4
5
6
7
8
9
10
# 测试样例
parameterst = initialize(50, 10)
Zt = linear_combination(trainX, parameterst['W1'], parameterst['b1'])
Ht = ReLU(Zt)
Ot = linear_combination(Ht, parameterst['W2'], parameterst['b2'])
y_pred = softmax(Ot)

print(y_pred.shape) # (1078, 10)
print(Ot.mean()) # 0.000600192658464
print(y_pred.mean()) # 0.1

接下来是实现损失函数,交叉熵损失函数:

这里又会出一个问题,交叉熵损失函数中,我们需要对$\rm softmax$的激活值取对数,也就是$\log{\hat{y}}$,这就要求我们的激活值全都是大于0的数,不能等于0,但是我们实现的$\rm softmax$在有些时候确实会输出0,比如:

1
softmax(np.array([[1e32, 0, -1e32]]))

这就使得在计算loss的时候会出现问题,解决这个问题的方法是$\rm log \ softmax$。所谓$\rm log \ softmax$,就是将交叉熵中的对数运算与$\rm softmax$结合起来,避开为0的情况

这样我们再计算$\rm loss$的时候就可以把输出层的输出直接放到$\rm log \ softmax$中计算,不用先激活,再取对数了。

我们先编写log_softmax

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
def log_softmax(x):
'''
log softmax

Parameters
----------
x: np.ndarray,待激活的矩阵

Returns
----------
log_activations: np.ndarray, 激活后取了对数的矩阵

'''
# 获取每行的最大值
max_ = np.max(x, axis = 1, keepdims = True)

# 指数运算
exp_x = np.exp(x - max_)

# 每行求和
Z = np.sum(exp_x, axis = 1, keepdims = True)

# 求log softmax
log_activations = x - max_ - np.log(Z)

return log_activations
1
2
3
4
5
6
7
8
# 测试样例
parameterst = initialize(50, 10)
Zt = linear_combination(trainX, parameterst['W1'], parameterst['b1'])
Ht = ReLU(Zt)
Ot = linear_combination(Ht, parameterst['W2'], parameterst['b2'])
t = log_softmax(Ot)
print(t.shape) # (1078, 10)
print(t.mean()) # -2.30259148717

然后编写cross_entropy_with_softmax

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
def cross_entropy_with_softmax(y_true, O):
'''
求解交叉熵损失函数,这里需要使用log softmax,所以参数分别是真值和未经softmax激活的输出值

Parameters
----------
y_true: np.ndarray,shape = (n, K), 真值

O: np.ndarray, shape = (n, K),softmax激活前的输出层的输出值

Returns
----------
loss: float, 平均的交叉熵损失值

'''

# 平均交叉熵损失
loss = - np.sum(log_softmax(O) * y_true) / len(y_true)

return loss
1
2
3
4
5
6
7
# 测试样例
parameterst = initialize(50, 10)
Zt = linear_combination(trainX, parameterst['W1'], parameterst['b1'])
Ht = ReLU(Zt)
Ot = linear_combination(Ht, parameterst['W2'], parameterst['b2'])
losst = cross_entropy_with_softmax(trainY_mat, Ot)
print(losst.mean()) # 2.30266707958
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
def forward(X, parameters):
'''
前向传播,从输入一直到输出层softmax激活前的值

Parameters
----------
X: np.ndarray, shape = (n, m),输入的数据

parameters: dict,参数

Returns
----------
O: np.ndarray, shape = (n, K),softmax激活前的输出层的输出值

'''
# 输入层到隐藏层
Z = linear_combination(X, parameters['W1'], parameters['b1'])

# 隐藏层的激活
H = ReLU(Z)

# 隐藏层到输出层
O = linear_combination(H, parameters['W2'], parameters['b2'])

return O
1
2
3
4
# 测试样例
parameterst = initialize(50, 10)
Ot = forward(trainX, parameterst)
print(Ot.mean()) # 0.000600192658464

5. 反向传播

先计算梯度

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
def compute_gradient(y_true, y_pred, H, Z, X, parameters):
'''
计算梯度

Parameters
----------
y_true: np.ndarray,shape = (n, K), 真值

y_pred: np.ndarray, shape = (n, K),softmax激活后的输出层的输出值

H: np.ndarray, shape = (n, h),隐藏层激活后的值

Z: np.ndarray, shape = (n, h), 隐藏层激活前的值

X: np.ndarray, shape = (n, m),输入的原始数据

parameters: dict,参数

Returns
----------
grads: dict, 梯度

'''

# 计算W2的梯度
dW2 = np.dot(H.T, (y_pred - y_true)) / len(y_pred)

# 计算b2的梯度
db2 = np.sum(y_pred - y_true, axis = 0) / len(y_pred)

# 计算ReLU的梯度
relu_grad = Z.copy()
relu_grad[relu_grad < 0] = 0
relu_grad[relu_grad >= 0] = 1

# 计算W1的梯度
dW1 = np.dot(X.T, np.dot(y_pred - y_true, parameters['W2'].T) * relu_grad) / len(y_pred)

# 计算b1的梯度
db1 = np.sum((np.dot(y_pred - y_true, parameters['W2'].T) * relu_grad), axis = 0) / len(y_pred)

grads = {'dW2': dW2, 'db2': db2, 'dW1': dW1, 'db1': db1}

return grads
1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 测试样例
parameterst = initialize(50, 10)

Zt = linear_combination(trainX, parameterst['W1'], parameterst['b1'])
Ht = ReLU(Zt)
Ot = linear_combination(Ht, parameterst['W2'], parameterst['b2'])
y_predt = softmax(Ot)

gradst = compute_gradient(trainY_mat, y_predt, Ht, Zt, trainX, parameterst)

print(gradst['dW1'].sum()) # 0.0429186117668
print(gradst['db1'].sum()) # -5.05985151857e-05
print(gradst['dW2'].sum()) # -2.16840434497e-18
print(gradst['db2'].sum()) # -1.34441069388e-17

梯度下降,参数更新

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
def update(parameters, grads, learning_rate):
'''
参数更新

Parameters
----------
parameters: dict,参数

grads: dict, 梯度

learning_rate: float, 学习率

'''
parameters['W2'] -= learning_rate * grads['dW2']
parameters['b2'] -= learning_rate * grads['db2']
parameters['W1'] -= learning_rate * grads['dW1']
parameters['b1'] -= learning_rate * grads['db1']

反向传播,参数更新

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
# 测试样例
parameterst = initialize(50, 10)
print(parameterst['W1'].sum()) # 0.583495454481
print(parameterst['b1'].sum()) # 0.0
print(parameterst['W2'].sum()) # 0.1888716431
print(parameterst['b2'].sum()) # 0.0
print()

Zt = linear_combination(trainX, parameterst['W1'], parameterst['b1'])
Ht = ReLU(Zt)
Ot = linear_combination(Ht, parameterst['W2'], parameterst['b2'])
y_predt = softmax(Ot)

gradst = compute_gradient(trainY_mat, y_predt, Ht, Zt, trainX, parameterst)
update(parameterst, gradst, 0.1)

print(parameterst['W1'].sum()) # 0.579203593304
print(parameterst['b1'].sum()) # 5.05985151857e-06
print(parameterst['W2'].sum()) # 0.1888716431
print(parameterst['b2'].sum()) # 1.24683249836e-18
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
def backward(y_true, y_pred, H, Z, X, parameters, learning_rate):
'''
计算梯度,参数更新

Parameters
----------
y_true: np.ndarray,shape = (n, K), 真值

y_pred: np.ndarray, shape = (n, K),softmax激活后的输出层的输出值

H: np.ndarray, shape = (n, h),隐藏层激活后的值

Z: np.ndarray, shape = (n, h), 隐藏层激活前的值

X: np.ndarray, shape = (n, m),输入的原始数据

parameters: dict,参数

learning_rate: float, 学习率

'''
grads = compute_gradient(y_true, y_pred, H, Z, X, parameters)
update(parameters, grads, learning_rate)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
# 测试样例
parameterst = initialize(50, 10)
print(parameterst['W1'].sum()) # 0.583495454481
print(parameterst['b1'].sum()) # 0.0
print(parameterst['W2'].sum()) # 0.1888716431
print(parameterst['b2'].sum()) # 0.0
print()

Zt = linear_combination(trainX, parameterst['W1'], parameterst['b1'])
Ht = ReLU(Zt)
Ot = linear_combination(Ht, parameterst['W2'], parameterst['b2'])
y_predt = softmax(Ot)

backward(trainY_mat, y_predt, Ht, Zt, trainX, parameterst, 0.1)

print(parameterst['W1'].sum()) # 0.579203593304
print(parameterst['b1'].sum()) # 5.05985151857e-06
print(parameterst['W2'].sum()) # 0.1888716431
print(parameterst['b2'].sum()) # 1.24683249836e-18

6. 训练

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
def train(trainX, trainY, testX, testY, parameters, epochs, learning_rate = 0.01, verbose = False):
'''
训练

Parameters
----------
Parameters
----------
trainX: np.ndarray, shape = (n, m), 训练集

trainY: np.ndarray, shape = (n, K), 训练集标记

testX: np.ndarray, shape = (n_test, m),测试集

testY: np.ndarray, shape = (n_test, K),测试集的标记

parameters: dict,参数

epochs: int, 要迭代的轮数

learning_rate: float, default 0.01,学习率

verbose: boolean, default False,是否打印损失值

Returns
----------
training_loss_list: list(float),每迭代一次之后,训练集上的损失值

testing_loss_list: list(float),每迭代一次之后,测试集上的损失值

'''
# 存储损失值
training_loss_list = []
testing_loss_list = []

for i in range(epochs):
# 这里要计算出Z和H,因为后面反向传播计算梯度的时候需要这两个矩阵
Z = linear_combination(trainX, parameters['W1'], parameters['b1'])
H = ReLU(Z)
train_O = linear_combination(H, parameters['W2'], parameters['b2'])
train_y_pred = softmax(train_O)
training_loss = cross_entropy_with_softmax(trainY, train_O)

test_O = forward(testX, parameters)
testing_loss = cross_entropy_with_softmax(testY, test_O)
if verbose == True:
print('epoch %s, training loss:%s'%(i + 1, training_loss))
print('epoch %s, testing loss:%s'%(i + 1, testing_loss))
print()

training_loss_list.append(training_loss)
testing_loss_list.append(testing_loss)

backward(trainY, train_y_pred, H, Z, trainX, parameters, learning_rate)
return training_loss_list, testing_loss_list
1
2
3
4
5
6
7
8
9
10
11
12
13
14
# 测试样例
parameterst = initialize(50, 10)
print(parameterst['W1'].sum()) # 0.583495454481
print(parameterst['b1'].sum()) # 0.0
print(parameterst['W2'].sum()) # 0.1888716431
print(parameterst['b2'].sum()) # 0.0
print()

training_loss_list, testing_loss_list = train(trainX, trainY_mat, testX, testY_mat, parameterst, 1, 0.1, False)

print(parameterst['W1'].sum()) # 0.579203593304
print(parameterst['b1'].sum()) # 5.05985151857e-06
print(parameterst['W2'].sum()) # 0.1888716431
print(parameterst['b2'].sum()) # 1.24683249836e-18

7. 绘制模型损失值变化曲线

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
def plot_loss_curve(training_loss_list, testing_loss_list):
'''
绘制损失值变化曲线

Parameters
----------
training_loss_list: list(float),每迭代一次之后,训练集上的损失值

testing_loss_list: list(float),每迭代一次之后,测试集上的损失值

'''
plt.figure(figsize = (10, 6))
plt.plot(training_loss_list, label = 'training loss')
plt.plot(testing_loss_list, label = 'testing loss')
plt.xlabel('epoch')
plt.ylabel('loss')
plt.legend()

8. 预测

模型训练完后,我们的就可以进行预测了,需要注意的是,我们的神经网络是针对每个样本,输出其分别属于$K$类的概率,我们要找最大的那个概率,对应的是哪个类。

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
def predict(X, parameters):
'''
预测,调用forward函数完成神经网络对输入X的计算,然后完成类别的划分,取每行最大的那个数的下标作为标记

Parameters
----------
X: np.ndarray, shape = (n, m), 训练集

parameters: dict,参数

Returns
----------
prediction: np.ndarray, shape = (n, 1),预测的标记

'''
# 用forward函数得到softmax激活前的值
O = forward(X, parameters)

# 计算softmax激活后的值
y_pred = softmax(O)

# 取每行最大的元素对应的下标
prediction = np.argmax(y_pred, axis = 1)

return prediction
1
2
3
4
5
6
7
8
# 测试样例
from sklearn.metrics import accuracy_score

parameterst = initialize(50, 10)
training_loss_list, testing_loss_list = train(trainX, trainY_mat, testX, testY_mat, parameterst, 1, 0.1, False)

predictiont = predict(testX, parameterst)
accuracy_score(predictiont, testY) # 0.15994436717663421

9. 训练一个三层感知机

隐藏层单元数设置为50,输出层单元数为10,我们设置学习率为0.03,迭代轮数为1000轮

1
2
3
4
5
6
7
8
9
start_time = time()

h = 50
K = 10
parameters = initialize(h, K)
training_loss_list, testing_loss_list = train(trainX, trainY_mat, testX, testY_mat, parameters, 1000, 0.03, False)

end_time = time()
print('training time: %s s'%(end_time - start_time))

计算测试集精度

1
2
prediction = predict(testX, parameters)
accuracy_score(prediction, testY)

绘制损失值变化曲线

1
plot_loss_curve(training_loss_list, testing_loss_list)

更换数据集

我们换一个数据集,使用MNIST手写数字数据集,我们使用的是kaggle提供的MNIST手写数字识别比赛的训练集。这个数据集还是手写数字的图片,只不过像素变成了 $28 \times 28$,图片的尺寸变大了,而且数据集的样本量也大了。我们取30%为测试集,70%为训练集。训练集样本数有29400个,测试集12600个。

1
2
3
4
5
6
7
8
9
10
11
12
13
import pandas as pd

data = pd.read_csv('data/kaggle_mnist/mnist_train.csv')
X = data.values[:, 1:].astype('float32')
Y = data.values[:, 0]

trainX, testX, trainY, testY = train_test_split(X, Y, test_size = 0.3, random_state = 32)

trainY_mat = np.zeros((len(trainY), 10))
trainY_mat[np.arange(0, len(trainY), 1), trainY] = 1

testY_mat = np.zeros((len(testY), 10))
testY_mat[np.arange(0, len(testY), 1), testY] = 1
1
trainX.shape, trainY.shape, trainY_mat.shape, testX.shape, testY.shape, testY_mat.shape

绘制训练集前10个图像

1
2
3
4
5
6
_, figs = plt.subplots(1, 10, figsize=(8, 4))
for f, img, lbl in zip(figs, trainX[:10], trainY[:10]):
f.imshow(img.reshape((28, 28)), cmap = 'gray')
f.set_title(lbl)
f.axes.get_xaxis().set_visible(False)
f.axes.get_yaxis().set_visible(False)

test:请你使用kaggle MNIST数据集,根据下表设定各个超参数,计算测试集上的精度,绘制损失值变化曲线,填写下表

任务流程:

  1. 对数据集进行标准化处理
  2. 设定学习率和迭代轮数进行训练
  3. 计算测试集精度
  4. 绘制曲线
双击此处填写

精度保留4位小数;训练时间单位为秒,保留两位小数。

隐藏层单元数 学习率 迭代轮数 测试集精度 训练时间(秒)
100 0.1 50 0.8071 50.47
100 0.1 100 0.8877 100.99
100 0.1 150 0.8998 147.81
100 0.1 500 0.9197 497.71
100 0.01 500 0.8125 484.04
1
2
3
stand = StandardScaler()
trainX_normalized = stand.fit_transform(trainX)
testX_normalized = stand.fit_transform(testX)
1
2
3
4
5
6
7
8
9
10
11
12
start_time = time()

h = 100
K = 10
parameters = initialize(h, K)
training_loss_list, testing_loss_list = train(trainX_normalized, trainY_mat, testX_normalized, testY_mat, parameters, 50, 0.1)
prediction = predict(testX_normalized, parameters)
print('testing accuracy:', accuracy_score(prediction, testY))
plot_loss_curve(training_loss_list, testing_loss_list)

end_time = time()
print('training time: %s s'%(end_time - start_time))
1
2
3
4
5
6
7
8
9
10
11
12
start_time = time()

h = 100
K = 10
parameters = initialize(h, K)
training_loss_list, testing_loss_list = train(trainX_normalized, trainY_mat, testX_normalized, testY_mat, parameters, 100, 0.1)
prediction = predict(testX_normalized, parameters)
print('testing accuracy:', accuracy_score(prediction, testY))
plot_loss_curve(training_loss_list, testing_loss_list)

end_time = time()
print('training time: %s s'%(end_time - start_time))
1
2
3
4
5
6
7
8
9
10
11
12
start_time = time()

h = 100
K = 10
parameters = initialize(h, K)
training_loss_list, testing_loss_list = train(trainX_normalized, trainY_mat, testX_normalized, testY_mat, parameters, 150, 0.1)
prediction = predict(testX_normalized, parameters)
print('testing accuracy:', accuracy_score(prediction, testY))
plot_loss_curve(training_loss_list, testing_loss_list)

end_time = time()
print('training time: %s s'%(end_time - start_time))
1
2
3
4
5
6
7
8
9
10
11
12
start_time = time()

h = 100
K = 10
parameters = initialize(h, K)
training_loss_list, testing_loss_list = train(trainX_normalized, trainY_mat, testX_normalized, testY_mat, parameters, 500, 0.1)
prediction = predict(testX_normalized, parameters)
print('testing accuracy:', accuracy_score(prediction, testY))
plot_loss_curve(training_loss_list, testing_loss_list)

end_time = time()
print('training time: %s s'%(end_time - start_time))
1
2
3
4
5
6
7
8
9
10
11
12
start_time = time()

h = 100
K = 10
parameters = initialize(h, K)
training_loss_list, testing_loss_list = train(trainX_normalized, trainY_mat, testX_normalized, testY_mat, parameters, 500, 0.01)
prediction = predict(testX_normalized, parameters)
print('testing accuracy:', accuracy_score(prediction, testY))
plot_loss_curve(training_loss_list, testing_loss_list)

end_time = time()
print('training time: %s s'%(end_time - start_time))