从零手写卷积神经网络CNN:NumPy实现、im2col与反向传播全解析

发布时间:2026/10/14 7:52:38
从零手写卷积神经网络CNN:NumPy实现、im2col与反向传播全解析 把书翻到第7章的时候我已经在编辑器里和卷积神经网络缠斗了三个晚上。说“手磕”一点都不夸张——没有框架顺手调用的Conv2d和MaxPool2d只有NumPy和自己的一双手。这期是《深度学习入门基于Python的理论与实现》共读系列的第三篇目标非常明确不用任何深度学习框架从零实现一个可以训练的CNN卷积神经网络并且手写数字识别任务上达到98%以上的测试准确率。这篇分享不是书上的章节复述而是我实际敲完代码之后沉淀下来的理解为什么作者要用im2col把卷积转成矩阵乘法反向传播里dW、dx是怎么一步步求出来的手写Conv层和Pooling层时最容易在哪几个地方翻车。如果你已经能手写全连接网络和反向传播想进一步把卷积神经网络也吃透这篇应该能给你省下一大笔Debug时间。1. 为什么放着框架不用偏要手磕CNN1.1 从两层全连接网络到卷积的跨越系列前两篇里我已经把Affine层、ReLU层、SoftmaxWithLoss层像搭积木一样搭了起来也手写了一个能在MNIST上跑出92%左右准确率的两层网络。当时最大的成就感来自“形状流转”——输入(batch, 784)经过矩阵乘法变成(batch, 100)再过一次ReLU最后变成(batch, 10)每一层的forward和backward都清清楚楚。但拿全连接网络直接处理图片有一个很明显的问题一张28x28的灰度图把它拉平成784维向量之后像素之间的空间位置关系就全被打散了。左上角的像素和右下角的像素在全连接层眼里没有任何区别都是输入向量里互不相干的两个分量。这显然不符合人类理解图像的方式——我们识别一张猫的图片靠的是边缘、纹理、局部形状的层层组合而不是把每个像素孤立地看待。卷积神经网络就是为这个而生的。它用卷积核在图像上滑动每次只看一个局部区域再把所有局部响应组合起来形成更高层的特征。这个“局部连接参数共享”的设计让网络天然适合图像数据参数量也比全连接小得多。1.2 手写一遍能收获什么有人问我PyTorch一行nn.Conv2d(1, 32, 3)不就好了何必自己造轮子我的回答是调框架只能让你知道卷积是什么手写才能让你明白卷积为什么是卷积。手写过程中你必须回答这些平时完全不会想的问题卷积核的输出尺寸是怎么算的padding加在哪个维度上反向传播的时候一个像素被多个卷积窗口覆盖它的梯度是直接覆盖还是累加池化层没有参数它的梯度往哪儿传这些问题每一个都是真刀真枪的形状运算答不上来代码就报错。还有一个很实际的好处手写层之后再看框架的官方文档就能秒懂。框架文档里的in_channels、out_channels、kernel_size这些参数本质上就是你在构造函数里要保存的那些变量。有了底层实现经验你在调参的时候脑子里的画面是数据矩阵在层与层之间流动的动图而不是一个个黑盒。1.3 动手前先建立三条心智模型在写第一行代码之前我建议你先在脑子里种下三句话第一卷积层和池化层在结构上只是“新积木”。它们要做到和之前的Affine、ReLU完全一样的接口——forward(x)返回输出backward(dout)返回梯度。只要接口统一后面的组装、训练、梯度检查全都复用老代码。第二卷积运算的核心是“滑窗内积”。滤波器大小就是窗口大小窗口每移动一个步幅就对窗口内的像素和滤波器做一次点乘。Padding是为了控制输出尺寸和保留边缘信息。第三反向传播永远遵循“链式求导梯度累加”。如果一个输入元素参与了多个前向计算路径那它接收的梯度就是这些路径的梯度之和。理解了这个col2im为什么是“加回去”而不是“赋值回去”就顺理成章了。这三条里面最后一条最容易忽略我见过不少人在写Conv反向传播时栽在这上面。2. im2col整个CNN实现的钥匙2.1 卷积运算怎么被“掰”成矩阵乘法的直接按定义实现卷积并不难遍历每个输出位置、每个输入通道嵌套好几层for循环做内积就行。但这种写法在Python里慢到让人崩溃——我们之前手写Softmax、Affine的时候已经体会过NumPy向量化带来的速度差异卷积要是写成多重for循环一个批次的训练都要等半天。im2col的思路非常巧妙既然for循环慢那就想办法把“滑动窗口”的操作也变成矩阵运算。具体做法是把输入数据重新排列成一个二维矩阵。对每个卷积窗口覆盖的区域把它“拉平”成一行。假设有个滤波器大小是3x3、输入通道数1那每个窗口就展成9个数的一行所有窗口的行拼在一起。这样卷积核也展平成一列两者一乘就同时得到了所有窗口的内积结果。原来的循环遍历变成了一个np.dot速度直接起飞。比如输入是(N, C, H, W) (1, 1, 5, 5)滤波器3x3、stride1、pad0那输出尺寸就是(1, 1, 3, 3)一共9个窗口。im2col会把它变成(9, 9)的矩阵9个窗口、每个窗口9个数。2.2 im2col的代码实现这是整个手写CNN里最核心的工具函数我把注释写在代码里一起看import numpy as np def im2col(input_data, filter_h, filter_w, stride1, pad0): 把输入图像按卷积窗口展开成二维矩阵 input_data : 四维数组 (N, C, H, W) filter_h, filter_w : 滤波器窗口尺寸 return col : (N*out_h*out_w, C*filter_h*filter_w) N, C, H, W input_data.shape out_h (H 2 * pad - filter_h) // stride 1 out_w (W 2 * pad - filter_w) // stride 1 # 先paddingpad参数在四个方向都有效 img np.pad(input_data, [(0, 0), (0, 0), (pad, pad), (pad, pad)], constant) # 这个col张量的形状是 # (N, C, filter_h, filter_w, out_h, out_w) col np.zeros((N, C, filter_h, filter_w, out_h, out_w)) for y in range(filter_h): y_max y stride * out_h for x in range(filter_w): x_max x stride * out_w # 把某个通道上所有窗口的第(y,x)个位置的像素抠出来 col[:, :, y, x, :, :] img[:, :, y:y_max:stride, x:x_max:stride] # 重排成行优先的大矩阵 # 每行是一个窗口每列是窗口内某个位置的像素 col col.transpose(0, 4, 5, 1, 2, 3).reshape(N * out_h * out_w, -1) return col这里最常用的写法是保留两层for循环但只遍历滤波器内部的坐标而不是遍历整个图像。因为滤波器尺寸通常很小3x3、5x5这个循环开销可以忽略不计。输出尺寸的公式需要特别注意它和池化层的尺寸公式不同卷积层的输出通常希望保持空间分辨率通过加padding而池化层则是有意缩小分辨率。2.3 col2im反向传播时把梯度“撒回去”既然im2col是把图像变成大矩阵反向传播的时候就要把梯度矩阵还原成图像形状。这个函数叫col2im是im2col的逆操作。def col2im(col, input_shape, filter_h, filter_w, stride1, pad0): N, C, H, W input_shape out_h (H 2 * pad - filter_h) // stride 1 out_w (W 2 * pad - filter_w) // stride 1 # 先把二维矩阵reshape回六维张量 col col.reshape(N, out_h, out_w, C, filter_h, filter_w) col col.transpose(0, 3, 4, 5, 1, 2) # 注意这里创建img时多加了stride-1的余量防止切片越界 img np.zeros((N, C, H 2 * pad stride - 1, W 2 * pad stride - 1)) for y in range(filter_h): y_max y stride * out_h for x in range(filter_w): x_max x stride * out_w # 这里用 而不是 一个像素被多个窗口覆盖时梯度要累加 img[:, :, y:y_max:stride, x:x_max:stride] col[:, :, y, x, :, :] # 去掉padding部分 return img[:, :, pad:pad H, pad:pad W]为什么反向传播必须是而不是赋值因为im2col展开的时候一个输入像素会被多个窗口重复使用。比如stride1且滤波器大于1时边缘附近的输入像素几乎出现在每个窗口里。前向它是被“复制”到多个位置的反向这些位置的梯度就应该累加回原像素。这就是链式法则里“一条路径一个梯度多条路径梯度求和”的直接体现。2.4 内存换时间到底值不值im2col最大的缺点就是内存膨胀。一个28x28的灰度图3x3滤波器、stride1、pad1展开之后单张图的列数是784行乘以9列比原始图像多了近10倍的数据。放大到大一点的输入特征图内存占用会非常可观。但为了速度这笔账是划算的。矩阵乘法的实现高度优化比Python层级的for循环快出几个数量级。实际训练中MNIST上的SimpleCNN用CPU跑一个epoch也就十几秒量级完全可以接受。如果输入再大还可以考虑把im2col改成按批次分块处理或者直接上框架。3. 从零写Conv层和Pooling层3.1 Convolution类初始化与正向传播有了im2col卷积层的前向就变得出奇简单。构造函数保存滤波器权重W、偏置b、步幅和padding。这里的W形状是(FN, C, FH, FW)也就是输出通道数、输入通道数、滤波器高、滤波器宽。class Convolution: def __init__(self, W, b, stride1, pad0): self.W W self.b b self.stride stride self.pad pad self.x None self.col None self.col_W None self.dW None self.db None def forward(self, x): FN, C, FH, FW self.W.shape N, C, H, W x.shape out_h (H 2 * self.pad - FH) // self.stride 1 out_w (W 2 * self.pad - FW) // self.stride 1 col im2col(x, FH, FW, self.stride, self.pad) col_W self.W.reshape(FN, -1).T # 每个滤波器展成一列 out np.dot(col, col_W) self.b # 先排成 (N, out_h, out_w, FN) 再转成 (N, FN, out_h, out_w) out out.reshape(N, out_h, out_w, -1).transpose(0, 3, 1, 2) self.x x self.col col self.col_W col_W return out这里最容易看晕的一步是col_W self.W.reshape(FN, -1).T。self.W.reshape(FN, -1)把每个滤波器拉平得到形状(FN, C*FH*FW)转置后变成(C*FH*FW, FN)。这样组合成矩阵后col是行数为窗口数、列数为窗口元素数的矩阵col_W的行数与col的列数正好匹配结果就是每个窗口和每个滤波器的内积。最后reshape和transpose的顺序很关键。np.dot(col, col_W)的结果形状是(窗口数, FN)必须先reshape(N, out_h, out_w, FN)再transpose(0, 3, 1, 2)最终的输出才是我们习惯的(N, FN, out_h, out_w)。如果先transpose再reshape得到的就是错乱的数据。3.2 反向传播dW、db、dx是怎么算出来的反向传播是手写网络的重头戏。它要回答三个问题权重怎么更新偏置怎么更新输入梯度怎么继续往回调def backward(self, dout): FN, C, FH, FW self.W.shape # dout形状是 (N, FN, out_h, out_w) # 把它变成 (N*out_h*out_w, FN)与forward对称 dout dout.transpose(0, 2, 3, 1).reshape(-1, FN) # 偏置梯度对样本数和窗口数求和即可 self.db np.sum(dout, axis0) # 权重梯度col.T dout 把“哪个窗口用了哪个像素”的信息和梯度相乘 self.dW np.dot(self.col.T, dout) self.dW self.dW.transpose(1, 0).reshape(FN, C, FH, FW) # 输入梯度dout 乘 col_W 的转置把梯度映射回窗口元素 dcol np.dot(dout, self.col_W.T) dx col2im(dcol, self.x.shape, FH, FW, self.stride, self.pad) return dxdb为什么是简单求个和因为偏置是加在卷积结果上的它对所有输出位置共享。根据链式法则某个偏置的梯度就是对所有输出位置的梯度求和。dW的理解是“当前这个滤波器哪些窗口用了它那些窗口的输入是什么梯度是多少。”col.T dout本质上是一次把所有窗口的贡献一次性聚合的运算。窗口数可能是几千几万但NumPy一个矩阵乘法就搞定了。dx是最有意思的。它先把输出梯度映射回窗口元素的维度得到dcol再用col2im把梯度撒回原始图像的每个像素。因为一个像素被多个窗口共用所以这一步必须累加。3.3 Pooling层没有参数反而更需小心池化层没有需要学习的参数它做的只是在每个窗口内取最大值Max Pooling。前向看起来简单反向传播却有个微妙的点最大值的位置要记住反向传播时把梯度只传给那个位置窗口内其他位置梯度为0。class Pooling: def __init__(self, pool_h, pool_w, stride1, pad0): self.pool_h pool_h self.pool_w pool_w self.stride stride self.pad pad self.x None self.arg_max None def forward(self, x): N, C, H, W x.shape out_h (H - self.pool_h) // self.stride 1 out_w (W - self.pool_w) // self.stride 1 col im2col(x, self.pool_h, self.pool_w, self.stride, self.pad) col col.reshape(-1, self.pool_h * self.pool_w) self.arg_max np.argmax(col, axis1) out np.max(col, axis1) out out.reshape(N, out_h, out_w, C).transpose(0, 3, 1, 2) self.x x return out def backward(self, dout): N, C, H, W self.x.shape dout dout.transpose(0, 2, 3, 1) pool_size self.pool_h * self.pool_w # 先全零再把最大值位置的梯度填进去 dmax np.zeros((dout.size, pool_size)) dmax[np.arange(self.arg_max.size), self.arg_max.flatten()] dout.flatten() dmax dmax.reshape(dout.shape (pool_size,)) dcol dmax.reshape(-1, pool_size) dx col2im(dcol, self.x.shape, self.pool_h, self.pool_w, self.stride, self.pad) return dx注意池化层的输出尺寸公式没有2*pad项因为实际使用中池化层几乎不做padding。如果给Pooling传入pad参数也可以支持但尺寸公式要相应调整最好保持“池化从不padding”的约定这样不容易出错。4. 组装SimpleCNN在MNIST上跑起来4.1 网络结构长什么样真正把层搭成网络的时候我之前写的那些工具类全都能复用只需要像拼积木一样按顺序组装。这里的思路是通过组合通用层类来构建一个端到端模型。class SimpleCNN: def __init__(self): # 第一层卷积1通道 - 32通道3x3pad1 self.conv1 Convolution(he_init(1, 32, 3), np.zeros(32), stride1, pad1) self.relu1 ReLU() self.pool1 Pooling(2, 2, stride2) # 第二层卷积32通道 - 64通道 self.conv2 Convolution(he_init(32, 64, 3), np.zeros(64), stride1, pad1) self.relu2 ReLU() self.pool2 Pooling(2, 2, stride2) # 全连接部分 self.affine1 Affine(he_init(64 * 7 * 7, 100), np.zeros(100)) self.relu3 ReLU() self.affine2 Affine(he_init(100, 10), np.zeros(10)) self.layers [ self.conv1, self.relu1, self.pool1, self.conv2, self.relu2, self.pool2, self.affine1, self.relu3, self.affine2 ] self.last_layer SoftmaxWithLoss() def forward(self, x): for layer in self.layers: x layer.forward(x) return x def predict(self, x): # 预测时不能过SoftmaxWithLoss它把loss也算进来了 y self.forward(x) return y def loss(self, x, t): y self.forward(x) return self.last_layer.forward(y, t) def backward(self, dout1): dout self.last_layer.backward(dout) for layer in reversed(self.layers): dout layer.backward(dout) return dout这里的he_init是权重的初始化函数。全连接网络里我们一般用较小的随机数初始化但卷积层和较深的网络需要更讲究。我用的是一种常见的均匀初始化方式标准差按输入尺寸倒数的平方根缩放。实际效果比固定小标准差要稳定很多尤其是网络层数一多浅了显不出来深了立刻见分晓。4.2 各层尺寸的变化轨迹手写代码调试最痛苦的就是“形状对不上”。我在设计网络时会把每一层的输出尺寸提前算清楚写下来层输入形状输出形状关键参数Conv1(N, 1, 28, 28)(N, 32, 28, 28)3x3, pad1, stride1ReLU1(N, 32, 28, 28)不变-Pool1(N, 32, 28, 28)(N, 32, 14, 14)2x2, stride2Conv2(N, 32, 14, 14)(N, 64, 14, 14)3x3, pad1, stride1ReLU2(N, 64, 14, 14)不变-Pool2(N, 64, 14, 14)(N, 64, 7, 7)2x2, stride2Affine1(N, 6477)(N, 100)-Affine2(N, 100)(N, 10)-这张表在写代码之前就应该画好。尤其是池化层把28x28变成14x14、再把14x14变成7x7这两步特征图尺寸的折半意味着后续全连接层的输入维度从64*14*14变成64*7*7差了好几倍。网络设计时如果不先算这些数字写Affine层初始化的时候一定会出错。4.3 训练循环与实测表现训练循环和之前全连接网络几乎一样我只需要替换数据集读取和网络实例化的部分def train(net, x_train, t_train, x_val, t_val, batch_size256, epochs20, lr0.1): optimizer SGD(lr) # 之前实现过的简单更新器 train_size x_train.shape[0] for epoch in range(epochs): # 打乱顺序 indices np.random.permutation(train_size) for i in range(0, train_size, batch_size): batch_idx indices[i:i batch_size] x_batch x_train[batch_idx] t_batch t_train[batch_idx] loss net.loss(x_batch, t_batch) net.backward() optimizer.update(net) # 对每一层执行 W - lr * dW val_acc evaluate(net, x_val, t_val) print(fepoch {epoch1}: val_acc {val_acc:.4f})实测下来这个结构在第4到第6个epoch时验证准确率就能过97%到第12个epoch左右稳定在98%以上。相比之前只有两层Affine的全连接网络CNN在同样的训练轮数下准确率高出近6个百分点这就是卷积结构对图像数据的天然优势。我用的优化器是最朴素的SGD学习率0.1。没有用动量也没有用Adam。这个网络结构在这样的超参数下非常稳定loss曲线平滑下降没有出现震荡。你如果自己实现的时候发现准确率卡在某个值上不去优先检查反向传播的梯度是否正确而不是调整学习率。5. 手写过程中踩过的坑与排查链路5.1 第一个坑形状错乱全在reshape和transpose的顺序手写卷积层的第一天我遇到最多的问题就是ValueError: shapes ... not aligned。排查下来十有八九是reshape和transpose顺序搞反了。前面说过np.dot的结果是(窗口数, FN)它对应的物理含义是“每一行是一个窗口每一列是该窗口在第几个滤波器上的响应”。要把这个矩阵变成(N, FN, out_h, out_w)必须先恢复出空间结构再调通道维正确out.reshape(N, out_h, out_w, -1).transpose(0, 3, 1, 2)错误out.transpose(0, 3, 1, 2).reshape(N, FN, out_h, out_w)——先转置再reshape数据顺序全乱了。这种错误不会报错只会让网络训练时loss死活不降。我的排查方法是写一个很小的“冒烟测试”构造随机输入手算一遍某个固定位置的卷积响应值和代码输出对比。这个测试只要跑过一次就能帮你确认整个前向链路是通的。5.2 第二个坑梯度检查是唯一的救命稻草手写反向传播光靠“我觉得没错”是不行的。我强烈建议你给每一层单独做一次数值梯度验证def check_grad(net, x, t): # 求解析梯度 net.loss(x, t) net.backward() # 对每个带参数的层计算数值梯度 for layer in net.layers: if hasattr(layer, W): dw_numerical compute_numerical_gradient( lambda W: net.loss(x, t), layer.W) diff np.abs(dw_numerical - layer.dW).max() assert diff 1e-7, f{type(layer).__name__} gradient error: {diff}数值梯度的原理是直接按导数的定义来(f(Wh) - f(W-h)) / (2h)h取极小值比如1e-6。这个方法很慢只在调试时对小规模输入跑一次。compute_numerical_gradient可以先对单个层、单个随机输入验证再对整个组合网络验证。分层验证的好处是定位精准——哪一层梯度有问题立刻就能知道。全部通过之后网络训练不收敛的概率就非常低了。5.3 第三个坑池化层反向传播的“稀疏梯度”不好想刚开始写Pooling backward的时候我一度试图把平均池化的梯度回传逻辑套到Max Pooling上导致loss在几百轮迭代后出现NaN。后来才明白Max Pooling的反向传播本质上是一个“选座位”的游戏。前向的时候记住每个窗口最大值的下标反向的时候构造一个和窗口一样大的向量只有最大值下标那个位置是梯度值其余全是0。因为最大值只有一个梯度被“路由”到了那一个输入元素上其他没有选中的元素在本次前向中对loss没有贡献梯度自然为0。这个逻辑本身不难但写代码时容易在一个地方翻车dmax构造完之后必须reshape成(N, out_h, out_w, C, pool_h, pool_w)这样的六维结构才能正确调用col2im。维度顺序一旦弄错col2im返回的dx和输入x形状不一致又是一通好找。5.4 第四个坑初始化不好训练直接死掉如果在SimpleCNN上把he_init换成固定标准差0.01的随机初始化你可能会看到loss在最初几百步疯狂下跌后突然变成NaN。原因很简单卷积层数量虽然不多但每个卷积核要处理的局部区域是3x39个数的内积一组参数不合适经过两层卷积再叠加ReLU数值很容易爆炸。我给Conv层和Affine层统一使用按输入维度缩放的初始化效果非常稳定。如果你用的是ReLU这类激活函数初始化标准差取“输入通道数的平方根倒数”是一个很稳妥的选择。偏置则一律初始化为0不用耍花样。5.5 第五个坑手写代码运行慢先从im2col找原因第一个完整版本跑起来的时候一个epoch要将近一分钟比预想慢不少。看profile发现瓶颈完全不在卷积计算而在im2col的padding操作上。我最初用的是np.zeros新开数组再手动拷贝改成np.pad之后速度提升非常明显。另外col矩阵本身占用的内存也比想象中大。如果网络输入尺寸再翻倍内存直接吃不消。我的处理方案是训练时保持原始MNIST图片尺寸不放大如果你要自己跑更大分辨率的图片建议在im2col内部按批次切片或者把展开后的矩阵分块做矩阵乘法别一次性全放进来。6. 手磕完之后的下一步玩法网络在MNIST上稳定98%之后我做了件很解压的事把第一层32个卷积核打印出来可视化。你会发现它们学到的全都是边缘、斜线、角点之类的基层纹理就像把图像的基本笔画拆了出来一样。我强烈建议你也看看——那是你写的代码自己生长出来的“视觉词汇”比任何理论描述都直观。再往下还可以尝试加BatchNorm层、调整网络深度、引入dropout防止过拟合。手写版本的好处是你对每个改动的影响都有体感——加了BatchNorm之后训练初期loss下降速度明显变化你会真实看到“为什么现代网络几乎都标配这个层”。我现在已经把手写CNN的整个模块化结构保留下来了后面继续看更深的网络时只需要在积木盒里再添新层就行。最后说一点个人体会手磕CNN这件事真正难的从来不是把前向和反向的公式背下来而是当你面对一堆ValueError和死活不掉的loss时愿意一层一层拆开去验证。我的排查顺序永远是先分层的数值梯度验证再跑小数据集最后才上全量数据训练。这套顺序看着笨但能帮你把所有可能的锅都提前排掉。如果你也想彻底弄懂卷积神经网络别急着上框架照着这个思路自己写一遍哪怕只是MNIST级别的小网络收获也会比看十遍教程大得多。