手写BP网络拟合一元非线性函数(附数学原理公式推导)

说明：本文采用了一个三层 BP 网络来拟合 $sin(x)$ 和其他的一些简单函数，并且训练集和测试集都是生成的随机数，在拟合更加复杂的模型时需要进行相应的数据预处理、正则化或归一化操作，同时调整网络的超参。

一、BP 网络原理简介

什么是 BP 网络？

​ 人工神经网络可以通过自身的训练，寻找并学习训练集中给定输入值和期望输出值之间的关系，从而总结出某些规则，并以此在应用过程中根据用户的输入得到一个合适的输出值。这些规则可以是为我们所知的，也有可能是未知的，神经网络的训练过程是一个黑盒操作，我们通过损失函数的变化、网络在测试集上的表现等等来判断一个网络的优劣。

​ BP 网络通过前向计算、反向传播两个简单的过程，利用梯度下降方法，逐步更新网络中的参数矩阵，从而实现对函数的拟合等其他任务。下图是一个三层神经网络的基本结构。

• input layer：输入层，可以输入一个向量。
• hidden layer：隐藏层/中间层，用来扩展输入向量的维度，对输入向量中蕴含的信息进行重排与组合，并在后面的反向转播过程中改进自己的组合方式以学习到我们期望的规律。hidden layer 可以不止一层，但是层数过多可能会出现梯度消失等问题。
• output layer：输出层。

BP 网络数学原理推导

前向计算过程：

​ 在前向传播过程中，我们需要将输入的向量乘以一个变换矩阵来实现线性组合，该矩阵也可以称为权重矩阵，因为该矩阵中的数字代表了输入层中相应神经元蕴含的信息在变化过程中的重要程度。

​ 以 h₁ 为例，$h_1=x_1{w}_{11}+x_2{w}_{21}+x_3{w}_{31}+x_4{w}_{41}$，从这个式子中很容易看出 w 的权重作用。

​ 应用到所有神经元，我们有：

​ $\mathbf{h}=\mathbf{x}\cdot {\mathbf{w}}_1$

​ 由于该组合是线性的，这意味着神经网络中的数据与输入始终都是线性关系，这就无法实现非线性函数的拟合。因此我们需要加一个激活函数 sigmoid，用来将线性数据转换为非线性。同时在做 sigmoid 操作之前需要加上一个偏置向量 b1 来提高拟合模型的灵活性。

​ $sigmoid：\sigma (x)=\frac{1}{1+e^{-x}}$​

​ sigmoid 函数也具有一个很好的特性：${\sigma }^{{}^{\prime }}(x)=\sigma (x)[1-\sigma (x)]$，这一点将会为后续梯度下降过程中的计算带来便利。

​ 回到前向传播过程上来，我们有：

​ $\mathbf{s}=sigmoid(\mathbf{h}+{\mathbf{b}}_1)$

​ 对 hidden layer 做相似的操作，只不过这次由升维改为了降维：

​ $\mathbf{y}=\mathbf{s}\cdot {\mathbf{w}}_2+{\mathbf{b}}_2$

​ 注意到输出层不要 sigmoid，因为 sigmoid 会把所有的负值转换成正值，显然没法拟合 $sin(x)$ 这种有负值的函数。

​ 这样我们就得到了输出的向量。接下来就是根据输出与期望输出之间的差值来调整网络中的权重矩阵。

反向传播过程：

​ 如何量化网络的输出与我们期望输出之间的差值？我们可以用 MSE(均方误差) 来定义损失函数。

​ MSE：$loss=\Sigma \frac{1}{2}(d_i-y_i{)}^2$，其中 $d$ 为期望输出，$i$ 为输出层的一个神经元索引。

​ 故 ：

​ $\frac{\partial (loss)}{\partial \mathbf{y}}=\mathbf{y}-\mathbf{d}$

​ $\frac{\partial (loss)}{\partial {\mathbf{w}}_2}=$ $\frac{\partial (loss)}{\partial \mathbf{y}}\cdot \frac{\partial \mathbf{y}}{\partial {\mathbf{w}}_2}$

​ 又 $\frac{\partial \mathbf{y}}{\partial {\mathbf{w}}_2}=\frac{\partial }{\partial {\mathbf{w}}_2}({\mathbf{w}}_2\cdot \mathbf{s}+{\mathbf{b}}_2)=\mathbf{s}$

​ 所以 $\frac{\partial (loss)}{\partial {\mathbf{w}}_2}=\mathbf{s}\cdot (\mathbf{y}-\mathbf{d})$

​ 同理，

​ $\frac{\partial (loss)}{\partial {\mathbf{b}}_2}=\frac{\partial (loss)}{\partial \mathbf{y}}\cdot \frac{\partial \mathbf{y}}{\partial {\mathbf{b}}_2}=\mathbf{y}-\mathbf{d}$

​ $\frac{\partial (loss)}{\partial {\mathbf{w}}_1}=\frac{\partial (loss)}{\partial \mathbf{y}}\cdot \frac{\partial \mathbf{y}}{\partial \mathbf{s}}\cdot \frac{\partial \mathbf{s}}{\partial {\mathbf{w}}_1}=(\mathbf{y}-\mathbf{d})\cdot {\mathbf{w}}_2\cdot \mathbf{x}\cdot {\sigma }^{\prime }(\mathbf{s})=(\mathbf{y}-\mathbf{d})\cdot {\mathbf{w}}_2\cdot \mathbf{x}\cdot \sigma (\mathbf{s})(1-\sigma (\mathbf{s}))$

​ $\frac{\partial (loss)}{\partial {\mathbf{b}}_1}=\frac{\partial (loss)}{\partial \mathbf{y}}\cdot \frac{\partial \mathbf{y}}{\partial \mathbf{s}}\cdot \frac{\partial \mathbf{s}}{\partial {\mathbf{b}}_1}=(\mathbf{y}-\mathbf{d})\cdot {\mathbf{w}}_2$

由此我们得到了损失函数关于网络中几个超参的偏导，组合起来也就成为了所谓的梯度。现在我们要应用梯度下降法调整超参，使得损失函数朝着减小的方向进行。注意减小的方向是梯度的反向。

​ $\Delta {\mathbf{w}}_2=r\cdot \mathbf{s}\cdot (\mathbf{y}-\mathbf{d})$

​ $\Delta {\mathbf{b}}_2=r\cdot (\mathbf{y}-\mathbf{d})$​

​ $\Delta {\mathbf{w}}_1=r\cdot (\mathbf{y}-\mathbf{d})\cdot {\mathbf{w}}_2\cdot \mathbf{x}\cdot \sigma (\mathbf{s})(1-\sigma (\mathbf{s}))$

​ $\Delta {\mathbf{b}}_1=r\cdot (\mathbf{y}-\mathbf{d})\cdot {\mathbf{w}}_2$

其中 $r$ 为 learning_rate 。

二、应用 BP 网络拟合一元非线性函数

在代码我以 $sin(x)$​ 为例，拟合其他函数可以修改 target_function() 实现。

代码如下：

import numpy as np
from matplotlib import pyplot as plt


# 定义 sigmoid 函数
def sigmoid(x):
    return 1 / (1 + np.exp(-x))


# 定义 sigmoid 函数的导数
def sigmoid_derivative(x):
    return x * (1 - x)


# 定义目标函数（sin 函数）
def target_function(x):
    return np.sin(x)


# 定义神经网络参数
input_size = 1
hidden_size = 12
output_size = 1
learning_rate = 0.1


# 定义函数拟合过程
def funcFitting():
    np.random.seed(0) # 固定 random.seed() 可以使得每次生成的随机数都一样，方便在出错时寻找问题
    train_set = np.zeros((100, input_size))
    # 生成训练数据
    for i in range(len(train_set)):
        train_set[i] = np.random.uniform(-np.pi, np.pi)
    # 初始化权重和偏置
    w1 = np.random.uniform(-10, 10, size=(input_size, hidden_size))
    b1 = np.random.uniform(-10, 0, size=(1, hidden_size))
    w2 = np.random.uniform(-10, 10, size=(hidden_size, output_size))
    b2 = np.random.uniform(-10, 0, size=(1, output_size))

    for epoch in range(1000):
        for i in range(len(train_set)):
            x = train_set[i]
            x = np.array(x).reshape(1, input_size)
            d = target_function(x)
            h = x @ w1 + b1
            s = sigmoid(h)
            y = s @ w2 + b2

            e = y - d
            d_b2 = e
            d_w2 = s.T * e
            d_b = e @ w2.T * sigmoid_derivative(s)
            d_w1 = x.T * d_b

            b2 -= learning_rate * d_b2
            w2 -= learning_rate * d_w2
            b1 -= learning_rate * d_b
            w1 -= learning_rate * d_w1

        loss = np.mean(np.square(y - d))
        print("Epoch:", epoch, "Loss:", loss)

    test_set = np.zeros((10000, 1))
    correct_num = 0
    for i in range(len(test_set)):
        test_set[i] = np.random.uniform(-np.pi, np.pi)
        x_test = test_set[i]
        x_test = np.array(x_test).reshape(1, input_size)
        d_test = target_function(x_test)
        h_test = x_test @ w1 + b1
        s_test = sigmoid(h_test)
        y_test = s_test @ w2 + b2
        if np.mean(np.square(y_test - d_test)) < 0.01:
            correct_num += 1
    print("Accuracy:", correct_num / 100, "%")

    # 绘制拟合结果
    x_plot = np.linspace(-np.pi, np.pi, 100).reshape(-1, 1)
    h_plot = x_plot @ w1 + b1
    s_plot = sigmoid(h_plot)
    z_plot = s_plot @ w2 + b2
    d_plot = target_function(x_plot)
    plt.plot(x_plot, d_plot, label="Actual")
    plt.plot(x_plot, z_plot, label="Fitted")
    plt.legend()
    plt.show()


if __name__ == "__main__":
    funcFitting()

下面是一些函数的拟合结果以及对应的超参(由于是一元函数，input_size 和output_size 都为 1，范围均为$(-\pi ,\pi )$)：

• $sin(x)$：hidden_size：12，epoch：1000，learn_rate：0.1

  测试集正确率：89.59%，允许误差：$\pm 0.001$

• $x^2$：hidden_size：12，epoch：1000，learn_rate：0.1

  测试集正确率：99.67%，允许误差：$\pm 0.01$

  

• $sinh(x)$​：

  hidden_size：12，epoch：1000，learn_rate：0.1

  测试集正确率：96.53%，允许误差：$\pm 0.01$

  

• $e^x$:

  hidden_size：12，epoch：1000，learn_rate：0.1

  测试集正确率：93.33%，允许误差：$\pm 0.01$

  

• $x^4+2x^3-11x^2+18x-6$：

  hidden_size：12，epoch：1000，learn_rate：0.1

  测试集正确率：2.19%，允许误差：$\pm 2$​

  

  这种是落入了局部最小值，调参的时候应该减小学习率相应调整 epoch。

  hidden_size：12，epoch：2000，learn_rate：0.01

  测试集正确率：89.52%，允许误差：$\pm 2$​