SageでTheanoの手書き数字認識（MNIST）を試す

参考サイト

ここでは、 人工知能に関する断創録 のTheanoに関連する記事をSageのノートブックで実装し、Thenoの修得を試みます。

今回は、TheanoのTutorialからMNISTの手書き数字認識を以下のページを参考にSageのノートブックで試してみます。

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

処理系をSageからPythonに変更

SageでTheanoのtutorialのMNISTを実行すると、以下のy_predがベクトルではなくスカラーになり、うまく処理できません。

	self.y_pred = T.argmax(self.p_y_given_x, axis=1)			
		

そこで、ノートブックの処理系をSageからPythonに切り替えます。上部の左から4つめのプルダウンメニューから 「python」を選択してください。

必要なライブラリをインポート

今回は、ちょっと多くのライブラリを使用します。

データのセットアップ

データのダウンロード

TheanoのTutorialで使用しているMNISTは手書き数字の画像データセットは、 cPickleモジュールでロードできる形式に圧縮したmnist.pkl.gzが以下のサイトから ダウンロードできます。

• http://deeplearning.net/data/mnist/mnist.pkl.gz

MNISTのデータには、70000の手書き数字データが収録されおり、それを以下の様に分割して使用します。

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

TensorFlowのMNIST for ML BeginnersのMNISTデータの説明が分かりやすいので、引用します。

1/10のサブセットを作成

Sageサーバのメモリは1Gと少ないため、MNISTのデータをそのまま使用するとメモリ不足になります。 そこで、1/10のサイズのサブセットデータmini_mnist.pkl.gzを作成し、これを使用して数字を認識することにします。

サンプルデータについて

どのような画像データが入っているのか、訓練用データ（train_set）の最初の100個を表示してみましょう。

データをTheanoの共有変数に格納

Theanoの共有変数を使用すると、GPUのメモリ領域に保存され、学習時（update）に高速に読み書きできます。 共有変数の値を取り出すときに、get_value()を使用するのは、GPUにセットされたデータをCPUのメモリにコピー するためと考えられます。

また、GPUで使用するデータは必ずfloat型で格納しなければなりません。 そのため、ラベルshared_yは、T.cast(shared_y, 'int32')でキャストして返しています。

これでdatasets変数に訓練用、検証用、テスト用のデータがセットされました。

多クラス分類

ニューラルネットの入力$x_i$と出力$u_i$の関係は、重み$W_i$とバイアス$b_i$を使って以下の関係になります。 $$ u_i = \sum_j W_{i,j} x_j + b_i $$

これを分かりやすく説明した図を再度、TensorFlowのMNIST for ML Beginnersから引用します。

この重みWとバイアスbをTheanoの共有変数でセットしているのが、以下の箇所です。

        # 重み行列を初期化
        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)			
		

ソフトマック関数

多クラス分類の活性化関数として使用されるのが、ソフトマック関数（Softmax）で、出力層のクラスiに分類される確率は、以下の様に表されます。 $$ y_i = \frac{exp(u_i)}{\sum_j exp(u_j)} $$

この部分をTheanoのシンボルで表現している箇所が以下の部分です。

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

人工知能に関する断創録の貴重なコメントがあります。

• T.dot(self.W, input)+ self.bは、サンプルxを1つだけ与えたときの確率を計算している
• T.dot(input, self.W)+ self.bは、複数のサンプル（ミニバッチ単位）を行列化してinputに与え、それらのサンプルの確率をまとめて計算している

確率が求まったら最終的に一番高い確率が得られるクラスをy_predにセットします。 $$ y_{pred} = argmax(y) $$ TheanoではT.argmax()にaxis=1を指定することでp_y_given_xの各行（サンプルに相当）において一番確率が高いインデックス（クラスに相当）がまとめて取得できるそうです。

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

コスト関数

分類されたクラスkだけ１で、他は0のベクトルを$d_n$とすると、事後分布は、以下の様になります。 $$ p(d | x) = \prod_{k=1}^K p(C_k | x)^{d_k} $$ 訓練データ${ (x_n, d_n)}(n=1,...,N)$に対するwの尤度は、以下の様になります。 $$ L(W) = \prod_{n=1}^N p(d_n | x_n; W) = \prod_{n=1}^N \prod_{k=1}^K p(c_k|x_n)^{d_{nk}} = \prod_{n=1}^N \prod_{k=1}^K (y_k(x; W))^{d_{nk}} $$

負の対数尤度を誤差関数とすると $$ E(W) = - \sum_{n=1}^N \sum_{k=1}^K d_{nk} log y_k(x_n; W) $$

人工知能に関する断創録では、Sumの代わりに平均meanを使っていることに注意！

        # 式通りに計算するとsumだがmeanの方がよい
        return -T.mean(T.log(self.p_y_given_x)[T.arange(y.shape[0]), y])			
		

モデルの訓練

モデルの訓練には、ミニバッチ確率的勾配降下法（MSGD）を使用しています。 確率的勾配降下法（SGD）はただ1つのサンプルで1回だけパラメータを更新するのに対し、 少数のサンプルをひとまとめにしてその単位で重みを更新します。 このひとまとめにしたサンプル集合をミニバッチ（minbatch）と呼びます。

何番目のミニバッチを使用するかを示すのがシンボルindexです。 コードのfunction()のgivensで定義されている部分がミニバッチの設定箇所です。

        givens={
            x: train_set_x[index * batch_size: (index + 1) * batch_size],
            y: train_set_y[index * batch_size: (index + 1) * batch_size]
        }			
		

確率的勾配降下法で使用するコスト関数とその微分は、負の対数尤度を計算する negative_log_likelihood関数とT.gradの箇所で計算しています。

    # 誤差（コスト）を計算 => 最小化したい
    cost = classifier.negative_log_likelihood(y)
    # コスト関数のtheta = (W,b)の微分を計算
    g_W = T.grad(cost=cost, wrt=classifier.W)
    g_b = T.grad(cost=cost, wrt=classifier.b)			
		

パラメータの更新式は、Wとbの2個をタプルとして指定します。

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

関数の定義は以下の通りです。

• inputs：何番目のミニバッチかを示すインデックスをindexで与え
• outputs：負の対数尤度を計算するcostを渡す
• updates：update式を渡す
• givens：関数を実行する際にシンボルを置き替えるオプションで、ここではじめてxとyに値をセットしています

モデルの評価

モデルの当てはまりの良さをエラー率で評価しています。LogisticRegressionのerrors関数で計算します。

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

Early-Stoppingによる収束判定

検証用データvalid_setのエラー率と訓練用データtrain_setのエラー率から適合（overfitting）を防ぐために、 検証用データのエラー率が増加した時点で学習を打ち切る早期終了（Early-Stopping）というテクニックを使用します。

人工知能に関する断創録に解説されているEarly-Stoppingの説明を引用します。

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

学習結果を使った予測

テスト用データtest_setの最初の25個に対して認識を行い、予測結果と画像を合わせて表示してみます。

1/10の学習データでもそこそこ良い結果が得られています。

重みWのパターン

最後に各クラスを識別した重みWがどのように形に求まったのかみてみましょう。

赤い部分で正で、葵部分が負の値です。各数字の特徴的な部分に赤の部分が見られることが分かります。