手写BP神经网络逼近二元函数：NumPy实现与调参实践

简介：一份基于反向传播算法的BP神经网络逼近二元函数的MATLAB源码，主要面向机器学习初学者、神经网络课程学生以及希望脱离工具箱手动实现算法的开发者。该资源完全依靠自编代码完成网络训练，不依赖任何神经网络工具箱，从网络结构定义、随机初始化、前向计算、误差反向传播到梯度下降更新权重均逐行实现，能够帮助读者看清BP算法每一步的数学本质。压缩包采用zip格式打包，文件总数仅1个m脚本文件，大小约1KB，代码紧凑，适合在MATLAB中直接打开、研读和调试。目前已有188人学习浏览，受到一定关注。运行该程序后，可以直观观察网络如何逐步逼近目标二元函数，同时掌握激活函数选择、均方误差计算、学习率调整等关键环节，非常适合作为入门深度学习底层原理的实操练习。

1. 为什么手写BP网络逼近二元函数：不依赖工具箱的场景与价值

用BP神经网络逼近二元函数，听起来是MATLAB里几行代码就能结束的任务，newff、train、sim三个函数一调就出结果。但真到了项目现场，经常会遇上两种情况：一是MATLAB环境里没有可用的神经网络工具箱授权；二是训练好的网络要移植到C/C++、嵌入式设备或Web后端，那套依赖工具箱的代码完全带不走。手写BP网络，就是不使用任何神经网络工具库，把前向传播、反向传播、梯度更新每一行都自己写出来的实现方式。这既是经典的算法基本功，也是做算子级调优和跨平台部署的必要前提。今天的实现完全基于纯NumPy，从零逼近z = f(x, y)这样的二元函数，并把采样、训练、验证和调参中容易踩的坑一并讲清楚。

2. BP网络逼近二元函数的数学基础：结构、激活与梯度推导

2.1 万能逼近定理给到我们的结构自由度

二元函数逼近本质上是回归问题：输入(x, y)这个二维向量，输出连续变化的z。BP网络能胜任这件事，理论基础是万能逼近定理——只要隐藏层使用非线性激活函数，单隐藏层前馈网络就能以任意精度逼近定义在紧集上的连续函数。这个定理保证的是“存在性”，说明从原理上这条路走得通，但它并不告诉我们隐藏层到底要放多少个节点。

网络结构在输入和输出两端是固定的：输入层2个节点对应x和y，输出层1个节点对应z。可变的部分全在中间。对z = sin(x)·cos(y)这类平滑函数，单隐藏层8到16个节点就够用；如果目标函数带尖角、突变或高频振荡，就需要更多节点甚至加深层数。常见做法是先锁住单隐藏层，把节点数从4、8、16、32往上试，同时观察验证误差是否持续改善。一旦验证误差不再下降，就说明结构容量已经够用，继续加节点只是在拟合训练数据里的噪声。

单看BP神经网络结构图，一个2-16-1的回归预测结构非常典型：输入层两个节点，隐藏层16个tanh神经元，输出层1个线性神经元。这张图虽然简单，却撑起了浅层函数逼近任务里绝大多数场景。

提示：万能逼近定理只给出存在性，节点数与精度的定量关系没有显式公式，工程上只能靠试。

2.2 激活函数选择与导数表达式

隐藏层激活函数常用sigmoid和tanh。sigmoid输出在(0, 1)之间，优点是平滑可导，缺点是均值不为0，梯度更新时会引入系统性偏移，收敛偏慢。tanh输出在(-1, 1)之间，均值为0，回归场景下的收敛速度通常明显快于sigmoid。输出层激活函数必须是线性的，因为我们要逼近任意实数值的z，任何有界激活都会把输出压在某个区间里，函数的峰值和谷值都表达不出来。

手写网络时必须同时准备激活函数和它的导数，因为反向传播需要用到：

import numpy as np def sigmoid(x): return 1 / (1 + np.exp(-x)) def sigmoid_derivative(a): return a * (1 - a) def tanh(x): return np.tanh(x) def tanh_derivative(a): return 1 - a ** 2

这里有一个非常容易写错的地方：sigmoid_derivative和tanh_derivative接收的参数，是前向传播缓存下来的激活值a，而不是原始输入x。链式法则要求计算的是dσ(z)/dz，而σ(z)恰好等于a，所以a * (1 - a)和1 - a²才是正确的。如果误把x传进去，导数曲线在数值上对不上，梯度方向会逐渐偏离，而且损失还在下降，属于最难排查的隐性问题。

2.3 反向传播的矩阵形式推导

损失函数取均方误差L = (1/(2m))·Σ(z_pred - z_true)²，多写一个1/2是为了让求导后消掉系数。设输入矩阵X为m×2，隐藏层权重W1为2×h，偏置b1为1×h，输出层权重W2为h×1，偏置b2为1×1。前向传播是：

Z1 = X·W1 + b1，A1 = tanh(Z1) Z2 = A1·W2 + b2，A2 = Z2

反向传播按链式法则展开，先计算输出层误差δ2 = (A2 - Y)/m，再逐层回传：

∂L/∂W2 = A1ᵀ·δ2 ∂L/∂b2 = Σδ2（按batch求和） δ1 = (δ2·W2ᵀ) ⊙ (1 - A1²) ∂L/∂W1 = Xᵀ·δ1 ∂L/∂b1 = Σδ1

这里⊙表示逐元素乘法。这套公式看上去不复杂，真正实现时最容易在δ1的计算上出错：δ2的维度是m×1，W2的维度是h×1，要得到m×h的中间矩阵，必须用δ2·W2ᵀ而不是W2·δ2ᵀ，方向错了维度立刻抛异常，方向对了但顺序错了则会出现静默的错误结果。理解到这一层，后面的代码就是顺理成章的事。

3. 纯NumPy实现BP网络：前向传播、反向传播与最小训练框架

3.1 网络类定义与权重初始化

把网络封装成一个类，权重、偏置全部作为成员变量保存。这样训练结束后可以单独提取W1、W2做推理，也可以把参数序列化保存到文件，供其他程序加载。初始化权重用随机数，但幅度必须控制：随机数太大会让隐藏层神经元一进入训练就饱和，导数接近零，梯度消失；太小则梯度信号微弱，收敛极慢。我习惯用np.random.randn生成标准正态随机数，再乘一个缩放因子。

import numpy as np class BPNN: def __init__(self, input_size=2, hidden_size=16, output_size=1, lr=0.1): self.lr = lr self.W1 = np.random.randn(input_size, hidden_size) * np.sqrt(2.0 / input_size) self.b1 = np.zeros((1, hidden_size)) self.W2 = np.random.randn(hidden_size, output_size) * np.sqrt(2.0 / hidden_size) self.b2 = np.zeros((1, output_size))

np.sqrt(2.0 / input_size)是Xavier/He初始化的简化版本。input_size为2时，W1的标准差约1.0；W2的缩放因子是np.sqrt(2.0 / hidden_size)，hidden_size为16时约0.35。这样逐层传播后，输出方差既不会爆炸也不会消失，训练起点比较健康。偏置初始化为零是通用做法，因为偏置的梯度更新本身就能让网络学到合理的平移量。

3.2 前向传播与反向传播实现

前向传播要做的事很直接：线性变换、激活、输出层线性变换，同时把中间量缓存下来。反向传播则根据缓存计算各层梯度和更新量，最后统一减去学习率乘以梯度的值。

def forward(self, X): self.X = X self.Z1 = np.dot(X, self.W1) + self.b1 self.A1 = np.tanh(self.Z1) self.Z2 = np.dot(self.A1, self.W2) + self.b2 self.A2 = self.Z2 return self.A2 def backward(self, y): m = y.shape[0] dZ2 = (self.A2 - y) / m dW2 = np.dot(self.A1.T, dZ2) db2 = np.sum(dZ2, axis=0, keepdims=True) dA1 = np.dot(dZ2, self.W2.T) dZ1 = dA1 * (1 - self.A1 ** 2) dW1 = np.dot(self.X.T, dZ1) db1 = np.sum(dZ1, axis=0, keepdims=True) self.W1 -= self.lr * dW1 self.b1 -= self.lr * db1 self.W2 -= self.lr * dW2 self.b2 -= self.lr * db2

三个细节需要单独说明。第一，dZ2 = (self.A2 - y) / m把损失对输出层加权输入的导数整体算好，后面dW2和db2直接使用dZ2，而不需要再除一次m。第二，tanh导数1 - A1²使用的是self.A1，这个值来自最近一次forward的缓存。如果backward之前忘了调用forward，self.A1是上一次训练迭代的残留，梯度计算会错乱且没有任何报错提示。第三，db1和db2用np.sum对m个样本求和，因为偏置对每个样本共享，梯度是各样本梯度的和，和函数前面的1/m配合正好是均值。

3.3 训练循环与超参数速查

有了forward和backward，训练循环就非常简单：每次迭代先forward得到预测值，再计算MSE损失，然后backward更新权重，同时记录损失值用于画收敛曲线。

def train(self, X, y, epochs=10000, verbose=1000): losses = [] for epoch in range(epochs): pred = self.forward(X) loss = np.mean((pred - y) ** 2) losses.append(loss) self.backward(y) if epoch % verbose == 0: print(f"epoch {epoch}, loss {loss:.6f}") return losses

上面这个循环是标准的全量梯度下降：每一轮把所有样本都算一遍，取平均梯度更新一次。样本量在几千个以内时完全可行，梯度无偏，收敛路径稳定。如果样本量到几万以上，就要考虑mini-batch，具体改法放在最后一章。

3.4 训练后的推理与权重导出

训练完成后，网络已经是一组确定的W1、b1、W2、b2。推理时只需要调用forward，不需要再计算梯度。参数保存用NumPy自带的np.savez，加载后重新建一个网络对象并赋值即可：

np.savez("bp_params.npz", W1=nn.W1, b1=nn.b1, W2=nn.W2, b2=nn.b2) loaded = np.load("bp_params.npz") nn2 = BPNN() nn2.W1, nn2.b1 = loaded["W1"], loaded["b1"] nn2.W2, nn2.b2 = loaded["W2"], loaded["b2"]

这样的设计为后续跨语言部署留好了接口，C++或Java端只需要按同样格式读取npz文件里的矩阵，就能脱离Python环境执行推理。

4. 训练与验证：数据采样、损失收敛与逼近效果评估

4.1 二元函数采样与数据预处理

以z = sin(x) · cos(y)作为待逼近的二元函数，定义域取[-2, 2] × [-2, 2]。为了让网络有足够样本学习这个光滑但非线性的函数面，用np.linspace生成40个等间距点，再用np.meshgrid构成网格，最后展平成坐标对：

n = 40 x_lin = np.linspace(-2, 2, n) y_lin = np.linspace(-2, 2, n) X_grid, Y_grid = np.meshgrid(x_lin, y_lin) X = np.stack([X_grid.ravel(), Y_grid.ravel()], axis=1) y = (np.sin(X_grid) * np.cos(Y_grid)).ravel().reshape(-1, 1)

这里有一个新手必踩的坑：np.meshgrid返回的X_grid和Y_grid各是40×40的网格矩阵，必须用ravel()展平成一位数组，再按列堆叠，才能得到形状为(1600, 2)的样本矩阵。如果直接把二维网格矩阵喂给网络，形状必然对不上。y用reshape(-1, 1)把一维数组变成列向量，因为网络输出是m×1，MSE计算和反向传播都要求y和A2保持相同维度。

数据准备好以后，可以在这时做一次输入归一化。当前定义域[-2, 2]并不极端，tanh在输入±2时还未饱和，但如果目标函数定义域扩大到[-10, 10]，不做归一化的网络几乎必然出现收敛缓慢或震荡。标准做法是计算训练集的均值和标准差，把输入标准化为均值为0、方差为1的分布。

4.2 训练参数设置与损失曲线解读

设置lr=0.1、hidden_size=16、epochs=10000，开始训练。观察损失曲线会看到典型的两个阶段：前几百次迭代损失从0.2量级快速下降到0.01量级，这是网络在拟合函数面的低频主体；之后进入缓慢下降阶段，逐步打磨峰谷细节。

nn = BPNN(input_size=2, hidden_size=16, output_size=1, lr=0.1) losses = nn.train(X, y, epochs=10000, verbose=2000) pred = nn.forward(X) mae = np.mean(np.abs(pred - y)) max_err = np.max(np.abs(pred - y)) print(f"MAE={mae:.4f}, MaxErr={max_err:.4f}")

用16个隐藏节点训练10000次，MAE通常在0.01到0.03之间，最大误差出现在函数值变化最剧烈的区域。sin(x)cos(y)在定义域角落附近变化平缓、精度高，在中心交叉处变化快、误差偏大。这种误差分布是逼近任务本身的特性，不是实现bug。

训练过程中如果看到损失在第几千次迭代后开始周期性震荡，第一反应应该是降低学习率，从0.1调到0.03再试。学习率过大时，每一步权重更新都跨过了损失函数的谷底，损失会在低数值区间反复横跳。若损失直接变成NaN，多半是学习率过大导致权重发散，或者输入数据里存在极端值。

4.3 从误差热力图判断逼近瓶颈

只看损失曲线无法知道逼近质量在空间上的分布，必须把误差画成热力图。网格点上的误差定义为预测值减去真值，正负都有。如果误差在定义域边缘区域扩大，说明样本点稀疏，网络外推能力弱；如果误差在峰谷区域集中，说明隐藏层节点数不够，函数曲率变化快的地方拟合不过来。

把误差热力图和损失曲线放在一起，可以分辨三类典型问题。第一，所有位置误差都偏大，说明结构容量不足或训练未收敛，加节点和加训练轮次都值得试。第二，只有局部区域误差大，说明该区域需要加密采样，或者网络对高频成分的建模能力不够，考虑增加节点数。第三，误差在正负之间大幅度振荡，说明学习率偏高或网络开始过拟合训练数据。这三条判断规则，同样适用于后续要逼近的其他二元函数。

5. 调参策略与复杂二元函数逼近的进阶技巧

5.1 从平滑函数到驼峰函数：节点数与学习率的匹配

换一个目标函数，比如z = exp(-(x² + y²))，在原点有一个高耸的峰值，曲率变化远大于sin(x)cos(y)。沿用原来的参数训练会发现最大误差明显增大，因为峰值附近曲率太陡，16个tanh节点刻画起来很吃力。把隐藏层增大到32，峰值区域误差出现肉眼可见的改善，训练时间也同步翻倍。

这张表能看出两点：节点数翻倍带来的精度改善在递减，从16到32，最大误差只降低0.02，训练开销却翻倍；学习率过大时整个训练直接发散，说明这个函数的损失面对学习率更敏感，节点数增加后尤其如此。

5.2 输入归一化与反归一化：标准做法

当定义域扩展到[-10, 10]，直接训练很难收敛。标准做法是保存训练集的均值和标准差，把输入标准化后送入网络，预测出结果以后再按输出层的均值和标准差还原：

X_mean, X_std = X.mean(axis=0), X.std(axis=0) y_mean, y_std = y.mean(), y.std() X_norm = (X - X_mean) / X_std y_norm = (y - y_mean) / y_std # 还原预测值 # pred_original = pred * y_std + y_mean

这里需要特别注意：预测阶段对新样本归一化时，必须使用训练集保存的X_mean和X_std，而不是新样本自己的均值。很多人在训练完接入在线服务时漏掉这一步，导致预测结果与训练阶段不在同一个坐标下，误差完全不可控。

5.3 加深层数与mini-batch训练的切入点

单隐藏层逼近复杂二元函数时，需要大量节点才能刻画高频细节。改成双隐藏层后，网络在参数量相近的情况下往往效果更好，因为第一层可以先提取局部结构，第二层再组合成完整曲面。实现双隐藏层时，反向传播只需要在原有基础上多回传一层delta，代码量增加很少，但初始学习率建议从0.1降到0.05，否则容易震荡。

另一个重要升级是把全量梯度下降改成mini-batch随机梯度下降。常见做法是先把样本随机打乱，按batch_size=32切分，每批计算一次梯度并更新权重，跑完一轮后再打乱。与全量梯度相比，mini-batch的梯度带有随机性，能帮助网络跳出局部极小点；批次越小随机性越强，批次越大收敛越稳定。训练时从batch_size=32开始调，结合损失曲线决定是加大还是缩小，这是手写BP网络从教学代码走向实用工具的关键一步。

本文还有配套的精品资源，点击获取