BP神经网络代码实战：从数学推导到Python实现的关键细节

1次阅读

共计 3575 个字符，预计需要花费 9 分钟才能阅读完成。

BP 神经网络通过多层非线性变换，能够自动提取图像中的层次化特征（如边缘→纹理→部件）。相比传统算法，它在 MNIST 手写识别等任务中可实现 >98% 的准确率。其核心价值在于端到端学习能力，无需人工设计特征提取器。

当使用 Sigmoid 激活函数时，其导数最大值为 0.25（当输入为 0 时），经过多层连乘后梯度呈指数衰减。例如 5 层网络的反向传播梯度可能仅为初始值的 0.0009（0.25^5）。

数学表现为：
$$\frac{\partial L}{\partial w_1} = \frac{\partial L}{\partial a_5}\cdot\sigma'(z_5)\cdot w_5 \cdot … \cdot \sigma'(z_1)\cdot x$$

学习率 η 的选择尤为关键：
– η=0.01 时可能收敛过慢（需要 500+epoch）
– η=0.1 时容易在最优值附近震荡
– η=1.0 典型导致损失爆炸(NaN)

全连接层的计算复杂度为 O(n³)，当隐藏层神经元达到 1000+ 时，单次迭代可能需要数秒。在 CPU 上训练 ResNet 级别的网络可能需要数月。

以 3 层网络为例，损失函数 L 对第一层权重 $w_1$ 的梯度计算过程：

前向传播：
$$z_2 = w_1^T x + b_1$$
$$a_2 = \sigma(z_2)$$
$$z_3 = w_2^T a_2 + b_2$$
$$a_3 = \sigma(z_3)$$
反向传播：
$$\frac{\partial L}{\partial w_2} = \frac{\partial L}{\partial a_3} \cdot \sigma'(z_3) \cdot a_2$$
$$\frac{\partial L}{\partial w_1} = \frac{\partial L}{\partial a_3} \cdot \sigma'(z_3) \cdot w_2 \cdot \sigma'(z_2) \cdot x$$

import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-5,5,100)

def sigmoid_grad(x):
    s = 1/(1+np.exp(-x))
    return s*(1-s)

def relu_grad(x):
    return np.where(x>0, 1, 0)

plt.plot(x, sigmoid_grad(x), label='Sigmoid')
plt.plot(x, relu_grad(x), label='ReLU')
plt.title('Activation Function Gradients')
plt.legend()
plt.show()

import numpy as np
from tqdm import tqdm

class DenseLayer:
    def __init__(self, n_input, n_output, activation='relu'):
        # Xavier 初始化
        self.W = np.random.randn(n_input, n_output) * np.sqrt(2/(n_input+n_output))
        self.b = np.zeros((1, n_output))
        self.activation = activation

    def forward(self, X):
        self.z = np.dot(X, self.W) + self.b
        if self.activation == 'relu':
            self.a = np.maximum(0, self.z)
        elif self.activation == 'sigmoid':
            self.a = 1/(1+np.exp(-self.z))
        return self.a

    def backward(self, X, grad_output):
        if self.activation == 'relu':
            grad_z = grad_output * (self.z > 0)
        elif self.activation == 'sigmoid':
            grad_z = grad_output * self.a * (1-self.a)

        grad_W = np.dot(X.T, grad_z)
        grad_b = np.sum(grad_z, axis=0, keepdims=True)
        grad_input = np.dot(grad_z, self.W.T)

        return grad_input, grad_W, grad_b

class NeuralNetwork:
    def __init__(self, layers):
        self.layers = layers

    def train(self, X, y, epochs=100, lr=0.01, batch_size=32):
        n_samples = X.shape[0]

        for epoch in tqdm(range(epochs)):
            # 学习率衰减
            current_lr = lr * (0.1 ** (epoch // 30))

            # Mini-batch
            indices = np.random.permutation(n_samples)
            for i in range(0, n_samples, batch_size):
                batch_idx = indices[i:i+batch_size]
                X_batch, y_batch = X[batch_idx], y[batch_idx]

                # 前向传播
                output = X_batch
                for layer in self.layers:
                    output = layer.forward(output)

                # 计算损失（假设使用 MSE）loss = np.mean((output - y_batch)**2)

                # 反向传播
                grad = 2*(output - y_batch)/batch_size
                for i in reversed(range(len(self.layers))):
                    layer = self.layers[i]
                    prev_output = X_batch if i==0 else self.layers[i-1].a
                    grad, grad_W, grad_b = layer.backward(prev_output, grad)

                    # 参数更新
                    layer.W -= current_lr * grad_W
                    layer.b -= current_lr * grad_b

对于第 l 层，初始化权重应满足：
$$W_{ij} \sim U\left(-\sqrt{\frac{6}{n_{in}+n_{out}}}, \sqrt{\frac{6}{n_{in}+n_{out}}}\right)$$

class BatchNormLayer:
    def __init__(self, n_features, eps=1e-5):
        self.gamma = np.ones((1, n_features))
        self.beta = np.zeros((1, n_features))
        self.eps = eps

    def forward(self, X, training=True):
        if training:
            self.mu = np.mean(X, axis=0)
            self.sigma2 = np.var(X, axis=0)
            X_norm = (X - self.mu)/np.sqrt(self.sigma2 + self.eps)
            self.cache = (X, X_norm)
            return self.gamma * X_norm + self.beta
        # 测试时使用移动平均
        else: 
            return self.gamma * (X - self.running_mu)/np.sqrt(self.running_sigma2 + self.eps) + self.beta

best_loss = float('inf')
patience = 5
counter = 0

for epoch in range(100):
    val_loss = evaluate(model, val_data)

    if val_loss < best_loss:
        best_loss = val_loss
        counter = 0
        # 保存最佳模型
    else:
        counter += 1
        if counter >= patience:
            print(f'Early stopping at epoch {epoch}')
            break

动量法 (Momentum) 需要维护速度变量 $v$，更新规则为：
$$v = \beta v + (1-\beta)\nabla_\theta J(\theta)$$
$$\theta = \theta – \alpha v$$
如何在现有代码中加入动量项？
PyTorch 的自动求导通过构建计算图实现，对比我们手动实现的优缺点是什么？当网络达到 100+ 层时，哪种方式更具优势？

正文完