共役勾配法によるニューラルネットのパラメータ推定

Courseraの機械学習ネタの続き。前回は、ロジスティック回帰のパラメータ推定（2014/4/15）に共役勾配法（Conjugate Gradient: CG法）を使いました。今回はより複雑なニューラルネット（多層パーセプトロン）のパラメータ推定に共役勾配法を適用してみました。

以前、多層パーセプトロンで手書き数字認識の実験をしたとき（2014/2/1）は、共役勾配法ではなく、勾配降下法（Gradient Descent）を用いてパラメータの更新式を自分で書いていました。

    self.weight1 -= learning_rate * np.dot(delta1.T, x)
    self.weight2 -= learning_rate * np.dot(delta2.T, z)

勾配降下法は、学習率（learning rate）を適切な値に設定しないと収束が遅い、発散するなど欠点があります。今回用いる共役勾配法（2014/4/14）は学習率を自分で設定する必要がない上に勾配降下法より高速というアルゴリズムです。実装は、scipy.optimizeにあるのでそれを使います。今回もMNISTの手書き数字認識を例題にしています。

ニューラルネットのコスト関数

Pythonの共役勾配法のライブラリは、scipy.optimizeモジュールにあります。このモジュールは、最適化したい関数とその導関数を与える必要があります。天下りですが、ニューラルネットのコスト関数は下のようなちょっとやばい式になります。

ここで、mは訓練データ数、Kは出力ユニットの数、 はk番目の出力ユニットに対する教師信号（1-of-K表記）、 は仮説の関数でニューラルネットではフォーワードプロパゲーションによる予測値です。この中の は推定するパラメータで入力層と隠れ層の間の重み（）と隠れ層と出力層の間の重み（）をマージしたものです。

このニューラルネットでは、二乗誤差関数ではなく、交差エントロピー誤差関数を使っています。また、過学習を防ぐために正則化しています。第二項が正則化項で、ニューラルネットのすべての重みを二乗して足すことによって重みの値が巨大にならないように調整しています。400は入力層のユニット数、25は隠れ層のユニット数、10は出力層のユニット数です（下のスクリプトでは少し変えています）。Pythonで書くと下のようなコードになります。ニューラルネットの構造は、入力層、隠れ層、出力層の3層を想定しています。

Coursera Machine Learning Week5 から引用

    J = 0
    for i in range(m):
        xi = X[i, :]
        yi = y[i]
        # forward propagation
        a1 = xi
        z2 = np.dot(Theta1, a1)
        a2 = sigmoid(z2)
        a2 = np.hstack((1, a2))
        z3 = np.dot(Theta2, a2)
        a3 = sigmoid(z3)
        J += sum(-yi * safe_log(a3) - (1 - yi) * safe_log(1 - a3))
    J /= m

    # 正則化項
    temp = 0.0;
    for j in range(hid_size):
        for k in range(1, in_size + 1):  # バイアスに対応する重みは加えない
            temp += Theta1[j, k] ** 2
    for j in range(num_labels):
        for k in range(1, hid_size + 1): # バイアスに対応する重みは加えない
            temp += Theta2[j, k] ** 2
    J += lam / (2.0 * m) * temp;

ニューラルネットでは、このコスト関数を最小化するパラメータ（Theta1とTheta2）を推定します。Theta1は入力層と隠れ層の間の重み、Theta2は隠れ層と出力層の間の重みです。

コスト関数の導関数

scipyの共役勾配法を使うには、最適化したいコスト関数の他にその偏微分が必要になります。ニューラルネットでは、バックプロパゲーションを用いることで効率的にコスト関数の偏微分を計算できます。

ここで、大文字のΔは、各層の誤差δと出力値aから計算できる行列

を全訓練データで足し合わせた値になります。Pythonで書くと下のようになります。

    for i in range(m):
        # forward propagation
        （省略）

        # backpropagation
        delta3 = a3 - yi
        delta2 = np.dot(Theta2.T, delta3) * sigmoidGradient(np.hstack((1, z2)))
        delta2 = delta2[1:]  # バイアス項に対応する要素を除外
        # ベクトル x ベクトル = 行列の演算をしなければならないので
        # 縦ベクトルへのreshapeが必要
        # 行数に-1を指定すると自動的に入る
        delta2 = delta2.reshape((-1, 1))
        delta3 = delta3.reshape((-1, 1))
        a1 = a1.reshape((-1, 1))
        a2 = a2.reshape((-1, 1))
        # 正則化ありのときのデルタの演算
        Theta1_grad += np.dot(delta2, a1.T)
        Theta2_grad += np.dot(delta3, a2.T)

    # 偏微分の正則化項
    Theta1_grad /= m
    Theta1_grad[:, 1:] += (lam / m) * Theta1_grad[:, 1:]
    Theta2_grad /= m
    Theta2_grad[:, 1:] += (lam / m) * Theta2_grad[:, 1:]

コスト関数と偏微分を1つの関数で計算したい

前回、ロジスティック回帰で共役勾配法を使ったとき（2014/4/15）は、コスト関数 J() とその偏微分を返す関数 gradient() を別々に用意してfminc_cg()に渡していました。

    # Conjugate Gradientでパラメータ推定
    theta = optimize.fmin_cg(J, initial_theta, fprime=gradient, args=(X, y))

ニューラルネットの場合、コスト関数を計算するためにフォワードプロパゲーションが必要です。また、偏微分を計算するためにはバックプロパゲーションが必要ですが、その前にフォワードプロパゲーションを行う必要があります。そのため、コスト関数と偏微分を計算する関数を分けるとフォワードプロパゲーションが重複して効率が悪くなります。このようなケースでは、fmin_cg()ではなく、minimize()というより低レベルな関数が使えます（scipy 0.11.0以上）。minimize()を使うとコスト関数と偏微分のリストを返す関数（サンプルではnnCostFunction）だけ渡せばよくなります。

    # Conjugate Gradientでパラメータ推定
    # NNはコスト関数と偏微分の計算が重複するため同じ関数（nnCostFunction）にまとめている
    # この場合、fmin_cgではなくminimizeを使用するとよい
    # minimize()はscipy 0.11.0以上が必要
    res = optimize.minimize(fun=nnCostFunction, x0=initial_nn_params, method="CG", jac=True,
                                  options={'maxiter':20, 'disp':True},
                                  args=(in_size, hid_size, num_labels, X, y, lam))

バッチ学習

前回（2014/2/1）は、確率的勾配降下法（Stochastic Gradient Descent）を用いていました。下のようにランダムに選択した訓練データを1つ使ってパラメータを1回更新する手法です。

    for k in range(epochs):
        # 訓練データからランダムに選択する
        i = np.random.randint(X.shape[0])
        x = X[i]
        (略）
        # パラメータを更新
        self.weight1 -= learning_rate * np.dot(delta1.T, x)
        self.weight2 -= learning_rate * np.dot(delta2.T, z)

今回は、確率的勾配降下法ではなく、バッチ学習を使っています。訓練データをすべて使って誤差を蓄積し、まとめてパラメータを一回更新する手法です。これを収束するまで繰り返し行います。

    for i in range(m):  # mは訓練データ数
        xi = X[i, :]
        yi = y[i]
        # 誤差を蓄積
        Theta1_grad += np.dot(delta2, a1.T)
        Theta2_grad += np.dot(delta3, a2.T)

この方法は最適解にまっすぐ近づきますが、訓練データ数が大きくなると一回の更新に時間がかかるという欠点があります。MNISTのデータ数は70000ですが、バッチ更新だとさすがにきついので訓練データ数は5000にしぼっています。確率的勾配法とバッチ学習の中間的な手法としてミニバッチ学習というのもよく使われるそうです。後で試してみたい。学習実行すると下のように更新を繰り返すことでコスト関数の値が減っていき最適化されることがわかります。

下のような出力が出ます。訓練データの精度は99%、テストデータの精度は92%という結果でした。

Warning: Maximum number of iterations has been exceeded.
         Current function value: 0.336109
         Iterations: 70
         Function evaluations: 201
         Gradient evaluations: 201

*** training set accuracy
[[498   0   0   1   1   0   1   0   0   0]
 [  0 584   0   0   1   0   0   1   0   1]
 [  1   0 473   1   4   1   4   3   1   0]
 [  0   1   1 480   0   3   0   1   3   0]
 [  1   1   0   0 446   0   0   1   0   0]
 [  1   0   1   2   4 465   2   0   3   0]
 [  0   0   0   0   1   0 459   0   0   0]
 [  0   0   3   0   1   1   0 525   1   2]
 [  0   1   1   2   0   0   0   0 505   0]
 [  0   0   0   4   2   0   0   2   0 498]]
             precision    recall  f1-score   support

        0.0       0.99      0.99      0.99       501
        1.0       0.99      0.99      0.99       587
        2.0       0.99      0.97      0.98       488
        3.0       0.98      0.98      0.98       489
        4.0       0.97      0.99      0.98       449
        5.0       0.99      0.97      0.98       478
        6.0       0.98      1.00      0.99       460
        7.0       0.98      0.98      0.98       533
        8.0       0.98      0.99      0.99       509
        9.0       0.99      0.98      0.99       506

avg / total       0.99      0.99      0.99      5000

*** test set accuracy

[[203   0   2   1   1   0   2   1   0   0]
 [  0 226   1   0   0   3   0   0   1   1]
 [  2   2 175   0   2   0   4   1   5   1]
 [  2   1   6 168   0  10   0   2   0   2]
 [  1   3   1   0 183   0   7   0   0   5]
 [  1   0   2   6   2 165   3   2   5   5]
 [  2   2   1   0   3   1 187   0   1   0]
 [  0   0   2   1   2   0   0 196   0   3]
 [  1   4   2   2   0   5   3   0 159   4]
 [  4   0   0   3   8   4   0   5   0 179]]
             precision    recall  f1-score   support

        0.0       0.94      0.97      0.95       210
        1.0       0.95      0.97      0.96       232
        2.0       0.91      0.91      0.91       192
        3.0       0.93      0.88      0.90       191
        4.0       0.91      0.92      0.91       200
        5.0       0.88      0.86      0.87       191
        6.0       0.91      0.95      0.93       197
        7.0       0.95      0.96      0.95       204
        8.0       0.93      0.88      0.91       180
        9.0       0.90      0.88      0.89       203

avg / total       0.92      0.92      0.92      2000

隠れ層の可視化

隠れ層の重みは下のコードで可視化することができます。

    # 隠れユニットを可視化
    displayData(Theta1[:, 1:])

ニューラルネットの隠れ層では数字画像の特徴をこのようにとらえているんですね。元の数字の面影はほとんどないけれど・・・顔画像を認識させたりするともっと面白い画像が得られるのかな？顔画像認識も今度試してみたい。

完成したスクリプト

https://github.com/sylvan5/PRML/blob/master/ch5/mlp_cg.py

#coding: utf-8
import numpy as np
from scipy import optimize
from matplotlib import pyplot
from sklearn.datasets import load_digits, fetch_mldata
from sklearn.preprocessing import LabelBinarizer
from sklearn.metrics import confusion_matrix, classification_report

"""
正則化多層パーセプトロンによる手書き文字認識
・バッチ処理
・共役勾配法（Conjugate Gradient Method）によるパラメータ最適化
"""

def displayData(X):
    """
    データからランダムに100サンプル選んで可視化
    画像データは8x8ピクセルを仮定
    """
    # ランダムに100サンプル選ぶ
    sel = np.random.permutation(X.shape[0])
    sel = sel[:100]
    X = X[sel, :]
    for index, data in enumerate(X):
        pyplot.subplot(10, 10, index + 1)
        pyplot.axis('off')
        image = data.reshape((28, 28))
        pyplot.imshow(image, cmap=pyplot.cm.gray_r,
                     interpolation='nearest')
    pyplot.show()

def randInitializeWeights(L_in, L_out):
    """
    (-epsilon_init, +epsilon_init) の範囲で
    重みをランダムに初期化した重み行列を返す
    """
    # 入力となる層にはバイアス項が入るので+1が必要なので注意
    epsilon_init = 0.12
    W = np.random.rand(L_out, 1 + L_in) * 2 * epsilon_init - epsilon_init
    return W

def sigmoid(z):
    return 1.0 / (1 + np.exp(-z))

def sigmoidGradient(z):
    return sigmoid(z) * (1 - sigmoid(z))

def safe_log(x, minval=0.0000000001):
    return np.log(x.clip(min=minval))

def nnCostFunction(nn_params, *args):
    """NNのコスト関数とその偏微分を求める"""
    in_size, hid_size, num_labels, X, y, lam = args

    # ニューラルネットの全パラメータを行列形式に復元
    Theta1 = nn_params[0:(in_size + 1) * hid_size].reshape((hid_size, in_size + 1))
    Theta2 = nn_params[(in_size + 1) * hid_size:].reshape((num_labels, hid_size + 1))

    # パラメータの偏微分
    Theta1_grad = np.zeros(Theta1.shape)
    Theta2_grad = np.zeros(Theta2.shape)

    # 訓練データ数
    m = X.shape[0]

    # 訓練データの1列目にバイアス項に対応する1を追加
    X = np.hstack((np.ones((m, 1)), X))

    # 教師ラベルを1-of-K表記に変換
    lb = LabelBinarizer()
    lb.fit(y)
    y = lb.transform(y)

    J = 0
    for i in range(m):
        xi = X[i, :]
        yi = y[i]
        # forward propagation
        a1 = xi
        z2 = np.dot(Theta1, a1)
        a2 = sigmoid(z2)
        a2 = np.hstack((1, a2))
        z3 = np.dot(Theta2, a2)
        a3 = sigmoid(z3)
        J += sum(-yi * safe_log(a3) - (1 - yi) * safe_log(1 - a3))
        # backpropagation
        delta3 = a3 - yi
        delta2 = np.dot(Theta2.T, delta3) * sigmoidGradient(np.hstack((1, z2)))
        delta2 = delta2[1:]  # バイアス項に対応する要素を除外
        # ベクトル x ベクトル = 行列の演算をしなければならないので
        # 縦ベクトルへのreshapeが必要
        # 行数に-1を指定すると自動的に入る
        delta2 = delta2.reshape((-1, 1))
        delta3 = delta3.reshape((-1, 1))
        a1 = a1.reshape((-1, 1))
        a2 = a2.reshape((-1, 1))
        # 正則化ありのときのデルタの演算
        Theta1_grad += np.dot(delta2, a1.T)
        Theta2_grad += np.dot(delta3, a2.T)
    J /= m

    # 正則化項
    temp = 0.0;
    for j in range(hid_size):
        for k in range(1, in_size + 1):  # バイアスに対応する重みは加えない
            temp += Theta1[j, k] ** 2
    for j in range(num_labels):
        for k in range(1, hid_size + 1): # バイアスに対応する重みは加えない
            temp += Theta2[j, k] ** 2
    J += lam / (2.0 * m) * temp;

    # 偏微分の正則化項
    Theta1_grad /= m
    Theta1_grad[:, 1:] += (lam / m) * Theta1_grad[:, 1:]
    Theta2_grad /= m
    Theta2_grad[:, 1:] += (lam / m) * Theta2_grad[:, 1:]

    # ベクトルに変換
    grad = np.hstack((np.ravel(Theta1_grad), np.ravel(Theta2_grad)))

    print "J =", J
    return J, grad

def predict(Theta1, Theta2, X):
    m = X.shape[0]
    num_labels = Theta2.shape[0]
    # forward propagation
    X = np.hstack((np.ones((m, 1)), X))
    h1 = sigmoid(np.dot(X, Theta1.T))

    h1 = np.hstack((np.ones((m, 1)), h1))
    h2 = sigmoid(np.dot(h1, Theta2.T))

    return np.argmax(h2, axis=1)

if __name__ == "__main__":
    # 入力層のユニット数
    in_size = 28 * 28

    # 隠れ層のユニット数
    hid_size = 25

    # 出力層のユニット数（＝ラベル数）
    num_labels = 10

    # 訓練データをロード
    digits = fetch_mldata('MNIST original', data_home=".")
    X = digits.data
    X = X.astype(np.float64)
    X /= X.max()
    y = digits.target

    # データをシャッフル
    p = np.random.permutation(len(X))
    X = X[p, :]
    y = y[p]

    # 最初の5000サンプルを選択
    X_train = X[:5000, :]
    y_train = y[:5000]

    # データを可視化
    displayData(X_train)

    # パラメータをランダムに初期化
    initial_Theta1 = randInitializeWeights(in_size, hid_size)
    initial_Theta2 = randInitializeWeights(hid_size, num_labels)

    # パラメータをベクトルにフラット化
    initial_nn_params = np.hstack((np.ravel(initial_Theta1), np.ravel(initial_Theta2)))

    # 正則化係数
    lam = 1.0

    # 初期状態のコストを計算
    J, grad = nnCostFunction(initial_nn_params, in_size, hid_size, num_labels, X_train, y_train, lam)

    # Conjugate Gradientでパラメータ推定
    # NNはコスト関数と偏微分の計算が重複するため同じ関数（nnCostFunction）にまとめている
    # この場合、fmin_cgではなくminimizeを使用するとよい
    # minimize()はscipy 0.11.0以上が必要
    res = optimize.minimize(fun=nnCostFunction, x0=initial_nn_params, method="CG", jac=True,
                                  options={'maxiter':70, 'disp':True},
                                  args=(in_size, hid_size, num_labels, X_train, y_train, lam))
    nn_params = res.x

    # パラメータを分解
    Theta1 = nn_params[0:(in_size + 1) * hid_size].reshape((hid_size, in_size + 1))
    Theta2 = nn_params[(in_size + 1) * hid_size:].reshape((num_labels, hid_size + 1))

    # 隠れユニットを可視化
    displayData(Theta1[:, 1:])

    # 訓練データのラベルを予測して精度を求める
    pred = predict(Theta1, Theta2, X_train)
    print "*** training set accuracy"
    print confusion_matrix(y_train, pred)
    print classification_report(y_train, pred)

    # 訓練データとかぶらない範囲からテストデータを選んで精度を求める
    print "*** test set accuracy"
    X_test = X[10000:12000, :]
    y_test = y[10000:12000]
    pred = predict(Theta1, Theta2, X_test)
    print len(pred)
    print confusion_matrix(y_test, pred)
    print classification_report(y_test, pred)