从头开始实现多层人工神经网络

从头开始实现多层人工神经网络

from IPython.display import Image
%matplotlib inline

1.使用人工神经网络建模复杂函数

神经网络中的Epoch、Iteration、Batchsize

神经网络中epoch与iteration是不相等的

• batchsize：中文翻译为批大小（批尺寸）。在深度学习中，一般采用SGD训练，即每次训练在训练集中取batchsize个样本训练；

• iteration：中文翻译为迭代，1个iteration等于使用batchsize个样本训练一次；一个迭代 = 一个正向通过+一个反向通过

• epoch：迭代次数，1个epoch等于使用训练集中的全部样本训练一次；一个epoch = 所有训练样本的一个正向传递和一个反向传递

• 举个例子，训练集有1000个样本，batchsize=10，那么：训练完整个样本集需要：100次iteration，1次epoch。

1.1单隐层神经网络综述

回顾前面的自适应神经元算法。在

前面，通过梯度下降优化算法来学习模型的权重系数。在每个epoch中(传递训练数据集)，权重更新过程如下：

$$\mathbf{w}:=\mathbf{w}+\Delta \mathbf{w},\text{ where }\Delta \mathbf{w}=-\eta \nabla J(\mathbf{w})$$

这里实际上使用的是批量梯度下降，即计算的梯度是基于整个训练集，同时基于梯度的负方向更新模型的权值。其中，定义目标

函数为误差平方和SSE，并记作$J(W)$。更进一步，通过将$-\nabla J(W)$乘上学习率$\eta$,用以控制

下降步伐，从而避免越过了代价函数的全局最小值。

基于上述优化方式，我们同时更新所有的权重系数，定义每个权重的偏导数如下：
$$\frac{\partial }{\partial w_j}J(\mathbf{w})=-\sum _i\left(y^{(i)}-a^{(i)}\right)x_j^{(i)}$$

这里的$y^{(i)}$代表的是特定样本$x^{(i)}$的类别标签，$a^{(i)}$代表的是神经元的激活函数，在自适应神经元中是一个线性

函数。形式如下：

$$\varphi (z)=z=a$$

其中，上式的$z$为连接输入层和输出层的权值线性组合：
$$z=\sum _j{w}_j{x}_j={\mathbf{w}}^T\mathbf{x}$$

使用上述激活函数计算梯度更新，进一步实现一个阈值函数，将连续值输出压缩成二进制类别标签：
$$\hat{y}=\{\begin{array}{rl}1& \text{ if }g(z)\ge 0\\ -1& \text{ otherwise }\end{array}$$

Image(filename='images/12_01.png', width=600) 

如上图，尽管Adaline由两个层组成，一个输入一个输出。但是由于输入层和输出增之间仅仅有一条链路，依然称其为单层网络。

另外一种加速模型学习的优化方式为随机梯度下降**Stochastic gradient descent SGD**，SGD 近似于单个训练样本(online learning),或者近似于使用一小部分训练样本（小批量学习）。

SGD相较于批量梯度下降，由于权重更新更加频繁，因此学习速度更快。同时噪声特性也使得其在训练具有非线性激活函数的

多层神经网络(其不具有凸代价函数)的时候，也很有益。

这里引入的噪声可以促进优化目标 避免陷入局部最小。

1.2介绍多层神经网络结构

将多个 单个的神经元连接到多层前馈神经网络，这种特殊类型的全连接网络也称为MLP。

# 三个层次的MLP几何表示
Image(filename='images/12_02.png', width=600) 

标记：第$L$层的第$i$个激活单元为$a_i^{(l)}$,则上图输入层的激活单元加上偏置单元后，表示如下：

$${\mathbf{a}}^{(in)}=\left[\begin{array}{c}{a}_0^{(in)}\\ a_1^{(in)}\\ ⋮\\ a_m^{(in)}\end{array}\right]=\left[\begin{array}{c}1\\ x_1^{(in)}\\ ⋮\\ x_m^{(in)}\end{array}\right]$$

对上述的网络结构，如果具有一个以上的隐含层，则称之为深层人工神经网络。

Image(filename='images/12_03.png', width=500) 

1.3通过前向传播激活神经网络

MLP学习过程总结：

1.从输入层开始，通过网络前向传播训练数据的模式，从而生成输出；

2.基于网络的输出，使用代价函数来计算我们想要最小化的误差；

3.对误差进行反向传播，求出误差相对于网络中每个权重的导数，同时对模型进行更新；

最后，在对多个epoch重复这三个步骤并学习MLP的权重之后，使用前向传播来计算网络输出，并应用阈值函数来获得在独热编码表示中预测的类别标签。

由于隐含层中的每个单元都与输入层中的所有单元相连，这里先计算隐含层的激活单元$a_1^{(h)}$:

$$\begin{array}{l}{z}_1^{(h)}=a_0^{(in)}{w}_{0,1}^{(h)}+a_1^{(in)}{w}_{1,1}^{(h)}+\cdots +a_m^{(in)}{w}_{m,1}^{(h)}\\ a_1^{(h)}=\varphi \left(z_1^{(h)}\right)\end{array}$$

其中，$z_1^{(h)}$是网络输入，$\varphi (\cdot )$是激活函数，其必须可微，这样才可以基于梯度的学习方法学习连接神经元的权重。

为了能够解决图像分类这样的复杂问题，MLP中需要使用非线性激活函数，例如sigmoid函数：

$$\varphi (z)=\frac{1}{1+e^{-z}}$$

# sigmoid函数几何表示如下
Image(filename='images/12_04.png', width=500) 

MLP是前馈人工神经网络的一个典型例子。前馈代表的是：每一层都作为下一层的输入，而没有循环，这与循环神经网络RNN不同。

multilayer perceptron多层感知机实际上有一定的混淆视听，其网络架构中的人工神经元是典型的sigmoid单元，而不是感知机。

我们可以把MLP中的神经元看做是逻辑回归单元，返回值在0到1之间的连续范围内。

为了表示简化，上式简写为：

$$\begin{array}{l}{\mathbf{z}}^{(h)}={\mathbf{a}}^{(in)}{\mathbf{W}}^{(h)}\\ {\mathbf{a}}^{(h)}=\varphi \left({\mathbf{z}}^{(h)}\right)\end{array}$$

这里的${\mathbf{a}}^{(in)}$为$1\times m$的特征向量，还加上了一个偏置单元，样本为${\mathbf{x}}^{in}$。

${\mathbf{W}}^{(h)}$为$m\times d$维的权重矩阵，其中$d$为隐含层神经元的个数。

经过矩阵-向量乘法之后，得到$1\times d$维的网络净输入(净活性值)${\mathbf{z}}^{(h)}$,以此来计算活性值${\mathbf{a}}^{(h)}$,其中${\mathbf{a}}^{(h)}\in {\mathbb{R}}^{1\times d}$。

净输入在经过非线性激活函数之后，得到网络的活性值。净输入为输入信息$x$的加权和，非线性函数也叫作激活函数。

更进一步，将上述计算推广到所有的样本上：
$${\mathbf{Z}}^{(h)}={\mathbf{A}}^{(in)}{\mathbf{W}}^{(h)}$$

这里的${\mathbf{A}}^{(in)}$是一个矩阵，形状为$n\times m$,则在矩阵乘法之后，得到网络的输入矩阵${\mathbf{Z}}^{(h)}$,形状为$n\times d$。

然后通过激活函数的计算得到网络的活性值，这里的输出层：
$${\mathbf{A}}^{(h)}=\varphi \left({\mathbf{Z}}^{(h)}\right)$$

类似的，将输出层的活性值写成向量形式如下：
$${\mathbf{Z}}^{(\text{out })}={\mathbf{A}}^{(h)}{\mathbf{W}}^{(\text{out })}$$

这里乘上$d\times t$维度的矩阵${\mathbf{W}}^{(\text{out })}$,$t$是输出单元的个数，${\mathbf{A}}^{(h)}$size为$n\times d$，最终得到$n\times t$维的${\mathbf{Z}}^{(\text{out })}$。

最后，添加上激活函数，得到网络的输出：
$${\mathbf{A}}^{(out)}=\varphi \left({\mathbf{Z}}^{(out)}\right),{\mathbf{A}}^{(out)}\in {\mathbb{R}}^{n\times t}$$

如果想通过sklearn API获取MNIST数据集，执行下述代码：

"""
from sklearn.datasets import fetch_openml
from sklearn.model_selection import train_test_split


X, y = fetch_openml('mnist_784', version=1, return_X_y=True)
y = y.astype(int)
X = ((X / 255.) - .5) * 2
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=10000, random_state=123, stratify=y)
"""

"\nfrom sklearn.datasets import fetch_openml\nfrom sklearn.model_selection import train_test_split\n\n\nX, y = fetch_openml('mnist_784', version=1, return_X_y=True)\ny = y.astype(int)\nX = ((X / 255.) - .5) * 2\nX_train, X_test, y_train, y_test = train_test_split(\n    X, y, test_size=10000, random_state=123, stratify=y)\n"

2.手写数字分类

…

2.1获取并准备MNIST数据集

MNIST数据集获取地址 ，其由以下几部分组成：

• Training set images: train-images-idx3-ubyte.gz (9.9 MB, 47 MB unzipped, 60,000 examples)
• Training set labels: train-labels-idx1-ubyte.gz (29 KB, 60 KB unzipped, 60,000 labels)
• Test set images: t10k-images-idx3-ubyte.gz (1.6 MB, 7.8 MB, 10,000 examples)
• Test set labels: t10k-labels-idx1-ubyte.gz (5 KB, 10 KB unzipped, 10,000 labels)

#  mnist数据解压代码

import sys
import gzip
import shutil
import os

if (sys.version_info > (3, 0)):
    writemode = 'wb'
else:
    writemode = 'w'

zipped_mnist = [f for f in os.listdir() if f.endswith('ubyte.gz')]
for z in zipped_mnist:
    with gzip.GzipFile(z, mode='rb') as decompressed, open(z[:-3], writemode) as outfile:
        outfile.write(decompressed.read()) 

---

如果没有使用gzip格式，需要确保文件命名形式如下：

• train-images-idx3-ubyte
• train-labels-idx1-ubyte
• t10k-images-idx3-ubyte
• t10k-labels-idx1-ubyte

---

# 图像数据以字节码的形式进行存储，将其读入为Numpy数组形式

import os
import struct
import numpy as np
 
def load_mnist(path, kind='train'):
    """Load MNIST data from `path`"""
    labels_path = os.path.join(path, 
                               '%s-labels-idx1-ubyte' % kind)
    images_path = os.path.join(path, 
                               '%s-images-idx3-ubyte' % kind)
        
    with open(labels_path, 'rb') as lbpath:
        magic, n = struct.unpack('>II', 
                                 lbpath.read(8))
        labels = np.fromfile(lbpath, 
                             dtype=np.uint8)

    with open(images_path, 'rb') as imgpath:
        magic, num, rows, cols = struct.unpack(">IIII", 
                                               imgpath.read(16))
        images = np.fromfile(imgpath, 
                             dtype=np.uint8).reshape(len(labels), 784)
        images = ((images / 255.) - .5) * 2  # 原始像素值在0-255之间，归一化之后缩放到正负一之间；
 
    return images, labels

# load_mnist函数返回两个数组，第一个为nxm维的NumPy数组(图像)，其中n为样本个数，m为特征个数。
# 训练数据集含60000个训练样本，测试数据含10000个样本。

# MNIST数据集中的图像为28x28维的像素点，每个像素点代表一个灰度强度值。
# 这里将28x28的图像拉伸为1维的行向量，其代表输入图像矩阵的行，对应784/行,或者/图像
# 上述函数返回的第二个数组为目标变量，即类别标签(这里为整数值0-9)

X_train, y_train = load_mnist('', kind='train')
print('Rows: %d, columns: %d' % (X_train.shape[0], X_train.shape[1]))

Rows: 60000, columns: 784

X_test, y_test = load_mnist('', kind='t10k')
print('Rows: %d, columns: %d' % (X_test.shape[0], X_test.shape[1]))

Rows: 10000, columns: 784

# 基于梯度的优化中，使用批量归一化是实现更好收敛性的一个技巧  Batch normalization

这里先通过imshow函数将特征矩阵从784像素的向量reshape成原始的28x28维度的图像，从而实现可视化

import matplotlib.pyplot as plt

fig, ax = plt.subplots(nrows=2, ncols=5, sharex=True, sharey=True)
ax = ax.flatten()

for i in range(10):
    img = X_train[y_train == i][0].reshape(28, 28)
    ax[i].imshow(img, cmap='Greys')

ax[0].set_xticks([])
ax[0].set_yticks([])
plt.tight_layout()
# plt.savefig('images/12_5.png', dpi=300)
plt.show()

看看同一个数字的不同写法

# 数字7的25中不同变体
fig, ax = plt.subplots(nrows=5, ncols=5, sharex=True, sharey=True,)
ax = ax.flatten()
for i in range(25):
    img = X_train[y_train == 7][i].reshape(28, 28)
    ax[i].imshow(img, cmap='Greys')

ax[0].set_xticks([])
ax[0].set_yticks([])
plt.tight_layout()
# plt.savefig('images/12_6.png', dpi=300)
plt.show()

在完成上述处理步骤之后，这里将处理过后的数据进行保存。通过使用NumPy的savez函数实现将多维数组保存到磁盘。

savez函数类似于python的pickle模块，其将创建压缩后的数据存档，生成包含.npy格式文件的.npz文件；此外，这里使用savez_compressed而不是

savez函数，因为它使用与savez相同的语法，但可以进一步将输出文件压缩到更小的文件体积。(从原本的400M到压缩之后的20M)。

# 将训练集和测试集保存到归档文件
import numpy as np

np.savez_compressed('mnist_scaled.npz', 
                    X_train=X_train,
                    y_train=y_train,
                    X_test=X_test,
                    y_test=y_test)

# 创建了压缩文件之后，可以使用Numpy的load方法进行加载
mnist = np.load('mnist_scaled.npz')
# mnist变量现在指向一个对象，该对象可以访问savez_compressed函数作为关键字参数提供的四个数据数组。通过mnist.files属性访问。
mnist.files

['X_train', 'y_train', 'X_test', 'y_test']

# 例如
# X_train = mnist['X_train']

# 使用列表生成式，可以检索所有四个数据数组
X_train, y_train, X_test, y_test = [mnist[f] for f in ['X_train', 'y_train', 
                                    'X_test', 'y_test']]

del mnist

X_train.shape

(60000, 784)

2.2实现多层感知机

import numpy as np
import sys


class NeuralNetMLP(object):
    """ Feedforward neural network / Multi-layer perceptron classifier.

    Parameters
    ------------
    n_hidden : int (default: 30)
        Number of hidden units.
    l2 : float (default: 0.)
        Lambda value for L2-regularization.
        No regularization if l2=0. (default)
    epochs : int (default: 100)
        Number of passes over the training set.
    eta : float (default: 0.001)
        Learning rate.
    shuffle : bool (default: True)
        Shuffles training data every epoch if True to prevent circles.
    minibatch_size : int (default: 1)
        Number of training examples per minibatch.
    seed : int (default: None)
        Random seed for initializing weights and shuffling.

    Attributes
    -----------
    eval_ : dict
      Dictionary collecting the cost, training accuracy,
      and validation accuracy for each epoch during training.

    """
    def __init__(self, n_hidden=30,
                 l2=0., epochs=100, eta=0.001,
                 shuffle=True, minibatch_size=1, seed=None):

        self.random = np.random.RandomState(seed)
        self.n_hidden = n_hidden
        self.l2 = l2
        self.epochs = epochs
        self.eta = eta
        self.shuffle = shuffle
        self.minibatch_size = minibatch_size

    def _onehot(self, y, n_classes):
        """Encode labels into one-hot representation

        Parameters
        ------------
        y : array, shape = [n_examples]
            Target values.
        n_classes : int
            Number of classes

        Returns
        -----------
        onehot : array, shape = (n_examples, n_labels)

        """
        onehot = np.zeros((n_classes, y.shape[0]))
        for idx, val in enumerate(y.astype(int)):
            onehot[val, idx] = 1.
        return onehot.T

    def _sigmoid(self, z):
        """Compute logistic function (sigmoid)"""
        return 1. / (1. + np.exp(-np.clip(z, -250, 250)))

    def _forward(self, X):
        """Compute forward propagation step"""

        # step 1: net input of hidden layer
        # [n_examples, n_features] dot [n_features, n_hidden]
        # -> [n_examples, n_hidden]
        z_h = np.dot(X, self.w_h) + self.b_h

        # step 2: activation of hidden layer
        a_h = self._sigmoid(z_h)

        # step 3: net input of output layer
        # [n_examples, n_hidden] dot [n_hidden, n_classlabels]
        # -> [n_examples, n_classlabels]

        z_out = np.dot(a_h, self.w_out) + self.b_out

        # step 4: activation output layer
        a_out = self._sigmoid(z_out)

        return z_h, a_h, z_out, a_out

    def _compute_cost(self, y_enc, output):
        """Compute cost function.

        Parameters
        ----------
        y_enc : array, shape = (n_examples, n_labels)
            one-hot encoded class labels.
        output : array, shape = [n_examples, n_output_units]
            Activation of the output layer (forward propagation)

        Returns
        ---------
        cost : float
            Regularized cost

        """
        L2_term = (self.l2 *
                   (np.sum(self.w_h ** 2.) +
                    np.sum(self.w_out ** 2.)))

        term1 = -y_enc * (np.log(output))
        term2 = (1. - y_enc) * np.log(1. - output)
        cost = np.sum(term1 - term2) + L2_term
        
        # If you are applying this cost function to other
        # datasets where activation
        # values maybe become more extreme (closer to zero or 1)
        # you may encounter "ZeroDivisionError"s due to numerical
        # instabilities in Python & NumPy for the current implementation.
        # I.e., the code tries to evaluate log(0), which is undefined.
        # To address this issue, you could add a small constant to the
        # activation values that are passed to the log function.
        #
        # For example:
        #
        # term1 = -y_enc * (np.log(output + 1e-5))
        # term2 = (1. - y_enc) * np.log(1. - output + 1e-5)
        
        return cost

    def predict(self, X):
        """Predict class labels

        Parameters
        -----------
        X : array, shape = [n_examples, n_features]
            Input layer with original features.

        Returns:
        ----------
        y_pred : array, shape = [n_examples]
            Predicted class labels.

        """
        z_h, a_h, z_out, a_out = self._forward(X)
        y_pred = np.argmax(z_out, axis=1)
        return y_pred

    def fit(self, X_train, y_train, X_valid, y_valid):
        """ Learn weights from training data.

        Parameters
        -----------
        X_train : array, shape = [n_examples, n_features]
            Input layer with original features.
        y_train : array, shape = [n_examples]
            Target class labels.
        X_valid : array, shape = [n_examples, n_features]
            Sample features for validation during training
        y_valid : array, shape = [n_examples]
            Sample labels for validation during training

        Returns:
        ----------
        self

        """
        n_output = np.unique(y_train).shape[0]  # number of class labels
        n_features = X_train.shape[1]

        ########################
        # Weight initialization
        ########################

        # weights for input -> hidden
        self.b_h = np.zeros(self.n_hidden)
        self.w_h = self.random.normal(loc=0.0, scale=0.1,
                                      size=(n_features, self.n_hidden))

        # weights for hidden -> output
        self.b_out = np.zeros(n_output)
        self.w_out = self.random.normal(loc=0.0, scale=0.1,
                                        size=(self.n_hidden, n_output))

        epoch_strlen = len(str(self.epochs))  # for progress formatting
        self.eval_ = {'cost': [], 'train_acc': [], 'valid_acc': []}

        y_train_enc = self._onehot(y_train, n_output)

        # iterate over training epochs
        for i in range(self.epochs):

            # iterate over minibatches
            indices = np.arange(X_train.shape[0])

            if self.shuffle:
                self.random.shuffle(indices)

            for start_idx in range(0, indices.shape[0] - self.minibatch_size +
                                   1, self.minibatch_size):
                batch_idx = indices[start_idx:start_idx + self.minibatch_size]

                # forward propagation
                z_h, a_h, z_out, a_out = self._forward(X_train[batch_idx])

                ##################
                # Backpropagation
                ##################

                # [n_examples, n_classlabels]
                delta_out = a_out - y_train_enc[batch_idx]

                # [n_examples, n_hidden]
                sigmoid_derivative_h = a_h * (1. - a_h)

                # [n_examples, n_classlabels] dot [n_classlabels, n_hidden]
                # -> [n_examples, n_hidden]
                delta_h = (np.dot(delta_out, self.w_out.T) *
                           sigmoid_derivative_h)

                # [n_features, n_examples] dot [n_examples, n_hidden]
                # -> [n_features, n_hidden]
                grad_w_h = np.dot(X_train[batch_idx].T, delta_h)
                grad_b_h = np.sum(delta_h, axis=0)

                # [n_hidden, n_examples] dot [n_examples, n_classlabels]
                # -> [n_hidden, n_classlabels]
                grad_w_out = np.dot(a_h.T, delta_out)
                grad_b_out = np.sum(delta_out, axis=0)

                # Regularization and weight updates
                delta_w_h = (grad_w_h + self.l2*self.w_h)
                delta_b_h = grad_b_h # bias is not regularized
                self.w_h -= self.eta * delta_w_h
                self.b_h -= self.eta * delta_b_h

                delta_w_out = (grad_w_out + self.l2*self.w_out)
                delta_b_out = grad_b_out  # bias is not regularized
                self.w_out -= self.eta * delta_w_out
                self.b_out -= self.eta * delta_b_out

            #############
            # Evaluation
            #############

            # Evaluation after each epoch during training
            z_h, a_h, z_out, a_out = self._forward(X_train)
            
            cost = self._compute_cost(y_enc=y_train_enc,
                                      output=a_out)

            y_train_pred = self.predict(X_train)
            y_valid_pred = self.predict(X_valid)

            train_acc = ((np.sum(y_train == y_train_pred)).astype(np.float) /
                         X_train.shape[0])
            valid_acc = ((np.sum(y_valid == y_valid_pred)).astype(np.float) /
                         X_valid.shape[0])

            sys.stderr.write('\r%0*d/%d | Cost: %.2f '
                             '| Train/Valid Acc.: %.2f%%/%.2f%% ' %
                             (epoch_strlen, i+1, self.epochs, cost,
                              train_acc*100, valid_acc*100))
            sys.stderr.flush()

            self.eval_['cost'].append(cost)
            self.eval_['train_acc'].append(train_acc)
            self.eval_['valid_acc'].append(valid_acc)

        return self

n_epochs = 200

## @Readers: PLEASE IGNORE IF-STATEMENT BELOW
##
## This cell is meant to run fewer epochs when
## the notebook is run on the Travis Continuous Integration
## platform to test the code on a smaller dataset
## to prevent timeout errors; it just serves a debugging tool

if 'TRAVIS' in os.environ:
    n_epochs = 20

# 执行代码之后，初始化了一个784-100-10的MLP。其具有784个输入单元，对应n_features为784,100个隐层单元(n_hidden),以及10个输出单元

使用55000个样本训练MLP，这些样本来自于已经打乱的MNIST训练数据集，并且在训练过程中使用剩下的5000个样本进行验证。

nn = NeuralNetMLP(n_hidden=100, 
                  l2=0.01, 
                  epochs=n_epochs, 
                  eta=0.0005,
                  minibatch_size=100, 
                  shuffle=True,
                  seed=1)

nn.fit(X_train=X_train[:55000], 
       y_train=y_train[:55000],
       X_valid=X_train[55000:],
       y_valid=y_train[55000:])

200/200 | Cost: 5065.78 | Train/Valid Acc.: 99.28%/97.98%  




<__main__.NeuralNetMLP at 0x1c20891a988>

import matplotlib.pyplot as plt


plt.plot(range(nn.epochs), nn.eval_['cost'])
plt.ylabel('Cost')
plt.xlabel('Epochs')
#plt.savefig('images/12_07.png', dpi=300)
plt.show()

可以看到，代价在前100个epoch期间大幅度下降，且貌似在100个epoch时缓慢收敛。但在175到200个epoch之间，下降缓慢，但依然呈现出下降趋势

# 可视化训练和验证的Accuracy变化情况
plt.plot(range(nn.epochs), nn.eval_['train_acc'], 
         label='Training')
plt.plot(range(nn.epochs), nn.eval_['valid_acc'], 
         label='Validation', linestyle='--')
plt.ylabel('Accuracy')
plt.xlabel('Epochs')
plt.legend(loc='lower right')
plt.savefig('images/12_08.png', dpi=300)
plt.show()

上图显示出，当epoch为50的时候，训练和验证accuracy相等，并在此之后，网络出现了过拟合情况。降低过拟合的有效措施可以是增加正则化强度，也可以使用丢弃法，即随机丢弃一部分神经元。and so on

# 通过计算模型在测试数据上的预测精度来评估模型的泛化性能
y_test_pred = nn.predict(X_test)
acc = (np.sum(y_test == y_test_pred)
       .astype(np.float) / X_test.shape[0])

print('Test accuracy: %.2f%%' % (acc * 100))

Test accuracy: 97.54%

尽管模型在训练数据集上具有轻微的过拟合情况，但相较于验证数据集的精度97.98%，模型依然在测试集上取得了不错的性能。

为了进一步微调模型，我们可以改变隐藏单元的数量，正则化参数的值，以及学习速率等；

自适应学习率、更加复杂的SGD、批量归一化、丢弃法

miscl_img = X_test[y_test != y_test_pred][:25]
correct_lab = y_test[y_test != y_test_pred][:25]
miscl_lab = y_test_pred[y_test != y_test_pred][:25]

fig, ax = plt.subplots(nrows=5, ncols=5, sharex=True, sharey=True)
ax = ax.flatten()
for i in range(25):
    img = miscl_img[i].reshape(28, 28)
    ax[i].imshow(img, cmap='Greys', interpolation='nearest')
    ax[i].set_title('%d) t: %d p: %d' % (i+1, correct_lab[i], miscl_lab[i]))

ax[0].set_xticks([])
ax[0].set_yticks([])
plt.tight_layout()
#plt.savefig('images/12_09.png', dpi=300)
plt.show()

3.训练一个人工神经网络

3.1计算Logistic损失函数

Logistic损失函数定义如下：
$$J(\mathbf{w})=-\sum _{i=1}^n{y}^{[i]}\log \left(a^{[i]}\right)+\left(1-y^{[i]}\right)\log \left(1-a^{[i]}\right)$$

这里的$a^{[i]}$是在前向传播过程中，计算的数据集中第$i$个样本的sigmoid激活：上标$[i]$代表的是样本索引。

$$a^{[i]}=\varphi \left(z^{[i]}\right)$$

通过在上述损失函数上增加一项损失函数，来降低模型过拟合程度，如L2正则化：
$$L2=\lambda \Vert w{\Vert }_2^2=\lambda \sum _{j=1}^m{w}_j^2$$

这样一来，就得到了如下的等式：
$$J(\mathbf{w})=-\left[\sum _{i=1}^n{y}^{[i]}\log \left(a^{[i]}\right)+\left(1-y^{[i]}\right)\log \left(1-a^{[i]}\right)\right]+\frac{\lambda }{2}\Vert \mathbf{w}{\Vert }_2^2$$

在此之前，我们为多类分类实现了一个mlp，它返回t个元素的输出向量，我们需要将其与独热码编码表示中的𝑡×1维目标向量进行比较。如果我们使用此MLP预测具有类别标签2的输入图像的类别标签，则第三层和目标的活性值可能如下所示：
$$a^{(out)}=\left[\begin{array}{c}0.1\\ 0.9\\ ⋮\\ 0.3\end{array}\right],y=\left[\begin{array}{c}0\\ 1\\ ⋮\\ 0\end{array}\right]$$

将上述损失函数推广到一般形式：
$$J(\mathbf{W})=-\sum _{i=1}^n\sum _{j=1}^t{y}_j^{[i]}\log \left(a_j^{[i]}\right)+\left(1-y_j^{[i]}\right)\log \left(1-a_j^{[i]}\right)$$

对应的，增加了正则项：
$$J(\mathbf{W})=-\left[\sum _{i=1}^n\sum _{j=1}^t{y}_j^{[i]}\log \left(a_j^{[i]}\right)+\left(1-y_j^{[i]}\right)\log \left(1-a_j^{[i]}\right)\right]+\frac{\lambda }{2}\sum _{l=1}^{L-1}\sum _{i=1}^{u_l}\sum _{j=1}^{u_{l+1}}{\left(w_{j,i}^{(l)}\right)}^2$$

其中，$u_l$代表的是第$l$层神经元个数，同时正则项表示如下：
$$\frac{\lambda }{2}\sum _{l=1}^{L-1}\sum _{i=1}^{u_l}\sum _{j=1}^{u_{l+1}}{\left(w_{j,i}^{(l)}\right)}^2$$

这里的目标是最小化$J(\mathbf{W})$，因此需要计算损失函数关于$\mathbf{W}$的偏导数：

$$\frac{\partial }{\partial w_{j,i}^{(l)}}J(\mathbf{W})$$

这里的$\mathbf{W}$实际上是由很多矩阵组成，譬如：${\mathbf{W}}^{(h)}$和${\mathbf{W}}^{(out)}$，对于具有一个隐层的MLP，

其三维的张量$\mathbf{W}$可视化即如果如下图所示：

# 几何表示形式
Image(filename='images/12_10.png', width=300) 

3.2反向传播算法backpropagation

在反向传播算法中，很常用的方法为链式求导法则，如对于嵌套函数$[f(g(x))]$：

$$\frac{d}{dx}[f(g(x))]=\frac{df}{dg}\cdot \frac{dg}{dx}$$

类似的，对于更长的复合函数$f(g(h(u(v(x)))))$：

$$\frac{dF}{dx}=\frac{d}{dx}F(x)=\frac{d}{dx}f(g(h(u(v(x)))))=\frac{df}{dg}\cdot \frac{dg}{dh}\cdot \frac{dh}{du}\cdot \frac{du}{dv}\cdot \frac{dv}{dx}$$

逆向自动微分的诀窍是从右往左进行，将一个矩阵乘以一个向量，将产生一个向量，该向量再乘以下一个矩阵，以此类推。

矩阵-向量乘法在计算上比矩阵乘法代价小很多，这就是反向传播算法在神经网络训练中最常用的原因之一。

3.3利用反向传播训练神经网络

使用误差反向传播算法，首先需要经过前向传播计算得到输出层的活性值。计算过程如下：
$$\begin{array}{l}{\mathbf{Z}}^{(h)}={\mathbf{A}}^{(in)}{\mathbf{W}}^{(h)}\text{ (net input of the hidden layer) }\\ {\mathbf{A}}^{(h)}=\varphi \left({\mathbf{Z}}^{(h)}\right)\text{ (activation of the hidden layer) }\\ {\mathbf{Z}}^{(\text{out })}={\mathbf{A}}^{(h)}{\mathbf{W}}^{(\text{out })}\text{ (net input of the output layer) }\\ {\mathbf{A}}^{(\text{out })}=\varphi \left({\mathbf{Z}}^{(\text{out })}\right)\text{ (activation of the output layer) }\end{array}$$

更具体地，通过网络中的连接前向传播输入特征，几何表示如下：

Image(filename='./images/12_11.png', width=400) 

在反向传播中，我们从右到左传播错误。首先计算输出层的误差向量：
$${\mathbf{\delta }}^{(out)}={\mathbf{a}}^{(out)}-\mathbf{y}$$

这里的$\mathbf{y}$为真实类别标签的向量。

接下来，计算隐藏层的误差项：

$${\mathbf{\delta }}^{(h)}={\mathbf{\delta }}^{(\text{out })}{\left({\mathbf{W}}^{(out)}\right)}^T\odot \frac{\partial \varphi \left(z^{(h)}\right)}{\partial z^{(h)}}$$

这里的$\frac{\partial \varphi \left(z^{(h)}\right)}{\partial z^{(h)}}$简单来说，就是激活函数的导数：

$$\frac{\partial \varphi (z)}{\partial z}=\left(a^{(h)}\odot \left(1-a^{(h)}\right)\right)$$

这里的$\odot$代表的是元素级别的乘法。

激活函数的导数：
$$\begin{array}{rl}{\varphi }^{\prime }(z)& =\frac{\partial }{\partial z}\left(\frac{1}{1+e^{-z}}\right)\\ & =\frac{e^{-z}}{{\left(1+e^{-z}\right)}^2}\\ & =\frac{1+e^{-z}}{{\left(1+e^{-z}\right)}^2}-{\left(\frac{1}{1+e^{-z}}\right)}^2\\ & =\frac{1}{\left(1+e^{-z}\right)}-{\left(\frac{1}{1+e^{-z}}\right)}^2\\ & =\varphi (z)-(\varphi (z){)}^2\\ & =\varphi (z)(1-\varphi (z))\\ & =a(1-a)\end{array}$$

接下来，计算$\mathbf{\delta }$层的误差矩阵：

$${\delta }^{(h)}={\delta }^{(out)}{\left({\mathbf{W}}^{(out)}\right)}^T\odot \left(a^{(h)}\odot \left(1-a^{(h)}\right)\right)$$

上面的等式，使用的矩阵转置${\left({\mathbf{W}}^{(out)}\right)}^T$,其中$\left({\mathbf{W}}^{(out)}\right)$为$h\times t$的矩阵，

这里的$t$为输出类别标签的个数，$h$为隐藏层单元的个数

最终在得到$\delta$之后，可以将损失函数的导数写为如下形式：
$$\begin{array}{rl}\frac{\partial }{\partial w_{i,j}^{(out)}}J(\mathbf{W})& =a_j^{(h)}{\delta }_i^{(\text{out })}\\ \frac{\partial }{\partial w_{i,j}^{(h)}}J(\mathbf{W})& =a_j^{(in)}{\delta }_i^{(h)}\end{array}$$

$$\begin{array}{l}{\Delta }^{(h)}={\left({\mathbf{A}}^{(in)}\right)}^T{\delta }^{(h)}\\ {\Delta }^{(out)}={\left({\mathbf{A}}^{(h)}\right)}^T{\delta }^{(out)}\end{array}$$

继而添加正则项：
$${\Delta }^{(l)}:={\Delta }^{(l)}+{\lambda }^{(l)}{\mathbf{W}}^{(l)}$$

最后，权重更新过程为：

$${\mathbf{W}}^{(l)}:={\mathbf{W}}^{(l)}-\eta {\Delta }^{(l)}$$

几何形式表示如下：

Image(filename='images/12_12.png', width=500) 

4.神经网络的收敛性

使用小批量梯度下降，这是一种关于随机梯度下降和批量梯度下降的一种中和方法，

虽然其依然是一种随机方法，但它通常可以得到非常精确的解，而且比常规梯度下降算法收敛更快；

如下图所示，优化算法很容易陷入局部极小；

通过提高学习率，可以较容易地摆脱局部极小，但也增加了超过全局最优值的概率；

Image(filename='images/12_13.png', width=500)