这篇记录的内容来自于Andrew Ng教授在coursera网站上的授课。
1.为什么使用神经网络：样本的特征数太多了，使用前面的算法都几乎不可能达到一个理想的结果。
2.神经网络(neural networks):
受到脑神经的启发，我们有以下结构：
术语：
输入层（input layer）：输入样本的特征。
输出层（output layer）：输出预测。
隐藏层（hidden layer）：因为“看不到”，所以叫隐藏层。（除非debug）
激活单元（activation units）：即图中的$a_1^{(2)},a_2^{(2)}$，通过前面一层的输出作为输入，经过激活函数来计算输出值。表示方法为$a_i^{j}$，含义为第j层（输入层是第一层）的第i个激活单元。
激活函数：这里用sigmoid函数。所以每一个输出的值都在(0,1)中。
z代表这层的输入，a代表这层的输出。
权重（weights）：层与层之间的映射。记为$Theta^{(j)}$，含义为第j层到第j+1层的权重。$Theta^{(j)}_{x,y}$表示第j层的第y个(不要搞反了)激活单元与第j+1层的第x个激活单元的权重。权重矩阵的大小是$s_{j+1}*(s_j+1)$的。+1是因为多了一个偏置(bias)项，老传统了。
前向传播（forward propagation）：根据输入，计算输出的方法。对于第j层的输出，有：$$a^{j}=g(Theta^{j-1}*a^{j-1})$$，其中g是sigmoid函数。
我们的假设h函数就等于它，它也是个向量而不止是一个数字了。
1.记号：
L：总层数。
$s_i$：第i层的节点数（不包括偏置项）。
K：分类数。
那么，整个网络的代价即为：
$$J(Theta)=-frac{1}{m}sum_{i=1}^{m}sum_{k=1}^{K}(y_k^{(i)}log((h_{ heta}(x^{(i)}))_k)+(1-y_k^{i})log(1-(h_{ heta}(x^{(i)}))_k))+frac{lambda}{2m}sum_{l=1}^{L-1}sum_{i=1}^{s_l}sum_{j=1}^{s_{l+1}}(Theta_{j,i})^2$$
2.反向传播(backforward propagation)：神经网络的梯度下降。
$delta_{j}^{(l)}$：第l层第j个激活单元的误差（bias从0开始）。对于输出层我们有：
$$delta^{(l)}=frac{1}{m}sum_{i=1}^{m}(h_{ heta}{(x^{(i)})}-y^{(i)})=frac{1}{m}sum_{i=1}^{m}(g(x^{(i)})-y^{(i)})$$
接着，
$$delta^{(l-1)}=(Theta^{(l-1)})^Tdelta^{(l)}.*g'(z^{(l-1)})$$其中$z^{(l)}$指的是第l层在前向传播时获得的输入，.*是点乘。
注意到，
$$g'(z^{(l)})=a^{(l)}.*(1-a^{(l)})$$
并且，
$$frac{partial}{partialTheta^{(l)}_{i,j}}J(Theta)=a^{(l)}_jdelta_i^{(l+1)}$$
具体推导此处省略。你可以观看3B1B的相关视频，来进行感性的理解。
神经网络的实现还是有点麻烦的。
我们先实现只有一层hidden layer的神经网络。
例子
对某一个人的手写数字进行识别。
样本：200张20*20的灰度图片。方便起见，我们事先将它们变成了用逗号隔开的长度为401的txt文件，对于前400个数字，由0/1代表对应位置是黑色还是白色，最后一位是它对应的数字。
样本下载：https://pan.baidu.com/s/1dEL27nMZzMbT3iAPTiJKqA
提取码：nq18
最后实现代码：
1 import numpy as np 2 import pandas as pd 3 import matplotlib.pyplot as plt 4 import matplotlib 5 from scipy.io import loadmat 6 from sklearn.preprocessing import OneHotEncoder 7 from scipy.optimize import minimize 8 from sklearn.metrics import classification_report###评价报告 9 """ 10 输入数据为20*20的灰度图像，标签为1~10 11 只有一层hidden layer 12 x:样本 13 y:标签 14 m:样本数 15 th:参数 16 """ 17 size_input=20*20 18 size_hidden=30 19 size_labels=10 20 alpha=3###就是lambda 21 22 #######################################sigmoid函数 23 def sigmoid(x): 24 return 1/(1+np.exp(-x)) 25 #######################################fp 26 def forprop(x,th1,th2): 27 m=x.shape[0] 28 a1=np.insert(x,0,values=np.ones(m),axis=1)###添加偏置项 29 z2=a1*th1.T 30 a2=np.insert(sigmoid(z2),0,values=np.ones(m),axis=1)###添加偏置项 31 z3=a2*th2.T 32 h=sigmoid(z3) 33 return a1,z2,a2,z3,h###现在的h的列数等于size_labels 34 #######################################cost函数 35 def cost(th1,th2,size_input,size_hidden,size_labels,x,y,alpha): 36 m=x.shape[0] 37 x=np.matrix(x) 38 y=np.matrix(y) 39 a1,z2,a2,z3,h=forprop(x,th1,th2) 40 A=np.multiply((-y),np.log(h))###数乘 41 B=np.multiply((1-y),np.log(1-h)) 42 reg=(np.sum(np.power(th1[:,1:],2))+np.sum(np.power(th2[:,1:],2)))*alpha/(2*m) 43 return np.sum(A-B)/m+reg 44 #######################################sigmoid导数 45 def sigmoidGradient(x): 46 return np.multiply(sigmoid(x),1-sigmoid(x)) 47 #######################################随机初始化 48 def random(): 49 return (np.random.random(size=size_hidden*(size_input+1)+size_labels*(size_hidden+1))-0.5)*0.24 50 #######################################bp 51 52 TOT=0 53 def backprop(th,size_input,size_hidden,size_labels,x,y,alpha): 54 global TOT 55 m=x.shape[0] 56 print(TOT) 57 TOT=TOT+1 58 x=np.matrix(x) 59 y=np.matrix(y) 60 th1=np.matrix(np.reshape(th[:size_hidden*(size_input+1)],(size_hidden,(size_input+1))))###请注意 61 th2=np.matrix(np.reshape(th[size_hidden*(size_input+1):],(size_labels,(size_hidden+1)))) 62 a1,z2,a2,z3,h=forprop(x,th1,th2) 63 A=np.multiply((-y),np.log(h)) 64 B=np.multiply((1-y),np.log(1-h)) 65 reg=(np.sum(np.power(th1[:,1:],2))+np.sum(np.power(th2[:,1:],2)))*alpha/(2*m) 66 J=np.sum(A-B)/m+reg 67 68 del1=np.zeros(th1.shape) 69 del2=np.zeros(th2.shape) 70 71 for t in range(m): 72 a1t=a1[t,:]###(1,401) 73 z2t=z2[t,:]###(1,25) 74 a2t=a2[t,:]###(1,26) 75 ht=h[t,:] ###(1,10) 76 yt=y[t,:] ###(1,10) 77 78 d3t=ht-yt 79 z2t=np.insert(z2t,0,values=np.ones(1))###(1,26)，老传统 80 d2t=np.multiply((th2.T*d3t.T).T,sigmoidGradient(z2t)) 81 82 del1=del1+(d2t[:,1:]).T*a1t 83 del2=del2+d3t.T*a2t 84 del1=del1/m 85 del2=del2/m 86 87 del1[:,1:]=del1[:,1:]+(th1[:,1:]*alpha)/m 88 del2[:,1:]=del2[:,1:]+(th2[:,1:]*alpha)/m 89 90 tmp=np.concatenate((np.ravel(del1),np.ravel(del2))) 91 92 return J,tmp 93 #######################################获得数据 94 maxlen=6 95 def getStr(x): 96 y=maxlen-len(str(x)) 97 return "0"*y+str(x) 98 99 def getdata(): 100 file=open("dig/NUM.txt") 101 m=int(file.read()) 102 x=np.zeros((m,size_input)) 103 y=np.zeros((m,size_labels)) 104 for i in range(m): 105 path="dig/"+getStr(i)+".txt" 106 file=open(path) 107 Q=file.read() 108 A=Q.split(",") 109 for j in range(size_input): 110 x[i,j]=int(A[j]) 111 y[i,int(A[size_input])]=1 112 print(y) 113 return x,y 114 #######################################预测 115 def predict(x,th1,th2): 116 a1,z2,a2,z3,h=forprop(x,th1,th2) 117 print(np.matrix(np.argmax(h,axis=1))) 118 #######################################主函数 119 def main(): 120 x,y=getdata() 121 th=random() 122 fmin=minimize(fun=backprop,x0=th,args=(size_input,size_hidden,size_labels,x,y,alpha),method='TNC',jac=True,options={'maxiter':1000}) 123 G=fmin.x 124 file=open("save++.txt","w") 125 for i in range(G.shape[0]): 126 file.write(str(G[i])+" ") 127 file.close() 128 x = np.matrix(x) 129 th1 = np.matrix(np.reshape(fmin.x[:size_hidden* (size_input + 1)], (size_hidden, (size_input + 1)))) 130 th2 = np.matrix(np.reshape(fmin.x[size_hidden * (size_input + 1):], (size_labels, (size_hidden + 1)))) 131 a1, z2, a2, z3, h = forprop(x, th1, th2 ) 132 y_pred = np.matrix(np.argmax(h, axis=1)) 133 m=x.shape[0] 134 Y=np.zeros(m) 135 for i in range(m): 136 for j in range(size_labels): 137 if(y[i,j]): 138 Y[i]=j 139 print(classification_report(Y,y_pred)) 140 141 if(__name__=="__main__"): 142 main()
但很显然，没有向量化运行速度太慢，没有结构化修改太复杂！因此博主重构了代码。
接下来的代码中，一个样本是按列排序，而不是上面按行排序（因为博主从头到尾推了一遍）。
这份代码是向量化之后的神经网络，只有一层hidden layer，没有正则项：
1 import numpy as np 2 #import pandas as pd 3 #import matplotlib.pyplot as plt 4 #import matplotlib 5 #from scipy.io import loadmat 6 #from sklearn.preprocessing import OneHotEncoder 7 from scipy.optimize import minimize 8 #from sklearn.metrics import classification_report 9 10 size_input=20*20 11 size_output=10 12 size=[size_input,30,size_output] 13 theta=[] 14 maxlen=6 15 16 def init():###随机初始化 17 global theta 18 tot=0 19 for i in range(len(size)-1): 20 tot+=size[i+1]*(size[i]+1) 21 theta=(np.array(np.random.random(size=tot))-0.5)*0.24 22 23 def sigmoid(x): 24 return 1/(1+np.exp(-x)) 25 26 def getWeight(theta,size):###获得矩阵 27 weight=[] 28 tot=0 29 for i in range(len(size)-1): 30 weight.append(np.matrix(np.array(theta[tot:tot+size[i+1]*(size[i]+1)])).reshape((size[i+1],size[i]+1))) 31 tot+=size[i+1]*(size[i]+1) 32 return weight 33 34 def forwardProp(theta,size,x):###前向传播 35 weight=getWeight(theta,size) 36 m=x.shape[1] 37 layer=len(size) 38 for i in range(layer-1): 39 x=np.insert(x,0,values=np.ones(m),axis=0)###0表示行 40 x=sigmoid(weight[i]*x) 41 return x 42 43 def sigmoidGradient(x): 44 return np.multiply(sigmoid(x),1-sigmoid(x)) 45 46 def BP(theta,size,x,y,alpha): 47 weight=getWeight(theta,size) 48 m=x.shape[1] 49 50 th1=weight[0] 51 th2=weight[1] 52 a1=np.insert(x,0,values=np.ones(m),axis=0)###添加偏置项 53 z2=th1*a1 54 a2=np.insert(sigmoid(z2),0,values=np.ones(m),axis=0)###添加偏置项 55 z3=th2*a2 56 h=sigmoid(z3) 57 58 A=np.multiply((-y),np.log(h)) 59 B=np.multiply((1-y),np.log(1-h)) 60 J=np.sum(A-B)/m 61 62 del1=np.zeros(th1.shape) 63 del2=np.zeros(th2.shape) 64 print(J) 65 a1t=a1###(401,t) 66 z2t=z2###(30,t) 67 a2t=a2###(31,t) 68 ht=h ###(10,t) 69 yt=y ###(10,t) 70 71 d3t=ht-yt 72 z2t=np.insert(z2t,0,values=np.ones(m),axis=0)###(31,t) 73 d2t=np.multiply(th2.T*d3t,sigmoidGradient(z2t)) 74 75 th1=(d2t[1:,:])*a1t.T 76 th2=d3t*a2t.T 77 print(th2.shape,d3t.shape,z2t.shape) 78 """ 79 del1=del1+(d2t[1:,:])*a1t.T 80 del2=del2+d3t*a2t.T 81 82 del1=del1/m 83 del2=del2/m 84 85 del1[1:,:]=del1[1:,:]#+(th1[1:,:]*alpha)/m 86 del2[1:,:]=del2[1:,:]#+(th2[1:,:]*alpha)/m 87 """ 88 tmp=np.concatenate((np.ravel(th1),np.ravel(th2))) 89 90 return J,tmp 91 92 def getStr(x): 93 y=maxlen-len(str(x)) 94 return "0"*y+str(x) 95 96 def getdata(): 97 file=open("dig/NUM.txt") 98 m=int(file.read()) 99 print(m) 100 x=np.zeros((size_input,m)) 101 y=np.zeros((size_output,m)) 102 tot=0 103 for i in range(m): 104 path="dig/"+getStr(i)+".txt" 105 file=open(path) 106 Q=file.read() 107 A=Q.split(",") 108 for j in range(size_input): 109 x[j,tot]=int(A[j]) 110 y[int(A[size_input]),tot]=1 111 tot+=1 112 return np.matrix(x),np.matrix(y) 113 114 def main(): 115 global theta 116 init() 117 X,Y=getdata() 118 119 fmin=minimize(fun=BP,x0=theta,args=(size,X,Y,1),method='TNC',jac=True,options={'maxiter':300}) 120 theta=fmin.x 121 prediction=forwardProp(theta,size,X) 122 prediction=np.matrix(np.argmax(prediction,axis=0)) 123 print(prediction) 124 true=np.matrix(np.argmax(Y,axis=0)) 125 print(true) 126 127 if(__name__=="__main__"): 128 main()
运行速度起码快了十倍。
需要特别注意的是，代码中a2t与z2t的效果是不一样的！a2t是第二层进行sigmoid后的结果，而z2t是没有sigmoid的输入值。如果混为一谈，那么你就会得到一个虚假的神经网络。
结构化的神经网络，向量化，能自定义层数、节点数，有正则化：
1 import numpy as np 2 #import pandas as pd 3 #import matplotlib.pyplot as plt 4 #import matplotlib 5 #from scipy.io import loadmat 6 #from sklearn.preprocessing import OneHotEncoder 7 from scipy.optimize import minimize 8 #from sklearn.metrics import classification_report 9 10 size_input=20*20 11 size_output=10 12 size=[size_input,30,15,size_output] 13 theta=[] 14 maxlen=6 15 16 def init():###随机初始化 17 global theta 18 tot=0 19 for i in range(len(size)-1): 20 tot+=size[i+1]*(size[i]+1) 21 theta=(np.array(np.random.random(size=tot))-0.5)*0.24 22 23 def sigmoid(x): 24 return 1/(1+np.exp(-x)) 25 26 def getWeight(theta,size):###获得矩阵 27 weight=[] 28 tot=0 29 for i in range(len(size)-1): 30 weight.append(np.matrix(np.array(theta[tot:tot+size[i+1]*(size[i]+1)])).reshape((size[i+1],size[i]+1))) 31 tot+=size[i+1]*(size[i]+1) 32 return weight 33 34 def forwardProp(theta,size,x):###前向传播 35 weight=getWeight(theta,size) 36 m=x.shape[1] 37 layer=len(size) 38 for i in range(layer-1): 39 x=np.insert(x,0,values=np.ones(m),axis=0)###0表示行 40 x=sigmoid(weight[i]*x) 41 return x 42 43 def sigmoidGradient(x): 44 return np.multiply(sigmoid(x),sigmoid(1-x)) 45 46 def BP(theta,size,x,y,alpha): 47 weight=getWeight(theta,size) 48 m=x.shape[1] 49 layer=len(size) 50 ###Z：到达这一层的输入，A：sigmoid后的输入 51 tmpA=[]###注意，下标从0开始 52 tmpZ=[] 53 for i in range(layer-1): 54 x=np.insert(x,0,values=np.ones(m),axis=0) 55 ###0表示行，添加偏置项 56 tmpA.append(x) 57 x=weight[i]*x 58 tmpZ.append(x) 59 x=sigmoid(x) 60 61 A=np.multiply((-y),np.log(x)) 62 B=np.multiply((1-y),np.log(1-x)) 63 reg=0###正则化项 64 for i in range(layer-1): 65 reg+=np.sum(np.power(weight[i][:,1:],2))*alpha/(2*m) 66 J=np.sum(A-B)/m+reg###总代价 67 print(J) 68 69 delta=x-y 70 for I in range(layer-1): 71 i=layer-I-1 72 if(i!=1): 73 z=np.insert(tmpZ[i-2],0,values=np.ones(m),axis=0) 74 ###获得第i-1层的没sigmoid的输入 75 tmpW=delta*tmpA[i-1].T 76 ###得到cost关于w的偏导 77 delta=np.multiply(weight[i-1].T*delta,sigmoidGradient(z)) 78 ###得到新的delta 79 delta=np.delete(delta,0,axis=0) 80 ###删掉偏置项 81 weight[i-1][:,0]=0 82 weight[i-1]*=alpha/m 83 weight[i-1]+=tmpW/m 84 else: 85 tmpW=delta*tmpA[i-1].T 86 weight[i-1][:,0]=0 87 weight[i-1]*=alpha/m 88 weight[i-1]+=tmpW/m 89 90 wait=[] 91 for i in range(layer-1): 92 wait=np.concatenate((wait,np.ravel(weight[i]))) 93 return J,wait 94 95 def getStr(x): 96 y=maxlen-len(str(x)) 97 return "0"*y+str(x) 98 99 def getdata(): 100 file=open("dig/NUM.txt") 101 m=int(file.read()) 102 x=np.zeros((size_input,m)) 103 y=np.zeros((size_output,m)) 104 tot=0 105 for i in range(m): 106 path="dig/"+getStr(i)+".txt" 107 file=open(path) 108 Q=file.read() 109 A=Q.split(",") 110 for j in range(size_input): 111 x[j,tot]=int(A[j]) 112 y[int(A[size_input]),tot]=1 113 tot+=1 114 return np.matrix(x),np.matrix(y) 115 116 def main(): 117 global theta 118 init() 119 X,Y=getdata() 120 121 fmin=minimize(fun=BP,x0=theta,args=(size,X,Y,1),method='TNC',jac=True,options={'maxiter':1000}) 122 theta=fmin.x 123 print(theta) 124 125 prediction=forwardProp(theta,size,X) 126 prediction=np.matrix(np.argmax(prediction,axis=0)) 127 print(prediction) 128 true=np.matrix(np.argmax(Y,axis=0)) 129 print(true) 130 131 if(__name__=="__main__"): 132 main()
要注意的是，当BP算法到最后一层时，对tmpZ的操作是没有用的，因此要特判。
但我发现一个奇怪的问题，当alpha=1时（正则化系数）就难以很快地梯度下降了，目前原因不明。