MNISTをCLAFIC法で識別する - 機械学習に詳しくなりたいブログ

今回はCLAFIC法と呼ばれる、KL展開を用いた手法でMNISTを識別してみたいと思います。これまでに、KL展開は分散最大基準と平均二乗誤差最小基準の2種類を見てきました。（参考：KL展開平均二乗誤差最小基準、KL展開分散最大基準）目的が分類の場合は平均二乗誤差最小基準が用いられるようです。

分散最大基準と平均二乗誤差基準のKL展開

それぞれの手法については、KL展開平均二乗誤差最小基準、KL展開分散最大基準に書いてあります。MNISTの画像をKL展開し、最大固有値に対応する固有ベクトルを図示してみました。（次元を１つに削減するとしたときの基底に相当します）

分散最大基準

平均二乗誤差最小基準

それぞれの求め方で、輪郭は似通ってますがけっこう違いますね。この違いを理屈で説明したいんですが、、、わかるような、わからないような。

CLAFIC法

さて、射影後のデータを元の次元で見るで、データ $\mathbf x$ の次元を行列 $\mathbf A$ によって削減した後のベクトル $\mathbf{A}^{T} \mathbf{x}$ を、元の次元で表現すると $\mathbf{A}\mathbf{A}^{T} \mathbf{x}$ であることを確認しました。*1 　このことは、 $\mathbf x$ の次元 $D^{\prime}$ に対する直交射影行列が $\mathbf{A}\mathbf{A}^{T}$ であると捉えることができ、これを

$\mathbf{P}=\mathbf{A}\mathbf{A}^{T} \tag{1}$

とします。また、射影行列の性質として以下の等式が成り立ちます。

$\begin{eqnarray*} \mathbf{P}\mathbf{P}&=&\mathbf{P} \tag{2} \\ \mathbf{P}^{T}&=&\mathbf{P} \tag{3} \end{eqnarray*}$

さて、データをある部分空間へ正射影することを考えます。このとき、その正射影ベクトルのノルムが大きければ、その部分空間が持つ特徴に似ていると捉えられます。よって識別対象全ての部分空間へ射影し、そのノルムが最大であるものが識別結果と考えることができます。このように分類する手法をCLAFIC法を呼びます。

$\mathbf P$ は直交射影行列ですから、データ $\mathbf x$ の $D$ 次元空間→ $D^{\prime}$ 次元空間の正射影は $\mathbf P \mathbf x$ となり、このノルムの二乗は、式(3)を使えば

$\begin{eqnarray*} \| \mathbf{P} \mathbf{x} \|^{2} &=& (\mathbf{P}\mathbf{x})^{T}\mathbf{P}\mathbf{x}\tag{4} \\ &=&\mathbf{x}^{T} \mathbf{P}\mathbf{P}\mathbf{x}\tag{5} \\ &=& \mathbf{x}^{T} \mathbf{P}\mathbf{x}\tag{6} \end{eqnarray*}$

と書けます。（参考：ベクトルのノルムの式(5)、転置行列の定理）この式(6)をMNISTの部分空間全て対して計算し、最大値を識別結果とします。

以上のことは下図ようなイメージです。紫と緑はそれぞれ違うクラスの部分空間を表します。いま、点 $\mathbf x$ をそれぞれの部分空間へ正射影したとき、そのノルム（矢印の長さ）が長いほうへ属するクラスだと識別されます。

CLAFIC法によるMNISTの識別実験結果

CLAFIC法で実際にMNISTを識別してみました。今回は平均二乗誤差最小基準、分散最大基準、それぞれで計算した部分空間を用います。また、部分空間の次元は適当にふって傾向を見てみました。結果は以下のグラフです。横軸は部分空間の次元数です。

平均二乗誤差最小基準のほうが精度が高く、95.7%くらいの正解率が出ています。*2 部分空間の次元を増やすにつれて、どちらの基準でも精度は悪くなっていきます。次元を増やしていくと、どの部分空間でも元の画像が表現できてしまいますから、各部分空間で差が表れにくいのでしょう。（参考：KL展開で得られた基底で画像を表現する） CLAFIC法は50年近く前に提案された手法だそうですが、この時代の手法でこれだけの精度が出るのは驚きました。でもMNISTの識別精度をいろいろネットで見て回ると、畳込みニューラルネットワークで実験した結果はよく目に付き、99.7%とかそれくらいの精度が出るみたいです。この数%の差が大きいんでしょうね。ただ問題の状況と求める精度によっては、ニューラルネットワークなんか使わなくても50年前の手法でも十分、ということもあるかもしれません。

今回のコードです。 MNISTのデータは、http://yann.lecun.com/exdb/mnist/に、訓練データ用、テスト用と分けられたものが公開されていましたのでこちらを使いました。また、式(6)は、以下のように変形して計算したほうが効率が良いそうなのでこちらを使っています。

$\begin{eqnarray*} \displaystyle \mathbf{x}^{T} \mathbf{P}\mathbf{x} &=&\mathbf{x}^{T} \left( \sum_{j=1}^{d} (\mathbf{u}_j\mathbf{u}_j^{T}) \right)\mathbf{x}\tag{7} \\ \displaystyle &=& \sum_{j=1}^{d} (\mathbf{x}^{T}\mathbf{u}_j)^2 \tag{8} \\ \end{eqnarray*}$

なお $\mathbf{u}$ などの文字の意味はKL展開分散最大基準　解の導出で定義してあるものと同じです。

# CLAFIC法でMNISTを識別
import numpy as np

# テストデータを読み込み
data = np.load("./mnist/test_data.npy")
label = np.load("./mnist/test_label.npy")

# データ数
N = data.shape[0]

# ラベル数
C = 10

# 次元数
D = data.shape[1]

# 各クラスの固有ベクトルを読み込み（固有値降順にソート済み）
A = np.empty([C, D, D])
for i in range(C):
    A[i, :, :] = np.load("./mnist/mnist_eig_{0}.npy".format(i))

# CLAFIC法（式(8)の計算）
out = np.empty(C)
for d in range(0, 200, 20):
    count_ok = 0
    for i in range(N):
        for j in range(C):
            norm = 0
            for k in range(d):
                norm += np.dot(data[i, :], A[j, :, k])**2
            out[j] = norm
        if np.argmax(out) == label[i]:
            count_ok += 1

    # 正解率
    print(count_ok/N * 100, d)

参考　MNISTデータをnumpyで扱える形式に変換

# http://yann.lecun.com/exdb/mnist/ でダウンロードしたMNISTをndarrayに変換して保存
import numpy as np
import gzip

with gzip.open("./mnist/t10k-images-idx3-ubyte.gz", "rb") as f:
    # 最初の16バイトは画像数などの情報
    x = np.frombuffer(f.read(), dtype=np.uint8, offset=16)
    x = x.reshape(-1,784)

np.save("./mnist/test_data", x)

with gzip.open("./mnist/train-images-idx3-ubyte.gz", "rb") as f:
    x = np.frombuffer(f.read(), dtype=np.uint8, offset=16)
    x = x.reshape(-1,784)

np.save("./mnist/train_data", x)

with gzip.open("./mnist/t10k-labels-idx1-ubyte.gz", "rb") as f:
    # 最初の8バイトはラベル数などの情報
    x = np.frombuffer(f.read(), dtype=np.uint8, offset=8)

np.save("./mnist/test_label", x)

with gzip.open("./mnist/train-labels-idx1-ubyte.gz", "rb") as f:
    x = np.frombuffer(f.read(), dtype=np.uint8, offset=8)

np.save("./mnist/train_label", x)

test_data.shapeは(10000, 784)であり、各行に文字データのグレースケール値が入ります。

test_label.shapeは(10000,)であり、それぞれ正解ラベルの0～9の値が入ります。

*1:変換行列 $\mathbf A$ の定義はこちら：KL展開分散最大基準　解の導出

*2:CLAFIC法でMNISTを識別した結果があまりネットになく、結果が合っているのかちょっと不安