人工知能に関する断創録

このブログでは人工知能のさまざまな分野について調査したことをまとめています（更新停止: 2019年12月31日）

Theanoによる畳み込みニューラルネットワークの実装 (1)

Deep Learning 機械学習 コンピュータビジョン Theano

Theanoによる多層パーセプトロンの実装（2015/6/18）のつづき。今回は、Deep Learning Tutorialの畳み込みニューラルネットワーク（Convolutional Neural Network: CNN, ConvNet）を実装してみる。

CNNは人間の視覚野を参考にした手法であり、画像認識に特化したDeep Learningアルゴリズムである。ImageNetの物体認識コンテストでぶっちぎりの成果を上げた手法はさまざまな工夫があるもののこのCNNをベースにしている。

本当は一般物体認識の実験をやりたいところだけどお楽しみは後に残しておいて、まずはMNISTの手書き数字認識を追試して感触をつかみたい。

ソースコード全体はここに置いた。

畳み込みニューラルネットワーク

まず今回実装する畳み込みニューラルネットワーク（CNN）の構成を図でまとめてみた（Deep Learning Tutorialの図は実装と関係ないじゃん・・・）

f:id:aidiary:20150626203849p:plain

一番左の入力から分類結果を出力する右側まで畳み込み層（convolution layer）とプーリング層（pooling layer）を何回か繰り返したあと最後に全結合した多層パーセプトロンが配置される構成になっている。前回実験した（2015/6/18）多層パーセプトロンは一番右側の部分にあたる。なんか一気に複雑になった・・・理解できんのこれ？

というわけでチュートリアルの実装に即して自分なりに理解したことをまとめてみたい。

畳み込み層

畳み込み層は入力画像に対してフィルタをかける（畳み込む）層である。画像の畳み込みによって画像内のパターンが検出できるようになる。畳み込みの数式は参考文献に書いてあるがここでは省略。Theanoではtheano.tensor.nnet.conv.conv2d()という関数が提供されているので実装する上ではあまり問題にならない。

先の図では28x28ピクセルの「7」という入力画像に対して5x5のそれぞれ異なるフィルタを20個かけて24x24の20枚の画像を出力している。フィルタをかけると画像サイズが少し小さくなる。入力画像サイズがW \times WでフィルタがH \times Hだと出力画像サイズは W - 2 \lfloor H/2 \rfloor \times W - 2 \lfloor H/2 \rfloorになる。ここで\lfloor \cdot \rfloorは小数点以下切り下げて整数化する演算。今回はW=28, H=5なので上式で計算すると24になる。フィルタをかけた出力画像は特徴マップと呼ばれることが多いようだ。

フィルタは4次元テンソルW[20,1,5,5]で表せる（テンソルの扱いnumpyのndarrayとほとんど同じ）この意味は5x5のフィルタで入力画像の枚数（チャネル数）が1で出力画像の枚数が20ということを意味している。MNISTのデータは白黒画像なので入力画像は1チャネルになる。もし入力がカラー画像のときはRGBの3チャネルである。今回の実装では1チャネルのフィルタが20枚使われる。フィルタを1つでなく、複数使うのがポイント。フィルタを増やすことで入力画像のさまざまな特徴を捉えられるようになる。

CNNがすごいのはこのフィルタを開発者が手動で設計するのではなく、学習によって自動獲得できるという点にある。最初、Wはランダムな値が入っていてよくわからない特徴にしか反応しないが、学習が進むと縦線や横線など画像認識に重要な特徴に強く反応するようになっていく。これがDeep Learningの力の源で表現学習と呼ばれる。これまでの説明からわかるようにCNNで更新（学習）するパラメータはフィルタWになる。

学習する前のランダムな値が入ったフィルタを可視化すると下のようになる。5x5のサイズのフィルタが20個ある。ランダムなパターンでなんかよくわからない（笑）

f:id:aidiary:20150626220554p:plain

CNNの学習が完了した後の同じ場所のフィルタを可視化すると下のようになる。

f:id:aidiary:20150626220604p:plain

少しわかりにくいが縦線や斜めの線が強調されるフィルタパターンが学習されたのががわかる。この20個のフィルタを数字の「7」の画像に畳み込むと下のような畳み込み層の出力画像が得られる。先の図で左から2つめの24x24ピクセルの20枚の画像がこれである。

f:id:aidiary:20150626220932p:plain

フィルタや畳み込み層の出力を可視化するコードはこの記事の最後に載せた。

プーリング層

プーリング層は畳み込み層の直後に置かれ、抽出された特徴の位置感度を低下させる働きがある。つまり、数字の「7」が画像の真ん中に書かれていても少し左や右にずれていても同じように「7」と判定するために必要ってこと。すごく難しそうだが処理自体は畳み込みよりずっと単純で2x2（poolsizeで指定）の矩形フィルタを入力画像内でずらしていき矩形内の最大の値を取り出して新しい画像を出力するような処理を行う。最大値に置き換える方法はmaxpoolingと呼ばれる。このプーリング層でも入力画像に比べて出力画像サイズは小さくなる。先の例では24x24の20枚の画像が12x12の20枚の画像になっている。画像のサイズは小さくなるが枚数は変わらない。Theanoではtheano.tensor.signal.downsample.max_pool_2d()という関数が用意されている。数式や動作の詳細は参考文献参照。なんでプーリングって名前が付いたのかよくわからないな。

LeNetConvPoolLayerクラス

Deep Learning Tutorialでは、畳み込み層とプーリング層のペアをLeNetConvPoolLayerというクラスで表現している。実際はペアである必要はなく、畳み込み層を3つ続けてプーリング層1つとかでもいいみたい。その場合はこのクラスは使えない。分けたほうがよいかもね。

LeNetというのは、畳み込みニューラルネットを発明したLeCunの最初のニューラルネットの名前に基づいているようだ。

class LeNetConvPoolLayer(object):
    """畳み込みニューラルネットの畳み込み層＋プーリング層"""
    def __init__(self, rng, input, image_shape, filter_shape, poolsize=(2, 2)):
        # 入力の特徴マップ数は一致する必要がある
        assert image_shape[1] == filter_shape[1]

        fan_in = np.prod(filter_shape[1:])
        fan_out = filter_shape[0] * np.prod(filter_shape[2:]) / np.prod(poolsize)

        W_bound = np.sqrt(6.0 / (fan_in + fan_out))
        self.W = theano.shared(
            np.asarray(rng.uniform(low=-W_bound, high=W_bound, size=filter_shape),
                       dtype=theano.config.floatX),  # @UndefinedVariable
            borrow=True)

        b_values = np.zeros((filter_shape[0],), dtype=theano.config.floatX)  # @UndefinedVariable
        self.b = theano.shared(value=b_values, borrow=T)

        # 入力の特徴マップとフィルタの畳み込み
        conv_out = conv.conv2d(
            input=input,
            filters=self.W,
            filter_shape=filter_shape,
            image_shape=image_shape)

        # Max-poolingを用いて各特徴マップをダウンサンプリング
        pooled_out = downsample.max_pool_2d(
            input=conv_out,
            ds=poolsize,
            ignore_border=True)

        # バイアスを加える
        self.output = T.tanh(pooled_out + self.b.dimshuffle('x', 0, 'x', 'x'))

        self.params = [self.W, self.b]

入力される画像がinputでフィルターはWで表現している。フィルタの値は多層パーセプトロンと同じくランダムな値で初期化される。bはバイアスでプーリングの後に各画像に対して加えられる。bはフィルタ数と同じ数の1次元配列なのでpooled_outの4Dテンソルと足し合わせるようにdimshuffleで次元を調整している。pooled_outの2つめの次元がbの次元数と一致するのがポイント。バイアスを加えたあとに活性化関数tanhを適用している。活性化関数にはtanhより高速で同じくらい高性能なReLU（Rectified Linear Unit）というのが最近ではよく使われるらしい。あとで置き換えて結果を比較してみたい。

Deep Learning Tutorialではプーリング層の出力に対してバイアスを加えて活性化関数を通すような実装になっているが、論文によっては畳み込み層の出力に対してバイアスを加えて活性化関数を通した後にプーリング層へというように説明されていることもある。どちらでもよいのかな？

このクラスは下のように使う。

    # 入力
    # 入力のサイズを4Dテンソルに変換
    # batch_sizeは訓練画像の枚数
    # チャンネル数は1
    # (28, 28)はMNISTの画像サイズ
    layer0_input = x.reshape((batch_size, 1, 28, 28))

    # 最初の畳み込み層+プーリング層
    # 畳み込みに使用するフィルタサイズは5x5ピクセル
    # 畳み込みによって画像サイズは28x28ピクセルから24x24ピクセルに落ちる
    # プーリングによって画像サイズはさらに12x12ピクセルに落ちる
    # 特徴マップ数は20枚でそれぞれの特徴マップのサイズは12x12ピクセル
    # 最終的にこの層の出力のサイズは (batch_size, 20, 12, 12) になる
    layer0 = LeNetConvPoolLayer(rng,
                input=layer0_input,
                image_shape=(batch_size, 1, 28, 28),  # 入力画像のサイズを4Dテンソルで指定
                filter_shape=(20, 1, 5, 5),           # フィルタのサイズを4Dテンソルで指定
                poolsize=(2, 2))

先の図では入力画像は「7」だけでたった1つのように描いたが、学習するときはミニバッチ単位でbatch_size=500枚の画像をまとめて渡す。このミニバッチ単位でフィルタのパラメータWとバイアスbを更新するためだ。

ここまでで先の図のlayer0がやっと終わった。

2つめの畳み込み層とプーリング層

layer0の出力は12x12ピクセルの画像（特徴マップ）が20枚になる。次はこの画像を入力としてlayer1でさらに畳み込み層とプーリング層に通す。layer0と異なり、複数チャネルの画像が入力される。

入力は20枚の画像から成るので20チャネルの画像と考えられる。この場合、畳み込みに使用するフィルタも20チャネル分用意する必要がある。複数チャネルの画像に複数チャネルのフィルタを畳み込むと1つの画像が出力される。はじめここを勘違いして理解するのに混乱していた・・・チュートリアルでは20チャネルのフィルタを50個用意しているので最終的に畳み込み層の出力として50個の画像（特徴マップ）が出力される。

フィルタは4次元テンソルW[50,20,5,5]で表せる。この意味は5x5のフィルタで入力画像のチャネル数が20で出力画像のチャネル数が50ということを意味している。前と同様に畳み込むと画像のサイズは少し小さくなり、8x8ピクセルになる。

さらにこの50枚の画像をプーリング層に通すことでサイズが小さくなり、4x4ピクセルの画像が50枚出力される。

コードで書くと下のようになる。前と同じくLeNetConvPoolLayerを再利用できる。入力はlayer0の出力なのでlayer0.outputが指定されている。

    # layer0の出力がlayer1への入力となる
    # layer0の出力画像のサイズは (batch_size, 20, 12, 12)
    # 12x12ピクセルの画像が特徴マップ数分（20枚）ある
    # 畳み込みによって画像サイズは12x12ピクセルから8x8ピクセルに落ちる
    # プーリングによって画像サイズはさらに4x4ピクセルに落ちる
    # 特徴マップ数は50枚でそれぞれの特徴マップのサイズは4x4ピクセル
    # 最終的にこの層の出力のサイズは (batch_size, 50, 4, 4) になる
    layer1 = LeNetConvPoolLayer(rng,
                input=layer0.output,
                image_shape=(batch_size, 20, 12, 12), # 入力画像のサイズを4Dテンソルで指定
                filter_shape=(50, 20, 5, 5),          # フィルタのサイズを4Dテンソルで指定
                poolsize=(2, 2))

ここまでで先の図のlayer1が終わった。

最後の多層パーセプトロン

畳み込みニューラルネットワークでは最後に全結合した多層パーセプトロンを配置して認識を行う。layer1の出力は4x4ピクセルの画像が50枚である。2次元の画像のままでは多層パーセプトロンに入力できないので、4x4x50=800次元のベクトルに変換する。

    # 隠れ層への入力
    # 画像のピクセルをフラット化する
    # layer1の出力のサイズは (batch_size, 50, 4, 4) なのでflatten()によって
    # (batch_size, 50*4*4) = (batch_size, 800) になる
    layer2_input = layer1.output.flatten(2)

多層パーセプトロンは前回（2015/6/18）作成したHiddenLayerとLogisticRegressionクラスを組み合わせて実装する。チュートリアルでは入力層が800ユニット、隠れ層（layer2）が500ユニット、出力層（layer3）が10ユニット（MNISTの手書き数字認識なので0-9の10種類）の多層パーセプトロンを構成している。

    # 全結合された隠れ層
    # 入力が800ユニット、出力が500ユニット
    layer2 = HiddenLayer(rng,
        input=layer2_input,
        n_in=50 * 4 * 4,
        n_out=500,
        activation=T.tanh)

    # 最終的な数字分類を行うsoftmax層
    layer3 = LogisticRegression(input=layer2.output, n_in=500, n_out=10)

確率的勾配降下法によるモデル訓練

モデルの訓練は多層パーセプトロンと同じく確率的勾配降下法によって行う。畳み込みニューラルネットワークは深い層構造を持つがプレトレーニング（各層を順番に学習して積み上げる方法）を使わずに古典的な勾配降下法だけで問題なく学習が行えることがわかっている。

いくつ層構造が積み重なってもTheanoの実装は非常にシンプル。

    # コスト関数を計算するシンボル
    cost = layer3.negative_log_likelihood(y)

    # パラメータ
    params = layer3.params + layer2.params + layer1.params + layer0.params

    # コスト関数の微分
    grads = T.grad(cost, params)

    # パラメータ更新式
    updates = [(param_i, param_i - learning_rate * grad_i) for param_i, grad_i in zip(params, grads)]

    # index番目の訓練バッチを入力し、パラメータを更新する関数を定義
    train_model = theano.function(
        [index],
        cost,
        updates=updates,
        givens={
            x: train_set_x[index * batch_size: (index + 1) * batch_size],
            y: train_set_y[index * batch_size: (index + 1) * batch_size]
        })

更新すべきパラメータが増えているだけで他は多層パーセプトロンと同じ。Theanoの自動微分最高！Early-Stoppingによる収束判定はロジスティック回帰 （2015/5/26）のときとほとんど同じなので省略。

実行時間の計測

CPUとGPUで実行時間を比較してみよう。まず、CPUを使った場合。

Optimization complete.
Best validation score of 0.910000 % obtained at iteration 19500, with test performance 0.920000 %
The code for file convolutional_mlp.py ran for 552.67m

学習に552分かかった。約9時間！夜寝るまえに動かして朝起きたらまだ動いてた。その間、ずっとファンが回ってたので電気代が心配だよ・・・エラー率は0.92%なので正解率99.08%。多層パーセプトロンのエラー率が1.65%で正解率が98.35%なのでさらに超えてる。もうMNISTでは限界値に近いかな。

次にGPUを使って同じ学習をしてみる。

Optimization complete.
Best validation score of 0.910000 % obtained at iteration 18500, with test performance 0.930000 %
The code for file convolutional_mlp.py ran for 41.39m

学習は41分で終わった。CPUに比べて約13倍速い。畳み込みニューラルネットワークの実験にはGPUが必須だな。

フィルタと畳み込み層の出力の可視化

チュートリアルにはないけどおまけでフィルタと畳み込み層の出力を可視化するスクリプトを書いた。学習の最後でlayer0とlayer1のオブジェクトをcPickleでファイルにダンプしておくコードを追加した。

    import cPickle
    cPickle.dump(layer0, open("layer0.pkl", "wb"))
    cPickle.dump(layer1, open("layer1.pkl", "wb"))

こうして学習結果をファイルにダンプしておけばまた最初から学習しなおさなくてすむ。学習したフィルタの重みとバイアスはこのオブジェクトのインスタンス変数Wやbで取得できる。

汎用的ではないけどフィルタと畳み込み層（プーリング層ではない）の出力画像を可視化するスクリプト。学習したフィルタ重みとバイアスをもとに再度畳み込みニューラルネットを構築しなおしている。layer0.outputやlayer1.outputがそのまま使えると期待したけどTheanoのエラーが出てうまいくいかなかった。あとでもっとよい方法が思いついたら書き直すかも。

#coding: utf-8
import cPickle
import pylab
import matplotlib.pyplot as plt

import theano
import theano.tensor as T
from theano.tensor.nnet import conv
from theano.tensor.signal import downsample
from logistic_sgd import load_data
from convolutional_mlp import LeNetConvPoolLayer

def visualize_filter(layer):
    """引数に指定されたLeNetConvPoolLayerのフィルタを描画する"""
    W = layer.W.get_value()
    n_filters, n_channels, h, w = W.shape
    plt.figure()
    pos = 1
    for f in range(n_filters):
        for c in range(n_channels):
            plt.subplot(n_channels, n_filters, pos)
            plt.subplots_adjust(wspace=0.1, hspace=0.1)
            plt.imshow(W[f, c], cmap=pylab.cm.gray_r)
            plt.axis('off')
            pos += 1
    plt.show()

def feedforward(input, layer, image_shape):
    """inputをlayerに通した畳み込み層の結果（畳み込み層の結果を可視化するため）と
    最終的な出力の特徴マップ（次の層への入力のため）を返す"""
    conv_out = conv.conv2d(input,
                           filters=layer.W,
                           filter_shape=layer.W.get_value().shape,
                           image_shape=image_shape)
    pooled_out = downsample.max_pool_2d(input=conv_out,
                                        ds=(2, 2),
                                        ignore_border=True)
    output = T.tanh(pooled_out + layer.b.dimshuffle('x', 0, 'x', 'x'))
    return conv_out, output

if __name__ == "__main__":
    layer0 = cPickle.load(open("layer0.pkl", "rb"))
    layer1 = cPickle.load(open("layer1.pkl", "rb"))
    print layer0, layer1

    # 最初の畳み込み層のフィルタの可視化
    visualize_filter(layer0)
    visualize_filter(layer1)

    # 各層の出力の可視化
    input = T.tensor4()
    # 入力画像は1のためimage_shapeのbatch_size=1となる
    layer0_conv_out, layer0_out = feedforward(input, layer0, image_shape=(1, 1, 28, 28))
    layer1_conv_out, layer1_out = feedforward(layer0_out, layer1, image_shape=(1, 20, 12, 12))

    # 画像を1枚受けて、畳み込み層の出力を返す関数を定義
    f0 = theano.function([input], layer0_conv_out)
    f1 = theano.function([input], layer1_conv_out)

    # 入力はテストデータから適当な画像を一枚入れる
    datasets = load_data("mnist.pkl.gz")
    test_set_x, test_set_y = datasets[2]
    input_image = test_set_x.get_value()[0]
    layer0_conv_image = f0(input_image.reshape((1, 1, 28, 28)))
    layer1_conv_image = f1(input_image.reshape((1, 1, 28, 28)))

    # 畳み込み層の出力画像を可視化
    print layer0_conv_image.shape
    plt.figure()
    for i in range(20):
        plt.subplot(4, 5, i + 1)
        plt.axis('off')
        plt.imshow(layer0_conv_image[0, i], cmap=pylab.cm.gray_r)
    plt.show()

    print layer1_conv_image.shape
    plt.figure()
    for i in range(50):
        plt.subplot(5, 10, i + 1)
        plt.axis('off')
        plt.imshow(layer1_conv_image[0, i], cmap=pylab.cm.gray_r)
    plt.show()

参考文献

はっきり言ってDeep Learning TutorialのCNNの説明はわかりにくく、これだけでは理解できなかった。次のサイト、書籍、論文で何とか自分なりに理解できたと思う。比較して確認したから勘違いはしてないと思うのだけど・・・

Theanoによる多層パーセプトロンの実装

Deep Learning 機械学習 Theano

Theanoによるロジスティック回帰の実装（2015/5/26）のつづき。今回は、Deep Learning Tutorialの多層パーセプトロン（Multilayer Perceptron）を実装してみる。タスクは前回と同じMNISTの手書き数字認識。多層パーセプトロンはこれまでも何回か実装してきた（2010/8/29）けど今回はTheanoを使ってみようという趣旨。これまでよりずいぶん簡単に実装できることがわかる。

ソースコード全体はここに置いた。

多層パーセプトロン

f:id:aidiary:20150617204542p:plain

上の図は入力層、隠れ層、出力層3層パーセプトロンの模式図。入力層に入力されたベクトルxは、入力層と隠れ層間の重み行列 W^{(1)} とバイアスベクトル b^{(1)} によって

h(x) = {\rm sigmoid}(b^{(1)} + W^{(1)})

で変換され、隠れ層の出力となる。活性化関数として使われるロジスティック・シグモイド関数 sigmoid によって入力ベクトルは非線形変換される。つまりもとの空間で線形分離不可能でもこの非線形変換によって隠れ層の出力は線形分離可能になる。sigmoidの代わりにtanhを使ってもよい。

隠れ層の出力 h(x) は、今度は出力層の入力になり、隠れ層と出力層間の重み行列 W^{(2)} とバイアスベクトル b^{(2)} によって

f(x) = {\rm softmax}(b^{(2)} + W^{(2)} h(x))

で変換され、出力層の出力となる。パラメータは W^{(1)}, b^{(1)}, W^{(2)}, b^{(2)} の4つ。

上の図で書いたように出力層での計算はロジスティック回帰の式（2015/5/26）

\displaystyle P(Y=i|x, W, b) = {\rm softmax}_{i} (W x + b) = \frac{e^{W_i x + b_i}}{\sum_{j} e^{W_j x + b_j}}

と同じになる点がポイント。多層パーセプトロンは入力データをそのままロジスティック回帰に入れるのではなく、入力層と隠れ層間で非線形変換してからロジスティック回帰に入れるいう見方ができる。この非線形変換によってロジスティック回帰より表現力が高くなるってのがポイントのようだ。

Deep Learningのチュートリアルがロジスティック回帰から始まってるのが不思議だったんだけどこういう理由があったのかと納得した。

隠れ層の実装

出力層の部分は前回のロジスティック回帰の実装がそのまま使えるため隠れ層だけ実装すればよい。

class HiddenLayer(object):
    def __init__(self, rng, input, n_in, n_out, W=None, b=None, activation=T.tanh):
        """隠れ層の初期化
        rng: 乱数生成器（重みの初期化で使用）
        input: ミニバッチ単位のデータ行列（n_samples, n_in)
        n_in: 入力データの次元数
        n_out: 隠れ層のユニット数
        W: 隠れ層の重み
        b: 隠れ層のバイアス
        activation: 活性化関数
        """
        self.input = input

        # 隠れ層の重み（共有変数）を初期化（[Xavier10]による）
        if W is None:
            W_values = np.asarray(
                rng.uniform(low=-np.sqrt(6.0 / (n_in + n_out)),
                            high=np.sqrt(6.0 / (n_in + n_out)),
                            size=(n_in, n_out)),
                dtype=theano.config.floatX)
            if activation == theano.tensor.nnet.sigmoid:
                W_values *= 4
            W = theano.shared(value=W_values, name='W', borrow=True)

        # 隠れ層のバイアス（共有変数）を初期化
        if b is None:
            b_values = np.zeros((n_out, ), dtype=theano.config.floatX)
            b = theano.shared(value=b_values, name='b', borrow=True)

        self.W = W
        self.b = b

        # 隠れ層の出力を計算するシンボルを作成
        lin_output = T.dot(input, self.W) + self.b
        if activation is None:  # 線形素子の場合
            self.output = lin_output
        else:  # 非線形な活性化関数を通す場合
            self.output = activation(lin_output)

        # 隠れ層のパラメータ
        self.params = [self.W, self.b]

重み行列とバイアスベクトルは学習時の更新対象なので共有変数に格納している。バイアスベクトルは0で初期化すればよいが、重み行列の初期値の与え方が少し巧妙。活性化関数にtanhを使うときは

\displaystyle \biggl[ -\sqrt{\frac{6}{in + out}}, \sqrt{\frac{6}{in + out}} \biggr]

の間の一様乱数とし、活性化関数にsigmoidを使うときは

\displaystyle \biggl[ -4 \sqrt{\frac{6}{in + out}}, 4 \sqrt{\frac{6}{in + out}} \biggr]

の間の一様乱数とすると訓練初期で学習が効率的に進むことが保証されるとのこと。ここでinは入力層のユニット数、outは隠れ層のユニット数を意味する。証明はよくわからない・・・

多層パーセプトロンの実装

多層パーセプトロンは隠れ層とロジスティック回帰を組み合わせて実装できる。隠れ層の出力がロジスティック回帰への入力となるのがポイント。

class MLP(object):
    def __init__(self, rng, input, n_in, n_hidden, n_out):
        # 多層パーセプトロンは隠れ層とロジスティック回帰で表される出力層から成る
        # 隠れ層の出力がロジスティック回帰の入力になる点に注意
        self.hiddenLayer = HiddenLayer(rng=rng, input=input, n_in=n_in, n_out=n_hidden, activation=T.tanh)
        self.logRegressionLayer = LogisticRegression(input=self.hiddenLayer.output, n_in=n_hidden, n_out=n_out)

        # L1/L2正則化の正則化項を計算するシンボル
        self.L1 = abs(self.hiddenLayer.W).sum() + abs(self.logRegressionLayer.W).sum()
        self.L2_sqr = (self.hiddenLayer.W ** 2).sum() + (self.logRegressionLayer.W ** 2).sum()

        # MLPの誤差関数を計算するシンボル
        # 出力層にのみ依存するのでロジスティック回帰の実装と同じでよい
        self.negative_log_likelihood = self.logRegressionLayer.negative_log_likelihood

        # 誤差を計算するシンボル
        self.errors = self.logRegressionLayer.errors

        # 多層パーセプトロンのパラメータ
        self.params = self.hiddenLayer.params + self.logRegressionLayer.params

誤差関数は出力層のみに依存するのでロジスティック回帰の負の対数尤度をそのまま使う。ただし、次で見るように過学習を防止するためにL1正則化とL2正則化項を追加している。正則化項のシンボルもTheanoを使うと非常に簡単に書ける。

確率的勾配降下法によるモデル訓練

最後にこれまで定義したクラスを使って実際にモデルを訓練する関数を定義する。データのロードなどは前回とあまり変わらないので省略している。

def test_mlp(learning_rate=0.01, L1_reg=0.00, L2_reg=0.0001, n_epochs=1000, dataset='mnist.pkl.gz', batch_size=20, n_hidden=500):

    # データのロードなどは省略

    # MLPを構築
    classifier = MLP(rng=rng, input=x, n_in=28 * 28, n_hidden=n_hidden, n_out=10)

    # コスト関数を計算するシンボル
    cost = classifier.negative_log_likelihood(y) + L1_reg * classifier.L1 + L2_reg * classifier.L2_sqr

    # index番目のテスト用ミニバッチを入力してエラー率を返す関数を定義
    test_model = theano.function(
        inputs=[index],
        outputs=classifier.errors(y),
        givens={
            x: test_set_x[index * batch_size: (index + 1) * batch_size],
            y: test_set_y[index * batch_size: (index + 1) * batch_size]
        })

    # index番目のバリデーション用ミニバッチを入力してエラー率を返す関数を定義
    validate_model = theano.function(
        inputs=[index],
        outputs=classifier.errors(y),
        givens={
            x: valid_set_x[index * batch_size: (index + 1) * batch_size],
            y: valid_set_y[index * batch_size: (index + 1) * batch_size]
        })

    # コスト関数の各パラメータでの微分を計算
    gparams = [T.grad(cost, param) for param in classifier.params]

    # パラメータ更新式
    updates = [(param, param - learning_rate * gparam) for param, gparam in zip(classifier.params, gparams)]

    # index番目の訓練バッチを入力し、パラメータを更新する関数を定義
    # 戻り値としてコストが返される
    train_model = theano.function(
        inputs=[index],
        outputs=cost,
        updates=updates,
        givens={
            x: train_set_x[index * batch_size: (index + 1) * batch_size],
            y: train_set_y[index * batch_size: (index + 1) * batch_size]
        })

入力層のユニット数は28x28=784、中間層のユニット数は500、出力層のユニット数は10。コスト関数は負の対数尤度に加えてL1とL2の正則化項がついている。

多層パーセプトロンではコスト関数の偏微分を効率よく計算するためにバックプロパゲーションが使われる（2014/1/22）がTheanoは自動微分を使えるので面倒な実装を考えなくてすむ。T.grad()でコスト関数のパラメータごとの偏微分を求め、updatesでパラメータ更新式を定義するだけ。パラメータは W^{(1)}, b^{(1)}, W^{(2)}, b^{(2)} と4つあるのでforループで回しているが非常に直観的な書き方だ。これはすごい！

Early-Stoppingによる収束判定はロジスティック回帰 （2015/5/26）とほとんど同じなので省略。

実行時間の計測

最後にCPUとGPUで実行時間を比較してみよう。まず、CPUを使った場合。

Optimization complete. Best validation score of 1.690000 % obtained at iteration 2070000, with test performance 1.650000 %
The code for file mlp.py ran for 211.46m

収束まで207万回のパラメータ更新をして3時間半もかかった・・・エラー率は1.65%なので正解率は98.35%となる。ロジスティック回帰の正解率は92.51%なのでかなりよい結果が得られた。エラー率の推移をグラフで表示してみよう。

f:id:aidiary:20150618202220p:plain:w450

急激に落ちた後は少しずつエラー率が下がっていった。Early-Stoppingで収束判定するために使うpatienceの初期値は10000だけど最終的に473万近くまで上がる。途中でなかなかエラー率が減らず、スクリプトを終了したくなるけどEarly-Stoppingを信じてじっと耐えると徐々にではあるがエラー率が着実に減ることがわかった。非常に巧妙な収束判定のようだ。

最後にGPUを使って同じ学習をしてみる。

Using gpu device 0: GeForce GTX 760 Ti OEM
Optimization complete. Best validation score of 1.690000 % obtained at iteration 2367500, with test performance 1.650000 %
The code for file mlp.py ran for 66.09m

60分で終わった。速い。

参考

Theanoによるロジスティック回帰の実装

Deep Learning 機械学習 Theano

Theanoによる2クラスロジスティック回帰の実装（2015/5/19）のつづき。

今回からDeep Learning Tutorialの内容にそって実装していく。このチュートリアルのコードやデータはGithub上で公開されているのでcloneしておいた。

git clone https://github.com/lisa-lab/DeepLearningTutorials

さっそく、最初の例のロジスティック回帰（Logistic Regression）でMNISTの手書き数字認識を試してみよう。載せているソースコードは上のリポジトリとほとんど同じ（コメントを日本語で補足したくらい）。ここでは全部掲載したりせずに自分がポイントだと思ったところだけまとめておきたい。

ソースコード全体はここに置いた。

MNISTデータのロード

MNISTは手書き数字の画像データセット。前に多層パーセプトロンで使用したことがあった。

このチュートリアルでは、MNISTのデータをPythonのcPickleモジュールでロードできる形式に圧縮したmnist.pkl.gzが提供されているのでそちらを使う。このデータはMNISTの70000サンプルが

  • 50000の訓練セット
  • 10000のバリデーションセット
  • 10000のテストセット

にあらかじめ分割されている。各サンプルは28x28=784次元のデータで値は0.0から1.0に正規化されている。

f:id:aidiary:20150525203441p:plain

学習データをロードする関数

def load_data(dataset):
    """データセットをロードしてGPUの共有変数に格納"""
    f = gzip.open(dataset, 'rb')
    train_set, valid_set, test_set = cPickle.load(f)
    f.close()

    def shared_dataset(data_xy, borrow=True):
        data_x, data_y = data_xy

        # 共有変数には必ずfloat型で格納
        shared_x = theano.shared(np.asarray(data_x, dtype=theano.config.floatX), borrow=borrow)
        shared_y = theano.shared(np.asarray(data_y, dtype=theano.config.floatX), borrow=borrow)

        # ラベルはint型なのでキャストして返す
        return shared_x, T.cast(shared_y, 'int32')

    test_set_x, test_set_y = shared_dataset(test_set)
    valid_set_x, valid_set_y = shared_dataset(valid_set)
    train_set_x, train_set_y = shared_dataset(train_set)

    rval = [(train_set_x, train_set_y),
            (valid_set_x, valid_set_y),
            (test_set_x, test_set_y)]

    return rval

ファイルからcPickleでロードした後に共有変数に格納している。共有変数にするとGPUのメモリ領域に保存されるため学習時に高速に読み込めるという利点がある。共有変数の値は、get_value()で取得できる。また、shapeで配列のサイズが調べられる。

train_set_x.get_value().shape
(50000, 784)

行にサンプル、列に特徴量という一般的な形式で格納されていることがわかる。

訓練データを共有変数に格納するときは必ずfloatX型にする必要がある。本来は整数型のラベル（今回は数字認識なので0から9）も同様。ただし、ラベルを使うときは整数型の方が都合がよいためT.cast()で整数型にキャストしてから返している。キャストするとSharedVariableがTensorVariableになってしまい、get_value()で値が取得できなくなるみたい・・・

先の数字の画像を描画するには、shared_dataset()の前に下のコードを挿入すればOK。共有変数に入れてからget_value()で取り出して描画しようとするとなぜかMemoryErrorで描画できなかった。

    import pylab
    train_set_x, train_set_y = train_set
    pylab.figure()
    for index in range(100):
        pylab.subplot(10, 10, index + 1)
        pylab.axis('off')
        pylab.imshow(train_set_x[index, :].reshape(28, 28), cmap=pylab.cm.gray_r, interpolation='nearest')
    pylab.show()
    exit()

多クラスロジスティック回帰

ロジスティック回帰のパラメータは、重み行列Wとバイアスベクトルbの2つ*1。パラメータは適当な値（今回は0）で初期化し、あとで更新できるように共有変数で定義しておく。

        # 重み行列を初期化
        self.W = theano.shared(value=np.zeros((n_in, n_out), dtype=theano.config.floatX),
                               name='W',
                               borrow=True)

        # バイアスベクトルを初期化
        self.b = theano.shared(value=np.zeros((n_out,), dtype=theano.config.floatX),
                               name='b',
                               borrow=True)

ロジスティック回帰では、入力ベクトルxを与えたとき、それがクラスiに分類される確率は

\displaystyle P(Y=i|x, W, b) = {\rm softmax}_{i} (W x + b) = \frac{e^{W_i x + b_i}}{\sum_{j} e^{W_j x + b_j}}

で定義される。この数式をTheanoのシンボルで表すと

        # 各サンプルが各クラスに分類される確率を計算するシンボル
        # 全データを行列化してまとめて計算している
        # 出力は(n_samples, n_out)の行列
        self.p_y_given_x = T.nnet.softmax(T.dot(input, self.W) + self.b)

最初に見たとき、T.dot()の中身が数式と逆じゃないかと思ったのだがこれで正しかった。先の数式はサンプルxを1つだけ与えたときの確率を計算しているのに対し、このコードは複数のサンプル（ミニバッチ単位）を行列化してinputに与え、それらのサンプルの確率をまとめて計算している。

確率が求まったら最終的に一番高い確率が得られるクラスに分類する。

y_{pred} = {\rm argmax}_{i} P(Y=i|x,W,b)

TheanoではT.argmax()にaxis=1を指定することでp_y_given_xの各行（サンプルに相当）において一番確率が高いインデックス（クラスに相当）がまとめて取得できる。

        # 確率が最大のクラスのインデックスを計算
        # 出力は(n_samples,)のベクトル
        self.y_pred = T.argmax(self.p_y_given_x, axis=1)

self.y_predは現在のパラメータにおける各サンプルの分類クラスを意味している。

Theanoの厄介なところはself.p_y_given_xなどのシンボルをprintしてみても具体的な値が表示されないのでデバッグしにくいところだと思う。一応、デバッグのコツがあるのであとでTheanoの使い方としてまとめたい。

コスト関数の定義

2クラスロジスティック回帰のコスト関数は交差エントロピー誤差だったが、多クラスロジスティック回帰のコスト関数は負の対数尤度（NLL）を使う（追記）ここによるとどうやら同じものらしい・・・要確認。

\displaystyle NLL(\theta=\{W, b\}, D) = - \sum_{i=0}^{|D|} \log (P(Y=y^{(i)} | x^{(i)}, W, b))

TheanoでNLLを実装すると下のようになる。

    def negative_log_likelihood(self, y):
        """誤差関数である負の対数尤度を計算するシンボルを返す
        yにはinputに対応する正解ラベルを渡す"""
        # 式通りに計算するとsumだがmeanの方がよい
        return -T.mean(T.log(self.p_y_given_x)[T.arange(y.shape[0]), y])

NLLを計算するためには訓練データ集合 (x^{(i)}, y^{(i)})が必要なので正解クラスを表すシンボルyを引数で与えている。対応するxはNLLで使われているシンボルself.p_y_given_xを計算するときにすでにinputとして与えているためこの関数で渡す必要はない（ここら辺がややこしく感じた）。

行列をスライスしている[T.arange(y.shape[0]), y]の理解が少し難しい。self.p_y_given_xは各訓練データが各クラスに分類される確率を表した行列になる。そのT.logをとっても同じ形状の行列になる。 その行列から訓練データの正解のクラスy^{(i)}に対応する確率P(Y=y^{(i)})のみを集めているのが先のスライスの意味。

数式ではスライスした確率の和を取っているが、コードでは平均を取っている。チュートリアル資料によると平均にすると学習率の設定がミニバッチのサイズに依存しにくくなるという利点があるらしい。これは初めて知った。

この関数の戻り値はNLLを計算するためのシンボルであることに注意が必要。何か具体的な値が計算されて戻されるわけではない。具体的な計算をするためにはこのシンボルをtheano.function()で使う必要がある。

このNLLを最小化するパラメータWとbを求めるのがロジスティック回帰の学習に当たる。勾配降下法を使うためにはコスト関数NLLの各パラメータ（Wとb）での偏微分が必要になるが、Theanoの自動微分で求められるため自分で定義しなくてよい。

モデルの訓練

モデルの訓練にはミニバッチ確率的勾配降下法（MSGD）を使う。確率的勾配降下法（SGD）はただ1つのサンプルで1回だけパラメータを更新するのに対し、MGSDでは複数のサンプルをまとめて1回パラメータを更新する。

何番目のミニバッチを用いるかはシンボルindexで表す。たとえば、訓練データのindex番目のサンプル集合は下のコードで取得できる。

  train_set_x[index * batch_size: (index + 1) * batch_size]
  train_set_y[index * batch_size: (index + 1) * batch_size]

train_set_xとtrain_set_yもシンボルなのでprintしても値は表示できない。

確率的勾配降下法では、コスト関数とその微分が必要になる。コスト関数は先に述べたようにNLL（negative_log_likelihood）を使う。コスト関数の微分は、Theanoの自動微分を使えば下記のように簡単に書ける。これはすごい！

    cost = classifier.negative_log_likelihood(y)
    g_W = T.grad(cost=cost, wrt=classifier.W)
    g_b = T.grad(cost=cost, wrt=classifier.b)

パラメータの更新式は、

    # パラメータ更新
    updates = [(classifier.W, classifier.W - learning_rate * g_W),
               (classifier.b, classifier.b - learning_rate * g_b)]

となる。今回はパラメータがWとbの2つあるため更新式も2つ指定している。ここまでいろいろシンボルを定義してきてようやくミニバッチ確率的勾配降下法を行う関数が定義できる。

    # index番目の訓練バッチを入力し、パラメータを更新する関数を定義
    # 戻り値としてコストが返される
    # この関数の呼び出し時にindexに具体的な値が初めて渡される
    train_model = theano.function(
        inputs=[index],
        outputs=cost,
        updates=updates,
        givens={
            x: train_set_x[index * batch_size: (index + 1) * batch_size],
            y: train_set_y[index * batch_size: (index + 1) * batch_size]
        })

関数の入力には何番目のミニバッチを使って訓練するかを表すインデックスを与える。関数の出力は更新されたパラメータでのコストが返される。givensは関数を実行する際にシンボルを置き換えるオプションである。ここでは、xをindex番目のミニバッチの訓練データのシンボル、yをindex番目のミニバッチのラベルのシンボルに置き換えている。関数の引数のindexはgivensで使われるだけみたい。

yはシンボルcostの計算に使われていたのにすぐ気付いたが、xがどこで使われているかすぐにわからなかった。いろいろ調べていくとupdatesで使われるclassifierオブジェクトの引数として渡されていいたのに気付いた。

    # MNISTの手書き数字を分類するロジスティック回帰モデル
    # 入力は28ピクセルx28ピクセルの画像、出力は0から9のラベル
    # 入力はシンボルxを割り当てておいてあとで具体的なデータに置換する
    classifier = LogisticRegression(input=x, n_in=28*28, n_out=10)

結局のところxはgivensでミニバッチのデータを表すシンボルに置き換えられ、それが引数inputに渡されるという仕組みになっている。このinputは前に見たように出力確率（classifier.p_y_given_x）の計算で使われていた。これは慣れないときつい・・・

モデルの評価

モデルの良さはエラー率で評価する。まずはエラー率を計算するシンボルを返す関数を定義する。

    def errors(self, y):
        """分類の誤差率を計算するシンボルを返す
        yにはinputに対応する正解クラスを渡す"""
        return T.mean(T.neq(self.y_pred, y))

self.y_predは前に書いたように現在のパラメータでの各サンプルの予測クラスを表すシンボル。yは正解クラスを意味する。T.neq()で両者が等しくない要素だけ1になるベクトルが返る。これの平均を取ればエラー率になる。T.sum()で和を取れば分類エラーになったサンプル数（Zero-One Loss）になる。

このerrors()の戻り値もシンボルのため関数化しないと使えない。バリデーションデータとテストデータのそれぞれについてエラー率を計算する関数を定義する。バリデーションデータはEarly-stoppingによる収束判定で使われる。最終的なモデルの評価はテストデータのエラー率で評価する。

    # index番目のバリデーション用ミニバッチを入力してエラー率を返す関数を定義
    validate_model = theano.function(
        inputs=[index],
        outputs=classifier.errors(y),
        givens={
            x: valid_set_x[index * batch_size: (index + 1) * batch_size],
            y: valid_set_y[index * batch_size: (index + 1) * batch_size]
        })


    # index番目のテスト用ミニバッチを入力してエラー率を返す関数を定義
    test_model = theano.function(
        inputs=[index],
        outputs=classifier.errors(y),
        givens={  # ここで初めてシンボル x, y を具体的な値で置き換える
            x: test_set_x[index * batch_size: (index + 1) * batch_size],
            y: test_set_y[index * batch_size: (index + 1) * batch_size]
        })

Early-Stoppingによる収束判定

Early-Stopping（早期終了）は、バリデーションセットを用いてモデルのパフォーマンスをモニタリングすることで過学習を防ぐ手法とのこと。WikipediaのEarly-Stoppingの項目を読むといくつか実装の仕方があるようだがこのチュートリアルではpatienceを利用した収束判定をしている。

    # モデル訓練
    print 'training the model ...'

    # eary-stoppingのパラメータ
    patience = 5000
    patience_increase = 2
    improvement_threshold = 0.995
    validation_frequency = min(n_train_batches, patience / 2)

    best_validation_loss = np.inf
    test_score = 0
    start_time = time.clock()

    done_looping = False
    epoch = 0

    while (epoch < n_epochs) and (not done_looping):
        epoch = epoch + 1
        for minibatch_index in xrange(n_train_batches):
            # minibatch_index番目の訓練データのミニバッチを用いてパラメータ更新
            minibatch_avg_cost = train_model(minibatch_index)

            # validation_frequency回の更新ごとにバリデーションセットによるモデル検証が入る
            iter = (epoch - 1) * n_train_batches + minibatch_index
            if (iter + 1) % validation_frequency == 0:
                # バリデーションセットの平均エラー率を計算
                validation_losses = [validate_model(i) for i in xrange(n_valid_batches)]
                this_validation_loss = np.mean(validation_losses)
                print "epoch %i, minibatch %i/%i, validation error %f %%" % (epoch, minibatch_index + 1, n_train_batches, this_validation_loss * 100)

                # エラー率が十分改善したならまだモデル改善の余地があるためpatienceを上げてより多くループを回せるようにする
                if this_validation_loss < best_validation_loss:
                    if this_validation_loss < best_validation_loss * improvement_threshold:
                        patience = max(patience, iter * patience_increase)
                        print "*** iter %d / patience %d" % (iter, patience)

                    best_validation_loss = this_validation_loss

                    # テストセットを用いたエラー率も求めておく
                    test_losses = [test_model(i) for i in xrange(n_test_batches)]
                    test_score = np.mean(test_losses)
                    print "    epoch %i, minibatch %i/%i, test error of best model %f %%" % (epoch, minibatch_index + 1, n_train_batches, test_score * 100)

            # patienceを超えたらループを終了
            if patience <= iter:
                done_looping = True
                break

    end_time = time.clock()
    print "Optimization complete with best validation score of %f %%, with test performance %f %%" % (best_validation_loss * 100, test_score * 100)
    print "The code run for %d epochs, with %f epochs/sec" % (epoch, 1.0 * epoch / (end_time - start_time))

このコードはモデル訓練をしている一番重要なループだけれど読み解くのが少し大変。私は以下のように理解した。

  • モデルの訓練は訓練データのミニバッチ単位で行う。
  • エポック epoch は訓練データのミニバッチを使い切ったら+1される。訓練データを一通り使いきっても収束していなければまた最初のミニバッチに戻って同じデータを繰り返し使ってパラメータ更新する。
  • 最初の while ループは、最大n_epochs 回エポックを繰り返されたら更新を終えることを意味する。これによりたとえ収束しなくてもプログラムはいつかは終わることが保証される。
  • iterはこれまで訓練に使用したミニバッチの数が入る。つまり、パラメータ更新回数を意味する。
  • パラメータ更新回数がpatienceを超えたらもうこれ以上更新してもモデルは改善しないと判断されてループを抜ける。逆に言うと、patienceの値が上がるとまだまだパラメータ更新によるモデル改善の余地があると判断されてことになる。
  • patienceを上げるかどうかの判断はバリデーションセットを用いて決められる。
  • validation_frequency回更新されるたびにモデルの検証が行われる。
  • バリデーションセットを用いて平均エラー率が計算され、この平均エラー率がこれまでの最良の平均エラー率より十分な程度下がっていればまだまだモデル改善の余地ありと判断してpatienceを引き上げる。十分かどうかはimprovement_thresholdで調整できる。
  • この検証の際にテストセットを用いたモデルのエラー率の評価も合わせて行う。

デフォルトで patience は5000になっているけど経過をたどると6140まで引き上げられていた。

実行時間の計測

最後にCPUとGPUで実行時間を比較してみよう。まず、CPUを使った場合。

Optimization complete with best validation score of 7.500000 %, with test performance 7.489583 %
The code run for 74 epochs, with 5.222703 epochs/sec

real    0m17.656s
user    0m0.015s
sys     0m0.031s

学習に17秒かかった。CPUでもかなり早い。エラー率は7.49%なので正解率は92.51%。単純なロジスティック回帰でもかなりの精度が出ていることがわかる。次にGPUを使った場合。

Using gpu device 0: GeForce GTX 760 Ti OEM

Optimization complete with best validation score of 7.500000 %, with test performance 7.489583 %
The code run for 74 epochs, with 9.571053 epochs/sec

real    0m11.927s
user    0m0.000s
sys     0m0.031s

17秒が11秒に改善されている。ロジスティック回帰はパラメータの収束が早いため改善効果はあまり大きくない。1秒あたりの更新エポック数を比較すると5エポックから9エポックになっているためGPUを使うことで着実に速くなっていることが確認できる。

次回は多層パーセプトロン。

参考

うちと同じようにDeep Learning Tutorialに沿って実装を解説されている。Theanoの検証例が豊富で非常に参考になった。今後も参考にさせてもらう予定。

*1:入力データに1を追加することでバイアスベクトルbを重み行列Wに含めることもできるが今回はわけている