神经网络——吴恩达课后作业python代码实现

在上一次多类别分类器练习中，通过多类别逻辑回归算法实现了识别数字笔迹。
但是，逻辑回归算法不能得到一个更为复杂的假设函数，因为它仅仅只是一个线性分类，而神经网络能很好地实现这一点。
神经网路中最基本的成分是“神经元”模型，在生物神经网络中，每个神经元与其他神经元相连，当他“兴奋”时，就会向相连的神经元发送化学物质，从而改变这些神经元类的电位；如果某神经元的电位超过一个“阈值”，那么他就会被激活。神经元接收来自n个神经元传递过来的输入信号，这些输入信号通过带权重的连接进行传输，将神经元接收到的总输入值与神经元的自身的阈值进行比较，然后通过“激活函数”处理产生神经元的输出。理想中的激活函数是阶跃函数，但是它不连续，通常用sigmoid函数作为激活函数。

一，基于神经网络的多类别分类器

1.1模型建立

建立模型示意图：

上图所示神经网络一共有三个层，分别是输入层（input layer），隐藏层（hidden layer），输出层（output layer）。
特点如下：
1、每层由单元（units）组成
2、输入层是有训练集的实例特征向量传入
3、经过连接接点的权重（weight）传入下一层，一层的输出是下一层的输入
4、隐藏层的个数可以是任意的，输入层有一层，输出层有一层
5、每个单元也可以称之为神经结点，根据生物学来源定义
6、以上成为两层的神经网络，输入层是不算在里面的
7、一层中加权求和，然后根据非线性方程转化输出
8、作为多层向前神经网络，理论上，如果有足够的隐藏层，和足够的训练集，可以模拟出任何方程

1.2前馈传播及多类别分类器的实现

练习的数据集给出了实现数字字迹识别算法中需要的$\theta$值，这里只需要运用神经网络前馈算法，利用已知的$\theta$便可以轻松完成多类别分类器的实现。

上图的计算过程简写:
$a_i^{(j)}=$第J层第i单元的“激活项”
${\theta }^{(j)}=$控制函数从第J层映射到第J+1层的权重矩阵
$a_1^{(2)}=g({{\theta }_{10}}^{(1)}{x}_0+{{\theta }_{11}}^{(1)}{x}_1+{{\theta }_{12}}^{(1)}{x}_2+{{\theta }_{13}}^{(1)}{x}_3)$
$a_2^{(2)}=g({{\theta }_{20}}^{(1)}{x}_0+{{\theta }_{21}}^{(1)}{x}_1+{{\theta }_{22}}^{(1)}{x}_2+{{\theta }_{23}}^{(1)}{x}_3)$
$a_3^{(2)}=g({{\theta }_{30}}^{(1)}{x}_0+{{\theta }_{31}}^{(1)}{x}_1+{{\theta }_{32}}^{(1)}{x}_2+{{\theta }_{33}}^{(1)}{x}_3)$
…
$h_{\theta }(x)=a^{(3)}=g({\theta }^{(2)}{x}_0+{\theta }^{(2)}{x}_1+{\theta }^{(2)}{x}_2+{\theta }^{(2)}{x}_3)$

…
令$z_1^{(2)}={{\theta }_{10}}^{(1)}{x}_0+{{\theta }_{11}}^{(1)}{x}_1+{{\theta }_{12}}^{(1)}{x}_2+{{\theta }_{13}}^{(1)}{x}_3$
$z_2^{(2)}={{\theta }_{20}}^{(1)}{x}_0+{{\theta }_{21}}^{(1)}{x}_1+{{\theta }_{22}}^{(1)}{x}_2+{{\theta }_{23}}^{(1)}{x}_3$
…
其中g(z)是sigmoid函数

前向传播中“激活项”相关的代码实现

 def feed_forward(theta, X,):
    '''得到每层的输入和输出'''
    t1, t2 = deserialize(theta)
    
    a1 = X
    z2 = a1 @ t1.T
    a2 = np.insert(sigmoid(z2), 0, 1, axis=1)
    z3 = a2 @ t2.T
    a3 = sigmoid(z3)
    
    return a1, z2, a2, z3, a3 

通过$\theta$矩阵，可以由输入通过神经网络前向传播得到输出，即多类别分类器的代码实现：

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

def load_weight(path):
    data = loadmat(path)
    return data['Theta1'], data['Theta2']

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

def load_data(path):
    data = loadmat(path)
    X = data['X']
    y = data['y']
    return X,y

theta1, theta2 = load_weight('ex3weights.mat')

theta1.shape, theta2.shape

X, y = load_data('ex3data1.mat')
y = y.flatten()
X = np.insert(X, 0, values=np.ones(X.shape[0]), axis=1)  # intercept

X.shape, y.shape
a1 = X
z2 = a1 @ theta1.T
z2.shape
z2 = np.insert(z2, 0, 1, axis=1)
a2 = sigmoid(z2)
a2.shape
z3 = a2 @ theta2.T
z3.shape
a3 = sigmoid(z3)
a3.shape
y_pred = np.argmax(a3, axis=1) + 1
accuracy = np.mean(y_pred == y)
print ('accuracy = {0}%'.format(accuracy * 100))  

可以看出基于神经网络的多类别分类器的准确率为0.9752。与上一篇博客中运用逻辑回归算法的效果要好。（逻辑回归算法得到的准确率为0.9446。）

1.3代价函数

已知逻辑回归的代价函数为：
$J(\theta )=-\frac{1}{m}[{\sum }_{i=1}^m{y}^{(i)}\log (h_{\theta }(x^{(i)})-(1-y^{(i)})\log (1-h_{\theta }(x^{(i)}))]+\frac{\lambda }{2m}{\sum }_{j=1}^n{\theta }_j^2$
而神经网络的代价函数也是这个式子的一般形式，但这里不再仅有一个逻辑回归输出单元，取而代之的是存在K个。
所以，有神经网络的代价函数为:
$J(\theta )=\frac{1}{m}{\sum }_{i=1}^m{\sum }_{k=1}^K[-y_k^{(i)}\log (h_{\theta }(x_k^{(i)})-(1-y_k^{(i)})(\log (1-h_{\theta }(x^{(i)}){)}_k)]+\frac{\lambda }{2m}{\sum }_{l=1}^{L-1}{\sum }_{i=1}^{S_l}{\sum }_{j=1}^{S_{(l+1)}}({\theta }_{(ji)}^{(l)}{)}^2$
其中上图表示了$h_{\theta }(x^{(i)})$的计算过程，K=10是可能的标签总数。注意，$h_{\theta }(x^{(i)})=a_k^{(3)}$是第k个输出单元的激活(输出值)。此外，r 既然原来的标签(在变量y中)是1，2，…，10，为了训练神经网络，我们需要将标签重新编码为只包含值0或1的向量，即

例如，如果$x^{(i)}$是数字5的图像，最终得到的输出结果$y^{(i)}$应该是一个十维向量，其中元素$y_5=1$,并且其它元素的值为0。

def cost(theta, X, y):
    a1, z2, a2, z3, h = feed_forward(theta, X)
    J = 0
    for i in range(len(X)):
        first = - y[i] * np.log(h[i])
        second = (1 - y[i]) * np.log(1 - h[i])
        J = J + np.sum(first - second)
    J = J / len(X)
    return J
'''
     # or just use verctorization
     J = - y * np.log(h) - (1 - y) * np.log(1 - h)
     return J.sum() / len(X)
'''
def regularized_cost(theta, X, y, l=1):
    '''正则化时忽略每层的偏置项，也就是参数矩阵的第一列'''
    t1, t2 = deserialize(theta)
    reg = np.sum(t1[:,1:] ** 2) + np.sum(t2[:,1:] ** 2)  # or use np.power(a, 2)
    return l / (2 * len(X)) * reg + cost(theta, X, y)

二、神经网络反向传播

在多类别分类器中，权重$\theta$矩阵是给定的，能够轻松构建分类器。在线性回归，逻辑回归中可以用梯度下降法，或者其他优化算法得到代价函数收敛后的$\theta$值。而在神经网络中，得到前馈传播权重的方法是反向传播。

2.1sigmoid函数梯度

sigmoid函数：
$sigmoid(z)=g(z)=\frac{1}{1+e^{-z}}$
它的导数
$g^{\prime }(z)=\frac{d}{dz}g(z)=g(z)(1-g(z))$
公式用代码实现

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

2.2随机初始化

当训练神经网络时，初始化参数的设置很重要。当设置的权重相同时，每层传播后得到的各个单元的数值相同，最后得到的y没有任何意义。所以需要随机初始化参数，破坏对称。
随机初始化的一种有效策略是在$[-\epsilon ,\epsilon ]$范围内均匀地为${\Theta }^{(l)}$选择值。其中应该使用$\epsilon =0.12$。这个范围使得参数足够小，训练更有效率

def random_init(size):
    '''从服从的均匀分布的范围中随机返回size大小的值'''
    return np.random.uniform(-0.12, 0.12, size)

2.3反向传播算法

根据已知训练集$(x^{(t)},y^{(t)})$,我们将利用神经网络前向传播，计算出整个网络中的“激活项”，以及最后的$h_{\Theta }(x)$.然后对于第$l$层的第$j$个节点，将计算出一个“误差项”${\delta }_j^{(l)}$,这衡量了该节点对我们输出中的任何错误“负责”的程度。
对于输出节点，我们可以直接测量网络的激活与真实目标值之间的差异，并且使用该差值来定义一个J(3)(因为层3是输出层)。对于隐藏层的 N个单位，将根据层(L+1)中节点的误差项的加权平均值来计算J(L)。
具体地，这里有反向传播算法(也在上图中示出)。
步骤: 1，将输入层的值(a(1))设置为第t个训练例x(t)。进行前向传播，计算第二层，第三层的“激活项”$(z^{(2)},a^{(2)},z^{(3)},a^{(3)})$。需要注意的是，需要添加一个1项，以确保层a(1)和a(2)的激活向量也包括偏置单元。

步骤： 2，对每一个输出层（第三层）的输出单元，设置
${\delta }_k^{(3)}=(a_k^{(3)}-y_k)$,
其中，$y_k\in 0,1$,表示当前训练样本是否属于第K类（$y_k=1$），或者属于其他类（$y_k=0$）。

步骤： 3，对于隐藏层$（l=2）$,设置
${\delta }^{(2)}=({\Theta }^{(2)}{)}^T{\delta }^{(3)}.\ast g^{{}^{\prime }}(z^{(2)})$

步骤： 4，计算$\Delta$值
对于训练集{$(x^{(1)},y^{(1)}),...,(x^{(m)},y^{(m)})$}
初始化设置${\Delta }_{ij}=0$ $(foralll,i,j)$
${\Delta }_{ij}:={\Delta }_{ij}+a_j^{(l)}{\delta }_i^{(l+1)}$

步骤： 5，计算代价函数的梯度
$D_{ij}^{(l)}:=1/m{\Delta }_{ij}^{(l)}+\lambda {\Theta }_{i{j}^{(l)}}ifj̸=0$
$D_{ij}^{(l)}:=1/m{\Delta }_{ij}^{(l)}ifj=0$
$\frac{\partial }{\partial {\theta }_{ij}^{(l)}}J(\theta )=D_{ij}^{(l)}$

循环步骤1到4，最终得到代价函数的偏导数。
代码实现

def gradient(theta, X, y):
    '''
    unregularized gradient, notice no d1 since the input layer has no error 
    return 所有参数theta的梯度，故梯度D(i)和参数theta(i)同shape，重要。
    '''
    t1, t2 = deserialize(theta)
    a1, z2, a2, z3, h = feed_forward(theta, X)
    d3 = h - y # (5000, 10)
    d2 = d3 @ t2[:,1:] * sigmoid_gradient(z2)  # (5000, 25)
    D2 = d3.T @ a2  # (10, 26)
    D1 = d2.T @ a1 # (25, 401)
    D = (1 / len(X)) * serialize(D1, D2)  # (10285,)
    
    return D

正则化

def regularized_cost(theta, X, y, l=1):
    '''正则化时忽略每层的偏置项，也就是参数矩阵的第一列'''
    t1, t2 = deserialize(theta)
    reg = np.sum(t1[:,1:] ** 2) + np.sum(t2[:,1:] ** 2)  # or use np.power(a, 2)
    return l / (2 * len(X)) * reg + cost(theta, X, y)

其中，deserialize函数功能是提取参数
当我们使用高级优化方法来优化神经网络时，我们需要将多个参数矩阵展开，才能传入优化函数，然后再恢复形状。

def serialize(a, b):
    '''展开参数'''
    return np.r_[a.flatten(),b.flatten()]
  
def deserialize(seq):
    '''提取参数'''
    return seq[:25*401].reshape(25, 401), seq[25*401:].reshape(10, 26)

2.4梯度检测

在神经网络算法中，希望得到最小的代价函数$J(\theta )$。要对参数执行梯度检查，可以想象将参数Θ(1)、Θ(2)“展开”为长向量θ。通过这样做，可以考虑代价函数为J(θ)，并使用以下梯度检查过程。
假设有一个函数fi(θ)，该函数旨在计算$\frac{\partial }{\partial {\theta }_i}J(\theta )$；希望检查fi是否输出正确的导数值。

这两个值的近似程度取决于J的细节，通过上式得到的结果去检测通过反向传播算法得到梯度。
代码实现

def gradient_checking(theta, X, y, e):
    def a_numeric_grad(plus, minus):
        """
        对每个参数theta_i计算数值梯度，即理论梯度。
        """
        return (regularized_cost(plus, X, y) - regularized_cost(minus, X, y)) / (e * 2)
   
    numeric_grad = [] 
    for i in range(len(theta)):
        plus = theta.copy()  # deep copy otherwise you will change the raw theta
        minus = theta.copy()
        plus[i] = plus[i] + e
        minus[i] = minus[i] - e
        grad_i = a_numeric_grad(plus, minus)
        numeric_grad.append(grad_i)
    
    numeric_grad = np.array(numeric_grad)
    analytic_grad = regularized_gradient(theta, X, y)
    diff = np.linalg.norm(numeric_grad - analytic_grad) / np.linalg.norm(numeric_grad + analytic_grad)

    print('If your backpropagation implementation is correct,\nthe relative difference will be smaller than 10e-9 (assume epsilon=0.0001).\nRelative Difference: {}\n'.format(diff))

三、优化参数

使用minimize函数优化参数
当隐藏层变多的时候我们的损失函数就变成了一个非凸函数如下图所示

这是考虑使用minimize进行非线性规划

scipy.optimize.minimize(fun, x0, args=(), method=None, jac=None, hess=None, hessp=None, bounds=None, constraints=(), tol=None, callback=None, options=None)

fun: 求最小值的目标函数

x0: 变量的初始猜测值，如果有多个变量，需要给每个变量一个初始猜测值。minimize是局部最优的解法，所以

args: 常数值,后面demo会讲解，fun中没有数字，都以变量的形式表示，对于常数项，需要在这里给值

method: 求极值的方法，官方文档给了很多种。一般使用默认。每种方法我理解是计算误差，反向传播的方式不同而已，这块有很大理论研究空间

constraints: 约束条件，针对fun中为参数的部分进行约束限制

def nn_training(X, y):
    init_theta = random_init(10285)  # 25*401 + 10*26

    res = opt.minimize(fun=regularized_cost,
                       x0=init_theta,
                       args=(X, y, 1),
                       method='TNC',
                       jac=regularized_gradient,
                       options={'maxiter': 400})
    return res

res = nn_training(X, y)#慢
res
'''
     fun: 0.5156784004838036
     jac: array([-2.51032294e-04, -2.11248326e-12,  4.38829369e-13, ...,
        9.88299811e-05, -2.59923586e-03, -8.52351187e-04])
 message: 'Converged (|f_n-f_(n-1)| ~= 0)'
    nfev: 271
     nit: 17
  status: 1
 success: True
       x: array([ 0.58440213, -0.02013683,  0.1118854 , ..., -2.8959637 ,
        1.85893941, -2.78756836])
'''

def accuracy(theta, X, y):
    _, _, _, _, h = feed_forward(res.x, X)
    y_pred = np.argmax(h, axis=1) + 1
    print(classification_report(y, y_pred))

accuracy(res.x, X, raw_y)
'''
             precision    recall  f1-score   support

          1       0.97      0.99      0.98       500
          2       0.98      0.97      0.98       500
          3       0.98      0.95      0.96       500
          4       0.98      0.97      0.97       500
          5       0.97      0.98      0.97       500
          6       0.99      0.98      0.98       500
          7       0.99      0.97      0.98       500
          8       0.96      0.98      0.97       500
          9       0.97      0.98      0.97       500
         10       0.99      0.99      0.99       500

avg / total       0.98      0.98      0.98      5000
'''