用python实现手写数字识别

前言

在之前的学习中，已经对神经网络的算法具体进行了学习和了解。现在，我们可以用python通过两种方法来实现手写数字的识别。这两种方法分别是多元逻辑回归和神经网络方法。

用多元逻辑回归手写数字识别

训练样本简介
首先，可以用逻辑回归实现手写数字识别的多元分类，在ex3data.mat文件中保存着一个包含5000个培训考试训练数据集的手写数字，.mat格式表示这些训练数据集以矩阵格式保存，而不是文本格式。训练集中有5000个训练样本，每个训练样本是一个20X20像素的灰度图像，每一个像素用一个float型数据表示当前位置的灰度强度。20X20像素形成一个400维向量，每一个训练样本形成一个单行矩阵 $X$ ，5000个训练样本构成一个5000X400的矩阵，每一个行向量代表手写数字图像的训练样本。矩阵 $X$ 可以如下表示：

数据集第二部分包含一个5000维的矩阵 $y$ , 代表输出值的标签，即0,1,2,3,4,5,6,7,8,9.
如下代码所示，首先导入算法实现所需要的python库

import matplotlib.pyplot as plt
import numpy as np
import scipy.io as scio

初始化一些参数

input_layer_size = 400  # 20x20 像素矩阵
num_labels = 10     # 手写数字 0,1,2,3,4,5,6,7,8,9

数据可视化
随机选取

X

矩阵中的100行向量，并将这些函数传递到displayData函数中，通过调用此函数，我们可以显示出这些向量的灰度图像，具体实现如下所示：

载入数据

  data = scio.loadmat('ex3data1.mat')
  X = data['X']
  y = data['y'].flatten()
  m = y.size

绘图函数实现
绘图函数的实现如下所示，具体实现思路就是确定显示图像的像素大小，然后读取训练集数据，依次将其填充到设定好大小的像素矩阵中，具体可参考代码注释。

def display_data(x):
 (m, n) = x.shape
 # m = 100
 # n = 400

 # 设置每个数字的宽度与高度(像素)
 example_width = np.round(np.sqrt(n)).astype(int)# example_width=20
 example_height = (n / example_width).astype(int) #example_height=20

 # 计算显示的行数与列数
 display_rows = np.floor(np.sqrt(m)).astype(int) #display_rows=10
 display_cols = np.ceil(m / display_rows).astype(int)#display_rows=10

 # 单个图像之间的间隔
 pad = 1

 # 设置并初始化显示像素矩阵的大小211*211 ,1+(10*20+1)
 display_array = - np.ones((pad + display_rows * (example_height + pad),
                           pad + display_rows * (example_height + pad)))

 # 将每个训练样本显示在矩阵中
 curr_ex = 0 
 for j in range(display_rows):
     for i in range(display_cols):
         if curr_ex > m:
             break

         # 每次每行和每列读取一个20*20像素的数字，矩阵大小加21
         # 实际上矩阵形式可以认为 10*10*400(20*20像素)
         max_val = np.max(np.abs(x[curr_ex]))
         display_array[pad + j * (example_height + pad) + np.arange(example_height),
                       pad + i * (example_width + pad) + np.arange(example_width)[:, np.newaxis]] = \
                       x[curr_ex].reshape((example_height, example_width)) / max_val
         curr_ex += 1

     if curr_ex > m:
         break

 # Display image
 plt.figure()
 plt.imshow(display_array, cmap='gray', extent=[-1, 1, -1, 1])
 plt.axis('off')

最后，绘图实现效果如下

矢量化逻辑回归
参考逻辑回归与正则化的相关知识，矢量化逻辑回归的算法实现代码如下所示：

def lr_cost_function(theta, X, y, lmd):
    m = y.size
    cost = 0
    grad = np.zeros(theta.shape)
    reg_theta = theta[1:]
    hy = sigmoid(np.dot(X,theta))
    cost = np.sum(-y * np.log(hy) - np.subtract(1, y) * np.log(np.subtract(1, hy))) / m + (lmd / (2 * m)) * np.sum(reg_theta * reg_theta)
    
    grad = np.dot(X.T,  np.subtract(hy, y)) / m
    grad[1:] = grad[1:] + reg_theta * (lmd / m) #正则化
    return cost, grad
def sigmoid(z):
    return 1/(1+np.exp(-z))

人为构造一些参数值，计算结果如下所示：

theta_t = np.array([-2, -1, 1, 2])
X_t = np.c_[np.ones(5), np.arange(1, 16).reshape((3, 5)).T/10]
y_t = np.array([1, 0, 1, 0, 1])
lmda_t = 3
cost, grad = lr_cost_function(theta_t, X_t, y_t, lmda_t)
print(cost,grad)

参数训练
在之前的代码中，已经实现了逻辑回归的梯度下降算法，利用训练样本，现在可以通过梯度下降算法训练 $\theta$ 参数了，其中，与二元分类不同的是，在本次代码中需要实现的是0-9的数字分类，与用0,1实现二元分类相比，更加复杂，所以在多元分类中，我们用多元逻辑回归算法实现多元分类时，通过0和1来区分每个数字是否识别成功，简而言之，相当于进行了10次二元逻辑回归，只不过每次 $y$ 的值会发生变化而已。具体实现代码如下所示：

def one_vs_all(X, y, num_labels, lmd):
    
    (m, n) = X.shape
    all_theta = np.zeros((num_labels, n + 1))
    X = np.c_[np.ones(m), X]

    for i in range(num_labels):
        initial_theta = np.zeros((n + 1, 1))
        iclass = i if i else 10
        y_i = np.array([1 if x == iclass else 0 for x in y])
        print('Optimizing for handwritten number {}...'.format(i))

        def cost_func(t):
            return lr_cost_function(t, X, y_i, lmd)[0]

        def grad_func(t):
            return lr_cost_function(t, X, y_i, lmd)[1]

        theta, *unused = opt.fmin_cg(cost_func, fprime=grad_func, x0=initial_theta, maxiter=100, disp=False, 
                                     full_output=True)
        print('Done')
        all_theta[i] = theta

    return all_theta

调用one_vs_all()函数，可以看到通过梯度下降算法得到的
$\theta$ 参数矩阵如下所示，

用多元逻辑回归算法实现一对多分类预测
现在，我们已经实现了多元逻辑回归算法，现在可以利用梯度下降算法获得的参数实现一对多分类的具体应用了，具体实现代码可以如下表示：

def predict_one_vs_all(all_theta, X):
    m = X.shape[0]
    num_labels = all_theta.shape[0]
    p = np.zeros(m)
    X = np.c_[np.ones(m), X]
    result = np.dot(all_theta, X.T)
    result = np.roll(result, -1, axis=0)
    result = np.vstack([np.zeros(m), result])
    p = np.argmax(result, axis=0)
    return p
pred = predict_one_vs_all(all_theta, X)
print('Training set accuracy: {}'.format(np.mean(pred == y)*100))

经过计算，识别精确率如下所示：

用神经网路实现手写数字识别

初始化一些参数
根据算法要求，首先需要初始化一些参数，具体代码如下所示：



data = scio.loadmat('ex3weights.mat')  #加载权重矩阵
theta1 = data['Theta1']
theta2 = data['Theta2']

计算神经网络算法识别效果的精确度
利用神经网络实现多元分类时，不再计算损失函数和梯度，利用输入层，隐藏层和输出层，我们可以直接利用神经网络实现手写数字识别的多元分类。具体实现代码如下所示：

def predict(theta1, theta2, x):
    num_labels = theta2.shape[0]
    p = np.zeros(m)
    x = np.c_[np.ones(m), x]
    h1 = sigmoid(np.dot(x, theta1.T))
    h1 = np.c_[np.ones(h1.shape[0]), h1]
    h2 = sigmoid(np.dot(h1, theta2.T))
    p = np.argmax(h2, axis=1) + 1

    return p

分析以上代码，我们可以发现这是有一个输入层（ $X$ 矩阵），一个
隐藏层(由 $X$ 矩阵计算得到的 $h1$ ），由隐藏层计算得到的输出层（ $h2$ ）,最后，由输出层得到分类结果。
通过神经网络得到识别的精确率如下所示：