基于numpy的CNN实现 进行MNIST手写数字识别

主要框架来自于这篇文章：/qq_36393962/article/details/99354969，下面会以原文来代称这篇文章

本文在原文的基础上增加了交叉熵以及mnist数据集的处理，并对卷积层反向传播更新权重的代码进行了一些相应的改动。建议先按原文脉络学习一下相关理论知识（原文讲的非常完整详细，但是个人感觉卷积层反向传播更新参数那部分有问题，卷积层更新参数的原理可参考刘建平前辈的文章），再参考下述实现代码。

文章目录

零、MNIST数据集读取+layer模板代码一、卷积层代码二、激活层代码三、池化层代码四、flatten层代码五、全连接层代码六、神经网络代码七、训练与测试八、模型保存与加载总结

零、MNIST数据集读取+layer模板代码

数据集读取直接用的这篇文章的方法1

import numpy as npimport structdef load_images(file_name):## 在读取或写入一个文件之前，你必须使用 Python 内置open()函数来打开它。#### file object = open(file_name [, access_mode][, buffering])#### file_name是包含您要访问的文件名的字符串值。#### access_mode指定该文件已被打开，即读，写，追加等方式。#### 0表示不使用缓冲，1表示在访问一个文件时进行缓冲。#### 这里rb表示只能以二进制读取的方式打开一个文件 ##binfile = open(file_name, 'rb') ## 从一个打开的文件读取数据buffers = binfile.read()## 读取image文件前4个整型数字magic,num,rows,cols = struct.unpack_from('>IIII',buffers, 0)## 整个images数据大小为60000*28*28bits = num * rows * cols## 读取images数据images = struct.unpack_from('>' + str(bits) + 'B', buffers, struct.calcsize('>IIII'))## 关闭文件binfile.close()## 转换为[60000,784]型数组images = np.reshape(images, [num, rows , cols,1])images=images.astype(float)return imagesdef load_labels(file_name):## 打开文件binfile = open(file_name, 'rb')## 从一个打开的文件读取数据 buffers = binfile.read()## 读取label文件前2个整形数字，label的长度为nummagic,num = struct.unpack_from('>II', buffers, 0) ## 读取labels数据labels = struct.unpack_from('>' + str(num) + "B", buffers, struct.calcsize('>II'))## 关闭文件binfile.close()## 转换为一维数组labels = np.reshape(labels, [num])# labels=labels.astype(float)return labels filename_train_images = 'train-images-idx3-ubyte'filename_train_labels = 'train-labels-idx1-ubyte'filename_test_images = 't10k-images-idx3-ubyte'filename_test_labels = 't10k-labels-idx1-ubyte'train_images=load_images(filename_train_images)train_labels=load_labels(filename_train_labels)test_images=load_images(filename_test_images)test_labels=load_labels(filename_test_labels)max_pixel=255train_images=train_images/255test_images=test_images/255

下面是每一层网络的模板代码：

from abc import abstractmethod# 网络层，正向传播相应方法为__call__，反向传播相应方法为backwardclass Layer:@abstractmethoddef __call__(self,x:np.ndarray)->np.ndarray:pass@abstractmethoddef backward(self,d:np.ndarray)->np.ndarray:pass

一、卷积层代码

卷积层反向传播时更新卷积核的参数时，我没办法理解原文说的以误差为矩阵，转置后的输入矩阵为卷积核的卷积操作，而且我用原文的代码跑出来的准确率也很糟糕：同样是两个epoch，按原文的写法测试集准确率只有38%，按我下面的写法准确率可以达到70%。

所以我按博客园上的文章，按如下格式进行了改写：

因此，与原文相比，我的代码主要是下面这块改了改：

from typing import Optional,Iterable# 卷积层class ConvolutionLayer(Layer):# 初始化def __init__(self,in_channels,out_channels,kernel_size,stride):self.in_channels=in_channelsself.out_channels=out_channelsself.ks=kernel_sizeself.kernel_shape=(kernel_size,kernel_size)self.stride=strideself.x:Optional[np.ndarray]=None# 卷积核：row,col,channel# 共有out_channels个(row,col,in_channels)的kernelsself.kernel=np.random.normal(loc=0.0,scale=1.0,size=(kernel_size,kernel_size,in_channels,out_channels,))# 每个卷积核共用一个bias，总共有out_channels个biasself.bias=np.random.normal(loc=0.0,scale=1.0,size=(out_channels,))def check_x_mat_shape(self,x_mat):'''要求卷积核在卷积时可以把矩阵铺满，右侧不能有多余的列，下侧不能有多余的行如28×28不能用（5×5，stride=2）的卷积核，因为它只能覆盖（27×27）''' row,col=x_mat.shape[0],x_mat.shape[1]k,s=self.ks,self.strideassert (row-k) // s*s + k ==rowassert (col-k) // s*s + k ==coldef __call__(self,x_mat:np.ndarray)->np.ndarray:self.check_x_mat_shape(x_mat)self.x=x_mat# 做卷积运算return self.__conv(stride=self.stride,mat=x_mat,kernel=self.kernel,bias=None,einsum_formula="ijk,ijkl->l",out_ele_shape=[self.out_channels])def backward(self,dc_dz:np.ndarray)->np.ndarray:# 反向卷积的目标是dc_dz补0后的矩阵（张量）# （padding+dilation）# 补0规则为：边缘padding kernel_size-1 层0；间隔处补stride-1层0# 只看横向，如果dc_dz有c列，那该矩阵有2kernel_size+(m-1)stride-1列# 反向卷积的stride固定为1# print("backward:ConvLayer")(kr,kc,in_ch,out_ch),s=self.kernel.shape,self.stridedc_dz_with_zeros_shape=(2*kr+(dc_dz.shape[0]-1)*s-1,2*kc+(dc_dz.shape[1]-1)*s-1,dc_dz.shape[2]) # 用D表示补0后的张量D=np.zeros(dc_dz_with_zeros_shape)for i in range(dc_dz.shape[0]):for j in range(dc_dz.shape[1]):D[kr+i*s-1,kc+j*s-1]=dc_dz[i,j]# 求dc_da,a指该层的输入self.x# 此时需注意步长是1# kernel[i,j,k,l]在正向推导时i表示row,j表示col,k表示in_ch,l表示out_chdc_da=self.__conv(stride=1,mat=D,kernel=self.kernel[::-1,::-1],bias=None,einsum_formula="ijl,ijkl->k",out_ele_shape=[in_ch]) # 求dc_dwdc_dw=self.__conv(stride=1,mat=self.x,kernel=dc_dz,bias=None,einsum_formula="ijk,ijl->kl",out_ele_shape=[in_ch,out_ch])# 求dc_dbdc_db=np.einsum("ijk->k",dc_dz)# 更新w和b，并返回dc_daself.kernel+=dc_dwself.bias+=dc_dbreturn dc_da# 它的返回值就是卷积层前向传播时的返回值def __conv(self,stride:int,mat:np.ndarray,kernel:np.ndarray,bias:np.ndarray=None,einsum_formula:str="ijk,ijkl->l",out_ele_shape:Iterable[int]=None)->np.ndarray:# 卷积运算 sub_np_tensor*kernel_np_tensor+biasif bias is None:def f(m):return np.einsum(einsum_formula,m,kernel)else:def f(m):return np.einsum(einsum_formula,m,kernel)+bias row,col=mat.shape[0],mat.shape[1]s=stride(kr,kc,*omit),s=kernel.shape,stride# 正向推导if out_ele_shape is None:assert len(kernel.shape)==3out_ch=kernel.shape[-1]out_ele_shape=(out_ch,)target_shape=((row-kr)//s+1,(col-kc)//s+1,*out_ele_shape)target = np.zeros(target_shape)for i in range(target_shape[0]):for j in range(target_shape[1]):r,c=i*s,j*starget[i,j]=f(mat[r:r+kr,c:c+kc])return target

二、激活层代码

# 激活函数的模板类class Func:def __init__(self,f:types.FunctionType,f_derivate:types.FunctionType,jacobin=False):self.f=fself.f_derivate=f_derivateself.jacobin=jacobin # True表明f导数将使用雅克比矩阵进行表示# 激活函数正向传播def __call__(self,z):return self.f(z)# 激活函数反向传播（求导）def derivate(self,z):return self.f_derivate(z) # 激活函数层class ActivateLayer(Layer):def __init__(self,activate_fn:Func):self.f=activate_fnself.z:np.ndarray=Nonedef __call__(self,x:np.ndarray)->np.ndarray:self.z=xreturn self.f(x)def backward(self,dc_da:np.ndarray)->np.ndarray:# print("backward:ActivateLayer")# 作为最后一个激活函数，它接收的dc_da为(-1/n)*(label/pred)*(-lr)da_dz=self.f.derivate(self.z)if self.f.jacobin:# 如果求导结果只能表示成雅克比矩阵，得用矩阵乘法dc_dz=dc_da.dot(da_dz.T)else:# 如果求导结果为对角矩阵，可以采用逐值相乘来简化运算dc_dz=dc_da*da_dzreturn dc_dz# 定义relu激活函数def f_relu(z):return np.maximum(z,0)def f_relu_derivate(z):return np.heaviside(z,0.5)relu=Func(f_relu,f_relu_derivate)# 定义sigmoid激活函数def f_sigmoid(z):return 1.0/(1.0+np.exp(-z))def f_sigmoid_derivate(z):y=f_sigmoid(z)return y*(1-y)sigmoid=Func(f_sigmoid,f_sigmoid_derivate)# 定义softmax激活函数# 之所以要乘一个exp(-max(z))，是为了防止计算结果超出浮点数的范围，比如exp(500)def f_softmax(z):# 得到的z是个1*n的数组，因此要取z=z[0]# z=z[0]exps=np.exp(z-np.max(z))exp_sum=np.sum(exps)b=exps/exp_sumreturn b[0]def f_softmax_derivate(z):y=f_softmax(z).reshape((-1,))# 相当于Iy-y·y^Treturn np.diag(y)-y.reshape((-1,1)).dot(y.reshape(1,-1))softmax=Func(f_softmax,f_softmax_derivate,True)

三、池化层代码

# 平均池化层class MeanPoolingLayer(Layer):def __init__(self,kernel_size:int,stride:int):self.ks=kernel_sizeself.kernel_shape=(kernel_size,kernel_size)self.channels:int=Noneself.stride=strideself.input_shape:tuple=None # row_cnt,col_cnt,channelsself.target_shape:tuple=None # 目标的shape# 第3维为通道数def __call__(self,mat:np.ndarray):self.input_shape=mat.shapeself.channels=mat.shape[2] row,col=mat.shape[0],mat.shape[1](kr,kc),s=self.kernel_shape,self.strideself.target_shape=((row-kr)//s+1,(col-kc)//s+1,self.channels) target=np.zeros(self.target_shape)for i in range(self.target_shape[0]):for j in range(self.target_shape[1]):r,c=i*s,j*starget[i,j]=np.average(mat[r:r+kr,c:c+kc],axis=(0,1))return targetdef backward(self,d_out:np.ndarray)->np.ndarray:# print("backward:MeanPoolingLayer")d_input=np.zeros(self.input_shape)n=self.kernel_shape[0]*self.kernel_shape[1]d_mat=d_out/n(kr,kc),s=self.kernel_shape,self.stridefor i in range(self.target_shape[0]):for j in range(self.target_shape[1]):r,c=i*s,j*sd_input[r:r+kr,c:c+kc]+=d_mat[i,j]# 如果需要传入学习率，那应该用下面这行# d_input[r:r+kr,c:c+kc]+=lr*d_mat[i,j]return d_input

四、flatten层代码

# reshape层，主要用于卷积层和全连接层的衔接# from_shape可以为None，输入第一个矩阵后自动获取shape# to_shape为(1,-1)，因为全连接层的输入就是一个一维张量，用二维表示只是方便用array.dotclass ReshapeLayer(Layer):def __init__(self,from_shape,to_shape):self.from_shape=from_shapeself.to_shape=to_shapedef __call__(self,x:np.ndarray)->np.ndarray:if self.from_shape is None:self.from_shape=x.shapereturn x.reshape(self.to_shape)def backward(self,d:np.ndarray)->np.ndarray:# print("backward:reshapeLayer")return d.reshape(self.from_shape)

五、全连接层代码

# 全连接层class FullConnectedLayer(Layer):def __init__(self,input_size,output_size):self.i_size=input_sizeself.o_size=output_sizeif self.i_size is not None:self.__init(self.i_size)# 初始化权重，并为输入矩阵x留一个坑位def __init(self,input_size):self.i_size=input_sizeself.w=np.random.normal(loc=0.0,scale=1.0,size=(self.i_size,self.o_size))self.b=np.random.normal(loc=0.0,scale=1.0,size=(1,self.o_size))self.x: np.ndarray=None # inputdef __call__(self,x:np.ndarray)->np.ndarray: x=x.reshape(1,-1)if self.i_size is None:# x的第一维是1，第二维就是输入的尺寸self.__init(x.shape[1])self.x=xself.z=x.dot(self.w)+self.breturn self.z# 反向传播def backward(self,dc_dz:np.ndarray)->np.ndarray:'''输入是单个样本的误差，输出也是单个样本的误差，因此这一层反向传播时，需要把各个样本的误差重新组合起来，形成一个完整的误差矩阵，才可以继续反向传播到其它层'''# print("backward:FullConnectedLayer")dc_dx=dc_dz.dot(self.w.T) self.w+=self.x.T.dot(dc_dz.reshape((1,-1)))self.b+=dc_dzreturn dc_dx

六、神经网络代码

# 神经网络class NN:def __init__(self,input_shape=(1,-1),output_shape=(1,-1)):self.layers=list()self.input_shape=input_shapeself.output_shape=output_shapedef predict(self,test_X:np.ndarray):'''输入：test_X.shape=n_samples*height*width*channels输出：pred.shape=n_samples*class_num'''n_samples=test_X.shape[0]pred=[]for i in range(n_samples):pred.append(self.forward(test_X[i]))return pred def forward(self,a:np.ndarray)->np.ndarray:for layer in self.layers:a=layer(a) return a # 神经网络的反向传播def backward(self,dc_da_last:np.ndarray)->np.ndarray:# print("begin NN backward")# dc_da_last=loss*-lr，即dc_da_last=-lr*1/n*(label/pred),形状为n_samples*n_classd=dc_da_last # 这样只是方便表示而已for layer in self.layers[::-1]:d=layer.backward(d) return ddef train(self,input_vec,label,loss_func:LossFunc,lr):n_samples=input_vec.shape[0]for i in tqdm(range(n_samples)):y=self.forward(input_vec[i])label_dilation=np.zeros((y.shape))label_dilation[label[i]]=1loss=loss_func.derivate(label_dilation,y)self.backward(loss*-lr)def set_layers(self,layers):self.layers=layersdef append(self,layer):self.layers.append(layer)

七、训练与测试

# 定义训练方法，输入为整个数据集，输出为训练后的模型def train_model(train_X,train_Y,epochs,class_num,lr):my_nn=NN(train_X[0].shape,(1,1))my_nn.set_layers([# 神经网络的结构# conv1ConvolutionLayer(1,4,3,1),# relu1ActivateLayer(relu),# pool1MeanPoolingLayer(2,2),# relu1ActivateLayer(relu),# flatReshapeLayer(None,(1,-1)),# fullyconnected(fc)FullConnectedLayer(None,class_num),# softmax(SM)ActivateLayer(softmax),])for i in tqdm(range(epochs)):my_nn.train(train_X,train_Y,cee,lr)print(f"convolutionlayer.kernel={my_nn.layers[0].kernel}")return my_nn# 定义评估方法def evaluate(pred,label):n_samples=len(label)t=0for i in range(n_samples):if np.argmax(pred[i])==label[i]:t+=1return t*1.0/n_samples # 训练、预测与评估：model_2=train_model(train_images,train_labels,epochs=10,class_num=10,lr=0.001)pred=model_2.predict(test_images)evaluate(pred,test_labels)

对于训练集中的60000个样本，设定10个epoch，总训练时长为67mins，准确率为87.23%；设定70个epoch时，准确率为96.3%

八、模型保存与加载

import picklesave_model=open("model_2.pkl","wb")str_model_2=pickle.dumps(model_2)save_model.write(str_model_2)save_model.close()with open("model_2.pkl","rb") as file:m2=pickle.loads(file.read())

总结

1.这次实验从学习理论到完全弄懂CNN的过程以及跑通，一共花了整整一周

N的整个过程为，在每个epoch中，每个样本都会产生一个误差，这个误差反向传播，就得到了参数的一次更新梯度。比如一个数据集有n_samples=1000，设置了epochs=10，那么整个训练过程中，每层的参数一共会经历epochs×n_samples次更新。

如果觉得《基于numpy的CNN实现 进行MNIST手写数字识别》对你有帮助，请点赞、收藏，并留下你的观点哦！