9、反向传播

'''
识别书写数字例子：
反向传播，学习神经网络的参数。
'''
import numpy as np
import matplotlib.pyplot as plt
from scipy.io import loadmat
import scipy.optimize as opt
from sklearn.metrics import classification_report


# 加载数据集
def load_mat(path):
    data = loadmat(path)  # loadmat读取mat文件
    # 这里的数据为MATLAB的格式，所以要使用SciPy.io的loadmat函数
    X = data['X']  # (5000,400)
    y = data['y'].flatten()  # (5000, )
    return X, y


# 随机输出100个样本图片,参数：样本数据X
def plot_100_image(X):
    index = np.random.choice(range(5000), 100)  # 从[0,5000)区间内随即取出100个数字，并将这100个数，组成一个数组输出 (100, )
    images = X[index]  # 取出100个样本图片 (100, 400)
    fig, ax_array = plt.subplots(nrows=10, ncols=10, sharey=True, sharex=True, figsize=(8, 8))  # 所有子图共享x轴和y轴
    # 显示100个样本图
    for row in range(10):
        for col in range(10):
            ax_array[row, col].matshow(images[10 * row + col].reshape((20, 20)), cmap='gray_r')
    plt.yticks([])
    plt.xticks([])
    plt.show()


# 加载数据集
X, y = load_mat('data\ex4data1.mat')
# 随机输出100个样本图片
# plot_100_image(X)


# 将标签值(1,2,...,10)转化成非线性相关的向量，例如将y[0]=6转化为y[0]=[0,0,0,0,0,1,0,0,0,0]
# 方法一
# def expend_y(y):
#     result = []
#     # 将y中的每个类别转化为一个向量
#     for i in y:
#         y_array = np.zeros(10)
#         y_array[i - 1] = 1  # 因为数组下标从0开始，所以要减1
#         result.append(y_array)
#     return np.array(result)
# 方法二
from sklearn.preprocessing import OneHotEncoder


def expend_y(y):
    encoder = OneHotEncoder(sparse=False)
    '''
    初始化OneHotEncoder实例时，默认sparse参数为True，
    编码后返回的是一个稀疏矩阵的对象，如果要使用一般要调用toarray()方法转化成array对象。
    若将sparse参数设置为False，则直接生成array对象，可直接使用。
    '''
    y_onehot = encoder.fit_transform(y.reshape(-1, 1))
    '''
    fit_transform()的作用就是先拟合数据，然后转化它将其转化为标准形式，一般应用在训练集中。
    '''
    return y_onehot


# print(expend_y(y))
'''
[[0. 0. 0. ... 0. 0. 1.]
 [0. 0. 0. ... 0. 0. 1.]
 [0. 0. 0. ... 0. 0. 1.]
 ...
 [0. 0. 0. ... 0. 1. 0.]
 [0. 0. 0. ... 0. 1. 0.]
 [0. 0. 0. ... 0. 1. 0.]]
'''

# 获取训练数据集，并做相应处理
raw_X, raw_y = load_mat('data\ex4data1.mat')
X = np.insert(raw_X, 0, 1, axis=1)  # 加一个元素值全为1的列
y = expend_y(raw_y)  # 转化为非线性相关向量形式


# print(X.shape, y.shape)  # (5000, 401) (5000, 10)

# 读取权重
# 这些参数的维度由神经网络的大小决定，第二层有25个单元，输出层有10个单元。
def load_weight(path):
    data = loadmat(path)
    return data['Theta1'], data['Theta2']


# 读取权重
t1, t2 = load_weight('data\ex4weights.mat')
# print(t1.shape, t2.shape)  # (25, 401) (10, 26)

'''
当我们使用高级优化方法来优化神经网络时，我们需要将多个参数矩阵展开，
之后才能传入优化函数，然后再恢复形状。
'''


# 展开参数矩阵
def serialize(a, b):
    return np.r_[a.flatten(), b.flatten()]  # np.r_是按列连接两个矩阵


theta = serialize(t1, t2)  # 扁平化权重参数


# print(theta.shape)  # (10285, )其中25*401+10*26=10285


# 提取参数(恢复形状)
def deserialize(seq):
    return seq[:25 * 401].reshape(25, 401), seq[25 * 401:].reshape(10, 26)
    # seq[:25 * 401]取seq[0,1,...,25*401-1] ; seq[25 * 401:]取seq[25*401,...,10285-1]


# sigmoid函数
def sigmoid(z):
    return 1 / (1 + np.exp(-z))


# 正向传播
def feed_forward(theta, X):
    # 返回每一层神经元的输入和输出
    t1, t2 = deserialize(theta)
    a1 = X
    z2 = a1 @ t1.T
    a2 = np.insert(sigmoid(z2), 0, 1, axis=1)  # 插入偏置单元
    z3 = a2 @ t2.T
    a3 = sigmoid(z3)
    return a1, z2, a2, z3, a3


# 不带正则化项的代价函数
def cost(theta, X, y):
    a1, z2, a2, z3, h = feed_forward(theta, X)
    J = -y * np.log(h) - (1 - y) * np.log(1 - h)
    return J.sum() / len(X)


# 不带正则化项的代价函数
print('不带正则化项的代价值：', cost(theta, X, y))  # 0.2876291651613189


# 正则化代价函数
def regularized_cost(theta, X, y, l=1):
    # 正则化时忽略每层的偏置项，也就是参数矩阵的第一列
    t1, t2 = deserialize(theta)
    reg = np.sum(t1[:, 1:] ** 2) + np.sum(t2[:, 1:] ** 2)
    return l / (2 * len(X)) * reg + cost(theta, X, y)


print('加正则化项的代价值：', regularized_cost(theta, X, y, 1))  # 0.38376985909092365


# sigmoid函数的导数
def sigmoid_gradient(z):
    return sigmoid(z) * (1 - sigmoid(z))


'''
参数随机初始化：
当我们训练神经网络时，随机初始化参数是很重要的，可以打破数据的对称性。
一个有效的策略是在均匀分布(−e，e)中随机选择值，
我们可以选择 e = 0.12 这个范围的值来确保参数足够小，使得训练更有效率。
'''


# 参数随机初始化
def random_init(size):
    # 从一个均匀分布[low,high)中随机采样,数量为size。结果返回数组。
    return np.random.uniform(-0.12, 0.12, size)


# 各种参数在网络中的维度
a1, z2, a2, z3, h = feed_forward(theta, X)
print('a1', a1.shape, 't1', t1.shape)
print('z2', z2.shape, 'a2', a2.shape, 't2', t2.shape)
print('z3', z3.shape, 'a3', h.shape)
'''
a1 (5000, 401) t1 (25, 401)
z2 (5000, 25) a2 (5000, 26) t2 (10, 26)
z3 (5000, 10) a3 (5000, 10)
'''


# 无正则化项的梯度
def gradient(theta, X, y):
    t1, t2 = deserialize(theta)
    a1, z2, a2, z3, h = feed_forward(theta, X)
    # 计算误差
    d3 = h - y  # (5000,10)
    d2 = d3 @ t2[:, 1:] * sigmoid_gradient(z2)  # (5000,25)
    # 计算梯度
    D2 = d3.T @ a2  # (10,26)
    D1 = d2.T @ a1  # (25,401)
    # 总梯度
    D = (1 / len(X)) * serialize(D1, D2)  # (10285, )
    return D


# 带正则化项的梯度 , 不惩罚偏置单元
def regularized_gradient(theta, X, y, l=1):
    D1, D2 = deserialize(gradient(theta, X, y))
    # 将偏置单元的梯度设为0
    t1[:, 0] = 0
    t2[:, 0] = 0
    # 加入正则化项后的梯度
    reg_D1 = D1 + (l / len(X)) * t1
    reg_D2 = D2 + (l / len(X)) * t2
    return serialize(reg_D1, reg_D2)  # 返回正则化总梯度


# # 梯度检测(运行很慢，检测完要及时注释掉)
# def gradient_checking(theta, X, y, e):
#     # 计算每个参数的数值梯度（理论梯度）
#     def a_numeric_grad(plus, minus):
#         return (regularized_cost(plus, X, y) - regularized_cost(minus, X, y)) / (2 * e)
#
#     # 存储这些理论梯度值
#     numeric_grad = []
#     for i in range(len(theta)):
#         plus = theta.copy()  # 深拷贝 （分配新地址）
#         minus = theta.copy()
#         plus[i] = plus[i] + e
#         minus[i] = minus[i] - e
#         grad_i = a_numeric_grad(plus, minus)
#         numeric_grad.append(grad_i)
#
#     numeric_grad = np.array(numeric_grad)  # 理论梯度
#     analytic_grad = regularized_gradient(theta, X, y)  # 正则化梯度
#     diff = np.linalg.norm(numeric_grad - analytic_grad) / np.linalg.norm(numeric_grad + analytic_grad)
#     # np.linalg.norm(ord=None) 默认情况下是求整个矩阵元素平方和再开根号
#     print('(如果反向传播实现正确，diff将小于10e-9)diff={}'.format(diff))
#
#
# gradient_checking(theta, X, y, e=0.0001)  # 梯度检测 diff=3.1753258526793174e-09


# 优化参数
def nn_training(X, y):
    init_theta = random_init(10285)  # 初始化参数
    res = opt.minimize(fun=regularized_cost, x0=init_theta, args=(X, y, 1), method='TNC', jac=regularized_gradient,
                       options={'maxiter': 400})  # maxiter设置要执行的最大迭代数
    return res


res = nn_training(X, y)  # 优化参数
# print(res)  # 因为权重参数是随机初始化的，所以每次结果都是不同的
print('本次最优权重参数为:', res.x)


# 准确率报告
def accuracy(theta, X, y):
    _, _, _, _, h = feed_forward(theta, X)  # 输出层h
    y_pred = np.argmax(h, axis=1) + 1  # 预测值（按行取最大值的索引并加1）
    print(classification_report(y, y_pred))


print('本次准确率报告:\n')  # 打印准确率报告
accuracy(res.x, X, raw_y)

'''
可视化隐藏层:
对于我们所训练的网络，注意到θ1中每一行都是一个401维的向量，代表每个隐藏层单元的参数。
如果我们忽略偏置项，我们就能得到400维的向量，这个向量代表每个样本输入到每个隐层单元的像素的权重。
因此可视化的一个方法是，reshape这个400维的向量为（20，20）的图像然后输出。
'''


# 可视化隐藏层 (这里只可视化t1的前25个隐藏层单元)
def plot_hidden(theta):
    t1, _ = deserialize(theta)
    t1 = t1[:, 1:]
    fig, ax_array = plt.subplots(5, 5, sharex=True, sharey=True, figsize=(6, 6))
    for r in range(5):
        for c in range(5):
            ax_array[r, c].matshow(t1[r * 5 + c].reshape(20, 20), cmap='gray_r')
            plt.xticks([])
            plt.yticks([])
    plt.show()


plot_hidden(res.x)  # 可视化优化参数后的隐藏层

详解：

正向传播：

正则化代价函数：

反向传播：

梯度检测，计算数值梯度公式：

正则化项梯度：

1、
函数原型： numpy.random.uniform(low,high,size)

功能：

从一个均匀分布[low,high)中随机采样，注意定义域是左闭右开，即包含low，不包含high.

参数介绍:

low: 采样下界，float类型，默认值为0；
high: 采样上界，float类型，默认值为1；
size: 输出样本数目，为int或元组(tuple)类型，例如，size=(m,n,k), 则输出mnk个样本，缺省时输出1个值。

返回值：

ndarray类型，其形状和参数size中描述一致。

运行结果：