机器学习|吴恩达机器学习之神经网络反向传播

正向传播

和机器学习|吴恩达机器学习之神经网络中的内容差不多，都是在给出$\Theta_{(1)}$和$\Theta_{(2)}$的情况下通过正向传播求个代价值或是验证一下准确率。

导入数据

def load_mat(path):
    '''读取数据'''
    data = loadmat('ex4data1.mat')  # return a dict
    X = data['X']
    y = data['y'].flatten()
    
    return X, y  

def load_weight():
    '''读取权重'''
    weight = loadmat('ex4weights.mat')
    theta1 = weight['Theta1']
    theta2 = weight['Theta2']

    return theta1, theta2

这里需要对y向量做一个处理：

def expand_y(y):
    result = []
    for i in y:
        y_array = np.zeros(10)
        y_array[i - 1] = 1
        result.append(y_array)
    return np.array(result)

原来的y是用数字表示：

1	[10 10 10 ... 9 9 9]

转换为矩阵，用1的位置来表示：

[[0. 0. 0. ... 0. 0. 1.]
 [0. 0. 0. ... 0. 0. 1.]
 [0. 0. 0. ... 0. 0. 1.]
 ...
 [0. 0. 0. ... 0. 1. 0.]
 [0. 0. 0. ... 0. 1. 0.]
 [0. 0. 0. ... 0. 1. 0.]]

可视化数据

随机选100张图，可视化观察一波

def plot_data(X):
    '''随机画100个数字'''
    index = np.random.choice(range(5000),
                             100)  # np.random.choice(arrange,size),返回ndarray
    images = X[index]  # 随机选择100个样本
    fig, ax_array = plt.subplots(
        10, 10, sharex=True, sharey=True, figsize=(8, 8))  # ax_array为Axes对象
    for r in range(10):
        for c in range(10):
            ax_array[r, c].matshow(
                images[r * 10 + c].reshape(20, 20), cmap='gray_r'
            )  # matshow() 第一个参数为要显示的矩阵（Display an array as a matrix in a new 				 figure window）
    plt.yticks([])
    plt.xticks([])
    plt.show()

关于plt.subplots展开，参考subplots画图。

计算

因为反向传播和正向传播的公式，代价函数，正则化都是一样的，在后面反向传播再展开讲。

正向传播

1.公式：

2.代码：

def feed_forward(theta1, theta2, X):
    z2 = X @ theta1.T
    a2 = sg.sigmoid(z2)  #(5000,25)
    a2 = np.insert(a2, 0, 1, axis=1)  #(5000,26)
    z3 = a2 @ theta2.T
    a3 = sg.sigmoid(z3)

    return z2, a2, z3, a3

代价函数

1.公式：$J(\theta)=\frac{1}{m} \sum_{i=1}^{m} \sum_{k=1}^{K}\left[-y_{k}^{(i)} \log \left(\left(h_{\theta}\left(x^{(i)}\right)\right)_{k}\right)-\left(1-y_{k}^{(i)}\right) \log \left(1-\left(h_{\theta}\left(x^{(i)}\right)\right)_{k}\right)\right.$

2.代码：

def cost(theta1, theta2, X, y):
    z2, a2, z3, h = feed_forward(theta1, theta2, X)
    # 这里的y是矩阵而不是向量了
    first = -y * np.log(h)
    second = (1 - y) * np.log(1 - h)

    return (np.sum(first - second)) / len(X)  # 这里不能用np.mean()，否则会相差10倍

代价函数正则化

1.公式：

$\begin{aligned} J(\theta)=& \frac{1}{m} \sum_{i=1}^{m} \sum_{k=1}^{K}\left[-y_{k}^{(i)} \log \left(\left(h_{\theta}\left(x^{(i)}\right)\right)_{k}\right)-\left(1-y_{k}^{(i)}\right) \log \left(1-\left(h_{\theta}\left(x^{(i)}\right)\right)_{k}\right)\right]+\\ & \frac{\lambda}{2 m}\left[\sum_{j=1}^{25} \sum_{k=1}^{400}\left(\Theta_{j, k}^{(1)}\right)^{2}+\sum_{j=1}^{10} \sum_{k=1}^{25}\left(\Theta_{j, k}^{(2)}\right)^{2}\right] \end{aligned}$

2.代码：

def cost_reg(theta1, theta2, X, y, lmd):
    c = cost(theta1, theta2, X, y)
    reg = (lmd / (2 * len(X))) * (
        np.sum(theta1[:, 1:]**2) + np.sum(theta2[:, 1:]**2))
    return reg + c

主函数

def main():
    path = 'ex4data1.mat'
    X, y = load_mat(path)
    X = np.insert(X, 0, 1, axis=1)
    y = expand_y(y)
    theta1, theta2 = load_weight()
    print(cost_reg(theta1, theta2, X, y, 1))  #0.38376985909092354
 

if __name__ == "__main__":
    main()

反向传播

随机初始化

之前我们是把$\theta$初始化为0向量，但是在神经网络中如果把$\theta_{1}$初始化为0，那么$S_{2}$中的激活单元都为相同值。同理，只要初始化为相同的数，那么结果都一样。

因此我们通常随机初始化，即在（-$\varepsilon$，$\varepsilon$）之间随机取值，为了保证效率，需要取值足够小，所以选择$\varepsilon=0.12$

1 2	def random_init(size): return np.random.uniform(-0.12, 0.12, size)

处理参数

使用优化参数opt.minimize()，需要把参数展开。

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)

代价函数（带正则项）

以多分类为例

1.公式： AW7Z40.md.png

其中

L：神经网络的层数
$s_{l}$：$l$层中的神经元个数（不包括bias unit）
K：输出层中的神经元个数
m：样本个数

累加项中表示从第一项累加到第k项（why？）

正则项表示神经网络中所有权重的平方和。

梯度项

一般情况下，我们只知道Input Layer和Output Layer两层的神经元个数，中间的Hidden Layer很难确定，不过对于初学者而言，都是参考别人算法里的，所以这里也直接给出了Hidden Layer的层数（1层）以及$\theta_{1}$（25，401）和$\theta_{2}$（10，26）的维度，一个神经元为一列。

计算前馈(feedforward)

参数含义及传递过程如下

1.参数含义：

$\Theta^{i}$第$i$层的参数矩阵
$z^{(l)}$第$l$层的输入
$a^{(l)}$第$l$层的输出

2.传递过程：

$a^{(1)}=x$
$z^{(2)}=\Theta^{(1)} a^{(1)}$
$a^{(2)}=g\left(z^{(2)}\right)\left(a d d \ bias \ \ a_{0}^{(2)}\right)$
$z^{(3)}=\Theta^{(2)} a^{(2)}$
$h=a^{(3)}=g\left(z^{(3)}\right)$

3.前馈代码：

def feed_forward(theta1, theta2, X):
    z2 = X @ theta1.T
    a2 = sg.sigmoid(z2)  #(5000,25)
    a2 = np.insert(a2, 0, 1, axis=1)  #(5000,26)
    z3 = a2 @ theta2.T
    a3 = sg.sigmoid(z3) 
    
    return z2, a2, z3, a3

计算梯度项

计算梯度项，也就是代价函数的偏导数$\frac{\partial}{\partial\Theta^{(l)}_{ij}}J\left(\Theta\right)$。通过前馈的计算我们得到了$h$，接下来计算“误差”，这里之所以用引号，是因为误差的实质是$\delta^{(l)}=\frac{\partial J}{\partial z^{(l)}}$

$\delta^{(3)}=h-y$ ……(1)
$\delta^{(2)}=(\Theta^{(2)})^{T}\delta^{(3)}g^{\prime}\left(z^{(2)}\right)$……(2)

第一层没有误差，接下去计算每层参数矩阵的梯度值，用$\Delta^{(l)}$表示

$\Delta^{(2)}=a^{(2)} \delta^{(3)}$……(3)
$\Delta^{(1)}=a^{(1)} \delta^{(2)}$……(4)

最后网络的总梯度为：

（这里并不是真的相加，而是将$\Delta^{(1)}$和$\Delta^{(2)}$合成为一个向量，方便后面计算）

$D=\frac{1}{m}\left(\Delta^{(1)}+\Delta^{(2)}\right)$

求梯度项代码：

def gradient(theta1,theta2,X,y):
    z2,a2,z3,h=feed_forward(theta1,theta2,X)
    d3=h-y # (5000,10) 
    d2=d3@theta2[:,1:]*sg.sigmoid_gradient(z2) # (5000,25)
    D2=d3.T@a2 # (10,26)
    D1=d2.T@X # (25,401)
    
    # 这里合并成1-D array是为了方便后面用优化函数处理
    D=(1/len(X))*serialize(D1,D2) #(10285,)

    return D

正则化

1.原理

2.代码：

def regularized_gradient(theta,X,y,lmd=1):
    theta1,theta2=deserialize(theta)
    D1,D2=deserialize(gradient(theta1,theta2,X,y))
    theta1[:,0]=0
    theta2[:,0]=0
    reg_D1=D1+(lmd/len(X))*theta1
    reg_D2=D2+(lmd/len(X))*theta2

    return serialize(reg_D1,reg_D2)

*推导$\delta$和$\Delta$

从代价函数入手，假设我们只有一个输入样本，那么代价函数为：$J(\theta)=-y \operatorname{logh}(x)-(1-y) \log (1-h)$，所谓梯度项，就是将代价函数对参数求导，即$\frac{\partial}{\partial \Theta^{(2)}} J(\theta), \frac{\partial}{\partial \Theta^{(1)}} J(\theta)$。而由传递过程函数：

$a^{(1)}=x$
$z^{(2)}=\Theta^{(1)} a^{(1)}$
$a^{(2)}=g\left(z^{(2)}\right)\left(a d d \ bias \ \ a_{0}^{(2)}\right)$
$z^{(3)}=\Theta^{(2)} a^{(2)}$
$h=a^{(3)}=g\left(z^{(3)}\right)$

我们可以使用链式求导法则，因此有$\frac{\partial J}{\partial \Theta^{(2)}}=\frac{\partial J}{\partial a^{(3)}} \frac{\partial a^{(3)}}{\partial z^{(3)}} \frac{\partial z^{(3)}}{\partial \Theta^{(2)}}=(h-y)a^{(2)}$

其中令$\delta^{(3)}=h-y$得到公式(1)；令$\Delta^{(2)}=\frac{\partial J}{\partial \Theta^{(2)}}$则得到公式(3)。

接着求$\frac{\partial J}{\partial \Theta^{(1)}}$=$\frac{\partial J}{\partial a^{(3)}} \frac{\partial a^{(3)}}{\partial z^{(3)}} \frac{\partial z^{(3)}}{\partial a^{(2)}} \frac{\partial a^{(2)}}{\partial z^{(2)}} \frac{\partial z^{(2)}}{\partial \Theta^{(1)}}$

=$\delta^{(3)} \Theta^{(2)} g^{\prime}\left(Z^{(2)}\right) a^{(1)}$

=$\delta^{(2)} a^{(1)}$

同样可以看出，令$\delta^{(3)}=\frac{\partial J}{\partial a^{(3)}} \frac{\partial a^{(3)}}{\partial z^{(3)}}$ ，$\delta^{(2)}=\delta^{(3)} \Theta^{(2)} g^{\prime}\left(Z^{(2)}\right)$则得到公式(2)。令$\frac{\partial J}{\partial \Theta^{(1)}}=\Delta^{(1)}$得到公式(4)。

梯度检验

目的：在反向传播的过程中，因为需要计算的参数很多，因此容易导致误差，使得最终的结果并非最优解。因此为了确定反向传播计算的梯度是否正确，需要用到梯度检验(gradient check)。

原理：通过计算$\frac{\partial}{\partial \Theta}J(\Theta)=\lim _{\varepsilon \rightarrow 0} \frac{J(\theta+\varepsilon)-J(\theta-\varepsilon)}{2 \varepsilon}$，估计出$J(\theta)$在$\theta$的值，和反向传播计算的梯度值$\Delta$进行对比。具体来说，对于某个特定的 $\theta$，我们计算出在 $\theta$-$\varepsilon $ 处和 $\theta$+$\varepsilon $ 的代价值（$\varepsilon $是一个非常小的值，通常选取 0.001），然后求两个代价的平均，用以估计在 $\theta$ 处的代价值。

具体做法：先将$\Theta_{(1)}$和$\Theta_{(2)}$展开成一个向量$\theta$，维度（10285，）。然后循环，对每个$\theta_{j}$求$\frac{\partial}{\partial \theta_{j}}J(\theta)$，以$\theta_{1}$为例，$\frac{\partial}{\partial \theta_{1}}=\frac{J\left(\theta_{1}+\varepsilon_{1}, \theta_{2}, \theta_{3} \dots \theta_{n}\right)-J\left(\theta_{1}-\varepsilon_{1}, \theta_{2}, \theta_{3} \ldots \theta_{n}\right)}{2 \varepsilon}$，注意这里虽然只有$\theta_{j}$每次都要把整个$J(\theta)$带进去，因此需要每次都复制整个$\theta$。

def gradient_check(theta1,theta2,X,y,e):
    theta_temp=serialize(theta1,theta2) # (10285,)
    numeric_grad=[]
    for i in range(len(theta_temp)):
        plus=copy.copy(theta_temp) # 复制theta
        minus=copy.copy(theta_temp) # 复制theta
        plus[i]+=e # 只改变当前的一个值
        minus[i]-=e # 只改变当前的一个值
        
        grad_i=(cost_reg(plus,X,y,1)-cost_reg(minus,X,y,1))/(e*2)
        numeric_grad.append(grad_i)

再和$\Theta_{(1)}$和$\Theta_{(2)}$比较求得准确度。

这里简单介绍一下数值梯度(numerical gradient)和解析梯度(analytic gradient)，数值梯度的优点是编程可以直接实现，不要求函数可微，缺点是运行速度非常慢，也就是上面中numeric_grad，且只能求出近似解；解析梯度能求出近似解，也是我们通常使用的方法，即analytic_grad。

得到数值梯度和解析梯度之后，要求他们的进行相似性度量(Similarity Measurement)，这里用标准化欧氏距离（存疑？）$diff=\frac{ | \text { numeric_grad }-\text {analytic_grad}\left|_{2}\right.}{ | \text { numeric_grad }\left|_{2}+\right| \text { analytic_grad }\left|_{2}\right.}$，当距离(diff)小于10e-9时为计算正确。

reg_D1,reg_D2=regularized_gradient(theta_temp,X,y)
analytic_grad=serialize(reg_D1,reg_D2)
   
   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))

np.linalg.norm()：linalg=linear(线性)+algebra(代数)，norm()表示范数$^{[1]}$。

优化参数

def nn_trainning(X,y):

    init_theta = random_init(10285) #随机初始化theta1和theta2
    res = opt.minimize(
        fun=cost_reg,
        x0=init_theta,
        args=(X, y, 1),
        method='TNC',
        jac=regularized_gradient,
        options={'maxiter': 400})

    return res

可视化隐藏层

一种明白神经网络是如何学习的方法就是将隐藏层捕获的内容可视化，通俗来说就是输入一个x，激活这个隐藏层。在我们的这个训练样本中，$\theta_{1}$（25,401）有401个参数，去掉偏置单元，将剩下的400个参数reshape为(20,20)。

def plot_hidden(theta):
    t1,_=deserialize(theta)
    t1=t1[:,1:]
    fig,ax_array=plt.subplots(5,5,sharex=True,sharey=True,figsize=(6,6))
    for r in range(5):
        for c in range(5):
            ax_array[r,c].matshow(t1[r*5+c].reshape(20,20),cmap='gray_r')
    plt.xticks([])
    plt.yticks([])
    plt.show()

计算精确度

def accuracy(theta, X, y):
    theta1, theta2 = deserialize(theta)
    _, _, _, h = feed_forward(theta1, theta2, X)
    y_pred = np.argmax(h, axis=1) + 1

    print(classification_report(y, y_pred))

主函数

def main():
    path = 'ex4data1.mat'
    X, raw_y = ff.load_mat(path)
    X = np.insert(X, 0, 1, axis=1)
    y = ff.expand_y(raw_y)  # y的一行表示一个样本
    # theta1, theta2 = ff.load_weight()
    # gradient_check(theta1,theta2,X,y,0.0001)
    theta_unroll = nn_training(X, y)
    accuracy(theta_unroll.x, X, raw_y)


if __name__ == "__main__":
    main()

参考资料

[1] np.linalg.norm(求范数)

[2] 吴恩达机器学习作业Python实现(四)：神经网络(反向传播)