用python的numpy实现mnist手写数字识别

完整代码的文章底部（Optimization_mnist.py和lr_utils.py），原理和公式部分可以看前面文章，转载文章请附上本文链接

学完前面（1到6）文章就完成了吴恩达deeplearning ai 课程前面2门课程的内容了，可以写出下面的代码，可以去参加一些比赛，这里推荐一个kaggle上面的一个mnist手写数字识别的知识竞赛，在没有使用深度学习框架情况下他的评分达到了0.94914，代码我直接给出来，就不解释了。

kaggle比赛链接
这下面是我的评分，代码再文章底部，我就不解释了,训练集需要去kaggle平台下载，这个比赛是一直开放的（lr_mnist.py和test_mnist.py）

内容分为三部分：

1.数据准备和处理
2.训练过程
3.测试过程

一.数据准备和处理：

我数据集是用tensorflow2.0获取mnist上的7万张手写数字图片（也可以用其他方式获取，不用tensorflow，我是为了方便就用tensorflow获取了），需要先安装tensorflow2.0版本。安装教程可以看前面文章：在win10上安装tensorflow，并且配置编译器

import  tensorflow as tf
from    tensorflow.keras import datasets
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
def load_dataset():(x, y), (x_test, y_test) = datasets.mnist.load_data()#使用tensorflow获取mnist手写数字图片训练集和测试集#对图片进行降维，训练集合x有6万张图片，测试集x_test有一万张图片，都是28*28，通道数为1（也就是黑白色），这里就是一个三维元组了#标签y是0到9的数字，表示图片是数字几，可以输出看一下。由于x是元组类型，无法用shape（），只能用shape【】y = y.reshape(1,-1)y_test = y_test.reshape(1, -1)print("训练集维度：",x.shape[0],x.shape[1],x.shape[2])#打印训练集维度print("测试集维度：",x_test.shape[0],x_test.shape[1],x_test.shape[2])#打印测试集维度#用独热码处理标签,调用pandas库中Series和get_dummier函数y = np.squeeze(y)#使得标签是1维的y = pd.Series(y)y = pd.get_dummies(y)#print(y)这个可以查看矩阵是如何的，y是60000万行，10列y = np.array(y).T#y变成10行，60000列，60000列是样本个数y_test = np.squeeze(y_test)  # 使得标签是1维的y_test = pd.Series(y_test)y_test = pd.get_dummies(y_test)# print(y_test)这个可以查看矩阵是如何的，y_test是10000万行，10列y_test = np.array(y_test).T  # y_test变成10行，10000列，10000列是样本个数train_set_x_orig, train_set_y, test_set_x_orig, test_set_y = x,y,x_test,y_testreturn train_set_x_orig, train_set_y, test_set_x_orig, test_set_y
load_dataset()
'''
train_set_x_orig, train_set_y, test_set_x_orig, test_set_y=load_dataset()
#可以打印图片看一下，比如训练集中第100张的图片
plt.imshow(train_set_x_orig[99])  # 打印的图片
print(train_set_y[0][99])  # 输出图片对应标签
plt.show()  # 展示需要打印的图片
'''

获取数据集很简单只需要一行代码：(x, y), (x_test, y_test) = datasets.mnist.load_data()
由于标签y和y_test是一维数组，我们要做的是把标签化为softmax分类器的类型，需要用独热码（one-hot）处理。

用pandas内的函数实现（具体看注释）：

y = pd.Series(y)
y = pd.get_dummies(y)

不清楚的地方可以print打印看看，我们打印一下训练集第99张图片：

对于图片降维的处理我们在主函数里面进行。

主函数代码：

if __name__ == "__main__":L = 10#神经网络层数dim = 20#隐藏层节点个数learning_rate = 0.05#学习率loss= []#损失函数lambd = 0.01 # L2正则化参数beta = 0.9#β值；1/（1-β）表示平均前多少轮的指数加权平均数decay_rate = 0.0009#学习率衰减率mini_batch = 300#一次的训练量，60000万张图片，要带入200次，全部训练一遍称为一代sigma2, mu = 0, 0  # 用于决定将batch正则化拟合进神经网络的均值和方差train_set_x_orig, train_set_y, test_set_x_orig, test_set_y = load_dataset()train_set_x = train_set_x_orig.reshape((train_set_x_orig.shape[0], -1)).T / 255  # 降维,化为区间（0，1）内的数test_set_x = test_set_x_orig.reshape((test_set_x_orig.shape[0], -1)).T / 255  # 降维,化为区间（0，1）内的数print("训练集降维后的维度： " + str(train_set_x.shape))print("训练集_标签的维数 : " + str(train_set_y.shape))print("测试集降维后的维度: " + str(test_set_x.shape))print("测试集_标签的维数 : " + str(test_set_y.shape))print()w,b,Vdw,Vdb,Sdw,Sdb = ward(L,train_set_x.shape[1],train_set_x.shape[0],dim)#vdw表示momentum，Sdw表示RMSpropfor i in range(0,200):Sigma2, Mu,J_average = 0,0,0#用于决定将batch正则化拟合进神经网络的均值和方差for j in range(0,(train_set_x.shape[1]//mini_batch)):#w,b,J,Vdw,Vdb = back_momentum(w,b,train_set_x,train_set_y,learning_rate,train_set_x.shape[1],L,lambd,beita,Vdw,Vdb)#w,b,J,Vdw,Vdb = back_RMSprop(w, b, train_set_x,train_set_y,learning_rate,train_set_x.shape[1],L,lambd,beita, Sdw, Sdb)w,b,J,Vdw,Vdb,Sdw,Sdb,sigma2,mu = back_Adam(w,b,((train_set_x.T)[j*mini_batch:(j+1)*mini_batch]).T,((train_set_y.T)[j*mini_batch:(j+1)*mini_batch]).T,learning_rate,mini_batch,L,lambd,beta,Vdw,Vdb, Sdw, Sdb)#如果有多个mini_batch应该在此处对返回的sigam2和mu使用指数加权平均数，也是滑动平均数Sigma2 = np.multiply(beta,Sigma2) + np.multiply((1-beta),sigma2)Mu = np.multiply(beta,Mu) + np.multiply((1-beta),mu)J_average = np.multiply(beta,J) + np.multiply((1-beta),J)learning_rate = learning_rate * (1 / (1 + i*decay_rate) )#使用学习率衰减if i % 10 == 0:print("loss：",J_average)loss.append(J_average)plt.plot(loss)#打印损失函数plt.show()text(test_set_x,test_set_y,w,b,L,test_set_x.shape[1],lambd,Sigma2,Mu)

训练之前我们需要初始化参数w和b，使用ward函数：顺便把Adam梯度下降法的表示式也一起初始化了，w其实不是三维数组，可以理解为是把多个矩阵放进一个列表里面，当然也可以放字典里面。训练集X是784行，60000列的矩阵，每一列表示的是特征值,60000表示样本个数，所以w1维度：（第一层节点个数，一个样本特征个数784）；由于标签是0到9，有10种情况，所以wL维度：（10，L-1层节点个数），节点也就是隐藏单元。

def ward(L,n,m,dim):#对参数进行初始化np.random.seed(1)w = []b = []Vdw, Vdb, Sdw, Sdb = [],[],[],[]#优化算法的值for i in range(0, L):if i != 0 and i != L - 1:#p = np.random.randn(dim, dim) *0.001p = np.random.randn(dim, dim) * np.sqrt(2 / dim)elif i == 0:#p = np.random.randn(dim, m) * 0.001p = np.random.randn(dim, m) * np.sqrt(2 / m)else:#p = np.random.randn(1, dim) * 0.001p = np.random.randn(10, dim) * np.sqrt(2 / dim)w.append(p)b.append(1)Vdw.append(0)#初始化为0Vdb.append(0)Sdw.append(0)Sdb.append(0)return w, b,Vdw,Vdb,Sdw,Sdb

三.训练过程:

1.我们先看一下前向传播过程：使用的是forward_back函数，前向传播中我使用了归一化输入（σ^2 和μ要缓存下来，用于测试时候的归一化输入；γ和β我是设置成刚好可以和z_norm约了的特殊情况，也就是数据方差为1，平均值为0的情况，你也可以设置成其他情况，还可以对γ和β进行梯度下降，不断跟新γ和β），我还使用了L2正则化，我是将神经网络层数设置为10层（你也可以用dropout正则化）。注意此时的成本函数J要重新定义我们前面的文章中得到的最后一层输出aL（也就是y^）是一个一行多列矩阵，这里是10分类情况，得到的aL是10行60000列的矩阵,具体看代码，使用softmax分类器。不理解地方可以调试看看。

def forward_back(w,b,a,Y,L,m,lambd):#用于训练的归一化前向传播z = []add = 0sigma2,mu = [],[]for i in range(0, L):zl = np.dot(w[i], a) + b[i]#归一化输入zmuL = (1/m)*np.sum(zl,axis=1,keepdims=True)sigmaL = (1/m)*np.sum(np.power(zl-muL,2),axis=1,keepdims=True)z_norm = (zl-muL)/(np.sqrt(sigmaL+0.00000001))gamma,beta_1 = 1*np.sqrt(sigmaL+0.00000001),muL+0#此时z的方差为1，均值为0zl = np.multiply(z_norm,gamma) + beta_1mu.append(muL)sigma2.append(sigmaL)#可以发现，如果使用zl = z_norm此时平均值为0，方差为1，对于我们这个数据是没有影响的，因为我们图片像素点已经除以了255，也比较均匀，是基本处于这个范围内的数#因为z = γ*z_norm+β可以看成线性函数，类似于z = wx+b，也可以对γ和β进行梯度下降更新。add += np.sum((lambd / (2 * m)) * np.dot(w[i], w[i].T))#L2正则化项z.append(zl)a = relu(zl)#使用softmax回归t = np.exp(zl)ti = np.sum(t,axis = 0,keepdims=True)#axis=0表示对行求和，对输出zl的行求和，这样可以保证最后概率之和为1a = np.divide(t,ti)#矩阵除法函数，也可以用/。不过我这样用/有时候会提示一些奇怪错误，提示无法进行除法#损失函数也应该重新定义J = (-1/m)*np.sum(Y*np.log(a))+ add #注意Y*nu.log(a)计算得到是多行的矩阵，但只有一行是非0行，我们之前定义的损失函数相当于这个特殊情况# a = sigmoid(zl)#J = (-1/m)*np.sum(1 * Y * np.log(a) + (1 - Y) * np.log(1 - a)) + add  # 损失函数return z, a, J,sigma2,mu

2.前向传播完成了，我们来看看反向传（播使用的是backward函数和back_Adam函数）：反向求导运算，缓存dw和db，然后对w和b进行更新，我使用的是Adam优化器。

不过有3个细节部分我没有完成：

反向传播也可以用递归写，这样好像会方便，我这样写代码可能比较丑，…嗯就是比较丑。
momentum优化器和RMSprop优化器里面一些参数我还没加进去，不能直接调用，你需要稍微改改，我比较懒不想改了。
还有一个就是Adam通常都是要用偏差修正，虽然说不用好像没多大影响，但我用上就数值就变得很大，比较奇怪；所以我没用，注释了，没用的话多轮训练后是没有什么影响的，也没关系。

我们训练过程中缓存了σ^2和μ，这个在使用mini_batch时候，需要进行用指数加权平均数重新计算σ ^2和μ，如果mini_batch为整个训练集时候，就不要用指数加权平均数计算了，直接用这个σ ^2和μ就好了。我们还可以使用学习率衰减率，对学习率进行跟新。

注意:设置mini_batch很重要，我是设置成300，如果你设置成整个训练集的话，也就是60000，不仅是对cpu的考验，训练效果也是非常差。我设置成60000后准确率只有9.8%，和没训练一个样，改成300后准确率就达到了93%，当然你也可以使用其他大小，听说2的指数次方运算速度比较好。

def backward(w,b,X,Y,m,L,lambd):#反向传播z,a,J,sigma2,mu = forward_back(w,b,X,Y,L,m,lambd)dw,db = [],[]for i in range(L - 1, 0, -1):if i == L - 1:dz = a - Yelse:dz = np.dot(w[i + 1].T, dz) * relu_1(z[i])Dw = 1 / m * (np.dot(dz, relu(z[i - 1]).T)) + (lambd / m) * w[i]Db = 1 / m * np.sum(dz, axis=1, keepdims=True)dw.append(Dw)db.append(Db)dz = np.dot(w[1].T, dz) * relu_1(z[0])Dw = 1 / m * np.dot(dz, X.T) + (lambd / m) * w[0]Db = 1 / m * np.sum(dz, axis=1, keepdims=True)dw.append(Dw)db.append(Db)return dw, db, J,sigma2,mu

def back_Adam(w,b,X,Y,learning,m,L,lambd,beta,Vdw,Vdb,Sdw,Sdb):#使用Adam梯度下降的反向传播dw, db, J ,sigma2,mu= backward(w,b,X,Y,m,L,lambd)# 通常使用偏差修正for i in range(0,L):#注意dw和db是由后往前存的，位置与w和b相反Vdw[i] = beta * Vdw[i] + (1 - beta) * dw[L - i - 1]Vdb[i] = beta * Vdb[i] + (1 - beta) * db[L - i - 1]Sdw[i] = beta * Sdw[i] + (1 - beta) * np.power(dw[L-i-1],2)Sdb[i] = beta * Sdb[i] + (1 - beta) * np.power(db[L-i-1],2)'''Vdw[i] /= 1 - np.power(beta, t)#t表示代数（所有mini_batch训练一次称为一代）Vdb[i] /= 1 - np.power(beta, t)Sdw[i] /= 1 - np.power(beta, t)Sdb[i] /= 1 - np.power(beta, t)'''w[i] = w[i] - learning * Vdw[i] / (np.power(Sdw[i],1/2) + 0.00000001)#加上一个较小的值，防止整个值变成无穷大b[i] = b[i] - learning * Vdb[i] / (np.power(Sdb[i],1/2) + 0.00000001)#加上一个较小的值，防止整个值变成无穷大return w,b,J,Vdw,Vdb,Sdw,Sdb,sigma2,mu

三.测试过程：使用text函数，调用的前向传播函数是forward_test函数

def forward_test(w, b, a, Y, L, m, lambd,sigma2,mu):  # 用于测试的归一化前向传播for i in range(0, L):zl = np.dot(w[i], a) + b[i]# 归一化输入zz_norm = (zl - mu[i]) / (np.sqrt(sigma2[i] + 0.00000001))gamma, beta_1 = 1 * np.sqrt(sigma2[i] + 0.00000001), mu[i] + 0  # 此时z的方差为1，均值为0zl = np.multiply(z_norm, gamma) + beta_1# 可以发现，如果使用zl = z_norm此时平均值为0，方差为1，对于我们这个数据是没有影响的，因为我们图片像素点已经除以了255，也比较均匀，是基本处于这个范围内的数a = relu(zl)# 使用softmax回归t = np.exp(zl)ti = np.sum(t, axis=0, keepdims=True)  # axis=0表示对行求和，对输出zl的行求和，这样可以保证最后概率之和为1a = np.divide(t, ti)  # 矩阵除法函数，也可以用/。不过我这样用/有时候会提示一些奇怪错误，提示无法进行除法# a = sigmoid(zl)return a

这个前向传播中归一化使用的σ^2和μ是训练过程中缓存的，其他基本和训练的前向传播一个样。

def text(x,y,w,b,L,m,lambd,sigma2, mu):#查看测试集准确率a = forward_test(w, b, x, y, L, m,lambd,sigma2,mu)#得到的a是一个10行，10000列的矩阵，里面数表示概率a = a.T#先将a转置for i in range(x.shape[1]):a[i] = np.where(a[i].max() == a[i], 1, 0)#最大数位置改为存1，不是最大数改为存0，这样就处理好y^了a = a.T  # 将a转回去lop_1 = np.sum(np.abs(y-a),axis = 0,keepdims=True)/2#这个lop_1中存的是1和0，1表示预测错误，0表示正确lop_2 = 1-np.sum(lop_1,axis=1,keepdims=True)/10000#表示正确所占比例，也就是准确率lop_2=np.squeeze(lop_2)print("测试集准确性：{0}%".format(lop_2*100))return 0

这个标签处理我用的方式你可能觉得有点怪，嗯…希望你可以理解它。

我们用训练好的w和b来计算测试集的图片输入数据，得到一个y^，y ^里面是概率，之和为1。
y^是的维度是：（10，10000），我们将它转置后用np.where()把每一行的数最大的改为1，其他改为0，然后再转置回来。这样aL（也就是y ^）和Y，我们选其第一列来看（也就是第一个样本），aL-Y取绝对值后，由于aL和y ^每一列中只有一个1，其他全部为0，aL如果和Y相等那么相减后这一列应该全为0，否则这一列里面有2个1，经过lop_1 = np.sum(np.abs(y-a),axis = 0,keepdims=True)/2和lop_2 = 1-np.sum(lop_1,axis=1,keepdims=True)/10000运算后lop_2就是识别正确的图片所占的比例，将这个作为准确率，当然你也可以自己拍一张照片，替代测试集，将你照片（必须为黑白照，因为训练集就是黑白照，也就是单通道的）数据输入使用训练好的w和b计算，他输出的y^就是你识别的结果（y ^你可能也需要简单的处理一下，我就不写这部分代码了，大家加油）。

下面这张图片是我的训练结果：训练200轮后，准确率93%，你可以多调试一下，参数多调调，应该会有百分之95%以上准确率

完整代码:

lr_utils.py:

import  tensorflow as tf
from    tensorflow.keras import datasets
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
def load_dataset():(x, y), (x_test, y_test) = datasets.mnist.load_data()#使用tensorflow获取mnist手写数字图片训练集和测试集#对图片进行降维，训练集合x有6万张图片，测试集x_test有一万张图片，都是28*28，通道数为1（也就是黑白色），这里就是一个三维元组了#标签y是0到9的数字，表示图片是数字几，可以输出看一下。由于x是元组类型，无法用shape（），只能用shape【】y = y.reshape(1,-1)y_test = y_test.reshape(1, -1)print("训练集维度：",x.shape[0],x.shape[1],x.shape[2])#打印训练集维度print("测试集维度：",x_test.shape[0],x_test.shape[1],x_test.shape[2])#打印测试集维度#用独热码处理标签,调用pandas库中Series和get_dummier函数y = np.squeeze(y)#使得标签是1维的y = pd.Series(y)y = pd.get_dummies(y)#print(y)这个可以查看矩阵是如何的，y是60000万行，10列y = np.array(y).T#y变成10行，60000列，60000列是样本个数y_test = np.squeeze(y_test)  # 使得标签是1维的y_test = pd.Series(y_test)y_test = pd.get_dummies(y_test)# print(y_test)这个可以查看矩阵是如何的，y_test是10000万行，10列y_test = np.array(y_test).T  # y_test变成10行，10000列，10000列是样本个数train_set_x_orig, train_set_y, test_set_x_orig, test_set_y = x,y,x_test,y_testreturn train_set_x_orig, train_set_y, test_set_x_orig, test_set_y
load_dataset()train_set_x_orig, train_set_y, test_set_x_orig, test_set_y=load_dataset()
#可以打印图片看一下，比如训练集中第100张的图片
plt.imshow(train_set_x_orig[99])  # 打印的图片
print(train_set_y[0][99])  # 输出图片对应标签
plt.show()  # 展示需要打印的图片

Optimization_mnist.py

# conding:utf-8
import numpy as np
import matplotlib.pyplot as plt
from lr_utils import load_datasetdef tanh(z):#tanh函数return (np.exp(z)-np.exp(-z))/(np.exp(z)+np.exp(-z))
def relu(z):#relu函数return np.maximum(0,z)
def tanh_1(z):#tanh函数的导数return 1-tanh(z)**2
def relu_1(z):#relu函数的导数return np.maximum(0, z/np.abs(z))
def sigmoid(z):return 1/(1+np.exp(-z))def ward(L,n,m,dim):#对参数进行初始化np.random.seed(1)w = []b = []Vdw, Vdb, Sdw, Sdb = [],[],[],[]#优化算法的值for i in range(0, L):if i != 0 and i != L - 1:#p = np.random.randn(dim, dim) *0.001p = np.random.randn(dim, dim) * np.sqrt(2 / dim)elif i == 0:#p = np.random.randn(dim, m) * 0.001p = np.random.randn(dim, m) * np.sqrt(2 / m)else:#p = np.random.randn(1, dim) * 0.001p = np.random.randn(10, dim) * np.sqrt(2 / dim)w.append(p)b.append(1)Vdw.append(0)#初始化为0Vdb.append(0)Sdw.append(0)Sdb.append(0)return w, b,Vdw,Vdb,Sdw,Sdbdef forward_back(w,b,a,Y,L,m,lambd):#用于训练的归一化前向传播z = []add = 0sigma2,mu = [],[]for i in range(0, L):zl = np.dot(w[i], a) + b[i]#归一化输入zmuL = (1/m)*np.sum(zl,axis=1,keepdims=True)sigmaL = (1/m)*np.sum(np.power(zl-muL,2),axis=1,keepdims=True)z_norm = (zl-muL)/(np.sqrt(sigmaL+0.00000001))gamma,beta_1 = 1*np.sqrt(sigmaL+0.00000001),muL+0#此时z的方差为1，均值为0zl = np.multiply(z_norm,gamma) + beta_1mu.append(muL)sigma2.append(sigmaL)#可以发现，如果使用zl = z_norm此时平均值为0，方差为1，对于我们这个数据是没有影响的，因为我们图片像素点已经除以了255，也比较均匀，是基本处于这个范围内的数#因为z = γ*z_norm+β可以看成线性函数，类似于z = wx+b，也可以对γ和β进行梯度下降更新。add += np.sum((lambd / (2 * m)) * np.dot(w[i], w[i].T))#L2正则化项z.append(zl)a = relu(zl)#使用softmax回归t = np.exp(zl)ti = np.sum(t,axis = 0,keepdims=True)#axis=0表示对行求和，对输出zl的行求和，这样可以保证最后概率之和为1a = np.divide(t,ti)#矩阵除法函数，也可以用/。不过我这样用/有时候会提示一些奇怪错误，提示无法进行除法#损失函数也应该重新定义J = (-1/m)*np.sum(Y*np.log(a))+ add #注意Y*nu.log(a)计算得到是多行的矩阵，但只有一行是非0行，我们之前定义的损失函数相当于这个特殊情况# a = sigmoid(zl)#J = (-1/m)*np.sum(1 * Y * np.log(a) + (1 - Y) * np.log(1 - a)) + add  # 损失函数return z, a, J,sigma2,mudef forward_test(w, b, a, Y, L, m, lambd,sigma2,mu):  # 用于测试的归一化前向传播for i in range(0, L):zl = np.dot(w[i], a) + b[i]# 归一化输入zz_norm = (zl - mu[i]) / (np.sqrt(sigma2[i] + 0.00000001))gamma, beta_1 = 1 * np.sqrt(sigma2[i] + 0.00000001), mu[i] + 0  # 此时z的方差为1，均值为0zl = np.multiply(z_norm, gamma) + beta_1# 可以发现，如果使用zl = z_norm此时平均值为0，方差为1，对于我们这个数据是没有影响的，因为我们图片像素点已经除以了255，也比较均匀，是基本处于这个范围内的数a = relu(zl)# 使用softmax回归t = np.exp(zl)ti = np.sum(t, axis=0, keepdims=True)  # axis=0表示对行求和，对输出zl的行求和，这样可以保证最后概率之和为1a = np.divide(t, ti)  # 矩阵除法函数，也可以用/。不过我这样用/有时候会提示一些奇怪错误，提示无法进行除法# a = sigmoid(zl)return adef backward(w,b,X,Y,m,L,lambd):#反向传播z,a,J,sigma2,mu = forward_back(w,b,X,Y,L,m,lambd)dw,db = [],[]for i in range(L - 1, 0, -1):if i == L - 1:dz = a - Yelse:dz = np.dot(w[i + 1].T, dz) * relu_1(z[i])Dw = 1 / m * (np.dot(dz, relu(z[i - 1]).T)) + (lambd / m) * w[i]Db = 1 / m * np.sum(dz, axis=1, keepdims=True)dw.append(Dw)db.append(Db)dz = np.dot(w[1].T, dz) * relu_1(z[0])Dw = 1 / m * np.dot(dz, X.T) + (lambd / m) * w[0]Db = 1 / m * np.sum(dz, axis=1, keepdims=True)dw.append(Dw)db.append(Db)return dw, db, J,sigma2,mu
def back_momentum(w,b,X,Y,learning,m,L,lambd,beta,Vdw,Vdb):#使用momentum梯度下降的反向传播dw, db, J,sigma2,mu = backward(w,b,X,Y,m,L,lambd)# 不使用偏差修正不影响最后结果for i in range(0,L):#注意dw和db是由后往前存的，位置与w和b相反Vdw[i] = beta * Vdw[i] + (1 - beta) * dw[L-i-1]Vdb[i] = beta * Vdb[i] + (1 - beta) * db[L-i-1]w[i] = w[i] - learning * Vdw[i]b[i] = b[i] - learning * Vdb[i]return w,b,J,Vdw,Vdb,sigma2,mu
def back_RMSprop(w,b,X,Y,learning,m,L,lambd,beta,Sdw,Sdb):#使用RMSprop梯度下降的反向传播dw, db, J ,sigma2,mu= backward(w,b,X,Y,m,L,lambd)# 不使用偏差修正不影响最后结果for i in range(0,L):#注意dw和db是由后往前存的，位置与w和b相反Sdw[i] = beta * Sdw[i] + (1 - beta) * np.power(dw[L-i-1],2)Sdb[i] = beta * Sdb[i] + (1 - beta) * np.power(db[L-i-1],2)w[i] = w[i] - learning * dw[L-i-1] / (np.power(Sdw[i],1/2) + 0.00000001)#加上一个较小的值，防止整个值变成无穷大b[i] = b[i] - learning * db[L-i-1] / (np.power(Sdb[i],1/2) + 0.00000001)#加上一个较小的值，防止整个值变成无穷大return w,b,J,Sdw,Sdb,sigma2,mu
def back_Adam(w,b,X,Y,learning,m,L,lambd,beta,Vdw,Vdb,Sdw,Sdb):#使用Adam梯度下降的反向传播dw, db, J ,sigma2,mu= backward(w,b,X,Y,m,L,lambd)# 通常使用偏差修正for i in range(0,L):#注意dw和db是由后往前存的，位置与w和b相反Vdw[i] = beta * Vdw[i] + (1 - beta) * dw[L - i - 1]Vdb[i] = beta * Vdb[i] + (1 - beta) * db[L - i - 1]Sdw[i] = beta * Sdw[i] + (1 - beta) * np.power(dw[L-i-1],2)Sdb[i] = beta * Sdb[i] + (1 - beta) * np.power(db[L-i-1],2)'''Vdw[i] /= 1 - np.power(beta, t)#t表示代数（所有mini_batch训练一次称为一代）Vdb[i] /= 1 - np.power(beta, t)Sdw[i] /= 1 - np.power(beta, t)Sdb[i] /= 1 - np.power(beta, t)'''w[i] = w[i] - learning * Vdw[i] / (np.power(Sdw[i],1/2) + 0.00000001)#加上一个较小的值，防止整个值变成无穷大b[i] = b[i] - learning * Vdb[i] / (np.power(Sdb[i],1/2) + 0.00000001)#加上一个较小的值，防止整个值变成无穷大return w,b,J,Vdw,Vdb,Sdw,Sdb,sigma2,mudef text(x,y,w,b,L,m,lambd,sigma2, mu):#查看测试集准确率a = forward_test(w, b, x, y, L, m,lambd,sigma2,mu)#得到的a是一个10行，10000列的矩阵，里面数表示概率a = a.T#先将a转置for i in range(x.shape[1]):a[i] = np.where(a[i].max() == a[i], 1, 0)#最大数位置改为存1，不是最大数改为存0，这样就处理好y^了a = a.T  # 将a转回去lop_1 = np.sum(np.abs(y-a),axis = 0,keepdims=True)/2#这个lop_1中存的是1和0，1表示预测错误，0表示正确lop_2 = 1-np.sum(lop_1,axis=1,keepdims=True)/10000#表示正确所占比例，也就是准确率lop_2=np.squeeze(lop_2)print("测试集准确性：{0}%".format(lop_2*100))return 0if __name__ == "__main__":L = 10#神经网络层数dim = 20#隐藏层节点个数learning_rate = 0.05#学习率loss= []#损失函数lambd = 0.01 # L2正则化参数beta = 0.9#β值；1/（1-β）表示平均前多少轮的指数加权平均数decay_rate = 0.0009#学习率衰减率mini_batch = 300#一次的训练量，60000万张图片，要带入200次，全部训练一遍称为一代sigma2, mu = 0, 0  # 用于决定将batch正则化拟合进神经网络的均值和方差train_set_x_orig, train_set_y, test_set_x_orig, test_set_y = load_dataset()train_set_x = train_set_x_orig.reshape((train_set_x_orig.shape[0], -1)).T / 255  # 降维,化为区间（0，1）内的数test_set_x = test_set_x_orig.reshape((test_set_x_orig.shape[0], -1)).T / 255  # 降维,化为区间（0，1）内的数print("训练集降维后的维度： " + str(train_set_x.shape))print("训练集_标签的维数 : " + str(train_set_y.shape))print("测试集降维后的维度: " + str(test_set_x.shape))print("测试集_标签的维数 : " + str(test_set_y.shape))print()w,b,Vdw,Vdb,Sdw,Sdb = ward(L,train_set_x.shape[1],train_set_x.shape[0],dim)#vdw表示momentum，Sdw表示RMSpropfor i in range(0,200):Sigma2, Mu,J_average = 0,0,0#用于决定将batch正则化拟合进神经网络的均值和方差for j in range(0,(train_set_x.shape[1]//mini_batch)):#w,b,J,Vdw,Vdb = back_momentum(w,b,train_set_x,train_set_y,learning_rate,train_set_x.shape[1],L,lambd,beita,Vdw,Vdb)#w,b,J,Vdw,Vdb = back_RMSprop(w, b, train_set_x,train_set_y,learning_rate,train_set_x.shape[1],L,lambd,beita, Sdw, Sdb)w,b,J,Vdw,Vdb,Sdw,Sdb,sigma2,mu = back_Adam(w,b,((train_set_x.T)[j*mini_batch:(j+1)*mini_batch]).T,((train_set_y.T)[j*mini_batch:(j+1)*mini_batch]).T,learning_rate,mini_batch,L,lambd,beta,Vdw,Vdb, Sdw, Sdb)#如果有多个mini_batch应该在此处对返回的sigam2和mu使用指数加权平均数，也是滑动平均数Sigma2 = np.multiply(beta,Sigma2) + np.multiply((1-beta),sigma2)Mu = np.multiply(beta,Mu) + np.multiply((1-beta),mu)J_average = np.multiply(beta,J) + np.multiply((1-beta),J)learning_rate = learning_rate * (1 / (1 + i*decay_rate) )#使用学习率衰减if i % 10 == 0:print("loss：",J_average)loss.append(J_average)plt.plot(loss)#打印损失函数plt.show()text(test_set_x,test_set_y,w,b,L,test_set_x.shape[1],lambd,Sigma2,Mu)

lr_mnist.py:

import pandas as pd
import numpy as npdef load_dataset1():df = pd.read_csv('train.csv')  # 得到的是一个字典集f1 = [f"pixel{i}" for i in range(0, 28 * 28)]  # 产生字符串列表，从pixel0到pixel783f2 = 'label'train_x = np.array(df[f1].values)  # 通过键获取字典数据，并且转化为矩阵train_y = np.array(df[f2].values)train_y=pd.Series(train_y)train_y=np.array(pd.get_dummies(train_y))#独热码实现softmax#print(train_y[0:12])train_set_y, test_set_y=train_y[0:40000].T,train_y[40000:42000].T#print(train_x.shape[0], train_x.shape[1])  # 输出维度print(train_y.shape[0],train_y.shape[1])  # 输出维度dp = pd.read_csv('test.csv')  # 得到的是一个字典集f = [f"pixel{i}" for i in range(0, 28 * 28)]  # 产生字符串列表，从pixel0到pixel783test_x = np.array(dp[f].values)  # 通过键获取字典数据，并且转化为矩阵#print(test_x.shape[0], test_x.shape[1])  # 输出维度train_set_x_orig,  test_set_x_orig = train_x[0:40000],train_x[40000:42000]return train_set_x_orig, train_set_y, test_set_x_orig, test_set_ydef load_dataset2():df = pd.read_csv('train.csv')  # 得到的是一个字典集f1 = [f"pixel{i}" for i in range(0, 28 * 28)]  # 产生字符串列表，从pixel0到pixel783f2 = 'label'train_x = np.array(df[f1].values)  # 通过键获取字典数据，并且转化为矩阵train_y = np.array(df[f2].values)train_y=pd.Series(train_y)train_y=np.array(pd.get_dummies(train_y))#独热码实现softmax#print(train_y[0:12])train_set_y=train_y.T#print(train_x.shape[0], train_x.shape[1])  # 输出维度print(train_y.shape[0],train_y.shape[1])  # 输出维度dp = pd.read_csv('test.csv')  # 得到的是一个字典集f = [f"pixel{i}" for i in range(0, 28 * 28)]  # 产生字符串列表，从pixel0到pixel783test_x = np.array(dp[f].values)  # 通过键获取字典数据，并且转化为矩阵#print(test_x.shape[0], test_x.shape[1])  # 输出维度train_set_x_orig,  test_set_x_orig = train_x,test_xreturn train_set_x_orig, train_set_y, test_set_x_orig#train_set_x_orig, train_set_y, test_set_x_orig, test_set_y=load_dataset()
#print(train_set_y.shape[0],train_set_y.shape[1])
#print(train_set_x_orig.shape[0],train_set_x_orig.shape[1])

test_mnist.py:

# conding:utf-8
import numpy as np
import csv
import matplotlib.pyplot as plt
from lr_mnist import load_dataset2def tanh(z):#tanh函数return (np.exp(z)-np.exp(-z))/(np.exp(z)+np.exp(-z))
def relu(z):#relu函数return np.maximum(0,z)
def tanh_1(z):#tanh函数的导数return 1-tanh(z)**2
def relu_1(z):#relu函数的导数return np.maximum(0, z/np.abs(z))
def sigmoid(z):return 1/(1+np.exp(-z))def ward(L,n,m,dim):#对参数进行初始化np.random.seed(1)w = []b = []Vdw, Vdb, Sdw, Sdb = [],[],[],[]#优化算法的值for i in range(0, L):if i != 0 and i != L - 1:#p = np.random.randn(dim, dim) *0.001p = np.random.randn(dim, dim) * np.sqrt(2 / dim)elif i == 0:#p = np.random.randn(dim, m) * 0.001p = np.random.randn(dim, m) * np.sqrt(2 / m)else:#p = np.random.randn(1, dim) * 0.001p = np.random.randn(10, dim) * np.sqrt(2 / dim)w.append(p)b.append(1)Vdw.append(0)#初始化为0Vdb.append(0)Sdw.append(0)Sdb.append(0)return w, b,Vdw,Vdb,Sdw,Sdbdef forward_back(w,b,a,Y,L,m,lambd):#用于训练的归一化前向传播z = []add = 0sigma2,mu = [],[]for i in range(0, L):zl = np.dot(w[i], a) + b[i]#归一化输入zmuL = (1/m)*np.sum(zl,axis=1,keepdims=True)sigmaL = (1/m)*np.sum(np.power(zl-muL,2),axis=1,keepdims=True)z_norm = (zl-muL)/(np.sqrt(sigmaL+0.00000001))gamma,beta_1 = 1*np.sqrt(sigmaL+0.00000001),muL+0#此时z的方差为1，均值为0zl = np.multiply(z_norm,gamma) + beta_1mu.append(muL)sigma2.append(sigmaL)#可以发现，如果使用zl = z_norm此时平均值为0，方差为1，对于我们这个数据是没有影响的，因为我们图片像素点已经除以了255，也比较均匀，是基本处于这个范围内的数#因为z = γ*z_norm+β可以看成线性函数，类似于z = wx+b，也可以对γ和β进行梯度下降更新。add += np.sum((lambd / (2 * m)) * np.dot(w[i], w[i].T))#L2正则化项z.append(zl)a = relu(zl)#使用softmax回归t = np.exp(zl)ti = np.sum(t,axis = 0,keepdims=True)#axis=0表示对行求和，对输出zl的行求和，这样可以保证最后概率之和为1a = np.divide(t,ti)#矩阵除法函数，也可以用/。不过我这样用/有时候会提示一些奇怪错误，提示无法进行除法#损失函数也应该重新定义J = (-1/m)*np.sum(Y*np.log(a))+ add #注意Y*nu.log(a)计算得到是多行的矩阵，但只有一行是非0行，我们之前定义的损失函数相当于这个特殊情况# a = sigmoid(zl)#J = (-1/m)*np.sum(1 * Y * np.log(a) + (1 - Y) * np.log(1 - a)) + add  # 损失函数return z, a, J,sigma2,mudef forward_test(w, b, a,  L,sigma2,mu):  # 用于测试的归一化前向传播for i in range(0, L):zl = np.dot(w[i], a) + b[i]# 归一化输入zz_norm = (zl - mu[i]) / (np.sqrt(sigma2[i] + 0.00000001))gamma, beta_1 = 1 * np.sqrt(sigma2[i] + 0.00000001), mu[i] + 0  # 此时z的方差为1，均值为0zl = np.multiply(z_norm, gamma) + beta_1# 可以发现，如果使用zl = z_norm此时平均值为0，方差为1，对于我们这个数据是没有影响的，因为我们图片像素点已经除以了255，也比较均匀，是基本处于这个范围内的数a = relu(zl)# 使用softmax回归t = np.exp(zl)ti = np.sum(t, axis=0, keepdims=True)  # axis=0表示对行求和，对输出zl的行求和，这样可以保证最后概率之和为1a = np.divide(t, ti)  # 矩阵除法函数，也可以用/。不过我这样用/有时候会提示一些奇怪错误，提示无法进行除法# a = sigmoid(zl)return adef backward(w,b,X,Y,m,L,lambd):#反向传播z,a,J,sigma2,mu = forward_back(w,b,X,Y,L,m,lambd)dw,db = [],[]for i in range(L - 1, 0, -1):if i == L - 1:dz = a - Yelse:dz = np.dot(w[i + 1].T, dz) * relu_1(z[i])Dw = 1 / m * (np.dot(dz, relu(z[i - 1]).T)) + (lambd / m) * w[i]Db = 1 / m * np.sum(dz, axis=1, keepdims=True)dw.append(Dw)db.append(Db)dz = np.dot(w[1].T, dz) * relu_1(z[0])Dw = 1 / m * np.dot(dz, X.T) + (lambd / m) * w[0]Db = 1 / m * np.sum(dz, axis=1, keepdims=True)dw.append(Dw)db.append(Db)return dw, db, J,sigma2,mu
def back_momentum(w,b,X,Y,learning,m,L,lambd,beta,Vdw,Vdb):#使用momentum梯度下降的反向传播dw, db, J,sigma2,mu = backward(w,b,X,Y,m,L,lambd)# 不使用偏差修正不影响最后结果for i in range(0,L):#注意dw和db是由后往前存的，位置与w和b相反Vdw[i] = beta * Vdw[i] + (1 - beta) * dw[L-i-1]Vdb[i] = beta * Vdb[i] + (1 - beta) * db[L-i-1]w[i] = w[i] - learning * Vdw[i]b[i] = b[i] - learning * Vdb[i]return w,b,J,Vdw,Vdb,sigma2,mu
def back_RMSprop(w,b,X,Y,learning,m,L,lambd,beta,Sdw,Sdb):#使用RMSprop梯度下降的反向传播dw, db, J ,sigma2,mu= backward(w,b,X,Y,m,L,lambd)# 不使用偏差修正不影响最后结果for i in range(0,L):#注意dw和db是由后往前存的，位置与w和b相反Sdw[i] = beta * Sdw[i] + (1 - beta) * np.power(dw[L-i-1],2)Sdb[i] = beta * Sdb[i] + (1 - beta) * np.power(db[L-i-1],2)w[i] = w[i] - learning * dw[L-i-1] / (np.power(Sdw[i],1/2) + 0.00000001)#加上一个较小的值，防止整个值变成无穷大b[i] = b[i] - learning * db[L-i-1] / (np.power(Sdb[i],1/2) + 0.00000001)#加上一个较小的值，防止整个值变成无穷大return w,b,J,Sdw,Sdb,sigma2,mu
def back_Adam(w,b,X,Y,learning,m,L,lambd,beta,Vdw,Vdb,Sdw,Sdb):#使用Adam梯度下降的反向传播dw, db, J ,sigma2,mu= backward(w,b,X,Y,m,L,lambd)# 通常使用偏差修正for i in range(0,L):#注意dw和db是由后往前存的，位置与w和b相反Vdw[i] = beta * Vdw[i] + (1 - beta) * dw[L - i - 1]Vdb[i] = beta * Vdb[i] + (1 - beta) * db[L - i - 1]Sdw[i] = beta * Sdw[i] + (1 - beta) * np.power(dw[L-i-1],2)Sdb[i] = beta * Sdb[i] + (1 - beta) * np.power(db[L-i-1],2)'''Vdw[i] /= 1 - np.power(beta, t)#t表示代数（所有mini_batch训练一次称为一代）Vdb[i] /= 1 - np.power(beta, t)Sdw[i] /= 1 - np.power(beta, t)Sdb[i] /= 1 - np.power(beta, t)'''w[i] = w[i] - learning * Vdw[i] / (np.power(Sdw[i],1/2) + 0.00000001)#加上一个较小的值，防止整个值变成无穷大b[i] = b[i] - learning * Vdb[i] / (np.power(Sdb[i],1/2) + 0.00000001)#加上一个较小的值，防止整个值变成无穷大return w,b,J,Vdw,Vdb,Sdw,Sdb,sigma2,mudef text(x,w,b,L,m,lambd,sigma2, mu):#查看测试集准确率a = forward_test(w, b, x,  L,sigma2,mu)#得到的a是一个10行，10000列的矩阵，里面数表示概率a = a.T#先将a转置y_n=[]#就是预测结果的标签for i in range(x.shape[1]):y_n.append(np.argmax(a[i]))my = open("my.csv", "w")my.write("ImageId,Label\n")for j in range(0, x.shape[1]):my.write(str(j+1) + "," + str(y_n[j]) + "\n")return 0if __name__ == "__main__":L = 5#神经网络层数dim = 30#隐藏层节点个数learning_rate = 0.05#学习率loss= []#损失函数lambd = 0.1 # L2正则化参数beta = 0.9#β值；1/（1-β）表示平均前多少轮的指数加权平均数decay_rate = 0.0009#学习率衰减率mini_batch = 200#一次的训练量，40000万张图片，要带入200次，全部训练一遍称为一代sigma2, mu = 0, 0  # 用于决定将batch正则化拟合进神经网络的均值和方差train_set_x_orig, train_set_y, test_set_x_orig = load_dataset2()train_set_x = train_set_x_orig.reshape((train_set_x_orig.shape[0], -1)).T / 255  # 降维,化为区间（0，1）内的数test_set_x = test_set_x_orig.reshape((test_set_x_orig.shape[0], -1)).T / 255  # 降维,化为区间（0，1）内的数print("训练集降维后的维度： " + str(train_set_x.shape))print("训练集_标签的维数 : " + str(train_set_y.shape))print("测试集降维后的维度: " + str(test_set_x.shape))print()w,b,Vdw,Vdb,Sdw,Sdb = ward(L,train_set_x.shape[1],train_set_x.shape[0],dim)#vdw表示momentum，Sdw表示RMSpropfor i in range(0,250):Sigma2, Mu,J_average = 0,0,0#用于决定将batch正则化拟合进神经网络的均值和方差for j in range(0,(train_set_x.shape[1]//mini_batch)):#w,b,J,Vdw,Vdb = back_momentum(w,b,train_set_x,train_set_y,learning_rate,train_set_x.shape[1],L,lambd,beita,Vdw,Vdb)#w,b,J,Vdw,Vdb = back_RMSprop(w, b, train_set_x,train_set_y,learning_rate,train_set_x.shape[1],L,lambd,beita, Sdw, Sdb)w,b,J,Vdw,Vdb,Sdw,Sdb,sigma2,mu = back_Adam(w,b,((train_set_x.T)[j*mini_batch:(j+1)*mini_batch]).T,((train_set_y.T)[j*mini_batch:(j+1)*mini_batch]).T,learning_rate,mini_batch,L,lambd,beta,Vdw,Vdb, Sdw, Sdb)#如果有多个mini_batch应该在此处对返回的sigam2和mu使用指数加权平均数，也是滑动平均数Sigma2 = np.multiply(beta,Sigma2) + np.multiply((1-beta),sigma2)Mu = np.multiply(beta,Mu) + np.multiply((1-beta),mu)J_average = np.multiply(beta,J) + np.multiply((1-beta),J)learning_rate = learning_rate * (1 / (1 + i*decay_rate) )#使用学习率衰减if i % 10 == 0:print("loss：",J_average)loss.append(J_average)plt.plot(loss)#打印损失函数plt.show()text(test_set_x,w,b,L,test_set_x.shape[1],lambd,Sigma2,Mu)