手書き数字データベース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行を読み取った後、ポインタ位置（オフセット位置）は0016を指します。
    fmt_image = '>' + str(
        image_size) + 'B'  # 画像データのピクセル値の型はunsigned char型で、対応するフォーマットはBです。ここで784のB形式データを読み取るために画像サイズを加えています。そうでなければ、画像内の1つのピクセル値のみを読み取ります。
    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):
    """
    大きい数を1に、小さい数を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],
        'Predicted 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)