机器学习之用逻辑回归实现手写数字识别

问题简介

本问题的训练集为：X（5000，400）代表5000个有400个特征的图像样本；Y（5000，1）代表5000个图像分别对应的数字值。
本问题的训练目标是：实现图片中手写数字的自动识别归类。

输入矩阵X的直观展示

给定的X数组的每一行代表一张图片，400个特征代表这张图片的20*20个灰度值。现通过下面的几行程序以第一行为例将这400个灰度值拼接成原图形

    #####################################
    #        展示第一个样本的图像           #
    X1 = X[0,:]
    X1_show = X1.reshape(20,20)
    plt.imshow(X1_show.T)
    plt.show()
    #        展示第一个样本的图像           #
    #####################################

运行效果如下，该图就是将X第一行转化为20*20的灰度值矩阵后显示出来的图片。

下面将X中的5000行按同样的方法显示出来，拼接为一个大的图像：

    #####################################
    #        展示所有样本的图像           #
    for i in range(100):
        X_show_row = False
        temp = 1
        for j in range(50):
            if temp:
                X_show_row = X[j+50*i,:].reshape(20,20)
                temp = 0
            else:
                X_show_row = np.concatenate((X_show_row,X[j+50*i,:].reshape(20,20)),axis=1)
        if type(X_show) == bool:
            X_show = X_show_row
        else:
            X_show = np.concatenate((X_show,X_show_row))

    plt.imshow(X_show.T)
    plt.show()
    #        展示所有样本的图像           #
    #####################################

运行效果如下所示：

输出向量Y的处理

在多分类问题中，因为逻辑回归的输出为0-1之间的某个数，在这里，我们将0-9中的某个数字对应到一个由0和1组成的向量上。于是原来代表所有样本输出的Y（5000，1）向量变为了一个（5000，10）的矩阵Ymatrix。相应的代码如下:

    #####################################
    #        将Y转换为矩阵Ymatrix         #
    i_temp = np.zeros((1,10))
    Ymatrix = False
    for i in range(5000):
        i_temp[...,Y[i]-1] = 1
        if type(Ymatrix) == bool:
            Ymatrix = i_temp
        else:
            Ymatrix = np.concatenate((Ymatrix,i_temp),axis=0)
        i_temp[..., Y[i] - 1] = 0

    Y_temp = Ymatrix[...,9]
    Y_temp = np.expand_dims(Y_temp, axis=1)
    Y_temp1 = Ymatrix[...,0:9]
    Ymatrix = np.concatenate((Y_temp,Y_temp1),axis=1)
    print(Ymatrix)
    #        将Y转换为矩阵Ymatrix          #
    #####################################

梯度下降的目的分析

在正式进入梯度下降训练阶段之前，我们要先想明白为何要进行梯度下降训练，这个训练到底要通过什么方式达到什么目的。
先来看看如何从特征和权值得到预估输出向量。下图是使用逻辑回归解决多分类问题的训练模型。实际上，逻辑回归算法可以看成是仅有输入层和输出层的神经网络，其结构图如下所示：

每个样本在进行上图的操作时，400个特征与特征矩阵θ结合后，得到的10个输出写在一起就可以得到预估输出向量hθ(x)。
进行5000次上述操作，得到5000个hθ(x)，把它们写在一起后构成一个（5000，10）的预估输出矩阵。这个预估输出矩阵与Ymatrix矩阵的构成相同。
梯度下降的目的就是通过不断调整权值矩阵θ来让预估输出矩阵与Ymatrix不断接近。

对十个输出分别进行梯度下降

有了梯度下降训练的目标之后，我们就针对预估输出矩阵和Ymatrix的误差，经过求导数之后进行多次迭代，一步步的通过改变权值矩阵将两个的误差缩小，从而找到一个拟合较好的权值矩阵。
这里对每个输出分别进行梯度下降的计算，训练好第一个输出需要的的权值向量之后在训练第二个输出需要的权值向量。这样训练10次之后得到最终的较好的权值矩阵。

    #####################################
    #            梯度下降法迭代            #
    sample_num = X.shape[0]
    feature_num = X.shape[1]
    X = np.insert(X, 0, 1, axis=1)
    trained_theta = False
    alpha = 0.2
    iterations = 10000
    for i in range(10):
        theta = np.zeros((feature_num + 1,))
        for k in range(iterations):
            sigmoid = 1 / (1 + np.exp(-X.dot(theta)))
            gradient = X.T.dot(sigmoid - Ymatrix[...,i]) / X.shape[0]
            theta -= alpha * gradient
        if type(trained_theta) == bool:
            trained_theta = theta.reshape(feature_num+1,1)
        else:
            trained_theta = np.concatenate((trained_theta,theta.reshape(feature_num+1,1)),axis=1)
    #            梯度下降法迭代            #
    #####################################

最后训练出的权值矩阵为：

全部代码展示

import numpy as np
import matplotlib.pyplot as plt
import scipy.io as sio
from sklearn.metrics import classification_report, accuracy_score
#plt.switch_backend('agg')

if __name__ == '__main__':
    data = sio.loadmat('ex3data1.mat')
    X = data['X']
    Y = data['y']


    X_show = False
    #####################################
    #        展示所有样本的图像           #
    for i in range(100):
        X_show_row = False
        temp = 1
        for j in range(50):
            if temp:
                X_show_row = X[j+50*i,:].reshape(20,20)
                temp = 0
            else:
                X_show_row = np.concatenate((X_show_row,X[j+50*i,:].reshape(20,20)),axis=1)
        if type(X_show) == bool:
            X_show = X_show_row
        else:
            X_show = np.concatenate((X_show,X_show_row))

    plt.imshow(X_show.T)
    plt.show()
    #        展示所有样本的图像           #
    #####################################

    #####################################
    #        将Y转换为矩阵Ymatrix         #
    i_temp = np.zeros((1,10))
    Ymatrix = False
    for i in range(5000):
        i_temp[...,Y[i]-1] = 1
        if type(Ymatrix) == bool:
            Ymatrix = i_temp
        else:
            Ymatrix = np.concatenate((Ymatrix,i_temp),axis=0)
        i_temp[..., Y[i] - 1] = 0

    Y_temp = Ymatrix[...,9]
    Y_temp = np.expand_dims(Y_temp, axis=1)
    Y_temp1 = Ymatrix[...,0:9]
    Ymatrix = np.concatenate((Y_temp,Y_temp1),axis=1)
    print(Ymatrix)
    #        将Y转换为矩阵Ymatrix          #
    #####################################

    #####################################
    #            梯度下降法迭代            #
    sample_num = X.shape[0]
    feature_num = X.shape[1]
    X = np.insert(X, 0, 1, axis=1)
    trained_theta = False
    alpha = 0.2
    iterations = 10000
    for i in range(10):
        theta = np.zeros((feature_num + 1,))
        for k in range(iterations):
            sigmoid = 1 / (1 + np.exp(-X.dot(theta)))
            gradient = X.T.dot(sigmoid - Ymatrix[...,i]) / X.shape[0]
            theta -= alpha * gradient
        if type(trained_theta) == bool:
            trained_theta = theta.reshape(feature_num+1,1)
        else:
            trained_theta = np.concatenate((trained_theta,theta.reshape(feature_num+1,1)),axis=1)
    #            梯度下降法迭代            #
    #####################################

    print (trained_theta)