机器学习练习四：用Python实现BP神经网络

mmd
2,517 阅读7分钟

数据介绍：ex4data1中有5000个训练实例。其中每个训练实例是一个20像素乘20像素灰度图像的数字。每个像素由一个浮点数表示，表示该位置的灰度强度。这个20x20的像素网格被“展开”成一个400维的向量。这些训练示例中的每一个都变成了数据矩阵X中的单行。这就得到了一个5000×400矩阵X，其中每一行都是一个手写数字图像的训练示例。

数据集是在MATLAB的本机格式，所以要加载它在Python，我们需要使用一个SciPy库。

1.前向传播算法

1.1 展示数据

import numpy as np
import pandas as pd
import matplotlib
import matplotlib.pyplot as plt
from scipy.io import loadmat

# 加载数据
data = loadmat('ex4data1.mat')
X, y = data['X'], data['y']
X.shape,y.shape

将数据可视化：

def plot_100figs():
    """
    画100幅图
    """
    # 随机100个索引,不重复
    sample_list = np.random.choice(np.arange(X.shape[0]), 100, replace=False)
    # 从数据中取出这100个值
    sample_data = X[sample_list, :]
    # 画图, sharex,sharey:子图共享x轴,y轴, ax_array：一个10*10的矩阵
    fig, ax_array = plt.subplots(nrows=10, ncols=10, sharex=True, sharey=True, figsize=(8, 8))
    # 画子图
    for r in range(10):
        for c in range(10):
            ax_array[r, c].matshow(sample_data[10 * r + c, :].reshape((20, 20)).T,
                             cmap=matplotlib.cm.binary)
    # 去掉刻度线
    plt.xticks([])
    plt.yticks([])
    plt.show()

1.2 One-hot编码

为了训练一个神经网络，我们需要将标签重新编码为仅包含0或1的向量。神经网络本质上是由许多的单节点逻辑回归（二分类）组成的。

将y转换为如下形状： 

def onehot(y, class_num):
    """
    onehot编码
    """
    # (5000, 10)
    y_onehot = np.zeros((y.shape[0], class_num))
    # 取y的值
    for row in range(y.shape[0]):
        # row -> [0,5000)
        # 把相应的位置置1
        cols = y[row] - 1
        y_onehot[row, cols] = 1
    return y_onehot

由于从原目标值中取出的数字是1-10，而新标签中每行的索引为0-9，故需cols = y[row] - 1使其对应。

1.3 参数展开(Unrolling Parameters)

为了使用诸如“scipy.optimize.minimize”这样的优化函数，我们需要“展开”所有元素并将它们放入一个长向量中。

展开元素：

def serialize(a, b):
    """
    展开元素
    """
    return np.concatenate((np.ravel(a), np.ravel(b)))

恢复元素（恢复成(25, 401)和(10, 26)的形状）：

def deserialize(seq):
    """
    恢复元素
    """
    # Θ的形状是(25, 401)和(10, 26)
    return seq[:25 * 401].reshape(25, 401), seq[25 * 401:].reshape(10, 26)

1.4 数据读取与形状改变

将数据改变为算法需要的形状。

# 读取训练好的权重
weight = loadmat('ex4weights.mat')
theta1, theta2 = weight['Theta1'], weight['Theta2']
# 插入全1行x0
X = np.insert(X, 0, values=np.ones(X.shape[0]), axis=1)  
# 将Θ展开
theta = serialize(theta1, theta2)

1.5 前馈运算

模型图示：

Sigmoid函数：

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

前馈算法：

def feed_forward(theta, X, y):
    # 将Θ恢复
    theta1, theta2 = deserialize(theta)
    # 第一层 --> 第二层
    a1 = X
    z2 = X @ theta1.T  # (5000, 401) @ (401, 25) = (5000, 25)
    a2 = sigmoid(z2)
    a2 = np.insert(a2, 0, values=np.ones(a2.shape[0]), axis=1)  # intercept:(5000, 26)
    # 第二层 --> 第三层
    z3 = a2 @ theta2.T  # (5000, 26) @ (26, 10) = (5000, 10)
    h = sigmoid(z3)
    # 这些参数需要在反向传播时使用
    return a1, z2, a2, z3, h

1.6 损失函数

  • 普通的损失函数

两个求和将输出层中每个单元的逻辑回归损失相加，现在我们有y 和 h_{\theta} \in R^{5000 \times 10}。

忽略m和K的维度，the eqation = -y*log(h_{\theta}) - (1-y)*log(1-h_{\theta}) ，最后将所得到的2维数组求和再除以m即可。

def cost(theta, X, y):
    X = np.matrix(X)
    y = np.matrix(y)
    # 前馈运算获取hθ(x)
    _,_,_,_,h = feed_forward(theta, X, y)
    # 计算
    # np.multiply是对应位置相乘
    first = np.multiply(-y, np.log(h))
    second = np.multiply((1 - y), np.log(1 - h))
    front = np.sum(first - second) / (len(X))
    return front

结果：0.2876291651613189

注意np.multiply()和@的区别，前者是对应位置相乘，后者是矩阵乘法。

  • 正则化损失函数
    正则化时忽略\theta^{(1)}和\theta^{(2)}的第一列，偏置不需要优化。

显然前半部分与普通损失函数一致，只需加正则项即可。

def regular_cost(theta, X, y, lambd):
    front = cost(theta, X, y)
    # 函数后半部分,忽略θ(1)和θ(2)的第一列
    # 恢复θ
    theta1, theta2 = deserialize(theta)
    last = lambd / (2 * len(X)) * (np.sum(np.power(theta1[:, 1:], 2)) + np.sum(np.power(theta2[:, 1:], 2)))
    return front + last

结果：0.38376985909092365

2.反向传播算法

使用反向传播算法的目的仍然是求出损失函数J(\theta)的导数，然后将其用于梯度下降算法优化中。关于线性回归和逻辑回归算法梯度下降的信息请看我之前的文章：机器学习练习二：用Python实现逻辑回归。

2.1 Sigmoid gradient（梯度sigmoid函数，即sigmoid函数的导数）

def sigmoid_gra(z):
    return sigmoid(z) * (1 - sigmoid(z))

2.2 BP算法

上下标的意义：

l代表目前计算的是第几层，i代表目前计算的是第几个样本，j代表目前计算的是这一层中的第几个节点。

模型图示：

伪代码可以表示为：

1.数据集：

2.设置Δ^{l}_{ij}=0

3.for-loop:

涉及到的数学公式：

  1. Δ^{(l)}_{ij} := Δ^{(l)}_{ij} + a^{(l)}_jδ^{(l+1)}_i

代码表示：

此时不考虑正则化，即\lambda=0，

其中δ^{(1)}是不需要求解的，因为信号输入时不需要考虑损失。再次强调注意np.multiply()和@的区别，前者是对应位置相乘，后者是矩阵乘法。

def backpropa(theta, X, y):
    """
    反向传播算法
    """
    m = X.shape[0]
    # 恢复θ
    theta1, theta2 = deserialize(theta)
    # set Δ=0(for all i,j,l)
    # 与Θ的形状一致
    delta1 = np.zeros(theta1.shape)  # (25, 401)
    delta2 = np.zeros(theta2.shape)  # (10, 26)
    # 获取前向传播中生成的参数
    a1, z2, a2, z3, h = feed_forward(theta, X, y)
    # 设置for-loop
    for i in range(m):
        # 获取各层的参数
        a1i = a1[i, :]  # (401,)
        z2i = z2[i, :]  # (25,)
        a2i = a2[i, :]  # (26,)
        
        hi = h[i, :]  # (10,)
        yi = y[i, :]  # (10,)
        
        # 计算各层δ
        # 没有δ1，不需要对输入层考虑误差
        d3i = hi - yi  # (10,)
        z2i = np.insert(z2i, 0, np.ones(1))  # (26,)
        d2i = np.multiply(theta2.T @ d3i, sigmoid_gra(z2i))  # [(26, 10) @ (10,)] * (26,) = (26,)
        
        # 计算Δ
        delta2 += np.mat(d3i).T @ np.mat(a2i)  # (10, 1) @ (1, 26)
        # 去掉d2的第0列
        delta1 += np.mat(d2i[1:]).T @ np.mat(a1i)  # (25, 1) @ (1, 401)
    
    # 得到Dij
    delta1 = delta1 / m
    delta2 = delta2 / m
    # 返回长序列
    return serialize(delta1, delta2)

返回值为展开后的D_{ij}的值，亦即J(\theta)的导数值。

2.3 正则化BP算法

此时引入了超参数\lambda。

代码实现：

显然当j=0，即计算该层第一个节点(偏置项)的D_{ij}时，无需正则化。而当j>=1时，前半部分与普通bp算法相同，只需加入正则项即可。代码的实现中用将θ的第一列置零的方法达到上述目的。

参数Θ中的第一列为偏置项，不需要优化。

def regularized_bp(theta, X, y, lambd=1):
    """偏置项（j = 0时）不进行优化"""
    m = X.shape[0]
    # 前半部分
    delta = backpropa(theta, X, y)
    # 恢复θ
    theta1, theta2 = deserialize(theta)
    # 恢复Dij
    delta1, delta2 = deserialize(delta)
    theta1[:, 0] = 0  # 将j=0时的θ置零
    reg_term_d1 = (lambd / m) * theta1
    delta1 += reg_term_d1

    theta2[:, 0] = 0  # j=0时的θ置零
    reg_term_d2 = (lambd / m) * theta2
    delta2 += reg_term_d2
    # 需要将theta拼接起来
    return serialize(delta1, delta2)

3.模型训练

3.1 随机初始化变量

将所有的权重初始化为零将使神经网络失去作用（这种方法在逻辑回归中是可行的）。因为当计算反向传播算法时，所有节点上的权重都将更新为相同的值。同理，如果我们初始所有的参数都为一个非零的数，结果也是一样的。这也称为参数对称。

为了打破对称，我们需要随机初始化变量。

# 随机初始化变量，打破参数对称
def random_init(size):
    # 均匀分布抽取样本在[-0.12,0.12]之间
    return np.random.uniform(-0.12, 0.12, size)

3.2 模型训练

# 使用 scipy.optimize.minimize 去寻找参数
import scipy.optimize as opt
def nn_training(X, y):
    # 随机初始化变量theta
    init_theta = random_init(10285)  # 25*401 + 10*26
    # 
    res = opt.minimize(fun=regular_cost,
                       x0=init_theta,
                       args=(X, y, 1),
                       method='TNC',
                       jac=regularized_bp,
                       options={'maxiter': 1000})
return res

运行结果：

可见经过1000次迭代将损失降低到了0.305。

3.3 计算准确率

关于准确率函数的解释请看：机器学习练习二：用Python实现逻辑回归。

def show_acc(theta, X, y):
    # 取得结果概率
    _, _, _, _, h = feed_forward(theta, X, y)
    # 预测值矩阵
    y_predict = np.mat(np.argmax(h,axis=1) + 1)
    y_true = np.mat(data['y']).ravel()  # 真实值矩阵
    # 矩阵进行比较返回各元素比对布尔值矩阵，列表进行比较返回整个列表的比对布尔值
    accuracy = np.mean(y_predict == y_true)
    print('accuracy = {}%'.format(accuracy * 100))

结果：

3.4 隐层可视化

与前面画100图的函数类似。

在你的训练网络中，你应该发现隐藏的单元大致对应于在输入中寻找笔画和其他图案的探测器。

def plot_hidden(final_theta):
    """
    将隐层画出来
    """
    # 获取theta1
    theta1, _ = deserialize(final_theta)  # (25, 401)
    # 去掉偏置列
    theta1 = theta1[:, 1:]
    # 画图, sharex,sharey:子图共享x轴,y轴, ax_array：一个5*5的矩阵
    fig, ax_array = plt.subplots(nrows=5, ncols=5, sharex=True, sharey=True, figsize=(8, 8))
    # 画子图
    for r in range(5):
        for c in range(5):
            ax_array[r, c].matshow(theta1[5 * r + c, :].reshape((20, 20)).T,
                             cmap=matplotlib.cm.binary)
    # 去掉刻度线
    plt.xticks([])
    plt.yticks([])
    plt.show()

结果：

avatar
文章
阅读
粉丝