吴恩达机器学习第三次编程作业-多元分类与神经网络

机器学习

课时67编程作业：多元分类与神经网络

工具：Pycharm2019.01，Python3.5

参考资料：Python库之SciPy，逻辑回归解决多元分类，
Python3 求最大/小值及索引值 Numpy，scipy.optimize.fmin_ncg，scipy.optimize，numpy矩阵乘法运算

多元分类完整代码: fmin_ncg
多元分类完整代码: fmin_cg
神经网络前向传播算法完整代码

注：本章所用参考资料均已给出

本次练习实现一对多逻辑回归与神经网络来识别手写数字。

Multi-class Classification

数据集格式

本次练习数据集格式为.mat类型，故采用Scipy.io.loadmat方式读入数据，以字典dictionary形式存在，其中有两个键’X’和’y’分别代表训练样本输入和输出，对应值是以ndarray形式存在。
其中，输入是一个(5000,400)的数组，400表示数字是20*20像素的，被拉伸为400维一行。
获取数据：

1
2
3
4
5
6
def Get_Data():
    data=sio.loadmat('./ex3/ex3data1.mat') 
    # print(type(data)) # dict
    # for key in data: # 如果不清楚文件内容，可以打印出来看看
    #     print(key,':',data[key])
    return data['X'],data['y'],data['y'].shape[0] # 输入，输出，以及样本数

Visualizing the data

由于数据集过大，故随机选取100个样本输入打印出来。
每次从[0,5000)中随机取一个数：rand_choose=np.random.random_integers(0,m) 

1
2
3
4
5
6
7
8
9
10
def Plot_Data(X,y,m): # random choose 100 samples to plot it
    fig=plt.figure(figsize=(6,6)) # 新建画布并指定size
    fig.subplots_adjust(left=0, right=1, bottom=0, top=1, hspace=0.05, wspace=0.05) # 指定显示比例
    for i in range(1,101): # 区域编号从1开始
        ax = fig.add_subplot(10, 10, i, xticks=[], yticks=[])
        rand_choose=np.random.random_integers(0,m) 
     ax.imshow(X[rand_choose].reshape(20,20).T,cmap=plt.cm.binary,interpolation='nearest')
        # label the image with the target value
        ax.text(0,20, str(y[rand_choose])) # 显示真实数字
    plt.show()

Figure_11.png

Vectorizing Logistic Regression

主要是定义代价函数与梯度下降，然后调用scipy.optimize库方法跑出参数θ，最后利用参数进行预测。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
def Cost_Func(theta,X,Y,m): # 原来的代价函数
    hypothesis=1/(1+np.exp(-X@theta)) # 假设函数
    Cost_J=-1/m*((Y@np.log(hypothesis)+(1-Y)@np.log(1-hypothesis)))
    return Cost_J

def Cost_reg(theta,X,Y,m,lamda): # 正则化代价函数，加上正则项
    Cost_J=Cost_Func(theta,X,Y,m)
    theta_tmp=theta[1:] # 不惩罚theta[0]
    Cost_J=Cost_J+lamda/(2*m)*np.sum(theta_tmp*theta_tmp)
    return Cost_J

def Gredient_Dec(theta,X,Y,m): # 原来的梯度函数
    hypothesis=1/(1+np.exp(-X.dot(theta))) 
    partial_derivatives=1/m*np.dot(X.T,hypothesis-Y) # 偏导数
    return partial_derivatives

def Gredient_reg(theta,X,Y,m,lamda): # 正则化梯度函数
    partial_derivations=Gredient_Dec(theta,X,Y,m) 
    theta_tmp=lamda/m*theta
    theta_tmp[0]=0 # 不惩罚theta[0]
    return partial_derivations+theta_tmp

One-vs-all Classification
  • 将复杂的分类问题简化为若干个二分类任务，最经典的拆分策略有三种：一对一（OvO）、一对其余（OvR）和多对多（MvM）。
  • 这里的分类问题需要10个分类器，故最终目标需要’训练’出10个θ参数向量。而在训练时，使用逻辑回归方法，将待分类的类别当作1，其余为0。
    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    def One_vs_All_classes(X,Y,theta,labda,m,K):
        y=Y.ravel()
        classes_y = np.zeros(m).astype(int)
        for i in range(K):
            classes_y[y!=i]=0 # 不属于此类设为0，否则为1
            classes_y[y==i]=1 # 变成了简单的逻辑回归问题
            # print(i,sum(classes_y==1),sum(y==i))
            result=opt.fmin_ncg(f=Cost_reg,fprime=Gredient_reg,x0=theta[:,i:i+1],args=(X,classes_y,m,lamda),maxiter=400) # 使用牛顿迭代的方法
            theta[:,i:i+1]=np.mat(result).T # 返回的result是一个向量(401,)
        return theta
One-vs-all Prediction

预测，使用参数θ向量评测每个样本，取Hθ(X)最大的那个作为预测值。

1
2
3
4
5
6
7
8
def One_vs_All_prediction(X,Y,fin_theta,m):
    # print(X.shape,fin_theta.shape)
    predict_matrix=1/(1+np.exp(-X.dot(fin_theta)))
    cnt=0
    for i in range(m):
        if np.argmax(predict_matrix[i])==y[i][0]:
            cnt=cnt+1
    print('%.2f%%' %(cnt/m*100))

在本机上跑了大概10分钟？
预测结果: 96.46%

不同的优化函数对比情况

上述所使用的fmin_ncg()牛顿法原理是利用损失函数二阶导数矩阵即海森矩阵来迭代优化损失函数。但在运行过程中一直打印警告信息提示hassian矩阵不是正数，而且运行时间太长。
在scipy.optimize库中使用其他优化方法试试？
参考资料见顶部开头。
由于一下三种方法的参数列表都相同

1
2
3
fmin_cg(f, x0[, fprime, args, gtol, norm, …]) Minimize a function using a nonlinear conjugate gradient algorithm.
fmin_bfgs(f, x0[, fprime, args, gtol, norm, …]) Minimize a function using the BFGS algorithm.
fmin_ncg(f, x0, fprime[, fhess_p, fhess, …]) Unconstrained minimization of a function using the Newton-CG method.

  1. 故先使用fmin_cg方法，原理上面已经给出：非线性动态梯度算法 .
    结果运行时长1min左右就出来了。
    预测结果：96.46%
    代码见开头

  2. 再把方法改为fmin_bfgs()会显示超过最大迭代次数警告，不过最终预测结果准确率竟然达到: 96.48%

  3. 用opt.minimize(method=method=’L-BFGS-B’)预测结果为：95.44%。
    完整代码为：

    1
    2
    result=opt.minimize(fun=Cost_reg,x0=theta[:,i:i+1],args=(X,classes_y,m),method='L-BFGS-B')
            theta[:,i:i+1]=np.mat(result.x).T

梯度下降不知道怎么嵌入。。。
路过大神望告知！

Neural Networks

目标：实现神经网络识别手写数字。
核心：实现前向传播算法，前向传播算法即神经网络从输入层到隐藏层到输出层的过程。

数据集

数字矩阵还是使用之前多元分类器的数据，这里再给出了第一层到第二层的θ1和第二层到第三层的θ2。其中，θ1是25 x 401的，即第二层有a(2)有25个单元；θ2是10 x 26，26是在第二层多加一个bias unit。

1
2
3
4
5
6
7
def Get_Data():
    data=sio.loadmat('./ex3/ex3data1.mat')
    theta=sio.loadmat('./ex3/ex3weights.mat')
    return data['X'],data['y'],theta['Theta1'],theta['Theta2'],data['y'].size

X,y,theta1,theta2,m=Get_Data()
X=np.hstack((np.ones((m,1)),X))

Model representation

2019-07-28 19-16-38 的屏幕截图.png
大致步骤：
不同的样本添加X0=1表示bias unit，分别通过第一层计算得到第二层
第二层得到的数据再添加一个bias unit，通过Theta2计算得到第三层，取第三层概率值Hθ(X)最大的就是最终的分类预测了。

Feedforward Propagation and Prediction

前向传播与预测。
前向传播就是信息从第一层到最后一层最终的出结果的过程。
预测则是类似多元分类，用的出的参数θ与样本或者其他数据进行运算，看计算结果属于哪一类。
这里可以直接用矩阵运算快速得出结果。
np.argmax()可得到ndarray类型最大值的索引，同理argmin则是取最小值的索引，可指定参数axis，默认为0表示列，1则表示行即每行最大值的索引

1
2
3
4
5
6
def prediction(X,y,theta1,theta2,m):
    # print(theta1.shape,X.shape)
    a_2=1/(1+np.exp(-np.dot(X,theta1.T))) # 直接计算第二层结果
    a_2=np.hstack((np.ones((m,1)),a_2)) # 添加一个bias unit
    a_3=1/(1+np.exp(-a_2.dot(theta2.T)))  # 用样的方法计算第三层
    print('%.2f%%' %(sum(np.argmax(a_3,axis=1)+1==y.ravel())/m*100)) # 加一是y的范围是[1,10]

预测结果：97.52%
完整代码见开头