神经网络基础

最近给本科生当机器学习课程的助教，给他们出的作业题需要看这些图，懒得放本地了，直接放博客里。发现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}$。
输出层没有激活函数，所以上面的神经网络的前向传播过程写为：

$$Z = XW + b$$

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

$$\mathrm{loss}(y, \hat{y}) = \frac{1}{n} \sum^n_{i=1}(y_i - \hat{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^n_{i = 1} (y_i - Z_i)^2$$

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

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

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

$$\begin{aligned} \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{aligned}$$

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

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

因为参数$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列作为特征

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%做训练集

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个特征

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

3. 参数初始化

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

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$ 的计算。

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

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

5. 损失函数

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

$$\mathrm{loss}(y, Z) = \frac{1}{n} \sum^n_{i = 1} (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 = ((y_true - y_pred) ** 2).sum() / len(y_true)
    
    return loss

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

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^n_{i=1}(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 = 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

# 测试样例
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. 梯度下降

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

$$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
    W -= learning_rate * dW
    
    # 更新b
    b -= learning_rate * db

# 测试样例
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组装在一起。

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)

# 测试样例
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. 训练

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

# 测试样例
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共同下降，这就是我们想要的结果，说明模型正在很好的学习。

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	绘制损失函数的变化曲线

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

# 特征数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)

绘制损失值的变化曲线

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) # -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. 标准化处理

标准化处理和第一题一样

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

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

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)

打印损失值变化曲线

plot_loss_curve(training_loss_list, testing_loss_list)

计算测试集上的MSE

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}$。

激活函数

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

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

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

前向传播

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

其中，$O \in \mathbb{R}^{n}$为输出层的结果，$\sigma$为$\rm 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^n_{i = 1}[y_i \log{\hat{y_i}} + (1 - y_i) \log{(1 - \hat{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。

参数更新

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

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

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

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

反向传播

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

$$\begin{aligned} \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{aligned}$$

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

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

第二项：

$$\begin{aligned} \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{aligned}$$

第三项：

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

综上：

$$\begin{aligned} \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^n_{i = 1} [\frac{y_i}{\hat{y_i}} - \frac{1 - y_i}{1 - \hat{y_i}}] [\sigma(Z_i)(1 - \sigma(Z_i))] {X_i}^{\mathrm{T}}\\ &= - \frac{1}{n} \sum^n_{i = 1} [\frac{y_i}{\hat{y_i}} - \frac{1 - y_i}{1 - \hat{y_i}}] [\hat{y_i}(1 - \hat{y_i})] {X_i}^{\mathrm{T}}\\ &= - \frac{1}{n} \sum^n_{i = 1} [y_i(1 - \hat{y_i}) - \hat{y_i}(1 - y_i)] {X_i}^{\mathrm{T}}\\ &= - \frac{1}{n} \sum^n_{i = 1} (y_i - y_i \hat{y_i} - \hat{y_i} + y_i \hat{y_i}) {X_i}^{\mathrm{T}}\\ &= - \frac{1}{n} \sum^n_{i = 1} (y_i - \hat{y_i}) {X_i}^{\mathrm{T}}\\ &= \frac{1}{n} [X^{\mathrm{T}}(\hat{y} - y)] \end{aligned}$$

同理，求$\rm loss$对$b$的偏导数：

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

$$\begin{aligned} \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^n_{i = 1} [\frac{y_i}{\hat{y_i}} - \frac{1 - y_i}{1 - \hat{y_i}}] [\sigma(Z_i)(1 - \sigma(Z_i))]\\ &= - \frac{1}{n} \sum^n_{i = 1} [\frac{y_i}{\hat{y_i}} - \frac{1 - y_i}{1 - \hat{y_i}}] [\hat{y_i}(1 - \hat{y_i})]\\ &= - \frac{1}{n} \sum^n_{i = 1} [y_i(1 - \hat{y_i}) - \hat{y_i}(1 - y_i)]\\ &= - \frac{1}{n} \sum^n_{i = 1} (y_i - y_i \hat{y_i} - \hat{y_i} + y_i \hat{y_i})\\ &= \frac{1}{n} \sum^n_{i = 1} (\hat{y_i} - y_i)\\ \end{aligned}$$

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

1. 导入数据集

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

我们生成半月形数据

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

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

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

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

2. 数据预处理

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

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。为什么要对神经网络的参数进行随机初始化，感兴趣的同学可以去查相关的资料。

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

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

3.2 前向传播

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

首先计算$Z = XW + 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                          # YOUR CODE HERE
    
    return Z

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

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

def my_sigmoid(x):
    '''
    simgoid 1 / (1 + exp(-x))
    
    Parameters
    ----------
    X: np.ndarray, 待激活的值
    
    '''
    # YOUR CODE HERE
    activations = 1 / (1 + np.exp(-x))
    
    return activations

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

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

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

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

np.exp(1e56)

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

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

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

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

其中，$\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$的计算。我们现在做的都是神经网络底层相关的运算，很容易出现数值不稳定性相关的问题，最好的办法就是使用别人已经实现好的函数，这样就能减少我们很多的工作量，同时又快速地完成任务。

from scipy.special import expit

def sigmoid(X):
    return expit(X)

# 测试样例
sigmoid(np.array([-1e56]))

接下来完成整个前向传播的函数，也就是 $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
    Z = linear_combination(X, W, b) # YOUR CODE HERE
    
    # 求激活后的预测值
    y_pred = sigmoid(Z)             # YOUR CODE HERE
    
    return y_pred

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

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

$$\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)
    
    # 求解对数损失
    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

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

3.3 反向传播

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

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

# 测试样例
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 参数更新

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

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

# 测试样例
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]

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

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)

# 测试样例
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. 训练函数的编写

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

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

# 测试样例
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. 绘制模型损失值变化曲线

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，这样输出的就是类标了。

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

# 测试样例
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轮

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

计算测试集精度

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

绘制损失值变化曲线

plot_loss_curve(training_loss_list, testing_loss_list)

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

双击此处填写

学习率	迭代轮数	测试集精度
0.0001	200	0.3325
0.1	1000	0.84

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

# 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

$$\mathrm{ReLU}(x) = \max (0, x)$$

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

$$Z = XW_1 + b_1\\ H_1 = \mathrm{ReLU}(Z)\\ O = H_1 W_2 + b_2$$

其中，$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函数激活，得到神经网络计算出的各类的概率值：

$$\begin{aligned} \hat{y_i} & = \mathrm{softmax}(O_i)\\ & = \frac{\exp{(O_i)}}{\sum^{K}_{k=1} \exp{(O_k)}} \end{aligned}$$

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

损失函数

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

$$\mathrm{cross\_entropy}(y, \hat{y}) = -\sum^{K}_{k=1}y_k \log{(\hat{y_k})}$$

这样，$n$个样本的平均损失为：

$$\mathrm{loss} = - \frac{1}{n} \sum_n \sum^{K}_{k=1} y_k \log{(\hat{y_k})}$$

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

反向传播

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

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

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

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

$$\frac{\partial \mathrm{loss}}{\partial O_i} = \frac{\partial \mathrm{loss}}{\partial \hat{y}} \frac{\partial \hat{y}}{\partial O_i}$$

首先求解第一项：

$$\frac{\partial \mathrm{loss}}{\partial \hat{y}} = - \frac{1}{n} \sum_n \sum^{K}_{k=1} y_k \frac{1}{\hat{y_k}}$$

然后求解第二项，因为$\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$时：

$$\begin{aligned} \frac{\partial \hat{y_i}}{\partial O_i} & = \frac{\partial \hat{y_i}}{\partial O_i}\\ & = \frac{\exp{(O_i)} (\sum^{K}_{k=1} \exp{(O_k)}) - (\exp{(O_i)})^2 }{(\sum^{K}_{k=1} \exp{(O_k)})^2}\\ & = \frac{\exp{(O_i)}}{\sum^{K}_{k=1} \exp{(O_k)}} \frac{\sum^{K}_{k=1} \exp{(O_k)} - \exp{(O_i)}}{\sum^{K}_{k=1} \exp{(O_k)}}\\ & = \hat{y_i} ( 1 - \hat{y_i} ) \end{aligned}$$

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

$$\begin{aligned} \frac{\partial \hat{y_j}}{\partial O_i} & = \frac{- \exp{(O_j)} \exp{(O_i)}}{(\sum^{K}_{k=1} \exp{(O_k)})^2}\\ & = - \hat{y_j} \hat{y_i} \end{aligned}$$

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

$$\begin{aligned} \frac{\partial \mathrm{loss}}{\partial O_i} & = \frac{\partial \mathrm{loss}}{\partial \hat{y}} \frac{\partial \hat{y}}{\partial O_i}\\ & = (- \frac{1}{n} \sum_n \sum^{K}_{k=1} y_k \frac{1}{\hat{y_k}}) \frac{\partial \hat{y}}{\partial O_i}\\ & = - \frac{1}{n} \sum_n (y_i \frac{1}{\hat{y_i}} \hat{y_i} ( 1 - \hat{y_i} ) + \sum^K_{k \not= i} y_k \frac{1}{\hat{y_k}}( - \hat{y_k} \hat{y_i}))\\ & = - \frac{1}{n} \sum_n ( y_i - y_i \hat{y_i} - \sum^K_{k \not= i} y_k \hat{y_i})\\ & = - \frac{1}{n} \sum_n ( y_i - \hat{y_i} \sum^K_{k = 1} y_k ) \end{aligned}$$

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

$$\begin{aligned} \frac{\partial \mathrm{loss}}{\partial O_i} &= - \frac{1}{n} \sum_n ( y_i - \hat{y_i})\\ & = \frac{1}{n} \sum_n (\hat{y_i} - y_i) \end{aligned}$$

将其写成矩阵表达式：

$$\begin{aligned} \frac{\partial \mathrm{loss}}{\partial O} &= \frac{1}{n} (\hat{y} - y) \end{aligned}$$

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

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

$$\begin{aligned} \frac{\partial loss}{\partial W_2} & = \frac{\partial \mathrm{loss}}{\partial \hat{y}} \frac{\partial \hat{y}}{\partial O} \frac{\partial O}{\partial W_2}\\ & = \frac{\partial loss}{\partial O} \frac{\partial O}{\partial W_2}\\ & = \frac{1}{n} (\hat{y} - y) \frac{\partial O}{\partial W_2}\\ & = \frac{1}{n} [{H_1}^\mathrm{T} (\hat{y} - y)] \end{aligned}$$ $$\begin{aligned} \frac{\partial loss}{\partial b_2} & = \frac{\partial \mathrm{loss}}{\partial \hat{y}} \frac{\partial \hat{y}}{\partial O} \frac{\partial O}{\partial b_2}\\ & = \frac{\partial loss}{\partial O} \frac{\partial O}{\partial b_2}\\ & = \frac{1}{n} (\hat{y} - y) \frac{\partial O}{\partial b_2}\\ & = \frac{1}{n} \sum^n_{i=1} (\hat{y_i} - y_i) \end{aligned}$$

其中，$\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$的偏导数：

$$\begin{aligned} \frac{\partial loss}{\partial W_1} & = \frac{\partial \mathrm{loss}}{\partial \hat{y}} \frac{\partial \hat{y}}{\partial O} \frac{\partial O}{\partial H_1} \frac{\partial H_1}{\partial Z} \frac{\partial Z}{\partial W_1}\\ & = \frac{\partial loss}{\partial O} \frac{\partial O}{\partial H_1} \frac{\partial H_1}{\partial Z} \frac{\partial Z}{\partial W_1}\\ & = \frac{1}{n} {X}^\mathrm{T} [(\hat{y} - y) {W_2}^\mathrm{T} \frac{\partial H_1}{\partial Z}]\\ \end{aligned}$$

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

$$\frac{\partial \mathrm{ReLU(x)}}{\partial x} = \begin{cases} 0 & \text{if } x < 0\\ 1 & \text{if } x \geq 0 \end{cases}$$

所以上式为：

$$\frac{\partial loss}{\partial {W_1}_{ij}} = \begin{cases} 0 & \text{if } {Z}_{ij} < 0\\ \frac{1}{n} {X}^\mathrm{T} (\hat{y} - y) {W_2}^\mathrm{T} & \text{if } {Z}_{ij} \geq 0 \end{cases}$$

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

$$\begin{aligned} \frac{\partial loss}{\partial b_1} & = \frac{\partial \mathrm{loss}}{\partial \hat{y}} \frac{\partial \hat{y}}{\partial O} \frac{\partial O}{\partial H_1} \frac{\partial H_1}{\partial Z} \frac{\partial Z}{\partial b_1}\\ & = \frac{\partial loss}{\partial O} \frac{\partial O}{\partial H_1} \frac{\partial H_1}{\partial Z} \frac{\partial Z}{\partial b_1}\\ & = \frac{1}{n} (\hat{y} - y) {W_2}^\mathrm{T} \frac{\partial H_1}{\partial Z}\\ & = \begin{cases} 0 &\text{if } {Z}_{ij} < 0\\ \frac{1}{n} \sum_n (\hat{y} - y) {W_2}^\mathrm{T} &\text{if } {Z}_{ij} \geq 0 \end{cases} \end{aligned}$$

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

参数更新

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

$$W_2 := W_2 - \alpha \frac{\partial \mathrm{loss}}{\partial W_2}\\ b_2 := b_2 - \alpha \frac{\partial \mathrm{loss}}{\partial b_2}\\ W_1 := W_1 - \alpha \frac{\partial \mathrm{loss}}{\partial W_1}\\ b_1 := b_1 - \alpha \frac{\partial \mathrm{loss}}{\partial b_1}\\$$

其中，$\alpha$是学习率

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

1. 导入数据

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

import matplotlib.pyplot as plt
%matplotlib inline

from time import time

import numpy as np
from sklearn.datasets import load_digits
from sklearn.model_selection import train_test_split

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

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

2. 数据预处理

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

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$的概率，那我们只要取最大的那个数的下标，就知道模型认为这个样本是哪类了。

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

trainY_mat.shape, testY_mat.shape

3. 参数初始化

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

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

# 测试样例
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的计算

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

# 测试样例
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$激活函数

def ReLU(X):
    '''
    ReLU激活函数
    
    Parameters
    ----------
    X: np.ndarray，待激活的矩阵
    
    Returns
    ----------
    activations: np.ndarray, 激活后的矩阵
    
    '''
    activations = X.copy()
    activations[activations < 0] = 0
    return activations

# 测试样例
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$激活

$$\mathrm{softmax}(O_i) = \frac{\exp{(O_i)}}{\sum^{K}_{k=1} \exp{(O_k)}}$$

def my_softmax(O):
    '''
    softmax激活
    '''
    denominator = np.exp(O).sum(axis = 1, keepdims = True)
    return np.exp(O) / denominator

# 测试样例
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$的向量，每个元素减去他们中的最大值，即

$$\mathrm{softmax}(O_i) = \mathrm{softmax}(O_i - \mathrm{max(O)})$$

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

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

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

# 测试样例
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

接下来是实现损失函数，交叉熵损失函数：

$$\mathrm{loss} = - \frac{1}{n} \sum_n \sum^{K}_{k=1} y_k \log{(\hat{y_k})}$$

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

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

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

$$\begin{aligned} \log{\frac{\exp{(O_i)}}{\sum_K \exp{(O_k)}}} &= \log{\frac{\exp{(O_i - \mathrm{max}(O))}}{\sum_K \exp{(O_k - \mathrm{max}(O))}}}\\ &= O_i - \mathrm{max}(O) - \log{\sum_K \exp{(O_k - \mathrm{max}(O))}} \end{aligned}$$

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

我们先编写log_softmax

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

# 测试样例
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

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

# 测试样例
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

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

# 测试样例
parameterst = initialize(50, 10)
Ot = forward(trainX, parameterst)
print(Ot.mean()) # 0.000600192658464

5. 反向传播

先计算梯度

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

# 测试样例
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

梯度下降，参数更新

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']

反向传播，参数更新

# 测试样例
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

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)

# 测试样例
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. 训练

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

# 测试样例
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. 绘制模型损失值变化曲线

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$类的概率，我们要找最大的那个概率，对应的是哪个类。

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

# 测试样例
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轮

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

计算测试集精度

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

绘制损失值变化曲线

plot_loss_curve(training_loss_list, testing_loss_list)

更换数据集

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

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

trainX.shape, trainY.shape, trainY_mat.shape, testX.shape, testY.shape, testY_mat.shape

绘制训练集前10个图像

_, 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

stand = StandardScaler()
trainX_normalized = stand.fit_transform(trainX)
testX_normalized = stand.fit_transform(testX)

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

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

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

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

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