[机器学习]第四、五周记录

这篇记录的内容来自于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_{\theta}(x^{(i)}))_k)+(1-y_k^{i})log(1-(h_{\theta}(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_{\theta}{(x^{(i)})}-y^{(i)})=\frac{1}{m}\sum_{i=1}^{m}(g(x^{(i)})-y^{(i)})$$

接着，

$$\delta^{(l-1)}=(\Theta^{(l-1)})^T\delta^{(l)}.*g'(z^{(l-1)})$$其中$z^{(l)}$指的是第l层在前向传播时获得的输入，.*是点乘。

注意到，

$$g'(z^{(l)})=a^{(l)}.*(1-a^{(l)})$$

并且，

$$\frac{\partial}{\partial\Theta^{(l)}_{i,j}}J(\Theta)=a^{(l)}_j\delta_i^{(l+1)}$$

具体推导此处省略。你可以观看3B1B的相关视频，来进行感性的理解。 

---

神经网络的实现还是有点麻烦的。

我们先实现只有一层hidden layer的神经网络。

例子

对某一个人的手写数字进行识别。

样本：200张20*20的灰度图片。方便起见，我们事先将它们变成了用逗号隔开的长度为401的txt文件，对于前400个数字，由0/1代表对应位置是黑色还是白色，最后一位是它对应的数字。

样本下载：https://pan.baidu.com/s/1dEL27nMZzMbT3iAPTiJKqA
提取码：nq18

最后实现代码：

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import matplotlib
from scipy.io import loadmat
from sklearn.preprocessing import OneHotEncoder
from scipy.optimize import minimize
from sklearn.metrics import classification_report###评价报告
"""
输入数据为20*20的灰度图像，标签为1~10
只有一层hidden layer
x:样本
y:标签
m:样本数
th:参数
"""
size_input=20*20
size_hidden=30
size_labels=10
alpha=3###就是lambda

#######################################sigmoid函数
def sigmoid(x):
    return 1/(1+np.exp(-x))
#######################################fp
def forprop(x,th1,th2):
    m=x.shape[0]
    a1=np.insert(x,0,values=np.ones(m),axis=1)###添加偏置项
    z2=a1*th1.T
    a2=np.insert(sigmoid(z2),0,values=np.ones(m),axis=1)###添加偏置项
    z3=a2*th2.T
    h=sigmoid(z3)
    return a1,z2,a2,z3,h###现在的h的列数等于size_labels
#######################################cost函数
def cost(th1,th2,size_input,size_hidden,size_labels,x,y,alpha):
    m=x.shape[0]
    x=np.matrix(x)
    y=np.matrix(y)
    a1,z2,a2,z3,h=forprop(x,th1,th2)
    A=np.multiply((-y),np.log(h))###数乘
    B=np.multiply((1-y),np.log(1-h))
    reg=(np.sum(np.power(th1[:,1:],2))+np.sum(np.power(th2[:,1:],2)))*alpha/(2*m)
    return np.sum(A-B)/m+reg
#######################################sigmoid导数
def sigmoidGradient(x):
    return np.multiply(sigmoid(x),1-sigmoid(x))
#######################################随机初始化
def random():
    return (np.random.random(size=size_hidden*(size_input+1)+size_labels*(size_hidden+1))-0.5)*0.24
#######################################bp

TOT=0
def backprop(th,size_input,size_hidden,size_labels,x,y,alpha):
    global TOT
    m=x.shape[0]
    print(TOT)
    TOT=TOT+1
    x=np.matrix(x)
    y=np.matrix(y)
    th1=np.matrix(np.reshape(th[:size_hidden*(size_input+1)],(size_hidden,(size_input+1))))###请注意
    th2=np.matrix(np.reshape(th[size_hidden*(size_input+1):],(size_labels,(size_hidden+1))))
    a1,z2,a2,z3,h=forprop(x,th1,th2)
    A=np.multiply((-y),np.log(h))
    B=np.multiply((1-y),np.log(1-h))
    reg=(np.sum(np.power(th1[:,1:],2))+np.sum(np.power(th2[:,1:],2)))*alpha/(2*m)
    J=np.sum(A-B)/m+reg
    
    del1=np.zeros(th1.shape)
    del2=np.zeros(th2.shape)
    
    for t in range(m):
        a1t=a1[t,:]###(1,401)
        z2t=z2[t,:]###(1,25)
        a2t=a2[t,:]###(1,26)
        ht=h[t,:]  ###(1,10)
        yt=y[t,:]  ###(1,10)
        
        d3t=ht-yt
        z2t=np.insert(z2t,0,values=np.ones(1))###(1,26)，老传统
        d2t=np.multiply((th2.T*d3t.T).T,sigmoidGradient(z2t))
        
        del1=del1+(d2t[:,1:]).T*a1t
        del2=del2+d3t.T*a2t
    del1=del1/m
    del2=del2/m
    
    del1[:,1:]=del1[:,1:]+(th1[:,1:]*alpha)/m
    del2[:,1:]=del2[:,1:]+(th2[:,1:]*alpha)/m
    
    tmp=np.concatenate((np.ravel(del1),np.ravel(del2)))
    
    return J,tmp
#######################################获得数据
maxlen=6
def getStr(x):
    y=maxlen-len(str(x))
    return "0"*y+str(x)
    
def getdata():
    file=open("dig/NUM.txt")
    m=int(file.read())
    x=np.zeros((m,size_input))
    y=np.zeros((m,size_labels))
    for i in range(m):
        path="dig/"+getStr(i)+".txt"
        file=open(path)
        Q=file.read()
        A=Q.split(",")
        for j in range(size_input):
            x[i,j]=int(A[j])
        y[i,int(A[size_input])]=1
    print(y)
    return x,y
#######################################预测
def predict(x,th1,th2):
    a1,z2,a2,z3,h=forprop(x,th1,th2)
    print(np.matrix(np.argmax(h,axis=1)))
#######################################主函数
def main():
    x,y=getdata()
    th=random()
    fmin=minimize(fun=backprop,x0=th,args=(size_input,size_hidden,size_labels,x,y,alpha),method='TNC',jac=True,options={'maxiter':1000})
    G=fmin.x
    file=open("save++.txt","w")
    for i in range(G.shape[0]):
        file.write(str(G[i])+"\n")
    file.close()
    x = np.matrix(x)
    th1 = np.matrix(np.reshape(fmin.x[:size_hidden* (size_input + 1)], (size_hidden, (size_input + 1))))
    th2 = np.matrix(np.reshape(fmin.x[size_hidden * (size_input + 1):], (size_labels, (size_hidden + 1))))
    a1, z2, a2, z3, h = forprop(x, th1, th2 )
    y_pred = np.matrix(np.argmax(h, axis=1))
    m=x.shape[0]
    Y=np.zeros(m)
    for i in range(m):
        for j in range(size_labels):
            if(y[i,j]):
                Y[i]=j
    print(classification_report(Y,y_pred))

if(__name__=="__main__"):
    main()

 但很显然，没有向量化运行速度太慢，没有结构化修改太复杂！因此博主重构了代码。

接下来的代码中，一个样本是按列排序，而不是上面按行排序（因为博主从头到尾推了一遍）。

这份代码是向量化之后的神经网络，只有一层hidden layer，没有正则项：

import numpy as np
#import pandas as pd
#import matplotlib.pyplot as plt
#import matplotlib
#from scipy.io import loadmat
#from sklearn.preprocessing import OneHotEncoder
from scipy.optimize import minimize
#from sklearn.metrics import classification_report

size_input=20*20
size_output=10
size=[size_input,30,size_output]
theta=[]
maxlen=6

def init():###随机初始化
    global theta
    tot=0
    for i in range(len(size)-1):
        tot+=size[i+1]*(size[i]+1)
    theta=(np.array(np.random.random(size=tot))-0.5)*0.24

def sigmoid(x):
    return 1/(1+np.exp(-x))

def getWeight(theta,size):###获得矩阵
    weight=[]
    tot=0
    for i in range(len(size)-1):
        weight.append(np.matrix(np.array(theta[tot:tot+size[i+1]*(size[i]+1)])).reshape((size[i+1],size[i]+1)))
        tot+=size[i+1]*(size[i]+1)
    return weight

def forwardProp(theta,size,x):###前向传播
    weight=getWeight(theta,size)
    m=x.shape[1]
    layer=len(size)
    for i in range(layer-1):
        x=np.insert(x,0,values=np.ones(m),axis=0)###0表示行
        x=sigmoid(weight[i]*x)
    return x

def sigmoidGradient(x):
    return np.multiply(sigmoid(x),1-sigmoid(x))

def BP(theta,size,x,y,alpha):
    weight=getWeight(theta,size)
    m=x.shape[1]
    
    th1=weight[0]
    th2=weight[1]
    a1=np.insert(x,0,values=np.ones(m),axis=0)###添加偏置项
    z2=th1*a1
    a2=np.insert(sigmoid(z2),0,values=np.ones(m),axis=0)###添加偏置项
    z3=th2*a2
    h=sigmoid(z3)
    
    A=np.multiply((-y),np.log(h))
    B=np.multiply((1-y),np.log(1-h))
    J=np.sum(A-B)/m
    
    del1=np.zeros(th1.shape)
    del2=np.zeros(th2.shape)
    print(J)
    a1t=a1###(401,t)
    z2t=z2###(30,t)
    a2t=a2###(31,t)
    ht=h  ###(10,t)
    yt=y  ###(10,t)
    
    d3t=ht-yt
    z2t=np.insert(z2t,0,values=np.ones(m),axis=0)###(31,t)
    d2t=np.multiply(th2.T*d3t,sigmoidGradient(z2t))
    
    th1=(d2t[1:,:])*a1t.T
    th2=d3t*a2t.T
    print(th2.shape,d3t.shape,z2t.shape)
    """
    del1=del1+(d2t[1:,:])*a1t.T
    del2=del2+d3t*a2t.T

    del1=del1/m
    del2=del2/m
    
    del1[1:,:]=del1[1:,:]#+(th1[1:,:]*alpha)/m
    del2[1:,:]=del2[1:,:]#+(th2[1:,:]*alpha)/m
    """
    tmp=np.concatenate((np.ravel(th1),np.ravel(th2)))
    
    return J,tmp

def getStr(x):
    y=maxlen-len(str(x))
    return "0"*y+str(x)

def getdata():
    file=open("dig/NUM.txt")
    m=int(file.read())
    print(m)
    x=np.zeros((size_input,m))
    y=np.zeros((size_output,m))
    tot=0
    for i in range(m):
        path="dig/"+getStr(i)+".txt"
        file=open(path)
        Q=file.read()
        A=Q.split(",")
        for j in range(size_input):
            x[j,tot]=int(A[j])
        y[int(A[size_input]),tot]=1
        tot+=1
    return np.matrix(x),np.matrix(y)

def main():
    global theta
    init()
    X,Y=getdata()
    
    fmin=minimize(fun=BP,x0=theta,args=(size,X,Y,1),method='TNC',jac=True,options={'maxiter':300})
    theta=fmin.x
    prediction=forwardProp(theta,size,X)
    prediction=np.matrix(np.argmax(prediction,axis=0))
    print(prediction)
    true=np.matrix(np.argmax(Y,axis=0))
    print(true)

if(__name__=="__main__"):
    main()

运行速度起码快了十倍。

需要特别注意的是，代码中a2t与z2t的效果是不一样的！a2t是第二层进行sigmoid后的结果，而z2t是没有sigmoid的输入值。如果混为一谈，那么你就会得到一个虚假的神经网络。

结构化的神经网络，向量化，能自定义层数、节点数，有正则化：

import numpy as np
#import pandas as pd
#import matplotlib.pyplot as plt
#import matplotlib
#from scipy.io import loadmat
#from sklearn.preprocessing import OneHotEncoder
from scipy.optimize import minimize
#from sklearn.metrics import classification_report

size_input=20*20
size_output=10
size=[size_input,30,15,size_output]
theta=[]
maxlen=6

def init():###随机初始化
    global theta
    tot=0
    for i in range(len(size)-1):
        tot+=size[i+1]*(size[i]+1)
    theta=(np.array(np.random.random(size=tot))-0.5)*0.24

def sigmoid(x):
    return 1/(1+np.exp(-x))

def getWeight(theta,size):###获得矩阵
    weight=[]
    tot=0
    for i in range(len(size)-1):
        weight.append(np.matrix(np.array(theta[tot:tot+size[i+1]*(size[i]+1)])).reshape((size[i+1],size[i]+1)))
        tot+=size[i+1]*(size[i]+1)
    return weight

def forwardProp(theta,size,x):###前向传播
    weight=getWeight(theta,size)
    m=x.shape[1]
    layer=len(size)
    for i in range(layer-1):
        x=np.insert(x,0,values=np.ones(m),axis=0)###0表示行
        x=sigmoid(weight[i]*x)
    return x

def sigmoidGradient(x):
    return np.multiply(sigmoid(x),sigmoid(1-x))

def BP(theta,size,x,y,alpha):
    weight=getWeight(theta,size)
    m=x.shape[1]
    layer=len(size)
    ###Z：到达这一层的输入，A：sigmoid后的输入
    tmpA=[]###注意，下标从0开始
    tmpZ=[]
    for i in range(layer-1):
        x=np.insert(x,0,values=np.ones(m),axis=0)
        ###0表示行，添加偏置项
        tmpA.append(x)
        x=weight[i]*x
        tmpZ.append(x)
        x=sigmoid(x)
    
    A=np.multiply((-y),np.log(x))
    B=np.multiply((1-y),np.log(1-x))
    reg=0###正则化项
    for i in range(layer-1):
        reg+=np.sum(np.power(weight[i][:,1:],2))*alpha/(2*m)
    J=np.sum(A-B)/m+reg###总代价
    print(J)
    
    delta=x-y
    for I in range(layer-1):
        i=layer-I-1
        if(i!=1):
            z=np.insert(tmpZ[i-2],0,values=np.ones(m),axis=0)
            ###获得第i-1层的没sigmoid的输入
            tmpW=delta*tmpA[i-1].T
            ###得到cost关于w的偏导
            delta=np.multiply(weight[i-1].T*delta,sigmoidGradient(z))
            ###得到新的delta
            delta=np.delete(delta,0,axis=0)
            ###删掉偏置项
            weight[i-1][:,0]=0
            weight[i-1]*=alpha/m
            weight[i-1]+=tmpW/m
        else:
            tmpW=delta*tmpA[i-1].T
            weight[i-1][:,0]=0
            weight[i-1]*=alpha/m
            weight[i-1]+=tmpW/m
    
    wait=[]
    for i in range(layer-1):
        wait=np.concatenate((wait,np.ravel(weight[i])))
    return J,wait

def getStr(x):
    y=maxlen-len(str(x))
    return "0"*y+str(x)

def getdata():
    file=open("dig/NUM.txt")
    m=int(file.read())
    x=np.zeros((size_input,m))
    y=np.zeros((size_output,m))
    tot=0
    for i in range(m):
        path="dig/"+getStr(i)+".txt"
        file=open(path)
        Q=file.read()
        A=Q.split(",")
        for j in range(size_input):
            x[j,tot]=int(A[j])
        y[int(A[size_input]),tot]=1
        tot+=1
    return np.matrix(x),np.matrix(y)

def main():
    global theta
    init()
    X,Y=getdata()
    
    fmin=minimize(fun=BP,x0=theta,args=(size,X,Y,1),method='TNC',jac=True,options={'maxiter':1000})
    theta=fmin.x
    print(theta)
    
    prediction=forwardProp(theta,size,X)
    prediction=np.matrix(np.argmax(prediction,axis=0))
    print(prediction)
    true=np.matrix(np.argmax(Y,axis=0))
    print(true)

if(__name__=="__main__"):
    main()

要注意的是，当BP算法到最后一层时，对tmpZ的操作是没有用的，因此要特判。

但我发现一个奇怪的问题，当alpha=1时（正则化系数）就难以很快地梯度下降了，目前原因不明。