机器学习反向传播算法实现分类 吴恩达作业题

import matplotlib.pyplot as plt
import numpy as np
import scipy.io as sio #此为载入matlab数据所需要的包
import matplotlib #作图函数
import scipy.optimize as opt #此为优化函数所需要的包
#提取数据函数
def load_data(path,transfer=True):
    data=sio.loadmat(path)
    y=data.get('y')
    print(y.shape)
    y=y.reshape(y.shape[0])
    print(y.shape)
    X=data.get('X')

    if transfer: #加上这个主要是为了让方向变得正确，原先像素为20*20的图像方向是错误的
        X=np.array([im.reshape((20,20)).T for im in X])
        X=np.array([im.reshape(400) for im in X])
    return X,y
X,y=load_data('ex4data1.mat')
#获取可视的图像
def plot_100_image(X):
    sample_number=np.random.choice(np.arange(X.shape[0]),100)#此为采样函数
    #np.randon.choice(array,size,p,replace) array为采样样本，size为采样的数量，p随机采样的概率分布，replace为是否放回
    sample_plot=X[sample_number,:]
    fig,ax_array=plt.subplots(nrows=10,ncols=10,sharey=True,sharex=True,figsize=(8,8))
    for r in range(10):
        for c in range(10):
            ax_array[r,c].matshow(sample_plot[r*10+c].reshape((20,20)),cmap='gray_r')
            plt.yticks([])#取消y的标度
            plt.xticks([])
plot_100_image(X)
plt.show()
#代价函数部分
X_raw,y_raw=load_data('ex4data1.mat',transfer=False)
X=np.insert(X_raw,0,values=1,axis=1)
#借鉴的向量化标签的说法，也就是y的原先shape为（5000,1），我们要将数字分类为1-10，所以我们将y的shape扩展为（5000,10），每一行的列的数字1就是
#y的值
def expand_f(y):
    res=[]#创建一个空的数组
    count=0
    for i in y:#for循环中若y为矩阵那么挨个提取y的列元素
        count=count+1
        y_array=np.zeros(10)#zeros创建0矩阵zero(1)创建一个零矩阵一维向量，内有一个0(一维指的是对列来说)
        y_array[i-1]=1#因为创建的是0-9而我们提取出来的数值为1-10所以我们存进去的时候要减1
        res.append(y_array)#在res的末尾加上我们新设立的一维矩阵总共5000个数值
    print(count)
    return np.array(res)
y=expand_f(y_raw)
#读取权重矩阵
data=sio.loadmat('ex4weights.mat')
Theta1=data['Theta1']
Theta2=data['Theta2']
print(Theta2.shape)
print(Theta1.shape)
def serialize(a,b):#这个地方为什么要将矩阵行列值拉升呢，因为之后函数优化时传入的数据需要用的拉伸
    return np.concatenate((np.ravel(a),np.ravel(b)))#np.ravel将数据变为一维数组
theta=serialize(Theta1,Theta2)
def deserialize(seq):#在代价函数和梯度下降算法中用到的权重矩阵还是需要正常的行列值，所以需要写回复函数
    return seq[:25*401].reshape(25,401),seq[25*401:].reshape(10,26)
def sigmod(z):#sigmod函数，分类里非常常用，线性回归里用的是普通的二次三次函数
    return 1/(1+np.exp(-z))
def feed_forward(theta,X):#前向传播算法
    theta1,theta2=deserialize(theta)#将一维矩阵恢复
    a1=X#输入层
    z2=a1@theta1.T#z2乘上权重矩阵
    a2=np.insert(sigmod(z2),0,values=1,axis=1)#查询到下一个权重矩阵的列26列所以这里需要加上一个数值为1的一列值
    z3=a2@theta2.T
    h=sigmod(z3)#输出层
    return a1,z2,a2,z3,h
def cost(theta,X,y):#代价函数（未正则化情况下的代价函数）
    _,_,_,_,h=feed_forward(theta, X)
    compute=-np.multiply(y,np.log(h))-np.multiply((1-y),np.log(1-h))#很正常的分类函数里的代价函数的公式
    return compute.sum()/X.shape[0]
cost(theta,X,y)
def regularised_cost(theta,X,y,l=1):#代价函数的正则化
    t1,t2=deserialize(theta)
    m=X.shape[0]
    reg1=(l/(2*m))*np.power(t1[:,1:],2).sum()#因为这里的权重矩阵是有两层的，而正则化的时候需要正则化所有的权重矩阵
    reg2=(l/(2*m))*np.power(t2[:,1:],2).sum()
    return cost(theta,X,y)+reg1+reg2
print(regularised_cost(theta,X,y))
def sigmod_grident(z):#这是一会反向传播中需要用到sigmod函数的求导之后的公式
    return np.multiply(sigmod(z),1-sigmod(z))
#重头戏来了 神经网络的梯度下降函数
def grident(theta,X,y):
    t1,t2=deserialize(theta)
    m=X.shape[0]
    delta1=np.zeros(t1.shape)#这是两个需要求的矩阵，因为系统给的权重矩阵是有误差的，我们需要求出最好的权重矩阵，代价最低的权重矩阵
    delta2=np.zeros(t2.shape)#由于又有隐藏层，所以我们在更新与输出层直接相连的权重矩阵时如果不使用反向传播算法是不能更新输入层与隐藏层直接的权重矩阵的
    a1,z2,a2,z3,h=feed_forward(theta,X)
    for i in range(m):#开始一个一个算各个数值的误差，总共五千个数值
        a1i=a1[i,:]#92-96行代表着我门提取第i行的各个数值[1,401]
        z2i=z2[i,:]#[1,25]
        a2i=a2[i,:]#[1,26]
        hi=h[i,:]
        yi=y[i,:]
        d3i=hi-yi#最末层的误差值
        z2i=np.insert(z2i,0,np.ones(1))#因为z2i的列是25，而权重矩阵d2i的列数为26，所以增加一层
        d2i=np.multiply(t2.T@d3i,sigmod_grident(z2i))
        delta2+=np.matrix(d3i).T@np.matrix(a2i)#反向传播算法中的更新权重矩阵的值
        delta1+=np.matrix(d2i[1:]).T@np.matrix(a1i)
    delta1=delta1/m
    delta2=delta2/m
    return serialize(delta1,delta2)
d1,d2=deserialize(grident(theta,X,y))
print(d1.shape,d2.shape)
def gridient_regularized(theta,X,y,l=1):#正则化反向传播算法
    d1,d2=deserialize(grident(theta,X,y))
    m=X.shape[0]
    t1,t2=deserialize(theta)
    t1[:,0]=0#因为第一列的不需要正则化，怕直接让整个式子没了
    reg1=(l/m)*t1
    d1+=reg1
    t2[:,0]=0
    reg2=(l/m)*t2
    d2+=reg2
    return serialize(d1,d2)
def init_theta(size):#随机化你的权重矩阵
    return np.random.uniform(-0.12,0.12,size)
def nn_train(X,y):#最后的最后就是运行你的优化函数了，slow
    theta=init_theta(10285)
    res=opt.minimize(fun=regularised_cost,x0=theta,args=(X,y,1),method='TNC',jac=gridient_regularized)
    return res