概要
MNISTの数字認識をPythonでフレームワーク等を使わずに行います。何番煎じかも分からないし、そもそもディープラーニングのスクラッチなら「ゼロから作るDeep Learning」[1]という名著があります。私も昔この本を少し読んでいましたし、この記事のコードも要所要所で影響を受けていると思います。
けれど、こういった文献のソースコードは優秀な方が(おそらく)何度も練って書いたもので、洗練されすぎているが故に頭に入ってこないです。オライリーの本の紙面上ではシンタックスハイライトが施されていないのも大きいと思います。※個人的感想です
そういった中で、先日統計検定®の準一級とデータサイエンスエキスパートを取得しまして、久々にディープラーニングを扱う機会があったので、復習を兼ねて、自分が自然な流れで書くなら...というコンセプトで、実行速度などの実用面を無視してMNISTに挑戦した記録がこの記事です。
おことわり
全般
・理論面の解説はほぼしません。私が引っかかっていた箇所について軽く書くことはあるかもしれません。私ははじめディープラーニングを[1]の文献で学んだので、考え方はこれに近いはずです。
・先述の経緯で作成したコードなので、多重ループなどが随所に存在します。効率化の余地は大いにあります。また、書いた後にコードを練る、いわゆるリファクタリングのようなことはしていません。逆に言えば、私はこの記事に載せた順でコードを書いていますし、それが私個人にとっては自然な流れだと考えています。
数学
行列なしとは書きましたが、流石にすべてfor文は逆に面倒なので、numpyの2次元配列に対しての、いわゆるアダマール積のみ使いました。以下のようなものです。同じ位置の要素同士の積を取ります。
\begin{bmatrix}
1 & 2 \\
3 & 4
\end{bmatrix}
\odot
\begin{bmatrix}
5 & 6 \\
7 & 8
\end{bmatrix}
=
\begin{bmatrix}
1 \times 5 & 2 \times 6 \\
3 \times 7 & 4 \times 8
\end{bmatrix}
=
\begin{bmatrix}
5 & 12 \\
21 & 32
\end{bmatrix}
演算子の*を使用します。
# 2つのnumpy2次元配列に対してアダマール積を考える
a = np.array([[1,2],[3,4]])
b = np.array([[5,6],[7,8]])
print(a*b)
[[ 5 12]
[21 32]]
依存するライブラリ
mathとnumpyを実装に、matplotlibを結果などの表示用に使います。
import math
import numpy as np
import matplotlib.pyplot as plt
ちなみに、今回のコーディングにおいては、学習目的であるため、本質的な箇所に生成AIは使っていません。(matplotlibやnumpyの書き方などの書式面においては一部使用しました)
下準備
まずは肝心のMNISTデータを用意します。こちらの記事を参考にしています。
mnist_data = np.load("mnist.npz") # パスは変えてくださいね
x_train = mnist_data['x_train']
y_train = mnist_data['y_train']
x_test = mnist_data['x_test']
y_test = mnist_data['y_test']
データサイエンスの第一歩はデータの確認からです。まずは各データの格納形式を確認しましょう。
print("x_trainのshape:" + str(x_train.shape))
print("y_trainのshape:" + str(y_train.shape))
print("x_testのshape:" + str(x_test.shape))
print("y_testのshape:" + str(y_test.shape))
x_trainのshape:(60000, 28, 28)
y_trainのshape:(60000,)
x_testのshape:(10000, 28, 28)
y_testのshape:(10000,)
x_trainとy_trainを詳しく見てみましょう。
x_trainの確認
print(x_train)
[[[0 0 0 ... 0 0 0]
[0 0 0 ... 0 0 0]
[0 0 0 ... 0 0 0]
...
[0 0 0 ... 0 0 0]
[0 0 0 ... 0 0 0]
[0 0 0 ... 0 0 0]]
<以下略>
28×28pxの画像データが60000枚入っているようですね。1枚当たりの中身も見てみましょう。
print(x_train[0])
[[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 3 18 18 18 126 136 175 26 166 255 247 127 0 0 0 0]
[ 0 0 0 0 0 0 0 0 30 36 94 154 170 253 253 253 253 253 225 172 253 242 195 64 0 0 0 0]
[ 0 0 0 0 0 0 0 49 238 253 253 253 253 253 253 253 253 251 93 82 82 56 39 0 0 0 0 0]
[ 0 0 0 0 0 0 0 18 219 253 253 253 253 253 198 182 247 241 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 80 156 107 253 253 205 11 0 43 154 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 14 1 154 253 90 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 139 253 190 2 0 0 0 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 11 190 253 70 0 0 0 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 35 241 225 160 108 1 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 81 240 253 253 119 25 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 45 186 253 253 150 27 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 16 93 252 253 187 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 249 253 249 64 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 46 130 183 253 253 207 2 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 39 148 229 253 253 253 250 182 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 24 114 221 253 253 253 253 201 78 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 23 66 213 253 253 253 253 198 81 2 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 18 171 219 253 253 253 253 195 80 9 0 0 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 55 172 226 253 253 253 253 244 133 11 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 136 253 253 253 212 135 132 16 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]
[ 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0]]
0-255の値のグレースケール画像のようです。matplotlibで表示してみると以下のような出力が得られます。
plt.imshow(x_train[0], cmap='gray')
plt.show()
y_trainの確認
y_trainはどうでしょうか。
print(y_train)
[5 0 4 ... 5 6 8]
0-9の正解ラベルのようですね。分布を確認してみましょう。
plt.hist(y_train)
全体が60000枚であり、計10種類のラベルがあることを考えると、おおよそ一様であると考えて問題なさそうです。
確認を踏まえたデータ加工
ニューラルネットワークの出力形式と合わせるべく、y_train、つまり正解ラベルはone-hot化してしまいましょう。
まずは1つの0-9の数値をone-hotラベルに変換する関数を定義します。
def make_one_hot_single(num):
res = np.zeros(10)
res[num] = 1
return res
# テスト
print(make_one_hot_single(8))
[0. 0. 0. 0. 0. 0. 0. 0. 1. 0.]
※値の型はfloatですが、ニューラルネットワークの出力も一般に小数(そのラベルが正解である確率)なので問題ありません。
これを用いて、y_trainのような配列を変換できるようにします。
def make_one_hot_multi(y_data):
res = np.array([])
for number in y_data:
res = np.append(res, make_one_hot_single(number))
res = np.reshape(res, (len(y_data), 10)) # こうしないと出力が1次元配列になってしまう
return res
# テスト
y_sample = np.array([1,2,3,5])
y_sample_one_hot = make_one_hot_multi(y_sample)
print(y_sample_one_hot)
print(y_sample_one_hot.shape)
[[0. 1. 0. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 1. 0. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 1. 0. 0. 0. 0. 0. 0.]
[0. 0. 0. 0. 0. 1. 0. 0. 0. 0.]]
(4, 10)
for文を使っています。お察しの通りとても遅いです。これでy_trainを変換しておきましょう。以下、one-hot化されたy_trainをy_train_one_hotとします。
y_train_one_hot = make_one_hot_multi(y_train)
print("y_train_one_hotのshape:" + str(y_train_one_hot.shape))
y_train_one_hotのshape:(60000, 10)
また、先ほどx_trainの値の範囲が0-255であると分かりましたが、これをそのままニューラルネットワークに入れると値が暴走してしまいます。そのため、値が0-1の範囲に収まるよう、255で割ります。
x_train = x_train / 255.0 # intへのキャスト回避のために、255.0と小数にする
ニューラルネットワークの設計
これを踏まえて、ニューラルネットワークの設計を行います。手作りする時点で各種フレームワーク・ライブラリより低速であるのはほぼ間違いないので、最低限の構成とします。...とはいえここを考えるのはこの記事の趣旨ではないので、Gemini君に聞いてみます。
結果、以下のようなCNNネットワークを作成することにしました。
※レイヤ→出力の形状
Input→(28,28,1)
Conv2D(Filter32,3×3,stride=1)→(26,26,32)
MaxPooling2D(2×2,stride=2)→(13,13,32)
Flatten→(5408)
Dense(ReLU)→(128)
Dense(Softmax)→(10)
損失関数は交差エントロピー誤差とし、出力層のsoftmaxと組み合わせることで誤差逆伝播法による実装を容易にします。
実装開始
実装の方針
どのように実装するか方針を立てます。
まず、各層をクラスとして設計し、インスタンス化して使うようにします。
この際、誤差逆伝播を行えるよう、各層は必ず順伝播用のforward()メソッドと逆伝播用のbackward()メソッドを持つようにします。
さらに、逆伝播に必要な情報(ex. ReLU関数のマスク情報)を保持できるよう、活性化関数も単なる関数ではなく、それをラップするクラスとして定義し、各層のインスタンスに活性化関数クラスのインスタンスを保持させます。活性化関数クラスもforward()メソッドとbackward()メソッドを持ちます。
学習率
この実装を通して使う学習率を定義しておきます。
learningRate = 0.03
損失関数
交差エントロピー誤差を計算する関数を定義します。
delta = 1e-7 # 真数条件でエラーを出さないために微小な値を加える
def crossEntropy(output, label):
return -np.sum(label*np.log(output+delta))
活性化関数
今回使う活性化関数は全結合層1層目のReLUと2層目にして出力層のsoftmaxのみです。それぞれ定義していきます。
ReLU
class HiReLU:
def __init__(self):
self.note = np.array([]) # マスク情報を保持する配列です
def __ReLU(self, x):
return max(0,x) # 値1つに対するReLU関数です
def __vectrizedReLU(self, x):
return np.vectorize(self.__ReLU)(x) # それをベクトルの入力に対応できるようにします
def forward(self,x):
self.note = (x <= 0) # xが0以下のインデックスをマスク情報として記録したうえで
return self.__vectrizedReLU(x) # ReLUを通したベクトルを返します
def backward(self,y):
y[self.note] = 0 # ReLU関数の微分値はxが0以下の場合0, そうでない場合下層から返ってきた勾配情報そのままです
return y
softmax
class HiSoftmax:
def __init__(self):
pass
def forward(self,x):
c = np.max(x) # 値の爆発を防ぐため、xの最大値が0となるようにします(結果は変化しません)
exp_x = np.exp(x - c)
sum_exp_x = np.sum(exp_x)
y = exp_x / sum_exp_x
return y
def backward(self,y):
return y # CrossEntropyとの併用前提
各層の定義
今回使う層は、平滑化のFlatten、畳み込みのConv2D、最大値プーリングのMaxPooling2D、全結合層のDenseです。
Flatten
配列をreshapeします。
class Flatten:
def __init__(self,k,x,y):
self.k = k # フィルタ数(Kernelのk)
self.x = x # 画像の縦幅
self.y = y # 画像の横幅
def forward(self, input):
return input.reshape([self.k*self.x*self.y]) # 3次元を1次元にreshape
def backward(self, gradient):
return gradient.reshape([self.k, self.x, self.y]) # 1次元を3次元にreshape
Conv2D
簡単のため、ストライドは1であることを前提とします。
class Conv2D:
def __init__(self, verticalKernelSize, horizontalKernelSize, filters):
self.verticalKernelSize = verticalKernelSize # フィルタの縦幅
self.horizontalKernelSize = horizontalKernelSize # フィルタの横幅
self.filters = np.random.normal(0, 0.5, (filters, horizontalKernelSize, verticalKernelSize)) # 引数のフィルタ数をもとに、フィルタの値を保持する配列を作成
self.data_recieved = np.array([]) # 逆伝播時のために、順伝播時の入力を記録する配列
def __lap(self, inputFragment, filter):
return np.sum(inputFragment * filter) # アダマール積を計算し、その合計を出す...つまり、画像の一部にフィルターを重ねた結果を返します
def forward(self, input):
self.data_recieved = input # 逆伝播時のため入力を記録する
result = np.zeros((self.filters.shape[0], input.shape[0] - self.verticalKernelSize + 1, input.shape[1] - self.horizontalKernelSize + 1)) # 次の層への出力を保持する配列
for i in range(input.shape[0] - self.verticalKernelSize + 1):
for j in range(input.shape[1] - self.horizontalKernelSize + 1): # 各重ね合わせについて
inputFragment = input[i:i+self.verticalKernelSize, j:j+self.horizontalKernelSize] # 重ね合わせる画像の断片を切り出し
for k in range(self.filters.shape[0]): # 各フィルタを重ね合わせ
result[k][i][j] = self.__lap(inputFragment, self.filters[k]) # その結果を記録する
return result
def backward(self, gradient):
# 今回は最初にのみConvを使うので、上に流す新しいgradientは計算する必要なしです。filterの値の更新のみ行います。
for k in range(self.filters.shape[0]): # フィルタごとに
errors = np.zeros([self.horizontalKernelSize, self.verticalKernelSize]) # 誤差を記録する配列を用意して
for i in range(self.verticalKernelSize):
for j in range(self.horizontalKernelSize):
errors[i][j] = np.sum(self.data_recieved[i:i+gradient.shape[1],j:j+gradient.shape[2]]*gradient[k]) # 各重ね合わせ時における誤差量を集計していきます。(逆畳み込み)
self.filters[k] = self.filters[k] - learningRate * errors # フィルタの値を更新
MaxPooling2D
簡単のため、ストライドはカーネルのサイズと等しいとします。また、カーネルは正方形とします。(Conv2Dでは縦横両方指定できたのにこちらは指定することができないのは、私の体力切れです(笑))。例えば、26×26の画像にこのクラスの処理をkernelSize=2で適用する場合、順伝播の出力は13×13になります。
class MaxPooling2D:
def __init__(self, kernelSize):
self.kernelSize = kernelSize # カーネルサイズ
self.input_shape = np.array([]) # 逆伝播時に元の配列サイズを復元できるよう記録しておく
self.note = np.array([]) # 各重ね合わせにおいて、どの位置の値が最大であったかを記録する
def forward(self, input):
self.input_shape = input.shape # 逆伝播に備え、入力サイズを記録
z = input.shape[0] # フィルタ数(奥行)
x = int(input.shape[1]/self.kernelSize) # 結果の縦幅
y = int(input.shape[2]/self.kernelSize) # 結果の横幅
result = np.zeros((z, x, y)) # プーリング結果を保持する配列
self.note = np.zeros((z, x, y, self.kernelSize, self.kernelSize)) # 3次元で、カーネルサイズと等しい2次元配列を保持し、最大値があった位置のみ1とし、逆伝播時に用いる各重ね合わせにおける最大値の位置の記録とする
for k in range(z):
for i in range(x):
for j in range(y):
area = input[k, i*self.kernelSize:(i+1)*self.kernelSize, j*self.kernelSize:(j+1)*self.kernelSize]
idxM = np.unravel_index(np.argmax(area), area.shape) # 最大値の位置を調べて
idxM = tuple(map(int, idxM))
result[k][i][j] = area[idxM]
self.note[k][i][j][idxM[0]][idxM[1]] = 1 # その位置に1を記録する
return result
def backward(self, gradient):
# 各重ね合わせにおいて最大値があった位置にのみ購買情報を反映
result = np.zeros(self.input_shape)
for k in range(self.input_shape[0]):
for i in range(gradient.shape[1]):
for j in range(gradient.shape[2]):
result[k, i*self.kernelSize:(i+1)*self.kernelSize, j*self.kernelSize:(j+1)*self.kernelSize] += gradient[k][i][j]*self.note[k][i][j]
return result
Dense
activationには任意の活性化関数クラスのインスタンスを格納します。
エッジの重みはInputに近い方のレイヤーの受け持ちとします。
エッジの重みは、former_nodes×nodesのnumpy2次元配列で管理します。
重みの初期化にはHeの初期化を用います。
class Dense():
def __init__(self, activation, nodes, former_nodes):
self.activation = activation # 活性化関数のインスタンスをセット
std = math.sqrt(2/former_nodes) # Heの初期化用の正規分布の標準偏差
self.weights = np.random.normal(0, std, (former_nodes, nodes)) # weightを初期化
self.data_recieved = np.zeros(former_nodes) # 重み更新時のため、入力値を記録しておく配列
self.biases = np.random.normal(0, std, nodes) # ノード内のバイアス
self.nodes = nodes # この層のノード数
self.former_nodes = former_nodes # 前の層のノード数
def forward(self, input):
# inputは長さformer_nodesのnumpy配列です
self.data_recieved = input # 逆伝播の時のために記録
inputToNodes = np.zeros(self.nodes) # 各ノードへの入力値の計算結果を保持する配列
for i in range(len(input)): # 前層の各ノードについて
weight = self.weights[i] # そのノードから生えてるエッジのウェイトの1次元配列をメモって
inputToNodes += weight*input[i] # 入力とウェイトを掛けて
inputToNodes += self.biases # バイアスを加える
inputToNodes = self.activation.forward(inputToNodes) # 最後に活性化関数に通す
return inputToNodes # 次の層へ
def backward(self, gradient):
# 自身のバイアスの更新と自分の上に生えてるエッジのウェイト更新
gradient = self.activation.backward(gradient) # 勾配を活性化関数に後ろから通して
self.biases = self.biases - learningRate * gradient # バイアスを更新して
new_gradient = np.zeros(self.former_nodes) # 前の層に渡すための新しい勾配の計算結果を保持する配列
for j in range(len(new_gradient)): # 前層の各ノードについて
new_gradient[j] = np.sum(self.weights[j]*gradient) # 勾配を更新して
for i in range(self.former_nodes): # 今度は前層の各ノードについて
for j in range(self.nodes): # 順伝播時の入力データをもとに、そのノードから生えている各エッジのウェイトを更新する
self.weights[i][j] = self.weights[i][j] - learningRate * self.data_recieved[i] * gradient[j]
return new_gradient # 前の層へ
実際に学習させてみる
10枚ずつ計300枚学習させて、1枚学習するたびに損失を計算し、また、10枚学習するたびに直近10枚の正答率を出すことにします。
imgToUse = 300
imgPerSet = 10
conv2D_1 = Conv2D(3,3,32)
maxPooling2D_1 = MaxPooling2D(2)
flatten_1 = Flatten(32,13,13)
dense_1 = Dense(HiReLU(), 128, 5408)
dense_2 = Dense(HiSoftmax(), 10, 128)
loss_list = np.array([])
acc_list = np.array([])
for i in range(int(imgToUse/imgPerSet)): # 1セットの処理
print("set:" + str(i))
correct = 0 # このセットにおける正解数を記録する変数
for j in range(imgPerSet): # 画像1枚の処理
imgIdx = i * imgPerSet + j # 読み込む画像のインデックスを計算
img = x_train[imgIdx] # 画像データ(28*28)
answer = y_train[imgIdx] # 正解(0-9)
label = y_train_one_hot[imgIdx] # one-hot化された正解
# 順伝播
f1 = conv2D_1.forward(img) # 畳み込み(→32*26*26)
f2 = maxPooling2D_1.forward(f1) # 最大値プーリング(→32*13*13)
f3 = flatten_1.forward(f2) # 平滑化(→5408)
f4 = dense_1.forward(f3) # 全結合層1(→128)
f5 = dense_2.forward(f4) # 全結合層2/出力層(→10)
loss = crossEntropy(f5, label) # 損失を計算
loss_list = np.append(loss_list, loss) # 記録
if(np.argmax(f5) == answer): # 正解している場合正解数を1増やす
correct += 1
# 逆伝播
error = f5 - label # softmax & crossentropyの場合の逆伝播後の勾配。これを踏まえsoftmaxクラスのbackwardは値をそのまま返す。
e5 = dense_2.backward(error) # 全結合層2のパラーメータ更新
e4 = dense_1.backward(e5) # 全結合層1のパラーメータ更新
e3 = flatten_1.backward(e4) # 平滑化層の上へ
e2 = maxPooling2D_1.backward(e3) # 最大値プーリング層の上へ
conv2D_1.backward(e2) # 畳み込み層のパラメータ更新
acc_list = np.append(acc_list, correct/imgPerSet) # 直近10枚について正答率を算出、記録
実行後のloss_list,acc_listは次のようになります。
plt.plot(loss_list)
plt.plot(acc_list)
300枚でも正しく学習できていることが分かります。
おまけ1:予測する
順伝播と正誤判定の処理を関数として切り出せば予測が可能です。
def predict(imgIdx):
img = x_test[imgIdx]
label = y_test[imgIdx]
print("answer is " + str(label))
f1 = conv2D_1.forward(img)
f2 = maxPooling2D_1.forward(f1)
f3 = flatten_1.forward(f2)
f4 = dense_1.forward(f3)
f5 = dense_2.forward(f4)
print("prediction is " + str(np.argmax(f5)))
# 使用する
predict(4000)
answer is 9
prediction is 9
おまけ2:畳み込み層の学習結果を見る
畳み込み層のfiltersの中身を見てみましょう。
# フィルタを取得
filters = conv2D_1.filters
# 描画用意
num_filters = filters.shape[0]
cols = 8
rows = math.ceil(num_filters / cols)
fig, axes = plt.subplots(rows, cols, figsize=(12, rows * 1.5))
axes = axes.flatten()
for i in range(num_filters):
# フィルタを表示
axes[i].imshow(filters[i], cmap='gray')
axes[i].axis('off') # 軸の数字を消す
axes[i].set_title(f"Filter {i}", fontsize=8)
plt.tight_layout()
plt.show()
...が、正直なところ3×3ではどれがどの特徴量かはあまり読み取れませんね。




