机器学习本科课程 实验2 神经网络

实验二 神经网络

• 1 sklearn多层感知机——手写数字分类任务
• 2 手动实现神经网络：线性回归
• 3 手动实现神经网络：对数几率回归
• 4 手动实现三层感知机——多分类

第一题sklearn多层感知机——手写数字分类任务

实验内容：

1. 使用sklearn.neural_network.MLPClassifier完成手写数字分类任务
2. 绘制学习率为3，1，0.1，0.01训练集损失函数的变化曲线

1. 读取数据集

from sklearn.datasets import load_digits
load_digits().keys()
print(load_digits()['DESCR'])

一共1797个样本，每个样本都是$8\times 8$的矩阵，因为是灰度的图像，所以只有一个通道
images对应的是原始的图像，data对应的是 $8\times 8$ reshape成 $1\times 64$ 的数据，target是标记，表示这张图片里面的数字是几

打印第一个样本

load_digits()['images'][0]

对数据集的前十张图片可视化

import matplotlib.pyplot as plt
%matplotlib inline
_, figs = plt.subplots(1, 10, figsize=(10, 4))
for f, img, lbl in zip(figs, load_digits()['images'][:10], load_digits()['target'][:10]):
    f.imshow(img, cmap = 'gray')
    f.set_title(lbl)
    f.axes.get_xaxis().set_visible(False)
    f.axes.get_yaxis().set_visible(False)

2. 划分数据集

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

3. 数据预处理

from sklearn.preprocessing import StandardScaler

神经网络的训练方法一般是基于梯度的优化算法，如梯度下降，为了让这类算法能更好的优化神经网络，我们往往需要对数据集进行归一化，这里我们选择对数据进行标准化

$$X^{\prime }=\frac{X-\bar{X}}{\mathrm{std}(X)}$$

其中，$\bar{X}$是均值，$\mathrm{std}$是标准差。减去均值可以让数据以0为中心，除以标准差可以让数据缩放到一个较小的范围内。这样可以使得梯度的下降方向更多样，同时缩小梯度的数量级，让学习变得稳定。
首先需要对训练集进行标准化，针对每个特征求出其均值和标准差，然后用训练集的每个样本减去均值除以标准差，就得到了新的训练集。然后用测试集的每个样本，减去训练集的均值，除以训练集的标准差，完成对测试集的标准化。

# 初始化一个标准化器的实例
standard = StandardScaler()

# 对训练集进行标准化，它会计算训练集的均值和标准差保存起来
trainX = standard.fit_transform(trainX)

# 使用标准化器在训练集上的均值和标准差，对测试集进行归一化
testX = standard.transform(testX)

4. 引入模型

from sklearn.neural_network import MLPClassifier
model = MLPClassifier(solver = 'sgd', learning_rate = 'constant', momentum = 0, learning_rate_init = 0.1, max_iter = 500)
model.fit(trainX, trainY)
prediction = model.predict(testX)

5. 预测与评估

from sklearn.metrics import accuracy_score
accuracy_score(prediction, testY)

6. 绘制训练集损失函数值的变化曲线

plt.figure(figsize = (10, 6))
plt.plot(model.loss_curve_)
plt.xlabel('epoch')
plt.ylabel('loss')

test：请你在一张图内，绘制出学习率为3，学习率为1，学习率为0.1，学习率为0.01，四个模型的损失函数变化曲线，最大迭代轮数为250轮。

提示：分别训练4个模型，然后在一张图中分别绘制4个模型的loss_curve_即可。

# your code
model_lr3 = MLPClassifier(solver = 'sgd', learning_rate = 'constant', momentum = 0, learning_rate_init = 3, max_iter = 250)
model_lr3.fit(trainX, trainY)
prediction_lr3 = model.predict(testX)
acc_lr3 = accuracy_score(prediction, testY)

model_lr1 = MLPClassifier(solver = 'sgd', learning_rate = 'constant', momentum = 0, learning_rate_init = 1, max_iter = 250)
model_lr1.fit(trainX, trainY)
prediction_lr1 = model.predict(testX)
acc_lr1 = accuracy_score(prediction, testY)

model_lr01 = MLPClassifier(solver = 'sgd', learning_rate = 'constant', momentum = 0, learning_rate_init = 0.1, max_iter = 250)
model_lr01.fit(trainX, trainY)
prediction_lr01 = model.predict(testX)
acc_lr01 = accuracy_score(prediction, testY)

model_lr001 = MLPClassifier(solver = 'sgd', learning_rate = 'constant', momentum = 0, learning_rate_init = 0.01, max_iter = 250)
model_lr001.fit(trainX, trainY)
prediction_lr001 = model.predict(testX)
acc_lr001 = accuracy_score(prediction, testY)

fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(nrows=2, ncols=2)
ax1.plot(model_lr3.loss_curve_)
ax1.set_xlabel('epoch')
ax1.set_ylabel('loss')
ax1.set_title('Learning Rate: 3')

ax2.plot(model_lr1.loss_curve_)
ax2.set_xlabel('epoch')
ax2.set_ylabel('loss')
ax2.set_title('Learning Rate: 1')

ax3.plot(model_lr01.loss_curve_)
ax3.set_xlabel('epoch')
ax3.set_ylabel('loss')
ax3.set_title('Learning Rate: 0.1')

ax4.plot(model_lr3.loss_curve_)
ax4.set_xlabel('epoch')
ax4.set_ylabel('loss')
ax4.set_title('Learning Rate: 0.01')

plt.tight_layout()
plt.show()

第二题手动实现神经网络：线性回归

实验内容：

使用Kaggle房价预测数据集，以LotArea, BsmtUnfSF, GarageArea三种特征作为模型的输入，SalePrice作为模型的输出

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}$。
输出层没有激活函数，所以上面的神经网络的前向传播过程写为：

$$Z=XW+b$$

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

$$\mathrm{loss}(y,\hat{y})=\frac{1}{n}\sum _{i=1}^n(y_i-\widehat{y_i}{)}^2$$

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

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

梯度下降

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

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

$$W:=W-\alpha \frac{\partial \mathrm{loss}}{\partial W}$$

$$b:=b-\alpha \frac{\partial \mathrm{loss}}{\partial b}$$

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

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

损失函数MSE记为：

$$\mathrm{loss}(y,Z)=\frac{1}{n}\sum _{i=1}^n(y_i-Z_i{)}^2$$

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

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

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

$$\begin{array}{rl}\frac{\partial \mathrm{loss}}{\partial W}& =\frac{\partial \mathrm{loss}}{\partial Z}\frac{\partial Z}{\partial W}\\ & =-\frac{2}{n}{X}^{\mathrm{T}}(y-Z)\\ & =\frac{2}{n}{X}^{\mathrm{T}}(Z-y)\end{array}$$

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

$$\begin{array}{rl}\frac{\partial \mathrm{loss}}{\partial b}& =\frac{\partial \mathrm{loss}}{\partial Z}\frac{\partial Z}{\partial b}\\ & =-\frac{2}{n}\sum _{i=1}^n(y_i-Z_i)\\ & =\frac{2}{n}\sum _{i=1}^n(Z_i-y_i)\end{array}$$

因为参数$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步，在反复迭代的过程中可以看到损失值不断减小的现象，如果没有下降说明出了问题

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

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline

# 读取数据
data = pd.read_csv('train.csv')

# 丢弃有缺失值的特征（列）
data.dropna(axis = 1, inplace = True)

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

2. 数据预处理

from sklearn.model_selection import train_test_split
# your code trainX, testX, trainY, testY = ...
num_of_samples = data.shape[0]
X = data.drop(columns=['SalePrice'])
y = data['SalePrice']
trainX, testX, trainY, testY = train_test_split(X, y, test_size=0.4, random_state=42)

3. 参数初始化

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. 前向传播

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
    # 可以用：np.dot()
    Z = np.dot(X, W) + b
    return Z
# 测试样例
Wt, bt = initialize(trainX.shape[1])
tmp = forward(trainX, Wt, bt)
print(tmp.mean()) 

5. 损失函数

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

$$\mathrm{loss}(y,Z)=\frac{1}{n}\sum _{i=1}^n(y_i-Z_i{)}^2$$

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 = np.sum((y_true - y_pred)**2)/y_true.shape[0]
    return loss
# 测试样例
Wt, bt = initialize(trainX.shape[1])
tmp = mse(trainY, forward(trainX, Wt, bt))
print(tmp) 

6. 反向传播

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

$$\frac{\partial \mathrm{loss}}{\partial W}=\frac{2}{n}{X}^{\mathrm{T}}(Z-y)$$

$$\frac{\partial \mathrm{loss}}{\partial b}=\frac{2}{n}\sum _{i=1}^n(Z_i-y_i)$$

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 = X.shape[0]
    X_T = np.transpose(X)
    # 计算W的梯度
    # YOUR CODE HERE   dW =...  
    dW = 2*np.dot(X_T, (Z - y_true))/n
    # 计算b的梯度
    # YOUR CODE HERE db = ... 
    db = 2*np.sum(Z - y_true)/n
    return dW, db

# 测试样例
Wt, bt = initialize(trainX.shape[1])
Zt = forward(trainX, Wt, bt)
dWt, dbt = compute_gradient(trainX, Zt, trainY)
print(dWt.shape) 
print(dWt.mean()) 
print(dbt) 

7. 梯度下降

这部分需要实现梯度下降的函数
$$W:=W-\alpha \frac{\partial \mathrm{loss}}{\partial W}$$

$$b:=b-\alpha \frac{\partial \mathrm{loss}}{\partial b}$$

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
    # your code
    W = W - learning_rate*dW
    # 更新b
    # your code
    b = b - learning_rate*db 
    # 经过测试, 实际上W和b并没有被修改, 需要提供返回值
    return W, b

# 测试样例
Wt, bt = initialize(trainX.shape[1])
print(Wt.mean()) 
print(bt.mean()) 
Zt = forward(trainX, Wt, bt)
dWt, dbt = compute_gradient(trainX, Zt, trainY)
Wt, bt = update(dWt, dbt, Wt, bt, 0.01)

print(Wt.shape) 
print(Wt.mean()) 
print(bt) 

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
    W, b = update(dW, db, W, b, learning_rate)
    return W, b

# 测试样例
Wt, bt = initialize(trainX.shape[1])
print(Wt.mean()) 
print(bt.mean()) 

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

print(Wt.shape) 
print(Wt.mean()) 
print(bt) 

8. 训练

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()
        
        # 反向传播，参数更新
        W, b = backward(trainX, Z, trainY, W, b, learning_rate)
        # print("Epoch:", epoch, "W,b", W, b)
    return training_loss_list, testing_loss_list, W, b

# 测试样例
Wt, bt = initialize(trainX.shape[1])
print(Wt.mean())          
print(bt.mean())          

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

print(training_loss_list) 
print(testing_loss_list)  
print(Wt.mean())          
print(bt)          

9. 检查

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

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

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

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

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()

# 特征数m
# m=...
m = 3
# 参数初始化
# W, b =...
W, b = initialize(m)
# 训练20轮，学习率为0.01
# training_loss_list, testing_loss_list = ...
training_loss_list, testing_loss_list, W, b= train(trainX, trainY, testX, testY, W, b, 20, learning_rate = 0.01, verbose = False)

plot_loss_curve(training_loss_list, testing_loss_list)

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

# 测试样例
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) # -364308.555764
print()

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

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

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

10. 标准化处理

from sklearn.preprocessing import StandardScaler
standard = StandardScaler()
trainX_norm = standard.fit_transform(trainX)
testX_norm = standard.fit_transform(testX)
W, b = initialize(m)
training_loss_list, testing_loss_list, W, b= train(trainX_norm, trainY, testX_norm, testY, W, b, 40, learning_rate = 0.1, verbose = False)
plot_loss_curve(training_loss_list, testing_loss_list)

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

实验内容：

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

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

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

激活函数

$$\mathrm{sigmoid}(x)=\frac{1}{1+e^{-x}}$$

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

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

前向传播

$$\begin{gathered} Z=XW+b \\ \hat{y}=\sigma (Z) \end{gathered}$$

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

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

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

$$\hat{y}=\frac{1}{1+e^{-XW+b}}$$

损失函数

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

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

$$\mathrm{loss}(y,\hat{y})=-y\log \hat{y}-(1-y)\log (1-\hat{y})$$

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

$$\mathrm{loss}(y,\hat{y})=-\frac{1}{n}\sum _{i=1}^n[y_i\log \widehat{y_i}+(1-y_i)\log (1-\widehat{y_i})]$$

注意，这里我们的提到的$\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。

参数更新

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

$$\begin{gathered} W:=W-\alpha \frac{\partial \mathrm{loss}}{\partial W} \\ b:=b-\alpha \frac{\partial \mathrm{loss}}{\partial b} \end{gathered}$$

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

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

反向传播

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

$$\begin{array}{rl}\frac{\partial \mathrm{loss}}{\partial W}& =\frac{\partial \mathrm{loss}}{\partial \hat{y}}\frac{\partial \hat{y}}{\partial Z}\frac{\partial Z}{\partial W}\end{array}$$

这里我们一项一项求，先求第一项：

$$\begin{array}{r}\frac{\partial \mathrm{loss}}{\partial \hat{y}}=-\frac{1}{n}\sum _{i=1}^n[\frac{y}{\hat{y}}-\frac{1-y}{1-\hat{y}}]\end{array}$$

第二项：

$$\begin{array}{rl}\frac{\partial \hat{y}}{\partial Z}& =\frac{\partial (\frac{1}{1+e^{-Z}})}{\partial Z}\\ & =\frac{e^{-Z}}{(1+e^{-Z}{)}^2}\\ & =\frac{e^{-Z}}{(1+e^{-Z})}\frac{1}{(1+e^{-Z})}\\ & =\frac{e^{-Z}}{(1+e^{-Z})}(1-\frac{e^{-Z}}{(1+e^{-Z})})\\ & =\sigma (Z)(1-\sigma (Z))\end{array}$$

第三项：

$$\frac{\partial Z}{\partial W}=X^{\mathrm{T}}$$

综上：

$$\begin{array}{rl}\frac{\partial \mathrm{loss}}{\partial W}& =\frac{\partial \mathrm{loss}}{\partial \hat{y}}\frac{\partial \hat{y}}{\partial Z}\frac{\partial Z}{\partial W}\\ & =-\frac{1}{n}\sum _{i=1}^n[\frac{y_i}{\widehat{y_i}}-\frac{1-y_i}{1-\widehat{y_i}}][\sigma (Z_i)(1-\sigma (Z_i))]{X_i}^{\mathrm{T}}\\ & =-\frac{1}{n}\sum _{i=1}^n[\frac{y_i}{\widehat{y_i}}-\frac{1-y_i}{1-\widehat{y_i}}][\widehat{y_i}(1-\widehat{y_i})]{X_i}^{\mathrm{T}}\\ & =-\frac{1}{n}\sum _{i=1}^n[y_i(1-\widehat{y_i})-\widehat{y_i}(1-y_i)]{X_i}^{\mathrm{T}}\\ & =-\frac{1}{n}\sum _{i=1}^n(y_i-y_i\widehat{y_i}-\widehat{y_i}+y_i\widehat{y_i}){X_i}^{\mathrm{T}}\\ & =-\frac{1}{n}\sum _{i=1}^n(y_i-\widehat{y_i}){X_i}^{\mathrm{T}}\\ & =\frac{1}{n}[X^{\mathrm{T}}(\hat{y}-y)]\end{array}$$

同理，求$\mathrm{loss}$对$b$的偏导数：

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

$$\begin{array}{rl}\frac{\partial \mathrm{loss}}{\partial b}& =\frac{\partial \mathrm{loss}}{\partial \hat{y}}\frac{\partial \hat{y}}{\partial Z}\frac{\partial Z}{\partial b}\\ & =-\frac{1}{n}\sum _{i=1}^n[\frac{y_i}{\widehat{y_i}}-\frac{1-y_i}{1-\widehat{y_i}}][\sigma (Z_i)(1-\sigma (Z_i))]\\ & =-\frac{1}{n}\sum _{i=1}^n[\frac{y_i}{\widehat{y_i}}-\frac{1-y_i}{1-\widehat{y_i}}][\widehat{y_i}(1-\widehat{y_i})]\\ & =-\frac{1}{n}\sum _{i=1}^n[y_i(1-\widehat{y_i})-\widehat{y_i}(1-y_i)]\\ & =-\frac{1}{n}\sum _{i=1}^n(y_i-y_i\widehat{y_i}-\widehat{y_i}+y_i\widehat{y_i})\\ & =\frac{1}{n}\sum _{i=1}^n(\widehat{y_i}-y_i)\end{array}$$

1. 导入数据集

import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from matplotlib.colors import ListedColormap

from sklearn.datasets import make_moons
X, y = make_moons(n_samples = 2000, noise = 0.3, random_state=0)
plt.figure(figsize = (6, 6))
cm_bright = ListedColormap(['#FF0000', '#0000FF'])
plt.scatter(X[:, 0], X[:, 1], c = y, cmap = cm_bright, edgecolors = 'k')

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

2. 数据预处理

from sklearn.preprocessing import StandardScaler
# 初始化一个标准化器的实例
standard = StandardScaler()

# 对训练集进行标准化，它会计算训练集的均值和标准差保存起来
trainX = standard.fit_transform(trainX)

# 使用标准化器在训练集上的均值和标准差，对测试集进行归一化
testX = standard.transform(testX)

3. 定义神经网络

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

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
    return Z
def my_sigmoid(x):
    '''
    simgoid 1 / (1 + exp(-x))
    
    Parameters
    ----------
    X: np.ndarray, 待激活的值
    
    '''
    # YOUR CODE HERE
    activations = 1 / (1 + np.exp(-x))
    return activations

在实现$\mathrm{sigmoid}$的时候，可能会遇到上溢(overflow)的问题，可以看到$\mathrm{sigmoid}$中有一个指数运算
$$\mathrm{sigmoid}(x)=\frac{1}{1+e^{-x}}$$
当$x$很大的时候，我们使用numpy.exp(x)会直接溢出

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

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

解决这种问题的方法有很多，比如，我们可以将$\mathrm{sigmoid}$进行变换：

$$\begin{array}{rl}\mathrm{sigmoid}(x)& =\frac{1}{1+e^{-x}}\\ & =\frac{e^x}{1+e^x}\\ & =\frac{1}{2}+\frac{1}{2}\tanh (\frac{x}{2})\end{array}$$

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

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

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

from scipy.special import expit
def sigmoid(X):
    return expit(X)

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

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
    # YOUR CODE HERE
    Z = linear_combination(X, W, b)
    # 求激活后的预测值
    # YOUR CODE HERE
    y_pred = sigmoid(Z)
    return y_pred

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

$$\mathrm{loss}(y,\hat{y})=-\frac{1}{n}[y\log \hat{y}+(1-y)\log (1-\hat{y})]$$

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

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

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)
    
    # 求解对数损失
    # YOUR CODE HERE
    loss = -(y_true*np.log(y_hat) + (1-y_true)*np.log(1-y_hat))
    return loss

我们接下来要完成损失函数对参数的偏导数的计算
$ \dfrac{\partial{loss}}{\partial{W}} = \dfrac{1}{n} X^T (\hat{y}-y)$

$ \dfrac{\partial{loss}}{\partial{b}} = \dfrac{1}{n} \sum{(\hat{y}-y)}$

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的偏导数
    
    '''
    n = X.shape[0]
    # 求损失函数对参数W的偏导数的导数值
    # YOUR CODE HERE
    dW = (1/n)*np.dot(X.T, (y_pred - y_true))
    # 求损失函数对参数b的偏导数的导数值
    # YOUR CODE HERE
    db = (1/n)*np.sum(y_pred-y_true)
    return dW, db

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 = W - learning_rate*dW
    # 对参数b进行更新
    b = b - learning_rate*db
    return W, b
    

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(trainY, forward(trainX, Wt, bt), trainX)
    
    # 梯度下降
    W, b = update(W, b, dW, db, learning_rate)
    return W, b

4. 训练函数的编写

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):
        
        # 计算训练集前向传播得到的预测值
        y_pred = forward(trainX, W, b)
        
        # 计算当前训练集的损失值
        training_loss = logloss(trainY, y_pred)
        
        # 计算测试集前向传播得到的预测值
        test_pred = forward(testX, W, b)
        testing_loss = logloss(testY, test_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.mean())
        testing_loss_list.append(testing_loss.mean())
        
        # 反向传播更新参数
        W, b = backward(trainY, y_pred, trainX, W, b, learning_rate)
        

    
    return training_loss_list, testing_loss_list, W, b

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

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. 预测

def predict(X, parameters):
    '''
    预测，调用forward函数完成神经网络对输入X的计算，然后完成类别的划分，取每行最大的那个数的下标作为标记
    
    Parameters
    ----------
    X: np.ndarray, shape = (n, m), 训练集
    
    parameters: dict，参数
    
    Returns
    ----------
    prediction: np.ndarray, shape = (n, 1)，预测的标记
    
    '''
    # 用forward函数得到softmax激活前的值
    # YOUR CODE HERE
    Z = forward(X, parameters)
    
    # 计算softmax激活后的值
    # YOUR CODE HERE
    H = softmax(Z)
    
    # 取每行最大的元素对应的下标
    # YOUR CODE HERE
    prediction = np.argmax(H, axis=1)
    
    return prediction

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, True)

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

prediction = predict(testX, parameters)
accuracy_score(prediction, testY)
plot_loss_curve(training_loss_list, testing_loss_list)

更换数据集

我们换一个数据集，使用MNIST手写数字数据集。

MNIST是最有名的手写数字数据集之一，主页：http://yann.lecun.com/exdb/mnist/

MNIST手写数字数据集有60000个样本组成的训练集，10000个样本组成的测试集，是NIST的子集。数字的尺寸都是归一化后的，且都在图像的中央。可以从上方的主页下载。

我们使用的数据集是kaggle手写数字识别比赛中的训练集。数据集一共60000行，785列，其中第1列是标记，第2列到第785列是图像从左上角到右下角的像素值。图像大小为28×28像素，单通道的灰度图像。

我们取30%为测试集，70%为训练集。

import pandas as pd

data = pd.read_csv('./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

_, 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. 绘制曲线

# YOUR CODE HERE
# 标准化操作
# 初始化一个标准化器的实例
standard = StandardScaler()

# 对训练集进行标准化，它会计算训练集的均值和标准差保存起来
trainX = standard.fit_transform(trainX)

# 使用标准化器在训练集上的均值和标准差，对测试集进行归一化
testX = standard.transform(testX)

# your code
parameters1 = initialize(100, 10)
print(parameters1['W1'].sum())  
print(parameters1['b1'].sum())  
print(parameters1['W2'].sum())  
print(parameters1['b2'].sum())  
print()

start_time = time()
training_loss_list, testing_loss_list = train(trainX, trainY_mat, testX, testY_mat, parameters1, 50, 0.1, True)
end_time = time()
print('training time: %s s'%(end_time - start_time))
print()

print(parameters1['W1'].sum())  
print(parameters1['b1'].sum())  
print(parameters1['W2'].sum())  
print(parameters1['b2'].sum())  

prediction = predict(testX, parameters1)
print("accuracy:", accuracy_score(prediction, testY))
plot_loss_curve(training_loss_list, testing_loss_list)

# your code
parameters1 = initialize(100, 10)
print(parameters1['W1'].sum())  
print(parameters1['b1'].sum())  
print(parameters1['W2'].sum())  
print(parameters1['b2'].sum())  
print()

start_time = time()
training_loss_list, testing_loss_list = train(trainX, trainY_mat, testX, testY_mat, parameters1, 100, 0.1, True)
end_time = time()
print('training time: %s s'%(end_time - start_time))
print()

print(parameters1['W1'].sum())  
print(parameters1['b1'].sum())  
print(parameters1['W2'].sum())  
print(parameters1['b2'].sum())  

prediction = predict(testX, parameters1)
print("accuracy:", accuracy_score(prediction, testY))
plot_loss_curve(training_loss_list, testing_loss_list)

# your code
parameters1 = initialize(100, 10)
print(parameters1['W1'].sum())  
print(parameters1['b1'].sum())  
print(parameters1['W2'].sum())  
print(parameters1['b2'].sum())  
print()

start_time = time()
training_loss_list, testing_loss_list = train(trainX, trainY_mat, testX, testY_mat, parameters1, 150, 0.1, True)
end_time = time()
print('training time: %s s'%(end_time - start_time))
print()

print(parameters1['W1'].sum())  
print(parameters1['b1'].sum())  
print(parameters1['W2'].sum())  
print(parameters1['b2'].sum())  

prediction = predict(testX, parameters1)
print("accuracy:", accuracy_score(prediction, testY))
plot_loss_curve(training_loss_list, testing_loss_list)

# your code
parameters1 = initialize(100, 10)
print(parameters1['W1'].sum())  
print(parameters1['b1'].sum())  
print(parameters1['W2'].sum())  
print(parameters1['b2'].sum())  
print()

start_time = time()
training_loss_list, testing_loss_list = train(trainX, trainY_mat, testX, testY_mat, parameters1, 500, 0.1, True)
end_time = time()
print('training time: %s s'%(end_time - start_time))
print()

print(parameters1['W1'].sum())  
print(parameters1['b1'].sum())  
print(parameters1['W2'].sum())  
print(parameters1['b2'].sum())  

prediction = predict(testX, parameters1)
print("accuracy:", accuracy_score(prediction, testY))
plot_loss_curve(training_loss_list, testing_loss_list)

parameters1 = initialize(100, 10)
print(parameters1['W1'].sum())  
print(parameters1['b1'].sum())  
print(parameters1['W2'].sum())  
print(parameters1['b2'].sum())  
print()

start_time = time()
training_loss_list, testing_loss_list = train(trainX, trainY_mat, testX, testY_mat, parameters1, 500, 0.01, True)
end_time = time()
print('training time: %s s'%(end_time - start_time))
print()

print(parameters1['W1'].sum())  
print(parameters1['b1'].sum())  
print(parameters1['W2'].sum())  
print(parameters1['b2'].sum())  

prediction = predict(testX, parameters1)
print("accuracy:", accuracy_score(prediction, testY))
plot_loss_curve(training_loss_list, testing_loss_list)