基于手写数字数据库MNIST,应用最小二乘法模型对数字$0$和其他数字进行分类

1. 目的与要求#

目的： 图像分类问题，基于手写数字数据库 MNIST, 应用最小二乘法模型对数字 $0$ 和其他数字进行分类。

MNIST 数据库介绍： MNIST 数据库来自美国国家标准与技术研究所，包含 70000 张 28×28（784 个像素）的手写数字灰度图像。图像数据已经被转化为 $28 × 28 = 784$ 维的向量形式存储（$784$ 个特征）；标签以 $1$ 维向量形式存储（每个元素表示图片对应的数字）。

要求： 首先用最小二乘模型将 MNIST 数据库中的数字 $0$ 和其他数字进行分类，给出模型在训练集和测试集上的准确率；然后添加 $5000$ 个随机特征重新分类，得出新的准确率；最后对比两次实验结果，给出实验结论。

2. 过程及结果#

2.1 应用最小二乘模型进行分类#

MNIST 数据库将数据分成 2 组，一组是包含 $60000$ 张图像的训练集，另一组是包 $10000$ 张图像的测试集。

特征选择：为了降低维度（维度太高，计算量太大）和噪声（很多像素点的灰度值都是 $0$ 或接近 $0$，没有太多信息），提高模型的效率和准确性，将特征矩阵设置为 $494$ 维的向量，第 $1$ 维是常量 $1$（表示截距或偏置）, 其余 $493$ 维是至少在 $600$ 张图片中像素值不为 $0$ 的像素。如果图像为数字 $0$, 则取 $y= 1$, 否则取 $y= −1$。
最下二乘法的标准模型为：

$Ax=b$，其最小二乘解为：$x=(A^TA)^{-1} A^Tb$
- $A$ 矩阵为特征矩阵，通过训练数据来构造
- $b$ 向量是对应的训练数据数字标签
这里不直接求解其最小二乘解，采用 QR 分解法求解最小二乘解，以提高效率。

2.1.1 具体步骤#

从 MINIST 数据库中加载已经划分好的训练数据集和测试数据集。

包括以下函数（列写函数名称，函数具体内容见文后完整代码）：

decode_idx3_ubyte(idx3_ubyte_file):解析idx3文件的通用函数

decode_idx1_ubyte(idx1_ubyte_file):解析idx1文件通用函数

load_train_images(idx_ubyte_file=train_images_idx3_ubyte_file):返回训练数据

load_train_labels(idx_ubyte_file=train_labels_idx1_ubyte_file):返回训练数据标签

load_test_images(idx_ubyte_file=test_images_idx3_ubyte_file):返回测试数据

load_test_labels(idx_ubyte_file=test_labels_idx1_ubyte_file):返回测试数据标签
数据预处理：

维度转换：将原始数据集的标签和图像数组进行 reshape 操作，将三维数组转换为二维数组。

标签转换：将标签转换为 $1$ 和 $-1$，如果数字为 $0$，则转换为 $1$，否则为 $-1$。

特征选择：找出至少在 $600$ 张图片中不为 0 的像素位置，提高模型的效率和准确度

将特征矩阵 $\mathcal {train_images}$ 设置为 $494$ 维的向量，构建新的数据集：

① 找出至少在 $600$ 张图片中不为 0 的像素位置，并记录哪些列被选做特征。

② 创建一个用于存储训练图像特征的数组 train_image_feature，并将具有足够非零元素的列复制到该数组的第 j 列。

③ 创建一个用于存储测试图像特征的数组 test_image_feature，并使用与 train_image_feature 相同的索引从 test_images 中选择相同的列作为特征。

④ 在训练集和测试集特征向量的开头添加常量 $1$ 作为偏置。

⑤ 返回处理后的训练和测试数据 train_labels, test_labels, train_images, test_images。

def pretreat(train_labels, test_labels, train_images, test_images):
  """
  数据预处理函数：特征选择
  :param 原始数据集
  :return: 处理好的新数据集
  """
  # 将三维数组转换为二维数组
  train_labels = train_labels.reshape(60000)
  test_labels = test_labels.reshape(10000)
  train_images = train_images.reshape(60000, 784)
  test_images = test_images.reshape(10000, 784)

  #  把标签转换为1和-1，数字为0则为1，否则为-1
  # np.where(condition,x,y):如果condition为真，返回x,否则返回y
  train_labels = np.where(train_labels == 0, 1, -1)
  test_labels = np.where(test_labels == 0, 1, -1)

  # 找出至少在600张图片中不为0的像素位置
  j = 0
  index = []
  train_image_feature = np.zeros([60000, 493])  # 用于存储训练图像的特征
  for i in range(784):
      non_zero = np.linalg.norm(train_images[:, i], ord=0)  # 计算每一列中非零元素的个数
      if non_zero >= 600:
          train_image_feature[:, j] = train_images[:, i]  # 将具有足够非零元素的列复制到train_image_feature的第j列
          j = j + 1
          index.append(i)  # 记录哪些列被选做特征
  test_image_feature = np.zeros([10000, 493])  # 创建一个零数组，用于存储测试图像的特征
  test_image_feature = test_images[:, index]  # 使用 index 从 test_images 中选择相同的列作为特征，并赋值给 test_image_feature

  # 在训练集和测试集特征向量开头添加常量1作为偏置
  train_images = np.insert(train_image_feature, 0, 1, axis=1)  
  test_images = np.insert(test_image_feature, 0, 1, axis=1)  

  return train_labels, test_labels, train_images, test_images

训练分类模型及模型评估

采用矩阵分解的最小二乘问题求解算法如下：

① 计算矩阵 $A$ 的 $QR$ 分解，$A=QR$

② 计算 $Q^Tb$

③ 求解 $Rx=Q^Tb$，求得 $x$

函数clf_1(X_train, y_train, X_test, y_test):训练分类模型及模型评估

函数tran1or_0(result):将结果中大于 $0$ 的数转换为 $1$，小于 $0$ 的转换为 $0$

函数show_result(labels, result, dataset):计算错误率及结果展示

def clf_1(X_train, y_train, X_test, y_test):
    """
    :param X_train: 训练数据
    :param y_train: 训练数据标签
    :param X_test: 测试数据
    :param y_test: 测试数据标签
    :return: 
    """
    # 开始训练
    start1 = time.time()
    A = X_train
    A_test = X_test
    b_train = y_train
    b_test = y_test
    print('进行QR分解中...')
    q, r = np.linalg.qr(A)
    print('已完成QR分解')
    print('\n')
    print('对x进行求解中...')
    x = np.linalg.pinv(r).dot(q.T.dot(b_train))
    print('已求得x的解')
    result = A.dot(x)
    print('已完成结果预测')
    tran1or_0(result)
    end1 = time.time()
    # 结果展示分析
    show_result(train_labels, result, 'Train_dataset')
    result_test = A_test.dot(x)
    tran1or_0(result_test)
    show_result(test_labels, result_test, 'Test_dataset')

    print("运行时间:%.2f秒\n" % (end1 - start1))

2.1.2 结果展示#

2.2 添加 5000 个随机特征重新分类#

添加随机特征

利用随机特征法构造 随机特征矩阵：随机矩阵 $R\in R_{5000\times 494}$，$R_{ij}$ 随机取值 $\pm1$，用如下方法计算作为新特征，与原来的特征结合后一共 $5494$ 个特征。

$\max(0,(R_{x})_j),\ j=1,2,3,\ldots,5000$

添加 $5000$ 个随机特征的实际意义是，通过 $\max$ 函数随机的保留特征，如果该特征比 $0$ 大，就保留原来的特征，如果随机的是 $-1$ 比 $0$ 小，$\max$ 取 $0$，不保留这个特征。

因此 $A$ 矩阵与 $R$ 矩阵的转置相乘，得到了 $60000\times5000$ 的矩阵，拼接到 $A (60000\times494)$ 的后面变得到了新的 $A (60000\times 5494)$ 矩阵，剩下的步骤和上述相同。

2.2.1 具体步骤#

加载数据集与数据预处理同上
训练分类模型及模型评估

函数clf_2(X_train, y_train, X_test, y_test):

def clf_2(X_train, y_train, X_test, y_test):
    """
    :param X_train: 训练数据
    :param y_train: 训练数据标签
    :param X_test: 测试数据
    :param y_test: 测试数据标签
    :return: 
    """
    Random_feature = np.random.choice([-1, 1], (5000, 494))
    # 开始训练
    start2 = time.time()
    A = X_train
    A_test = X_test
    A_random = A.dot(Random_feature.T)
    A_random[A_random < 0] = 0
    A_add = np.hstack([X_train, A_random])
    b_train = y_train
    b_test = y_test
    print('进行QR分解中...')
    q_add, r_add = np.linalg.qr(A_add)
    print('已完成QR分解')
    print('对x进行求解中...')
    x_add = np.linalg.pinv(r_add).dot(q_add.T.dot(b_train))
    print('已求得x的解')
    result_add = A_add.dot(x_add)
    print('已完成结果预测')
    tran1or_0(result_add)
    end2 = time.time()
    # 结果展示分析
    show_result(train_labels, result_add, 'Train_dataset_ADD')
    A_random_test = np.dot(A_test, Random_feature.T)  ##10000*50000
    A_random_test[A_random_test < 0] = 0
    A_add_test = np.hstack([A_test, A_random_test])  ## 10000*5494
    result_add_test = A_add_test.dot(x_add)
    tran1or_0(result_add_test)
    show_result(test_labels, result_add_test, 'Test_dataset_ADD')

    print("运行时间:%.2f秒\n" % (end2 - start2))

2.2.2 结果展示#

3. 结果及分析#

	训练集分类错误率	测试集分类错误率	运行时间 (s)
$494$ 个特征数据集	$1.553%$	$1.580%$	$1.38$
$5494$ 个特征数据集	$0.218%$	$0.220%$	$145.23$

准确率：增加 $5000$ 个随机矩阵后，增加了模型的复杂度，模型可以学习到更泛化的特征，提高模型的性能和泛化能力，从而在数据集上获得更好的性能。从实验结果可以看出，模型在训练集和测试集上准确率都得到了很大的提升。
运行时间分析：增加随机特征后，将 $A$ 矩阵扩充，错误率显著降低。但是运行时间很慢，原因是生成的随机矩阵过大，约 $60000\times 5000=3$ 亿个数据。

4. 完整代码#

import numpy as np
import pandas as pd
import struct
import prettytable as pt
import time

# 训练集文件
train_images_idx3_ubyte_file = r'C:\Users\路正桥\Documents\train-images.idx3-ubyte'
# 训练集标签文件
train_labels_idx1_ubyte_file = r'C:\Users\路正桥\Documents\train-labels.idx1-ubyte'

# 测试集文件
test_images_idx3_ubyte_file = r'C:\Users\路正桥\Documents\t10k-images.idx3-ubyte'
# 测试集标签文件
test_labels_idx1_ubyte_file = r'C:\Users\路正桥\Documents\t10k-labels.idx1-ubyte'


def decode_idx3_ubyte(idx3_ubyte_file):
    """
    解析idx3文件的通用函数
    :param idx3_ubyte_file: idx3文件路径
    :return: 数据集
    """
    # 读取二进制数据
    bin_data = open(idx3_ubyte_file, 'rb').read()

    # 解析文件头信息，依次为魔数、图片数量、每张图片高、每张图片宽
    offset = 0
    fmt_header = '>iiii'  # 因为数据结构中前4行的数据类型都是32位整型，所以采用i格式，但我们需要读取前4行数据，所以需要4个i。我们后面会看到标签集中，只使用2个ii。
    magic_number, num_images, num_rows, num_cols = struct.unpack_from(fmt_header, bin_data, offset)
    print('图片数量: %d张, 图片大小: %d*%d\n' % (num_images, num_rows, num_cols))

    # 解析数据集
    image_size = num_rows * num_cols
    offset += struct.calcsize(fmt_header)  # 获得数据在缓存中的指针位置，从前面介绍的数据结构可以看出，读取了前4行之后，指针位置（即偏移位置offset）指向0016。
    fmt_image = '>' + str(
        image_size) + 'B'  # 图像数据像素值的类型为unsigned char型，对应的format格式为B。这里还有加上图像大小784，是为了读取784个B格式数据，如果没有则只会读取一个值（即一副图像中的一个像素值）
    images = np.empty((num_images, num_rows, num_cols))
    for i in range(num_images):
        images[i] = np.array(struct.unpack_from(fmt_image, bin_data, offset)).reshape((num_rows, num_cols))
        # print(images[i])
        offset += struct.calcsize(fmt_image)

    return images


def decode_idx1_ubyte(idx1_ubyte_file):
    """
    解析idx1文件的通用函数
    :param idx1_ubyte_file: idx1文件路径
    :return: 数据集
    """
    # 读取二进制数据
    bin_data = open(idx1_ubyte_file, 'rb').read()

    # 解析文件头信息，依次为魔数和标签数
    offset = 0
    fmt_header = '>ii'
    magic_number, num_images = struct.unpack_from(fmt_header, bin_data, offset)
    print('标签数量: %d张\n' % (num_images))

    # 解析数据集
    offset += struct.calcsize(fmt_header)
    fmt_image = '>B'
    labels = np.empty(num_images)
    for i in range(num_images):
        labels[i] = struct.unpack_from(fmt_image, bin_data, offset)[0]
        offset += struct.calcsize(fmt_image)
    return labels


def load_train_images(idx_ubyte_file=train_images_idx3_ubyte_file):
    """
    :param idx_ubyte_file: idx文件路径
    :return: n*row*col维np.array对象，n为图片数量
    """
    print("训练图片数据已解析：")
    return decode_idx3_ubyte(idx_ubyte_file)


def load_train_labels(idx_ubyte_file=train_labels_idx1_ubyte_file):
    """
    :param idx_ubyte_file: idx文件路径
    :return: n*1维np.array对象，n为图片数量
    """
    print("训练标签数据已解析：")
    return decode_idx1_ubyte(idx_ubyte_file)


def load_test_images(idx_ubyte_file=test_images_idx3_ubyte_file):
    """
    :param idx_ubyte_file: idx文件路径
    :return: n*row*col维np.array对象，n为图片数量
    """
    print("测试图片数据已解析：")
    return decode_idx3_ubyte(idx_ubyte_file)


def load_test_labels(idx_ubyte_file=test_labels_idx1_ubyte_file):
    """
    :param idx_ubyte_file: idx文件路径
    :return: n*1维np.array对象，n为图片数量
    """
    print("测试标签数据已解析：")
    return decode_idx1_ubyte(idx_ubyte_file)


def pretreat(train_labels, test_labels, train_images, test_images):
    """
    数据预处理函数：特征选择
    :param 原始数据集
    :return: 处理好的新数据集
    """
    # 将三维数组转换为二维数组
    train_labels = train_labels.reshape(60000)
    test_labels = test_labels.reshape(10000)
    train_images = train_images.reshape(60000, 784)
    test_images = test_images.reshape(10000, 784)

    #  把标签转换为1和-1，数字为0则为1，否则为-1
    # np.where(condition,x,y):如果condition为真，返回x,否则返回y
    train_labels = np.where(train_labels == 0, 1, -1)
    test_labels = np.where(test_labels == 0, 1, -1)

    # 找出至少在600张图片中不为0的像素位置
    j = 0
    index = []
    train_image_feature = np.zeros([60000, 493])  # 用于存储训练图像的特征
    for i in range(784):
        non_zero = np.linalg.norm(train_images[:, i], ord=0)  # 计算每一列中非零元素的个数
        if non_zero >= 600:
            train_image_feature[:, j] = train_images[:, i]  # 将具有足够非零元素的列复制到train_image_feature的第j列
            j = j + 1
            index.append(i)  # 记录哪些列被选做特征
    test_image_feature = np.zeros([10000, 493])  # 创建一个零数组，用于存储测试图像的特征
    test_image_feature = test_images[:, index]  # 使用 index 从 test_images 中选择相同的列作为特征，并赋值给 test_image_feature

    # 在训练集和测试集特征向量开头添加常量1作为偏置
    train_images = np.insert(train_image_feature, 0, 1, axis=1)
    test_images = np.insert(test_image_feature, 0, 1, axis=1)

    return train_labels, test_labels, train_images, test_images


def clf_1(X_train, y_train, X_test, y_test):
    """
    :param X_train: 训练数据
    :param y_train: 训练数据标签
    :param X_test: 测试数据
    :param y_test: 测试数据标签
    :return: 
    """
    # 开始训练
    start1 = time.time()
    A = X_train
    A_test = X_test
    b_train = y_train
    b_test = y_test
    print('进行QR分解中...')
    q, r = np.linalg.qr(A)
    print('已完成QR分解')
    print('\n')
    print('对x进行求解中...')
    x = np.linalg.pinv(r).dot(q.T.dot(b_train))
    print('已求得x的解')
    result = A.dot(x)
    print('已完成结果预测')
    tran1or_0(result)
    end1 = time.time()
    # 结果展示分析
    show_result(train_labels, result, 'Train_dataset')
    result_test = A_test.dot(x)
    tran1or_0(result_test)
    show_result(test_labels, result_test, 'Test_dataset')

    print("运行时间:%.2f秒\n" % (end1 - start1))


def clf_2(X_train, y_train, X_test, y_test):
    """
    :param X_train: 训练数据
    :param y_train: 训练数据标签
    :param X_test: 测试数据
    :param y_test: 测试数据标签
    :return: 
    """
    Random_feature = np.random.choice([-1, 1], (5000, 494))
    # 开始训练
    start2 = time.time()
    A = X_train
    A_test = X_test
    A_random = A.dot(Random_feature.T)
    A_random[A_random < 0] = 0
    A_add = np.hstack([X_train, A_random])
    b_train = y_train
    b_test = y_test
    print('进行QR分解中...')
    q_add, r_add = np.linalg.qr(A_add)
    print('已完成QR分解')
    print('对x进行求解中...')
    x_add = np.linalg.pinv(r_add).dot(q_add.T.dot(b_train))
    print('已求得x的解')
    result_add = A_add.dot(x_add)
    print('已完成结果预测')
    tran1or_0(result_add)
    end2 = time.time()
    # 结果展示分析
    show_result(train_labels, result_add, 'Train_dataset_ADD')
    A_random_test = np.dot(A_test, Random_feature.T)  ##10000*50000
    A_random_test[A_random_test < 0] = 0
    A_add_test = np.hstack([A_test, A_random_test])  ## 10000*5494
    result_add_test = A_add_test.dot(x_add)
    tran1or_0(result_add_test)
    show_result(test_labels, result_add_test, 'Test_dataset_ADD')

    print("运行时间:%.2f秒\n" % (end2 - start2))


def tran1or_0(result):
    """
    将大于0的数转换为1，小于0的转换为0
    :param
    :return:
    """
    result[result > 0] = 1
    result[result < 0] = -1
    result.reshape(len(result), 1)
    return result


def show_result(labels, result, dataset):
    """
    计算错误率及结果展示
    :param labels: 真实数据标签
    :param result: 预测得到的数据标签
    :param dataset: 数据名称（train or test）
    :return: 
    """
    TP = 0  # 正类预测为正类
    FN = 0  # 正类预测为负类
    FP = 0  # 负类预测为正类
    TN = 0  # 负类预测为负类
    for i in range(len(labels)):
        if labels[i] == 1 and result[i] == 1:
            TP = TP + 1
        elif labels[i] == 1 and result[i] == -1:
            FN = FN + 1
        elif labels[i] == -1 and result[i] == 1:
            FP = FP + 1
        elif labels[i] == -1 and result[i] == -1:
            TN = TN + 1

    data = {
        dataset: ['y=+1', 'y=-1', 'All'],
        'Predicted y=+1': [TP, FN, TP + FN],
        'Prediected y=-1': [FP, TN, FP + TN],
        'Total': [TP + FP, FN + TN, TP + FP + FN + TN]
    }
    df = pd.DataFrame(data)
    # 设置表格样式
    styled_df = df.style \
        .set_properties(**{'text-align': 'center', 'border': '1px solid black'}) \
        .set_table_styles([{'selector': 'th', 'props': [('border', '1px solid black')]}]) \
        .set_caption(dataset + '预测结果：')
    # 输出表格
    display(styled_df)

    error = (FN + FP) / (TP + FP + FN + TN) * 100
    print('错误率为', "%.3f" % error, '%')
    print('\n')


if __name__ == '__main__':
    # 加载原始数据
    train_images = load_train_images()
    train_labels = load_train_labels()
    test_images = load_test_images()
    test_labels = load_test_labels()
    # 得到数组形式的数据，此时数组为三维的
    [train_labels, test_labels, train_images, test_images] = pretreat(train_labels, test_labels, train_images,
                                                                      test_images)
    # 训练模型得出结果
    clf_1(train_images, train_labels, test_images, test_labels)
    # 增加5000个随机特征后重新训练模型，评估模型
    print("增加5000个随机特征后的结果：\n")
    clf_2(train_images, train_labels, test_images, test_labels)