読者です 読者をやめる 読者になる 読者になる

俵言

しがない社会人が書く、勉強とかゆるふわプロコンのこと

【jupyterで学ぶ】 ゼロから作るDeep Learning - 第4回:3章(その3)ニューラルネットワークを試してみる

はじめに

今回も引き続き3章です。今回は その3:試してみる編 です。

この前作ったニューラルネットワークを試すだけなのでさらっと終わると思います。

その3) 実装したニューラルネットワークを試してみる

今回は前回作ったNNを用いて手書き数字(0~9)の認識を試みます。用いるデータセットMNISTというモノクロ画像のデータセットで、とても有名なやつらしいです。一枚の画像は 28 × 28 のピクセルで表現されていて、各ピクセルは 0~255 の値をとります。これらの数字を認識させて、正解ラベルとどれくらい一致するか試してみます。

なお、データを取得する手段やNNのパラメータ(重みとバイアス)は著者の方が用意してくれているため、本当に何もする必要がありません。

準備

実際に認識を行わせる前にデータの取り込みなどを行います。
ただその前に、現在作業しているディレクトリの構成を出しときます。現在の作業ディレクトリはDLFZという名前にしていて、その下の構成は以下の通りです。

  • DLFZ(←ここがカレントディレクトリ)
    • ch01
    • ch02
      .
      .
    • ch08
    • common
    • dataset
    • my_functions.py

ch01~datasetまでは著者の方が用意してくれているスクリプトなどが入っているディレクトリ、my_functions.pyは僕が今まで作った関数やクラスが入ってるpythonファイルです。今回は著者の方が作った物を使うことが多いので一応構成を出しときました。いきなりfrom dataset import xxって出てきても「えっ?」ってなるかなって思ったので。

ではまずはデータセットのダウンロードです。dataset下のmnist_for_py2.pyにある関数load_mnistを実行します。

## import
from dataset.mnist_for_py2 import load_mnist

## データをロード
(x_train, t_train), (x_test, t_test) = load_mnist(normalize=True, flatten=True)

内部で何やってるか簡単に説明すると、

  1. MNISTのページから画像データ(をあらわす28×28ピクセル)を、 訓練データ60000枚・テストデータ10000枚 ダウンロード
  2. 各画像データを一次元の配列に直し、0~1の値に正規化(255で割っている)
  3. データをpickleファイルとして保存
  4. 訓練データおよびそのラベル、テストデータおよびそのラベルを返す

pickleはオブジェクトをそのまま保存できるモジュールです(公式に書かれてる直列化の意味がよくわからない..)。こいつはかなり便利なので頻繁につかってます。

テストデータの形状を確認すると

# テストデータの形状の確認
print x_test.shape
< (10000L, 784L)

となっています。28×28=784を一次元に直しており、それが10000枚あるわけですね。
せっかくなのでテストデータ画像の0番目を表示させます。pythonの画像処理にはPIL(python Imaging Library)を使うとのこと(シラナカッタ..)。

from PIL import Image
def img_show(img):
    pil_img = Image.fromarray(np.uint8(img))
    return pil_img
## 正解ラベルを表示
print "正解のラベル: {0}".format(t_test[0])
## 一次元配列を元の二次元配列に戻し、正規化を解除(255をかける)
img_2darray = x_test[0].reshape(28,28)*255
## 二次元配列を画像に変換 (np.uint8()は各要素を符号なし8ビット整数に変換)
pil_img = img_show(np.uint8(img_2darray))
## 保存と表示
pil_img.save("data_sample.png")
pil_img.show()

< 正解のラベル: 7

f:id:Tawara:20161105092636p:plain

画像ちっさ!いや28×28だしそりゃそうか..

次に著者の方が用意してくれたNNの重み(重み行列とバイアスベクトル)をロードします。これはpickleファイルとして保存されているので、pickle モジュールを使ってオブジェクトに戻します。ディレクトリch03の中にあるものです。

import pickle as pic
## ファイルを開いてロード
with open("ch03\\sample_weight_for_py2.pkl","rb") as f:
    params = pic.load(f)

中身は辞書(dict)で、重み行列を表す二次元配列(W1,W2,W3)とバイアスベクトルを表す一次元配列(b1,b2,b3)が入ってるのが確認できます。

print "ファイルの中身の型: {0}".format(type(params))
print "辞書の中身"
for k,v in sorted(params.iteritems()):
    print "{0} {1: <12} {2}".format(k, v.shape, type(v))

< ファイルの中身の型: <type 'dict'>
< 辞書の中身
< W1 (784L, 50L)  <type 'numpy.ndarray'>
< W2 (50L, 100L)  <type 'numpy.ndarray'>
< W3 (100L, 10L)  <type 'numpy.ndarray'>
< b1 (50L,)       <type 'numpy.ndarray'>
< b2 (100L,)      <type 'numpy.ndarray'>
< b3 (10L,)       <type 'numpy.ndarray'>

実装したNNを試す

データの取り込みも終わったところで本題です。前回作ったNNにロードした重みを食わせて、これを使って手書き文字認識をします。実装は以下の通り。

## 多層ニューラルネットワークの実装
class MultiLayerNN:
    # 初期化で重みの行列とバイアス項のベクトルを渡す
    def __init__(self, weight, bias, act_func = sigmoid, out_func = indentify_function):
        ## 重み行列のリスト
        self.W = weight
        ## バイアスベクトルのリスト
        self.B = bias
        ## レイヤー(層)の数 ※ただし入力層を除いた数
        self.layer_num = len(self.B)
        ## 活性化関数
        self.act_func = act_func
        ## 活性化関数(出力時)
        self.out_func = out_func

    # NNによって入力xを出力yに変換
    def forward(self, x):
        z = x
        for i in xrange(self.layer_num-1):
            ## 第i層 -> 第(i+1)層
            a = self.W[i].dot(z) + self.B[i]
            z = self.act_func(a)

        ## 第(layer_num-1)層 -> 第layer_num層
        a = self.W[self.layer_num-1].dot(z) + self.B[self.layer_num-1]
        y = self.out_func(a)

        return y

前回作ったといいつつ実装変更しちゃいました。だって前回のやつ三層にしか対応してないんですもん。重み行列、バイアスベクトルをリストで受け取り、活性化関数も引数で受け取れるようにしています。

それでは自作NNを初期化して実行してみます。まずは重みを食わせましょう。

from my_functions import MultiLayerNN, softmax
## 重みの形式を自作関数用に変更
Weight = [params["W{0}".format(i+1,)].T for i in xrange(3)]
Bias = [params["b{0}".format(i+1,)] for i in xrange(3)]

## 初期化
mnn = mnn = MultiLayerNN(weight = Weight, bias = Bias, out_func = softmax)

ここで、重み行列を転置(.T)していることに注意です。これは、前回述べたように僕の実装では入力に対して重みを左からかける( W \cdot z)ためです。本では右からかける( z \cdot W)のが前提になってるので与えられた重みもそれに合わせて作られています。
この時点では、これが悲劇につながるとは思いもしませんでした..

さて、では実際にテストデータを認識させて見ます。ここでは、テストデータを一個ずつ認識させて正解ラベルと一致する数を数えます。

N = len(x_test)
match_cnt = 0
for i in xrange(N):
    ## forwardによって各数字(0-9)である確率が入った一次元配列が得られる(インデックスが数字に対応)
    ## この一次元配列から一番確率の高いインデックスを取得
    t_predicated = np.argmax(mnn.forward(x_test[i]))
    ## 正解ラベルと比較
    if t_predicated == t_test[i]:
        match_cnt += 1
print "Accuracy: {0} %".format(match_cnt * 100. / N)
< Accuracy: 93.52 %

93.52%、なかなか高いですね。この本ではこれから更に工夫をして、最終的に99%を超えさせるまで行くそうです。

バッチ処理

先ほどは、テストデータであるベクトルを一個ずつ認識させていました。でも認識の処理って実装してきて分かった通り実際は行列演算です。

f:id:Tawara:20161105123402p:plain

図のように、一個のベクトル  x に対して行列をかけていって(本当はバイアスを足す処理がありますが図が見にくくなるので省略)、最終的にベクトル  y が出てきています。
でもこれって入力を行列にしても一斉に処理できますよね?例えば100個のベクトルをまとめて行列にしてしまって一気に認識処理をさせることができるわけです。

f:id:Tawara:20161105124459p:plain

で、このようにデータのまとまり(バッチ : batch)を一度に認識処理させることをバッチ処理と言うそうです。多くの数値計算ライブラリでは大きな配列の計算を高度に最適化していて、小さな配列の計算を一個一個やるよりも効率的になるとのこと。なお、バッチ学習とは意味が違うことに注意!

さてこのバッチ処理を実際にやってみましょう。

## 入力を行列で与えて一気に計算する
N = len(x_test)
match_cnt = 0
batch_size = 100
for i in xrange(0,N,batch_size):
    ## t_predicted がベクトルになる
    t_predicted = np.argmax(mnn.forward(x_test[i:i+batch_size].T), axis = 0)
    match_cnt += np.sum(t_predicted == t_test[i:i+batch_size])
print "Accuracy: {0}%".format(match_cnt * 100. / N)

結果自体はさっきと同じになるはずですよね。

< ValueError: operands could not be broadcast together with shapes (50,100) (50,) 

エラー!? ..何故だ、計算は間違っていなかったはず...

と思っていたのですが、実は大きな過ちを犯していたことにここで気づきました。

ブロードキャストの罠

まあ罠って言っても僕が仕様をちゃんと読んでなかっただけですが。
※この節の話は本の通りにやってれば問題ないので読み飛ばしてOKです。やたらと長いので。

numpy.ndarrayを使う上でめっちゃ便利な機能がブロードキャストです。たとえば、配列に対して単一の値を足したりできるのはこの機能のおかげ。

A = np.array([1,2,3])
B = A + 1
print "A:\n",A
print "B:\n",B

< A:
< [1 2 3]
< B:
< [2 3 4]

そして何より重要なのが、形状(要素数や次元)の異なる配列同士で計算できることです。

A = np.array([[1],[3]]) # A.shape : (2,1)
B = np.array([2,5])      # B.shape : (2) 
C = A+B
print "A:\n",A
print "B:\n",B
print "C:\n",C

< A:
< [[1]
<  [3]]
< B:
< [2 5]
< C:
< [[6 8]
<  [8 10]]

ブロードキャストは要素数が少ない方を勝手に補完してくれる非常に便利な機能なのですが、実は制約があります。

Broadcasting — NumPy v1.11 Manual

When operating on two arrays, NumPy compares their shapes element-wise. It starts with the trailing dimensions, and works its way forward. Two dimensions are compatible when

  1. they are equal, or
  2. one of them is 1

形状を後ろから比較していって、要素数が等しい or 片方が1なら計算が可能です。この後ろからというのが重要。例えば以下の二つの配列を考えたとき、

WZ = np.array([[1,4,7],[2,5,8]]) # WZ.shape:(2,3)
b = np.array([1,2]) # b.shape(2,)
print "WZ:\n",WZ
print "b:\n",b

< WZ:
< [[1 4 7]
< [2 5 8]]
< b:
< [1 2]

こう b が勝手に90度回転して縦向きになって計算できそうじゃん!って思うんですよ。でも実際は

print WZ + b

< ValueError: operands could not be broadcast together with shapes (2,3) (2,) 

エラーが出ます。shapeの数値からわかるように、一番後ろが 2 != 3 だから計算できないんですよね。そして全く同じことが僕の実装で起きた訳です。

何でエラーが起きるのかというと、各層における伝達処理において重みを左からかけるようにした( z^{(i+1)} = W \cdot z^{(i)} + b)から。そう、素直に本の通りに  z^{(i+1)} = z^{(i)} \cdot W + b で実装していたらこのエラーは起きなかったんですよね。

左からかけるやり方の場合は、入力  z^{(i)} を縦ベクトルとして計算を行います。バッチ処理をしない場合は、 W \cdot z^{(i)} b も実装上は一次元配列なので、足し算が可能です。
でもバッチ処理を行う場合は入力がベクトル  z^{(i)} じゃなくて行列 Z^{(i)} になるのでうまくいきません。  Z^{(i+1)} = W \cdot Z^{(i)} + b を考えると、実装上は行列  W \cdot Z^{(i)} は二次元配列、 b は一次元配列です。ここでブロードキャストが働いてうまくいくって思ってたんですが、実際はさっきに示した例と同じ現象が起きてアウトとなるわけ。

本のやり方でやる場合は

WZ = np.array([[1,2],[4,5],[7,8]]) # WZ.shape:(3,2)
b = np.array([1,2]) # b.shape(2,)
print "WZ:\n",WZ
print "b:\n",b

print "\nWZ + b:\n",WZ + b # WZ + b.shape:(3,2)

< WZ:
< [[1 2]
<  [4 5]
<  [7 8]]
< b:
< [1 2]
< 
< WZ + b:
< [[ 2  4]
<  [ 5  7]
<  [ 8 10]]

このように一番後ろが揃うので問題なく計算できます。本の書き方が右からかけるやり方になってるのって、単に図の方向と合わせるためとか好みのためかと思ってたんですが、入力ベクトルを横ベクトルにすることでブロードキャストをスムーズに済ませるためだったのかもしれません。

なお、対処する方法としては以下が考えられます。
1. 本のやり方に屈する(素直に右から重みをかける)
2. バイアスの形を明確に与える(縦にする)

もちろん2を選びます。せっかくですし。さっきの例なら、bを二次元配列で表現した縦ベクトルにしてあげます。

WZ = np.array([[1,4,7],[2,5,8]]) # WZ.shape:(2,3)
b = np.array([[1],[2]]) # b.shape(2,1)
print "WZ:\n",WZ
print "b:\n",b

< WZ:
< [[1 4 7]
< [2 5 8]]
< b:
< [[1]
<  [2]]
< 
< WZ + b:
< [[ 2  5  8]
<  [ 4  7 10]]

ちゃんと計算できました!

本題に戻る

さてかなり寄り道しましたが、バッチ処理を再びやって見ましょう。さきほど述べたようにバイアスパラメータを明示的に縦ベクトルの形にすることでエラーが取れるはずです。

## もう一度重み、バイアス、出力関数をNNにわたして初期化し直す
Weight = [params["W{0}".format(i+1,)].T for i in xrange(3)]
### バイアスを二次元配列で渡すことで対処する
Bias = [np.array([params["b{0}".format(i+1,)]]).T for i in xrange(3)]
### NNを初期化し直す
mnn = MultiLayerNN(weight = Weight, bias = Bias, out_func = softmax)

てわけでもっかい実行.

N = len(x_test)
match_cnt = 0
batch_size = 100
for i in xrange(0,N,batch_size):
    t_predicted = np.argmax(mnn.forward(x_test[i:i+batch_size].T), axis = 0)
    match_cnt += np.sum(t_predicted == t_test[i:i+batch_size])
print "Accuracy: {0} %".format(match_cnt * 100. / N)

< Accuracy: 93.52 %

一個ずつやった場合と同じ結果だ!やったッ!第3章完!

バッチ処理、いかほどのものか

「完!」と言ったものの、バッチ処理でどれくらい効率化するのか気になりません?なので実行時間を比較してみます。(メモリ消費量は今回は無視で。)

さっきのバッチ処理のコードをほぼそのままbatch_sizeを引数とする関数batch_testとして定義し、実行時間を計測します。batch_sizeは、{1,10,25,50,100,250,500,1000,2500,5000,10000} の11種類です。
計測にはtimeitというモジュールを使用します。ちょっと引数の渡し方がややこしいのですが、複数のケースを試したい場合はこれがよさそう。計測については100回実行して最も速かった場合・遅かった場合を計測します。結構ばらつきがあるみたいなので、平均は見ません。

import timeit

## batch_testを定義(出力する部分は省略)
def batch_test(batch_size):
    for i in xrange(0,N,batch_size):
        match_cnt = 0
        t_predicted = np.argmax(mnn.forward(x_test[i:i+batch_size].T), axis = 0)
        match_cnt += np.sum(t_predicted == t_test[i:i+batch_size])

## バッチのサイズ
 size_list = [1, 5, 10, 25, 50, 100, 250, 500, 1000, 2500, 5000, 10000]

## 時間を計測
Time_list = []
for batch_size in size_list:
    s = "batch_test({0})".format(batch_size)
    res = timeit.repeat(setup = "from __main__ import batch_test", stmt = s, repeat = 100, number = 1)
    min_time = np.min(res); max_time = np.max(res)
    Time_list.append((batch_size, min_time, max_time))

表を出すだけでもいいんですが練習がてらグラフで可視化しましょう。色々楽なので一旦pandas.DataFrameに直してからプロットします。

## データフレーム化
from pandas import DataFrame
Time_df = DataFrame(Time_list,  columns = ["batch_size","min_time","max_time"])
print Time_df

<     batch_size  min_time  max_time
< 0            1  3.051640  6.661964
< 1            5  0.851746  1.937732
< 2           10  0.541192  1.394774
< 3           25  0.363268  0.797433
< 4           50  0.244394  0.466779
< 5          100  0.194512  0.466373
< 6          250  0.151969  0.376893
< 7          500  0.136464  0.372573
< 8         1000  0.116628  0.361746
< 9         2500  0.114554  0.312740
< 10        5000  0.128757  0.186853
< 11       10000  0.145491  0.324528

## プロット
from matplotlib import pyplot as plt
fig, axs = plt.subplots(1,1,figsize = (15,8), sharex = True)
Time_df.plot(kind = "bar", x = u"batch_size", ax = axs)
axs.set_ylabel("execution time (s)")

fig.savefig("batch_test.png")

f:id:Tawara:20161106103751p:plain

結果を見ると、max_time と min_time で結構な差があるものの、共通して言えるのはバッチ処理した方が明らかに速いってことですね。
今回はテストデータの数が10000枚なのでこれが限界ですが、テストデータ数が増えると、バッチのサイズが大き過ぎると逆に遅いってケースが出てきそうです。上の結果でもなんとなくそれっぽいですが、データ数もっと増えたほうがわかりやすそう。

なにはともあれ、バッチ処理で(少なくとも処理時間的には)効率的な処理ができました!

今回のまとめ

さらっと終わるはずだった、しかし現実は厳しかった...まあ100%自業自得ですが(処理の仕方を変えたり python2 使ったり)。その分勉強になることは色々あったのでまあいいかなと。

今回のポイントは以下。

ブロードキャストの都合上、重み行列を右からかける実装の方が楽なんですよね。まあ本の通りにするのもなんなので今後も今の方針で行こうと思います。

次回から 第4章 ニューラルネットワークの学習 です。遂に一番勉強したかった学習の部分に入れますね。ここはしっかり押さえつつ、第5章 誤差逆伝播法 につなげたいところ。
導入である2, 3章すら時間かかったのでこれから更に時間かかりそうですが、完走目指してがんばります!

ゼロから作るDeep Learning ―Pythonで学ぶディープラーニングの理論と実装

おまけ

今回の記事で、さりげなくmnist_for_py2とかsample_weight_for_py2.pklとか出てきました。 でもこの本は python3 環境を前提として書かれているので、こんなものは用意してくれてません。なので python2 使いは対処する必要がありました。この話をしようと思います。 「いや今だけでも python3 使えよ」って話ですが、python2使いたいんやん? 本当は移行した方がよさそうですけど、今のところそこまで困ってないし。

序盤で引っかかる

最初の方に、mnistのデータセットを取り込むフェーズがあったんですが、用意されているload_mnistを python2 環境でimportすると、

from dataset.mnist import load_mnist

以下のようなエラーが出ます。

ImportError: You should use Python 3.x

実はこのImportError著者が用意したエラーです。スクリプト覗くと、urllib.requestimport出来ないときに出るようになってます。

# coding: utf-8
try:
    import urllib.request
except ImportError:
    raise ImportError('You should use Python 3.x')
.
.
(以下コードが続く)
.
.

オラ、エラー出したったぞ。素直にpython3使えや..」って声が聞こえてきそうですが、
「嫌や!僕はpyhon2使うんや!」
てわけで対処するで。

なんでこのurllib.requestが必要なんかなーと思ってmnist.pyの中を見てみると、 urllib.request.urlretrieve (指定したURLからファイルを取得して指定した場所に置く関数)を使うためだと判明しました。幸いにして python2 のurllib.urlretrieveが対応していたので、 urllib.requestの代わりにurllibimportして、元々のurllib.request.urlretrieveurllib.urlretrieveに置き換えるだけでOK。

def _download(file_name):
    file_path = dataset_dir + "/" + file_name
    
    if os.path.exists(file_path):
        return

    print("Downloading " + file_name + " ... ")
    # urllib.request.urlretrieve(url_base + file_name, file_path) ## 元々のの記述
    urllib.urlretrieve(url_base + file_name, file_path) ## python2のために追加
    print("Done")

名前同じだけど python2 と python3 でモジュールの構成が異なっています。まあそもそもpython2ってurlliburllib2があってややこしいんですよね。python3だとこれらがurllibに一本化されていて(公式によるとurllib2の機能が urllib.request, urllib.errorに分割されて組み込まれたらしい)そこらへんいいですよね。..まあ python2 使いますけど。

ちなみに、さっきの公式ページでも紹介されていますがrequestsというモジュールがあります。こいつが相当便利みたいなのでHTTP関連の処理したかったらこっちの方がいいかもしれません。

途中でも引っかかる

データセットをダウンロードした後、今度は著者が用意してくれた重み(重み行列とバイアスベクトル)をロードして自作のNNに突っ込んであげるのですが..

import pickle as pic
with open("ch03\\sample_weight.pkl","rb") as f:
    parameters = pic.load(f)

結果がこれ↓。

< ValueError: unsupported pickle protocol: 3

pickleよく使ってるけど、何このエラー見たこと無い..と思って調べたら python3 の pickle の公式にも載ってますし、以下のような記事もありました。

stackoverflow.com

pickleで保存する際にいくつかプロトコルがあって、python2 が {0,1,2} しかないのに対し、python3 は {0,1,2,3,4} の 5 種類あります。で、python3 はデフォルトではプロトコルに 3 を選ぶので python2 でオブジェクトに戻そうとするとエラーが起きたという訳です。二年ぐらい使ってるのにprotocolって引数あるの初めて知ったわ..

対処法としては、python3 で一回 pickle ファイルをオブジェクトに戻して、プロトコル:2 で保存し直します。python2 を使うために python3 を使ってるので「もう python3 でやれよ」って幻聴が聞こえる..

python2 と python3 の違いの話って、rangeとかinputprintの仕様変更とか、文字コードの話とか、演算子の話とか、intの長さの話だけだた思ってたんですが pickle は知りませんでした。urllibも何となく変更されてるのは知ってたんですけど結構がっつりでしたね。そのうち移行するべきなんでしょうけどいつやるかってのはみんな悩みどころなのでしょうか。