【DL-CV】反向传播，（随机）梯度下降

摘要：下面介绍梯度下降算法以并用反向传播来求梯度梯度下降看名字就和梯度脱不了关系了。运用反向传播，我们可以计算损失函数关于每一层权重的梯度，然后实现每一层权重的训练。应用小批量梯度下降法的随机梯度下降法已经成为当前深度学习的主流算法。

【DL-CV】损失函数，SVM损失与交叉熵损失<前篇---后篇>【DL-CV】激活函数及其选择

有了损失函数L，我们能定量的评价模型的好坏。我们希望损失能最小化，或具体点，我们希望能找到使损失最小化的权重W。当然这个过程不是一步完成的，我们会使用梯度下降算法来一步步修改权重W，好让损失逐渐逼近最小值，这是一个收敛的过程。下面介绍梯度下降算法以并用反向传播来求梯度

看名字就和梯度脱不了关系了。其原理很简单，学过高数都知道，梯度是一个向量，方向指向函数增大最快的方向；那反过来梯度的负值指向函数衰减最快的方向。损失函数展开后是关于权重W的函数L(W)，那其梯度负值 -∇L 指向损失下降最快的方向，我们让权重W往该方向走一小步得到新的权重W_new（更新权重），它对应更低的损失。不停地重复计算梯度，更新权重这两步，权重/模型就会趋于完美（所谓迭代训练过程）

负梯度： $-∇L =- {∂L over ∂W}$, 权重更新： $W_{nwe}={W-alpha {∂L over ∂W}}$，从公式中也知道梯度$∂L over ∂W$$是和W形状一样的矩阵（标量对矩阵求导的特性）。
当然这里的权重W是一个矩阵而不是一个值，会涉及矩阵乘法求梯度。

这里的α叫学习率（通常是个很小的正数），用于控制权重变化的幅度，可以理解为步长。学习率的选择是神经网络训练中最重要的超参数设定之一，学习率太大，容易越过损失函数的最小值然后在最小值周围浮动无法收敛；学习率太小，收敛速度太慢，训练效率降低。关于学习率的设定我们延后详细讲解。

总的来说，这就是原味不加特效的梯度下降法，也就是完全跟着梯度走。随着知识的深入，后面会介绍更高级的更新法则（不再死跟梯度），以获得更高的性能。

原理就是高数中求复合函数偏导数/偏微分的链式法则，我还记得刚学的时候会画树状图，用着老师教的口诀“分叉相加，分层相乘”进行推导。不过那都是很简单的题目了，遇上这种层数多参数多涉及矩阵的神经网络，还是很容易犯错的，最好还是老老实实的看下反向传播吧。

在应用于神级网络模型之前，我们拿个简单的函数直观地展示反向传播原理。
设有函数$f(w,x)={1over 1+e^{w*x}}$,w是行向量，x是列向量。我们可以原函数的计算细分成很多小步，交给不同的子函数完成，这些子函数嵌套（复合函数）起来便能计算出原函数的值。把每个子函数看做一个节点，相连得下图

正向传播是计算函数值，输出结果给下一个函数计算（从左往右）。反向传播求全局梯度的原理就是先求局部梯度的公式，然后接受正向传播传来的结果带入局部梯度公式得到局部梯度，于是每个节点都有自己的局部梯度，最后我们在反向的过程中（从右往左）运用口诀“分叉相加，分层相乘”把它们拼起来得到全局梯度了。从图中可以看出对w的梯度大小是和w一样的。

有人可能会一口气算出复合函数的导数公式然后一个劲地带w和x进去得到结果，像这样简单函数当然是可以的。但这样对计算机来说运算量就加倍了：正向传播时计算了一遍，而用你的公式代w和x又会重复正向传播的计算。面对多层神经网络这是超耗时的，所以正向传播时要分层保存结果供每层局部梯度计算使用（实际上深度学习框架就是这样干的）

我们的多层神经网络的模型就相当于一个巨大的函数，它由多个线性分类器$f(x,W,b) = {W*x+b}$，激活函数，最后加个损失函数（这些函数都是连续可微的）复合而成。运用反向传播，我们可以计算损失函数关于每一层权重的梯度，然后实现每一层权重的训练。

上面反向传播的例子是对权重某个具体值进行求导的，实际使用上我们是对整个权重这个矩阵进行求导的，但概念是通用的，建议初学时写出一个很小很明确的向量化例子，在纸上演算梯度，然后对其一般化，得到一个高效的向量化操作形式。
但更多时候我们会用一个小技巧，这个技巧的关键是分析维度，我们可以通过维度的拼凑使得拼出来的${∂L over ∂W}与{W}形状相同$，总有一个方式是能够让维度之间能够对的上的

# 前向传播
W = np.random.randn(5, 10) # 假设W是 5x10 矩阵
X = np.random.randn(10, 3) # x是 10x3 矩阵
D = W.dot(X) # D=W*x 是 5x3 矩阵


# 假设我们得到了D的梯度∂L/∂D
dD = np.random.randn(*D.shape) # 和D一样的尺寸 5x3

#X.T指X的转置
dW = dD.dot(X.T) #∂L/∂W应该是 5x10 矩阵，由 dD*X.T拼出（[5x3]*[3x10]=[5x10]）
dX = W.T.dot(dD) #∂L/∂X应该是 10x3 矩阵，由 W.T*dD拼出（[10x5]*[5x3]=[10x3]）

有随机当然就有不随机，这种不随机的算法叫批量学习（batch learning）算法，为了引入主角“随机”我们先来聊聊其它。批量学习算法在进行迭代训练（计算损失，计算梯度，权重更新三循环）的时候要遍历全部训练样本，也因此，这种算法能够有效地抑制训练集内带噪声的样本所导致的剧烈变动，并且获得全局最优解；但同时也难免顾此失彼，由于每次更新权重所有样本都要参与训练，训练集一大起来非常耗时。

为了弥补面对大量数据时用时上的缺陷，就有了随机梯度下降法，这里介绍其中一种叫小批量梯度下降（mini-batch gradient descent）的算法：
每次迭代训练时从训练集中随机部分样本（也就是batch size，数量通常为2ⁿ,常用32，64，128，256，根据情况定）进行迭代更新权重。由于每次迭代只使用部分样本，所以和批量学习相比，能减少单次训练时间。它保持收敛性的同时还能减少了迭代结果陷入局部最优解的情况。应用小批量梯度下降法的随机梯度下降法已经成为当前深度学习的主流算法。

# 大概思路
while True:
    data_batch = sample_training_data(data, 256) # 从训练集中随机取256个样本用于训练
    weights_grad = evaluate_gradient(loss_fun, data_batch, weights) #获取梯度
    weights += - step_size * weights_grad # 更新权重

下面我们通过一个小网络例子并运用求导技巧来展示神经网络反向传播的实现。现在我们有一个使用softmax损失的双层网络,我们的目标是求$partial Loverpartial W_1$，$partial Loverpartial W_2$，$partial Loverpartial b_1$，$partial Loverpartial b_2$

输入（X）→→线性分类器+relu激活函数（h1 =max(0,X*W1 + b1) ）→→线性分类器（h2 = h1*W2 + b2）→→输出（S=h2）→→softmax损失函数并使用L2正则化
X是（N*D）的矩阵，包含N个样本数据，每行Xi是第i个样本的数据
以上的h1，h2/s，都是矩阵，每行是第i个样本的层激活值
以下的 j 代表类别的数字，假如这个网络输出的是10个类别的分数，则j=0,1,2,...9

从右至左，首先来推导softmax损失函数的梯度。我们有损失公式：
$$L = { frac{1}{N} sum_i^N L_i }+ { lambda R(W) }$$
共有N个样本，其中第i个样本带来的损失是：
$$L_i=log (sum_je^{s_{ij}})-s_{{iy_i}}$$
当 j!=y_i 时：
$${partial L_iover partial s_{ij}}={e^{s_{ij}}over sum_je^{s_{ij}}}$$
当 j==y_i 时：
$${partial L_iover partial s_{ij}}={e^{s_{ij}}over sum_je^{s_{ij}}}-1$$

据此容易求得${partial L_iover partial s_{i}}$，使用维度技巧开始拼凑，${partial L_iover partial W_2}(20*10)应该由{partial L_iover partial s_i}(1*10)和h^1_i.T(20*1)组合而成$，于是：
$${partial L_iover partial W_2}=h^1_i.T*{partial L_iover partial s_i}quad→(20*10)=(20*1)(1*10)$$

继续往上走会遇到ReLU激活函数逐元素运行max(0,x)，输入小于等于0，回传梯度为0；输入大于0,原封不动回传梯度。我们有$h^1_i =max(0,X_i*W_1+b_1)，另设f_i=X_i*W_1+b_1$，运用维度技巧，${partial L_iover partial h_i}(1*20)应该由{partial L_iover partial s_i}(1*10)和W_2.T(10*20)组合而成$，于是：
$${partial L_iover partial h^1_i}={partial L_iover partial s_i}*W_2.Tquad→(1*20)=(1*10)(10*20)$$
$${partial L_iover partial f_{ik}}={partial L_iover partial h^1_{ik}}quad (h^1_{ik}>0)$$
$${partial L_iover partial f_{ik}}=0quad (h^1_{ik}<=0)$$

继续使用技巧，剩下的${partial L_iover partial W_1}(3072*20)$ 应该有由${partial L_iover partial f_{i}}(1*20)$和$X_i.T(3072*1)$组合而成，于是：
$${partial L_iover partial W_1}=X_i.T*{partial L_iover partial f_{i}}quad→(3072*20)=(3072*1)(1*20)$$

至于偏转值b的梯度就好求多了，由于在公式中局部梯度为1，直接接受传来的梯度即可
$${partial L_ioverpartial b_2}={partial L_ioverpartial s_i}$$
$${partial L_ioverpartial b_1}={partial L_ioverpartial f_i}$$
以上的推导的都是${L_i}$对谁谁的梯度，那么${L}$对谁谁的梯度怎么求？难道求N次梯度求和取平均吗？——不用，矩阵的运算给了我们很好的性质，把上面提到的${X_i}$，${h^1_i}$，${f_i}$，${s_i}$这些向量全换成完整的矩阵${X}$，${h^1}$，${f}$，${s}$进行运算即可，矩阵乘法会执行求和操作，你只需最后求平均即可加上正则化损失的梯度即可。${L}$对$b$的梯度些许不同，需要对传来的梯度在列上求和取平均
$${partial Lover partial W_2}={1over N}h^1.T*{partial Lover partial s}+2lambda W_2quad→(20*10)=(20*N)(N*10)$$
$${partial Lover partial h^1}={partial Lover partial s}*W_2.Tquad→(N*20)=(N*10)(10*20)$$
$${partial Lover partial f_{ik}}={partial Lover partial h^1_{ik}}quad (h^1_{ik}>0)$$
$${partial Lover partial f_{ik}}=0quad (h^1_{ik}<=0)$$
$${partial Lover partial W_1}={1over N}X.T*{partial Lover partial f_{i}}+2lambda W_1quad→(3072*20)=(3072*N)(N*20)$$
${L}$对$b$的梯度些许不同，需要对传来的梯度在列上求和取平均。

import numpy as np
def loss_and_grads(X, y, reg):
    N, D = X.shape  #记录X尺寸，N个样本
    # 正向传播
    h1 = np.maximum(0, np.dot(X, W1) + b1)  #
    h2 = np.dot(h1, W2) + b2
    scores = h2

    scores -= np.repeat(np.max(scores, axis=1), self.num_classes).reshape(scores.shape)  # 预处理防溢出
    exp_class_scores = np.exp(scores)
    exp_corrext_class_scores = exp_class_scores[np.arange(N), y]
    loss = np.log(np.sum(exp_class_scores, axis=1)) - exp_corrext_class_scores
    loss = np.sum(loss)/N
    loss += reg*(np.sum(W2**2)+np.sum(W1**2))

    grads = {}
    # layer2
    dh2 = exp_class_scores / np.sum(exp_class_scores, axis=1, keepdims=True)  # 按行求和保持二维特性以广播
    dh2[np.arange(N), y] -= 1  #j==y_i 时要减一
    dh2 /= N

    dW2 = np.dot(h1.T, dh2)
    dW2 += 2*reg*W2  # 加上正则化损失的梯度

    db2 = np.sum(dh2, axis=0) / N

    # layer1
    dh1 = np.dot(dh2, W2.T)

    df = dh1
    df[h1 <= 0] = 0  # 布尔标记使激活值非正数的梯度为0

    dW1 = np.dot(X.T, df)
    dW1 += 2*reg*W1

    db1 = np.sum(df, axis=0) / N

    grads["W2"] = dW2
    grads["b2"] = db2
    grads["W1"] = dW1
    grads["b1"] = db1

    return loss, grads