从矩阵链式求导的角度来深入理解BP算法（原理+代码）

本文链接：https://blog.csdn.net/Cyril_KI/article/details/122440968

I. 前言

在反向传播算法(Backpropagation)----Gradient Descent的推导过程中我简单推导过BP算法，但是只是简单用符号来表示参数间的求导关系，并没有详细介绍求导的每一步具体步骤。因此，本篇文章着重推导反向传播时的求导过程。

系列文章：

II. 网络结构

本次搭建的神经网络参考了我的另一篇博文：手写神经网络识别MNIST数据集（pytorch读入数据）。神经网络由1个输入层、3个隐藏层以及1个输出层组成，激活函数全部采用Sigmoid函数。

III. 前向传播

网络各层间的运算关系，也就是前向传播过程如下所示：
$z_1=Iw_1,h_1=\sigma(z_1)\\ z_2=h_1w_2,h_2=\sigma(z_2)\\ z_3=h_2w_3,h_3=\sigma(z_3)\\ z_4=h_3w_4,O=\sigma(z_4)\\ loss=\frac{1}{2} (O-y)^2$

其中：

输入 $I$ 的shape为 $100 \times32$ ，每一次输入100个样本，每个样本的特征数为32
$w_1$ 的shape为 $32 \times 20$
$z_1$ 和 $h_1$ 的shape都为 $100 \times 20$
$w_2$ 的shape为 $20 \times 20$
$z_2$ 和 $h_2$ 的shape都为 $100 \times 20$
$w_3$ 的shape为 $20 \times 20$
$z_3$ 和 $h_3$ 的shape都为 $100 \times 20$
$w_4$ 的shape为 $20 \times 1$
$z_4$ 和 $O$ 的shape都为 $100 \times 1$ ，表示100个样本的输出
$l oss$ 为损失函数，shape和输出 $O$ 一致

代码实现：

初始化

class BP:
    def __init__(self):
        self.input = np.zeros((100, 32))  # 100 samples per round
        self.w1 = 2 * np.random.random((32, 20)) - 1  # limit to (-1, 1)
        self.z1 = 2 * np.random.random((100, 20)) - 1
        self.hidden_layer_1 = np.zeros((100, 20))
        self.w2 = 2 * np.random.random((20, 20)) - 1
        self.z2 = 2 * np.random.random((100, 20)) - 1
        self.hidden_layer_2 = np.zeros((100, 20))
        self.w3 = 2 * np.random.random((20, 20)) - 1
        self.z3 = 2 * np.random.random((100, 20)) - 1
        self.hidden_layer_3 = np.zeros((100, 20))
        self.w4 = 2 * np.random.random((20, 1)) - 1
        self.z4 = 2 * np.random.random((100, 1)) - 1
        self.output_layer = np.zeros((100, 1))
        self.loss = np.zeros((100, 1))
        self.learning_rate = 0.08

前向传播

def forward_prop(self, data, label):  # label:100 X 1,data: 100 X 32
    self.input = data
    self.z1 = np.dot(self.input, self.w1)
    self.hidden_layer_1 = self.sigmoid(self.z1)
    self.z2 = np.dot(self.hidden_layer_1, self.w2)
    self.hidden_layer_2 = self.sigmoid(self.z2)
    self.z3 = np.dot(self.hidden_layer_2, self.w3)
    self.hidden_layer_3 = self.sigmoid(self.z3)
    self.z4 = np.dot(self.hidden_layer_3, self.w4)
    self.output_layer = self.sigmoid(self.z4)
    # loss
    self.loss = 1 / 2 * (label - self.output_layer) ** 2
    
    return self.output_layer

IV. 反向传播

为了实现对4个 $w$ 进行更新，我们需要求出损失函数 $l$ 对4个 $w$ 的导数，然后再对4个 $w$ 分别进行更新：
$\frac{\partial l}{\partial w}$
因此，反向传播的过程中唯一任务就是找出损失函数 $l$ 对4个 $w$ 的导数。

首先分析一下 $w_4$ ：

由
$z_4=h_3w_4$

可知：
$\frac{\partial l}{\partial w_4}=h_3^T \frac{\partial l}{\partial z_4}$
在矩阵的链式求导过程中，若求导对象在前向传播式子中在后，则相应梯度在后；如果求导对象在前，则相应梯度在前。也就是说，如果要算损失函数对 $h_3$ 的导数， $h_3$ 在前，则有：
$\frac{\partial l}{\partial h_3}= \frac{\partial l}{\partial z_4}w_4^T$
因此，现在需要寻找 $\frac{\partial l}{\partial z_4}$ ，而因为又有：
$O=\sigma(z_4)$
和
$loss=\frac{1}{2} (y-O)^2$
所以我们有：
$\frac{\partial l}{\partial z_4} = \frac{\partial l}{\partial O}\frac{\partial O}{\partial z_4}$

故：
$\frac{\partial l}{\partial z_4}=(O-y)\odot \sigma^{'}(O)$

因此，最终求得的损失函数对 $w_4$ 的导数就为：
$\frac{\partial l}{\partial w_4} = h_3^T(\frac{\partial l}{\partial O} \odot \sigma^{'}(O))$
相应代码为：

l_deri_out = self.output_layer - label
l_deri_z4 = l_deri_out * self.sigmoid_deri(self.output_layer)
l_deri_w4 = np.dot(self.hidden_layer_3.T, l_deri_z4)

那么同理，对 $w_3$ 我们有：
$\frac{\partial l}{\partial w_4} = h_2^T(\frac{\partial l}{\partial h_3} \odot \sigma^{'}(h_3))$
相应代码为：

l_deri_h3 = np.dot(l_deri_z4, self.w4.T)
l_deri_z3 = l_deri_h3 * self.sigmoid_deri(self.hidden_layer_3)
l_deri_w3 = np.dot(self.hidden_layer_2.T, l_deri_z3)

对 $w_2$ 有：
$\frac{\partial l}{\partial w_2} = h_1^T(\frac{\partial l}{\partial h_2} \odot \sigma^{'}(h_2))$
相应代码为：

l_deri_h2 = np.dot(l_deri_z3, self.w3.T)
l_deri_z2 = l_deri_h2 * self.sigmoid_deri(self.hidden_layer_2)
l_deri_w2 = np.dot(self.hidden_layer_1.T, l_deri_z2)

对 $w_1$ 有：
$\frac{\partial l}{\partial w_1} = I^T(\frac{\partial l}{\partial h_1} \odot \sigma^{'}(h_1))$
相应代码为：

l_deri_h1 = np.dot(l_deri_z2, self.w2.T)
l_deri_z1 = l_deri_h1 * self.sigmoid_deri(self.hidden_layer_1)
l_deri_w1 = np.dot(self.input.T, l_deri_z1)

因此，反向传播的完整代码为：

def backward_prop(self, label):
    # w4
    l_deri_out = self.output_layer - label
    l_deri_z4 = l_deri_out * self.sigmoid_deri(self.output_layer)
    l_deri_w4 = np.dot(self.hidden_layer_3.T, l_deri_z4)
    # w3
    l_deri_h3 = np.dot(l_deri_z4, self.w4.T)
    l_deri_z3 = l_deri_h3 * self.sigmoid_deri(self.hidden_layer_3)
    l_deri_w3 = np.dot(self.hidden_layer_2.T, l_deri_z3)
    # w2
    l_deri_h2 = np.dot(l_deri_z3, self.w3.T)
    l_deri_z2 = l_deri_h2 * self.sigmoid_deri(self.hidden_layer_2)
    l_deri_w2 = np.dot(self.hidden_layer_1.T, l_deri_z2)
    # w1
    l_deri_h1 = np.dot(l_deri_z2, self.w2.T)
    l_deri_z1 = l_deri_h1 * self.sigmoid_deri(self.hidden_layer_1)
    l_deri_w1 = np.dot(self.input.T, l_deri_z1)
    # update
    self.w4 -= self.learning_rate * l_deri_w4
    self.w3 -= self.learning_rate * l_deri_w3
    self.w2 -= self.learning_rate * l_deri_w2
    self.w1 -= self.learning_rate * l_deri_w1

V. 实验

1.数据处理

本次实验拟对负荷进行预测：利用前24个小时的负荷+环境因素来预测下一时刻的负荷值。

def nn_seq(B):
    dataset = load_data()
    # split
    train = dataset[:int(len(dataset) * 0.6)]
    val = dataset[int(len(dataset) * 0.6):int(len(dataset) * 0.8)]
    test = dataset[int(len(dataset) * 0.8):len(dataset)]
    m, n = np.max(train[train.columns[1]]), np.min(train[train.columns[1]])

    def process(data):
        load = data[data.columns[1]]
        data = data.values.tolist()
        load = (load - n) / (m - n)
        load = load.tolist()
        X, Y = [], []
        for i in range(len(data) - 24):
            train_seq = []
            train_label = []
            for j in range(i, i + 24):
                x = [load[j]]
                for c in range(2, 8):
                    x.append(data[j][c])
                train_seq.append(x)
            train_label.append(load[i + 24])
            # 展平
            train_seq = [x for y in train_seq for x in y]
            X.append(train_seq)
            Y.append(train_label)

        X, Y = np.array(X), np.array(Y)
        s = int(len(X) / B) * B
        X, Y = X[:s], Y[:s]

        return X, Y

    train_x, train_y = process(train)
    val_x, val_y = process(val)
    test_x, test_y = process(test)

    return [train_x, train_y], [val_x, val_y], [test_x, test_y], m, n

2.训练

def train(nn, Dtr, Val):
    print('training...')
    train_x, train_y = Dtr[0], Dtr[1]
    batch_size = nn.B
    min_epochs = 0
    best_nn = None
    min_val_loss = 5
    epochs = 50
    batch = int(len(train_x) / batch_size)
    train_loss = []
    for epoch in tqdm(range(epochs)):
        for i in range(batch):
            start = i * batch_size
            end = start + batch_size
            nn.forward_prop(train_x[start:end], train_y[start:end])
            nn.backward_prop(train_y[start:end])
            train_loss.append(np.mean(nn.loss))

        # val
        val_loss = get_val_loss(nn, Val)
        if epoch + 1 >= min_epochs and val_loss < min_val_loss:
            min_val_loss = val_loss
            best_nn = copy.deepcopy(nn)

        print('epoch {:03d} train_loss {:.8f} val_loss {:.8f}'.format(epoch, np.mean(train_loss), val_loss))

    return best_nn

3.测试

def test(nn, Dte, m, n):
    test_x, test_y = Dte[0], Dte[1]
    pred = []
    batch = int(len(test_y) / nn.B)
    for i in range(batch):
        start = i * nn.B
        end = start + nn.B
        res = nn.forward_prop(test_x[start:end], test_y[start:end])
        res = res.tolist()
        res = list(chain.from_iterable(res))
        # print('res=', res)
        pred.extend(res)
    pred = np.array(pred)
    test_y = test_y.flatten()
    pred = pred * (m - n) + n
    test_y = test_y * (m - n) + n
    print('mape:', get_mape(test_y, pred))