吴恩达《机器学习》课后测试Ex4：神经网络反向传播（详细Python代码注解）

基于吴恩达《机器学习》课程
参考黄海广的笔记

---

import numpy as np
from scipy.io import loadmat
from sklearn.preprocessing import OneHotEncoder
from scipy.optimize import minimize

对于这个练习，我们将再次处理手写数字数据集，这次使用反向传播的前馈神经网络。 我们将通过反向传播算法实现神经网络成本函数和梯度计算的非正则化和正则化版本。 我们还将实现随机权重初始化和使用网络进行预测的方法。

sigmoid 函数

$h_{\theta }\left(x\right)=\frac{1}{1+e^{-{\theta }^{T}X}}$

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

前向传播函数

我们构建的神经网络输入层为400个单位（加偏置401），25个单位的隐藏层（加偏置26），以及10个单位输出层（对应我们的一个one-hot编码类标签）。

def forward_propagate(X, theta1, theta2):
    # INPUT：参数值theta，数据X  OUTPUT：当前参数值下前项传播结果  TODO：根据参数和输入的数据计算前项传播结果
    # STEP1：获取样本个数
    m = X.shape[0]
    # STEP2：实现神经网络正向传播  
    a1 =np.insert(X, 0, values=np.ones(m), axis=1) # 给X矩阵插入一列1元素，axis=1表列
    z2 =a1 * theta1.T
    a2 =np.insert(sigmoid(z2), 0, values=np.ones(m), axis=1)  # 插入一列1元素
    z3 = a2 * theta2.T
    h = sigmoid(z3)
    return a1, z2, a2, z3, h

代价函数

$J(\Theta )=-\frac{1}{m}\left[\sum _{i=1}^m\sum _{k=1}^K{y_k}^{(i)}\log {(h_{\Theta }(x^{(i)}))}_k+\left(1-y_k^{(i)}\right)\log \left(1-{\left(h_{\Theta }\left(x^{(i)}\right)\right)}_k\right)\right]+\frac{\lambda }{2m}\sum _{l=1}^{L-1}\sum _{i=1}^{s_l}\sum _{j=1}^{s_{l+1}}{\left({\Theta }_{ji}^{(l)}\right)}^2$

这里只有两层，所以偏置也可以写成：

$\frac{\lambda }{2m}\left[\sum _{i=1}^{s_l}\sum _{j=1}^{s_{l+1}}{\left({\Theta }_{ji}^{(1)}\right)}^2+\sum _{i=1}^{s_l}\sum _{j=1}^{s_{l+1}}{\left({\Theta }_{ji}^{(2)}\right)}^2\right]$

def cost(params, input_size, hidden_size, num_labels, X, y, lamda):
    # INPUT：神经网络参数，输入层维度，隐藏层维度，训练数据及标签，正则化参数
    # OUTPUT：当前参数值下的代价函数  TODO：根据公式计算代价函数
    # STEP1：获取样本个数,这里有5000
    m = X.shape[0]
    # STEP2：将矩阵X,y转换为numpy型矩阵
    X = np.matrix(X)
    y = np.matrix(y)
    # STEP3：从params中获取神经网络参数，并按照输入层维度和隐藏层维度重新定义reshape参数的维度，拆成两个theta
    theta1 = np.matrix(np.reshape(params[:hidden_size * (input_size + 1)], (hidden_size, (input_size + 1))))
    theta2 = np.matrix(np.reshape(params[hidden_size * (input_size + 1):], (num_labels, (hidden_size + 1))))
    # STEP4：调用前面写好的前项传播函数
    a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2)
    # STEP5：初始化代价函数
    J = 0
    # STEP6：根据公式计算代价函数
    for i in range(m):  # 遍历每个样本
        first_term = np.multiply(-y[i, :], np.log(h[i, :]))
        second_term = np.multiply((1 - y[i, :]), np.log(1 - h[i, :]))
        J += np.sum(first_term - second_term)
    J = J / m
    # STEP7：计算代价函数的正则化部分，不包括theta0
    J += (float(lamda) / (2 * m)) * (np.sum(np.power(theta1[:, 1:], 2)) + np.sum(np.power(theta2[:, 1:], 2)))
    return J

计算Sigmoid函数的倒数

接下来是反向传播算法。 反向传播参数更新计算将减少训练数据上的网络误差。 我们需要的第一件事是计算我们之前创建的Sigmoid函数的梯度的函数。

def sigmoid_gradient(z):
    return np.multiply(sigmoid(z), (1 - sigmoid(z)))

反向传播函数

现在我们准备好实施反向传播来计算梯度。

由于反向传播所需的计算是代价函数中所需的计算过程，我们实际上将扩展代价函数以执行反向传播并返回代价和梯度。

注意！这里极易混淆A * B与np.multiply(A，B)使用。 前者是矩阵乘法，后者是元素乘法（除非A或B是标量值，在这种情况下没关系）。

def backprop(params, input_size, hidden_size, num_labels, X, y, lamda):
    # INPUT：神经网络参数，输入层维度，隐藏层维度，训练数据及标签，正则化参数
    # OUTPUT：当前参数值下的代价函数、梯度
    # STEP1：获取样本个数
    m = X.shape[0]
    # STEP2：将矩阵X,y转换为numpy型矩阵
    X = np.matrix(X)
    y = np.matrix(y)
    # STEP3：从params中获取神经网络参数，并按照输入层维度和隐藏层维度重新定义参数的维度
    theta1 = np.matrix(np.reshape(params[:hidden_size * (input_size + 1)], (hidden_size, (input_size + 1))))
    theta2 = np.matrix(np.reshape(params[hidden_size * (input_size + 1):], (num_labels, (hidden_size + 1))))
    # STEP4：调用前面写好的前项传播函数
    a1, z2, a2, z3, h =forward_propagate(X, theta1, theta2)
    # STEP5：初始化
    J = 0
    delta1 = np.zeros(theta1.shape)  # (25, 401) delta用于存theta的梯度
    delta2 = np.zeros(theta2.shape)  # (10, 26)
    # STEP6：计算代价函数(调用函数)
    for i in range(m):
        first_term = np.multiply(-y[i,:], np.log(h[i,:]))
        second_term = np.multiply((1 - y[i,:]), np.log(1 - h[i,:]))
        J += np.sum(first_term - second_term)
    J = J / m
    # STEP7：实现反向传播
    for t in range(m):  #遍历每个样本
        a1t = a1[t,:]  # (1, 401)
        z2t = z2[t,:]  # (1, 25)
        a2t = a2[t,:]  # (1, 26)
        ht = h[t,:]  # (1, 10)
        yt = y[t,:]  # (1, 10)
        d3t = ht - yt
        z2t = np.insert(z2t, 0, values=np.ones(1)) # (1, 26) 插入偏置
        d2t = np.multiply((theta2.T * d3t.T).T, sigmoid_gradient(z2t))  # (1, 26)
        # 将所有样本的梯度加起来
        delta1 = delta1 + (d2t[:,1:]).T * a1t  # d2t[:,1:]去掉d2t的第一维偏置，这里得到代价函数对theta1的梯度
        delta2 = delta2 + d3t.T * a2t  # d3t.T * a2t即为代价函数对theta2的梯度
    # STEP8：加入正则化的梯度
    delta1[:,1:] = delta1[:,1:] + (theta1[:,1:] * lamda) / m  # 正则化的theta去掉第一项偏置
    delta2[:,1:] =  delta2[:,1:] + (theta2[:,1:] * lamda) / m
    # STEP9：将梯度矩阵转换为单个数组
    grad = np.concatenate((np.ravel(delta1), np.ravel(delta2)))  # np.concatenate拼接数组
    return J, grad

主程序代码

我们在练习3中使用的数据集是相同的，所以我们用同样方法来加载数据。

def ex4_1():
    data = loadmat('ex4data1.mat')
    X = data['X']
    y = data['y']
    X.shape, y.shape
    print(X.shape, y.shape)  # 看下维度
    print(np.unique(data['y']))  # 看下有几类标签

输出如下：

(5000, 400) (5000, 1)
[ 1  2  3  4  5  6  7  8  9 10]

可以看到y中有1-10共十类标签，我们需要对我们的y标签进行一次one-hot编码。

one-hot编码将k类标签转换为长度为k的向量，其中索引n为“hot”（1），而其余为0。例如便签2就是[0,1,0,0,0,0,0,0,0,0] 。这样y_onehot就变为了5000*10矩阵，一个训练数据标签为10维，对应着神经网络的十个分类输出。

Scikitlearn有一个内置的实用程序，我们可以使用这个。

# 代码接上
    encoder = OneHotEncoder(sparse=False)   # 默认sparse参数为True，编码后返回的是一个稀疏矩阵的对象
    y_onehot = encoder.fit_transform(y)
    print(y_onehot.shape)
    print(y[0], y_onehot[0, :])

这里sparse=False返回的是array对象，如果默认True，例如标签2会返回稀疏矩阵对象(2,1)表示第二位为1。此处输出如下：

(5000, 10)
[10] [0. 0. 0. 0. 0. 0. 0. 0. 0. 1.]

接下来就是初始化数据：

# 代码接上
    # 初始化设置
    input_size = 400    # 输入训练数据维数
    hidden_size = 25    # 隐藏层维数
    num_labels = 10     # 分类数
    lamda = 1   # 正则化参数
    # 随机初始化整个网络的参数数组，维度25*401+10*26，-0.5再乘0.25让其处于正负0.125之间
    params = (np.random.random(size=hidden_size * (input_size + 1) + num_labels * (hidden_size + 1)) - 0.5) * 0.25
    m = X.shape[0]
    X = np.matrix(X)
    y = np.matrix(y)
    # 将参数数组解开为每个层的参数矩阵
    theta1 = np.matrix(np.reshape(params[:hidden_size * (input_size + 1)], (hidden_size, (input_size + 1))))
    theta2 = np.matrix(np.reshape(params[hidden_size * (input_size + 1):], (num_labels, (hidden_size + 1))))
    print(theta1.shape, theta2.shape)

输出如下：

(25, 401) (10, 26)

然后调用前向传播函数，计算代价。由于随机初始值不同，每次计算代价也有区别。

# 代码接上  
    a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2)
    print(a1.shape, z2.shape, a2.shape, z3.shape, h.shape)
    print(cost(params, input_size, hidden_size, num_labels, X, y_onehot, lamda))

输出如下：

(5000, 401) (5000, 25) (5000, 26) (5000, 10) (5000, 10)
6.92922699465799

接下来就可以进行反向传播了，由于随机初始值不同，每次计算代价也有区别。

# 代码接上
    # 反向传播
    J, grad = backprop(params, input_size, hidden_size, num_labels, X, y_onehot, lamda)
    print(J, grad.shape)

输出如下：

7.0856460879022025 (10285,)

接下来调用minimize最小化参数：

# 代码接上
    # 最小化参数
    fmin = minimize(fun=backprop, x0=params, args=(input_size, hidden_size, num_labels, X, y_onehot, lamda),
                    method='TNC', jac=True, options={'maxiter': 250})
    # options：用来控制最大的迭代次数  jac=True则fun里要有两个返回值：代价函数、梯度
    print(fmin)

输出如下，可以看到，找到的参数在fmin.x中：

     fun: 0.23857388315436762
     jac: array([-4.51400702e+01,  2.42009028e-05, -1.50489283e-05, ...,
       -2.73699447e+00,  7.14064502e-01, -9.48905699e-01])
 message: 'Max. number of function evaluations reached'
    nfev: 251
     nit: 15
  status: 3
 success: False
       x: array([ 1.12513   ,  0.12100451, -0.07524464, ..., -0.90956441,
       -3.88185143, -0.57911355])

由于目标函数不太可能完全收敛，我们对迭代次数进行了限制。 我们的总代价已经下降到0.5以下，这是算法正常工作的一个很好的指标。

使用它找到的的参数，并通过网络前向传播以获得预测，然后计算准确度。

# 代码接上
    # 用找到的参数进行预测并计算准确度
    X = np.matrix(X)
    theta1 = np.matrix(np.reshape(fmin.x[:hidden_size * (input_size + 1)], (hidden_size, (input_size + 1))))
    theta2 = np.matrix(np.reshape(fmin.x[hidden_size * (input_size + 1):], (num_labels, (hidden_size + 1))))
    a1, z2, a2, z3, h = forward_propagate(X, theta1, theta2)
    y_pred = np.array(np.argmax(h, axis=1) + 1)
    print(y_pred)
    correct = [1 if a == b else 0 for (a, b) in zip(y_pred, y)]
    accuracy = (sum(map(int, correct)) / float(len(correct)))
    print('accuracy = {0}%'.format(accuracy * 100))

输出如下：

[[10]
 [10]
 [10]
 ...
 [ 9]
 [ 9]
 [ 9]]
accuracy = 93.10000000000001%