Упражнение по машинному обучению 4: Реализация нейронных сетей BP в Python

Введение данных: в ex4data1 есть 5000 обучающих примеров. где каждый обучающий экземпляр представляет собой изображение в градациях серого размером 20 на 20 пикселей. Каждый пиксель представлен числом с плавающей запятой, представляющим интенсивность оттенков серого в этом месте. Эта сетка размером 20x20 пикселей «развёрнута» в 400-мерный вектор. Каждый из этих обучающих примеров становится одной строкой в матрице данных X. В результате получается матрица X размером 5000×400, где каждая строка является обучающим примером изображения рукописной цифры.

Набор данных находится в собственном формате 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 Горячее кодирование

Чтобы обучить нейронную сеть, нам нужно перекодировать метки в векторы, содержащие только 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 Развертывание параметров

Чтобы использовать функцию оптимизации, такую как «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 Упреждающая операция

Иллюстрация модели:

Сигмовидная функция:

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 Функция потерь

Обыкновенная функция потерь

Две суммы суммируют потери логистической регрессии для каждой единицы в выходном слое, теперь мы имеем

и $h_{\theta} \in R^{5000 \times 10}$ .

Без учета размерностей m и K уравнение $= -y*log(h_{\theta}) - (1-y)*log(1-h_{\theta})$ , и, наконец, просуммируйте полученные двумерные массивы и разделите на 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)$ Производная от , которая затем используется при оптимизации алгоритма градиентного спуска. Для получения информации о градиентном спуске для алгоритмов линейной регрессии и логистической регрессии см. мою предыдущую статью:Упражнение по машинному обучению 2. Реализация логистической регрессии в Python.

2.1 Сигмовидный градиент (градиентная сигмовидная функция, то есть производная от сигмовидной функции)

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

2.2 Алгоритм БП

Значение нижних индексов:

Указывает, какой слой вычисляется в данный момент,Он представляет количество рассчитанных в данный момент выборок,Представляет количество узлов в этом слое, которые вычисляются в данный момент.

Иллюстрация модели:

Псевдокод может быть выражен как:

1. Набор данных:

2. Настройки $Δ^{l}_{ij}=0$

3.for-loop:

Математическая формула включала:

$Δ^{(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

результат операции:

Видно, что потери уменьшаются до 0,305 после 1000 итераций.

3.3 Точность расчета

Для объяснения функции точности см.:Упражнение по машинному обучению 2. Реализация логистической регрессии в 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()

результат: