机器学习——神经网络

神经网络(neural-networks)

准备数据集

神经网络一般用于比较复杂的问题，而且计算量较大，这次处理的数据集是长和宽都为20的图像(用一个长度为400的向量存储)，首先看一下数据集长什么样：

def plot_100_image(X):
    """
    :param X: (5000,400)
    """
    size = int(np.sqrt(X.shape[1]))  # 20
    sample_index = np.random.choice(np.arange(X.shape[0]), 100)  # 100 * 400, 随机选取100个样本
    sample_images = X[sample_index, :]
    fig, ax_array = plt.subplots(nrows=10, ncols=10, sharey='all', sharex='all', figsize=(8, 8))  # 10*10的图像, 共享y轴和x轴
    for r in range(10):
        for c in range(10):
            # 重塑为20*20的图像
            ax_array[r, c].matshow(sample_images[10 * r + c].reshape((size, size)), cmap=matplotlib.cm.binary)
            plt.xticks(np.array([]))  # 去掉x坐标轴
            plt.yticks(np.array([]))  # 去掉y坐标轴

这里展现了100个数据集，可以看出是手写数字的图像集，我们的模型要做的就是根据训练集来训练一个手写体数字识别的模型，我们进一步将问题进行转化：输入一个20*20像素的手写数字图像来进行0-9的多分类。

前向传播

我们先根据已有的模型参数来进行前向传播看看分类效果

首先我们需要了解到神经元的网络架构，每一层的输出变量都是下一层的输入变量。下图为一个3层的神经网络，第一层成为输入层（Input Layer），最后一层称为输出层（Output Layer），中间一层成为隐藏层（Hidden Layers）。我们为每一层都增加一个偏差单位（Bias unit）：

先说明一些符号描述：

当然我们可以使用向量的手段使得公式更加简洁：

注: g表示sigmoid激活函数，这样公式就相当简化了

下面给出前向传播代码：

# 前向传播
def feed_forward(theta, X):
    """
    :param theta: 两个权重矩阵扁平化合并后的向量，t1(25, 401), t2(10, 26)
    :param X: 5000 * 401
    """
    t1 = theta[:25 * 401].reshape(25, 401)
    t2 = theta[25 * 401:].reshape(10, 26)
    m = X.shape[0]  # 5000
    a1 = X  # 5000 * 401
    z2 = a1 @ t1.T  # 5000 * 25
    a2 = np.insert(sigmoid(z2), 0, values=np.ones(m), axis=1)  # 5000 * 26
    z3 = a2 @ t2.T  # 5000 * 10
    h = sigmoid(z3)  # 5000 * 10, 最终的输出
    return a1, z2, a2, z3, h	# 之后的反向传播也会用到这些量

代价函数

再进行反向传播前，我们先完成代价函数，此问题是多分类问题，不考虑正则化的代价函数如下(不考虑数据集的数量m的维度)：
$$J(\theta )=\frac{1}{m}\sum _{i=1}^m{y}^{(i)}\ast log(h_{\theta })-(1-y^{(i)})\ast log(1-h_{\theta })$$

# 代价函数
def cost(theta, X, y):
    """
    :param theta: 两个权重矩阵扁平化合并后的向量，t1(25, 401), t2(10, 26)
    :param X: 5000 * 401
    :param y: 5000 * 10
    """
    m = X.shape[0]  # 5000
    _, _, _, _, h = feed_forward(theta, X)  # 5000 * 10
    pair_computation = -y * np.log(h) - (1 - y) * np.log(1 - h)  # 5000 * 10
    return np.sum(pair_computation) / m

正则化代价函数

$$\begin{gathered} {\theta }^{(i)}表示第i层到第i+1层的权重,注意不惩罚偏置项,即\theta 的下标应该从1开始取 \\ 正则化项：\frac{\lambda }{2m}\sum _{i=1}^K({\theta }^{(i)}{)}^2 \end{gathered}$$

# 正则化代价函数
def regularized_cost(theta, X, y, l=1):
    """
    :param theta: 两个权重矩阵，t1(25, 401), t2(10, 26)
    :param X: 5000 * 401
    :param y: 5000 * 10
    :param l: lambda
    """
    m = X.shape[0]  # 5000
    t1 = theta[:25 * 401].reshape(25, 401)
    t2 = theta[25 * 401:].reshape(10, 26)
    # 不惩罚偏置项, 索引从1开始
    regularized_term = (l / (2 * m)) * (np.sum(np.power(t1[:, 1:], 2)) + np.sum(np.power(t2[:, 1:], 2)))
    return cost(theta, X, y) + regularized_term

反向传播

sigmoid梯度函数

$$g^{\prime }(z)=\frac{\partial }{\partial z}g(z)=g(z)(1-g(z))$$

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

# sigmoid梯度下降函数
def sigmoid_gradient(z):
    return np.multiply(sigmoid(z), (1 - sigmoid(z)))

θ的梯度

现在准备实施反向传播来计算梯度，由于反向传播所需的计算是代价函数中所需的计算过程，我们实际上将扩展代价函数以执行反向传播并返回代价和梯度，计算公式如下：
$$\begin{gathered} {\delta }^{(i)}表示第i层的误差 \\ {\delta }^{\left(3\right)}=a^{(3)}-y \\ {\delta }^{\left(2\right)}={\left({\Theta }^{(2)}\right)}^T{\delta }^{(3)}.\ast g^{\prime }(z^{(2)}) \\ {\Delta }^{\left(l\right)}={\Delta }^{(l)}+{\delta }^{(l+1)}{\left(a^{(l)}\right)}^T \\ \frac{\partial }{\partial {\Theta }_{ij}^{(l)}}J(\Theta )=\frac{1}{m}{\Delta }_{ij}^l \end{gathered}$$

# 梯度下降
def gradient(theta, X, y):
    m = X.shape[0]  # 5000
    t1 = theta[:25 * 401].reshape(25, 401)  # 25 * 401
    t2 = theta[25 * 401:].reshape(10, 26)  # 10 * 26

    delta1 = np.zeros(t1.shape)  # 25 * 401
    delta2 = np.zeros(t2.shape)  # 10 * 26

    # 前向传播
    a1, z2, a2, z3, h = feed_forward(theta, X)

    for i in range(m):
        a1_i = a1[i, :]  # 1 * 401，第i个样本的输入层
        z2_i = z2[i, :]  # 1 * 25，第i个样本的隐藏层
        a2_i = a2[i, :]  # 1 * 26，第i个样本的隐藏层
        h_i = h[i, :]  # 1 * 10，第i个样本的输出层
        y_i = y[i, :]  # 1 * 10，第i个样本的真实值
        # 第三层
        d3_i = h_i - y_i  # 1 * 10，第i个样本的输出层误差
        # 第二层
        d2_i = np.multiply(d3_i @ t2, sigmoid_gradient(np.insert(z2_i, 0, values=np.ones(1))))  # 1 * 26
        # 累加
        delta2 += np.matrix(d3_i).T @ np.matrix(a2_i)  # (1*10).T @ (1*26) = 10*26
        delta1 += np.matrix(d2_i[1:]).T @ np.matrix(a1_i)  # (1*25).T @ (1*401) = 25*401

    delta1 /= m
    delta2 /= m

    return np.concatenate((np.ravel(delta1), np.ravel(delta2)))  # 10285的向量，扁平化处理合并后的权重向量

正则化的梯度

加入正则化项后的梯度函数:
$$\begin{array}{rl}\frac{\partial }{\partial {\Theta }_{ij}^{(l)}}J(\Theta )& =\frac{1}{m}{\Delta }_{ij}^l\text{ , for }j=0\\ \frac{\partial }{\partial {\Theta }_{ij}^{(l)}}J(\Theta )& =\frac{1}{m}{\Delta }_{ij}^l+\frac{\lambda }{m}{\Theta }_{ij}^{(l)}\text{ , for }j\ge 1\end{array}$$

# 正则化梯度下降
def regularized_gradient(theta, X, y, l=1):
    m = X.shape[0]  # 5000
    t1 = theta[:25 * 401].reshape(25, 401)  # 25 * 401
    t2 = theta[25 * 401:].reshape(10, 26)  # 10 * 26

    delta = gradient(theta, X, y)
    delta1 = delta[:25 * 401].reshape(25, 401)
    delta2 = delta[25 * 401:].reshape(10, 26)

    t1[:, 0] = 0  # 不惩罚偏置项
    reg_term_d1 = (l / m) * t1  # 第一个参数的正则化项
    delta1 += reg_term_d1

    t2[:, 0] = 0  # 不惩罚偏置项
    reg_term_d2 = (l / m) * t2  # 第二个参数的正则化项
    delta2 += reg_term_d2

    return np.concatenate((np.ravel(delta1), np.ravel(delta2)))  # 10285

训练模型

首先神经网络的初始化参数不能全为0，我们需要随机化处理参数来打破对称性

import scipy.optimize as opt

# 随机初始化
def random_init(size):
    return np.random.uniform(-0.12, 0.12, size)  # uniform均匀分布

# 训练
def train(X, y):
    init_theta = random_init(10285)  # 10285 = 25 * 401 + 10 * 26
    res = opt.minimize(fun=regularized_cost,  # 代价函数
                       x0=init_theta,  # 参数
                       args=(X, y, 5),  # 代价函数的参数
                       method='TNC',  # 梯度下降算法
                       jac=regularized_gradient,  # 梯度函数
                       options={'maxfun': 400}  # 最大迭代次数
                       )
    # res.x为最终参数, 10285 = 25 * 401 + 10 * 26
    return res

我们来看看训练结果：

# 显示准确率
from sklearn.metrics import classification_report

def show_accuracy(theta, X, y):
    _, _, _, _, h = feed_forward(theta, X)
    y_pred = np.argmax(h, axis=1) + 1  # 5000 * 1, 每一行最大值的索引
    print(classification_report(y, y_pred))
    correct = [1 if a == b else 0 for (a, b) in zip(y_pred, y)]
    accuracy = sum(correct) / len(correct)
    print('total accuracy = {0}%'.format(accuracy * 100))

主函数相关代码：

# 梯度下降
    res = train(X, y)
    print('训练完成')
    print(res)  # res.x为最终参数, 10285 = 25 * 401 + 10 * 26
    # 训练后评价
    print('训练后评价为')
    show_accuracy(res.x, X, y_raw)

看看最终的训练效果：

参考链接: https://github.com/loveunk/machine-learning-excercise-notebook-python/tree/master/4.nurual_network_back_propagation