ニューラルネットワークの入門の内容です。復習も兼ねて演習問題をPythonでも再現しました。Scipyの最小化関数があまり頼もしくなく、Octaveと同じ結果を出すのに苦労しました。
前回まで1つ目を組み込みで実装→1つ目を自分で実装→2つ目の組み込み→2つ目を自分でとやっていましたが、見づらかったので組み込み(1)→組み込み(2)→セルフ(1)→セルフ(2)の順番で書きます。

追記 gistからソースファイルをダウンロードできます

多クラス、自分で実装　https://gist.github.com/koshian2/62514efd7c0d7b85f4a68fd5f935b230
多クラス、組み込み　https://gist.github.com/koshian2/4f3f78a3839c00ad93dee54005f82b28
ニューラルネットワーク、自分で実装　https://gist.github.com/koshian2/4d0daeddc85f2514083ac71a3affd4c6
ニューラルネットワーク、組み込み　https://gist.github.com/koshian2/8c558accf13aa03eccaa8481666026cf

これまでの目次

##◎組み込みの多クラス分類（隠れ層なしのニューラルネットワーク）
ロジスティック回帰の拡張で、MNISTの10個の手書き数字を判別します。10個のクラスに分ける多クラス分類ですが、言い方を変えれば隠れ層のないニューラルネットワークとなります。最も基本的なニューラルネットワークと言ってもいいでしょう。

まずは組み込みで実装してみます。多クラス分類になってもScikit-learnの組み込み```LogisticRegression``で対応できます。

import numpy as np
from scipy.io import loadmat
from sklearn.linear_model import LogisticRegression

# データの読み込み
def load_data1():
    data = loadmat("ex3data1")
    # yが元データだと5000x1の行列なので、ベクトルに変換する
    return np.array(data['X']), np.ravel(np.array(data['y']))

X_data, y = load_data1()

# ロジスティック回帰
regr = LogisticRegression(multi_class="ovr", solver="newton-cg") #multi_class="ovr"でOneVsRest（1対多クラス分類になる）
print("X : shape =", X_data.shape)
print("y : shape =", y.shape)
print()
regr.fit(X_data, y)
print("切片 : shape =", regr.intercept_.shape)
print(regr.intercept_)
print("係数 : shape =", regr.coef_.shape)
print(regr.coef_)
print("精度 =",regr.score(X_data, y) * 100)

出力

X : shape = (5000, 400)
y : shape = (5000,)

切片 : shape = (10,)
[-2.38368769 -3.18302619 -4.79746133 -2.01355322  0.17550751 -3.14856414
 -1.90394385 -7.98717568 -4.57254937 -5.40568903]
係数 : shape = (10, 400)
[[  0.00000000e+00   0.00000000e+00   3.53650804e-05 ...,   1.30431448e-03
   -7.29175668e-10   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00  -2.80610230e-05 ...,   4.46082357e-03
   -5.08589273e-04   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00  -8.08986150e-06 ...,  -2.87050197e-05
   -2.47483300e-07   0.00000000e+00]
 ...,
 [  0.00000000e+00   0.00000000e+00  -4.78980835e-06 ...,  -8.94785180e-05
    7.21469722e-06   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00  -5.50664592e-07 ...,  -1.33539649e-03
    9.98597489e-05   0.00000000e+00]
 [  0.00000000e+00   0.00000000e+00  -2.35768210e-09 ...,  -1.16635470e-04
    7.88299299e-06   0.00000000e+00]]
精度 = 94.46

数秒で出てきます。Octaveでも10秒ぐらいかかって精度95.0だったので、これはかなりコスパいいと思います。自分で実装した例では10分近くかかって94.7だったので、自分で実装するよりかは組み込み推奨です。出力の切片と係数もy=1…10順に10セットできているので想定通りです。loadmatについては後ほど。

##◎sklearn.neural_networkの組み込みでニューラルネットワーク（隠れ層あり）
scikit-learnで用意されているsklearn.neural_network.MLPClassifierを使うと、ニューラルネットワークによる分類器を自分で実装することなく使うことができます。自分で実装する場合は（問題の例では）、第1層、第2層の係数が既知のものとして実装していますが、今回はより実践的にこれらの係数が未知のものとして実装してみます。
与えるのは隠れ層のノードの数だけです。「こういうネットワーク作って」という設計だけ与えて、あとは勝手にやらせるイメージです。

import numpy as np

from scipy.io import loadmat
from sklearn.neural_network import MLPClassifier

# データの読み込み
def load_data1():
    data = loadmat("ex3data1")
    return np.array(data['X']), np.ravel(np.array(data['y']))

X_data, y = load_data1()
m = len(X_data[:, 1])

# 分類器
clf = MLPClassifier(hidden_layer_sizes=(26, ), solver="adam", random_state=114514, max_iter=10000)
clf.fit(X_data, y)

# 切片と係数のデータ型はlist
for idx, (intercept, coef) in enumerate(zip(clf.intercepts_, clf.coefs_)):
    print("第",idx,"層")
    print(" 切片 =", intercept.shape)
    print(" 係数 =", coef.shape)
    #print(coef)
print("")
print("訓練データの精度 :",clf.score(X_data, y)*100)

これだけで終わりです。solverはsgd（確率的勾配降下法）よりデフォルトのadamのほうが精度が良かったです。自分の環境では20秒ぐらいで答えが出ました。

出力

第 0 層
 切片 = (26,)
 係数 = (400, 26)
第 1 層
 切片 = (10,)
 係数 = (26, 10)

訓練データの精度 : 100.0

「100%当てることができますよ（ﾄﾞﾔｧ）」ということらしい。バイアスノードの扱いが変わるので若干次元は違いますが、切片・係数は期待通りの結果になりました。本当にできてるの？ということで係数を見てみます。

追記：正則化しないで学習させると、訓練データを100%近く当てられるサンプルのようです。過学習（オーバーフィッティング）の可能性あり。

第 0 層
 切片 = (26,)
 係数 = (400, 26)
[[  9.13465305e-136  -3.60318571e-136  -4.54186384e-140 ...,
    1.14727123e-121   8.78660405e-141   7.27236139e-141]
 [ -2.15214607e-123   4.47523327e-143  -1.70900500e-138 ...,
    6.17161518e-140   3.19286303e-128  -3.19701789e-142]
 [  2.31728319e-005  -9.11762043e-005   1.55968062e-004 ...,
   -8.32028382e-005  -1.13577123e-004  -8.33986978e-005]
 ...,
 [ -6.18642262e-003  -2.36230505e-002  -1.02453868e-002 ...,
   -4.61198105e-003  -1.04477710e-002  -5.23835175e-004]
 [  2.39883965e-005   1.91689394e-005   1.52643880e-004 ...,
    3.73256334e-005  -6.93261934e-005  -2.55899623e-005]
 [ -2.42048502e-123   8.73562267e-132  -2.33770448e-138 ...,
    5.94129597e-125  -1.51688638e-140   3.46603214e-132]]

これは入力層の係数です。個々の係数の数字の意味は、人間には正直よくわかりませんが、これらの係数から計算すると訓練データを100%当てることができます。ニューラルネットワークしゅごい。

##・多クラス分類を自分で実装する
ロジスティック回帰に戻ります。前回までデータをハードコーディングしていましたが、今回は7MBもあるので断念。ファイルを同じディレクトリに置いて読み込んでください。

###1.データの読み込み

import numpy as np
import matplotlib.pyplot as plt

from scipy.io import loadmat

# データの読み込み
def load_data1():
    data = loadmat("ex3data1")
    # yが元データだと5000x1の行列なので、ベクトルに変換する
    return np.array(data['X']), np.ravel(np.array(data['y']))

X_data, y = load_data1()
m = len(X_data[:, 1])

元がMatlabのデータファイルですが、scipy.io.loadmat()で読み込むことができます。Xは素直に行列なのですが、yはベクトルなのに5000x1の行列に変換されてしまいコスト関数でハマるので、ravel()でベクトルに変換しておきます（こういうのが面倒くさい）。

###2.データのプロット

# ランダムに100個を画像で表示
np.random.seed(114514)# ここをコメントアウトすると再現性はなくなる
sel = np.arange(m)
np.random.shuffle(sel)
sel = sel[:100]
fig = plt.figure(figsize = (10, 10))
fig.subplots_adjust(hspace=0.05, wspace=0.05)
for i in range(100):
    ax = fig.add_subplot(10, 10, i + 1, xticks=[], yticks=[])
    ax.imshow(X_data[sel[i]].reshape((20, 20)).T, cmap='gray')
plt.show()

ランダムでシャッフルしていますが、再現性をもたせるためにシードを固定しています。reshapeすると横に回転した数字が出てくるので、転置して戻しておきます。おなじみの数字が出てきます。

###3.コスト関数の定義

# シグモイド関数
def sigmoid(z):
    return 1 / (1 + np.exp(-z))

# ロジスティック回帰のコスト関数
def lr_cost_function(theta, X, y, lambda_):
    m = len(y)
    h_theta = sigmoid(np.dot(X, theta))
    J = np.sum(-y * np.log(h_theta) - (1 - y) * np.log(1 - h_theta)) / m + lambda_ / 2 / m * np.sum(theta[1:] ** 2)
    # θ0を正則化しないようにする
    grad = np.dot(X.T, h_theta - y) / m
    temp = theta[:]
    temp[0] = 0
    grad += lambda_ / m * temp
    return J, grad

$\theta_0$を正則化しないようにするのが面倒くさい。データのyのほうを(5000, 1)の行列から(5000, )のベクトルになおしておかないと、本番でハマリます。

###4.コスト関数のテスト

# コスト関数のテスト
print("Testing lrCostFunction() with regularization")
theta_t = np.array([-2, -1, 1, 2])
# orderの設定をしないとreshapeのデフォルト設定はOctaveと逆なので注意
X_t = np.c_[np.ones(5), np.arange(1, 16).reshape(5, 3, order='F') / 10] 
y_t = (np.array([1, 0, 1, 0, 1]) >= 0.5).astype(int)
lambda_t = 3
J, grad = lr_cost_function(theta_t, X_t, y_t, lambda_t)
print("Cost:", J)
print("Expected cost: 2.534819")
print("Gradients:")
print(grad)
print("Expected gradients:")
print(" 0.146561\n -0.548558\n 0.724722\n 1.398003\n")

ここにも罠があって、reshapeのデフォルトのオーダーはOctaveとPython/Numpyでは逆です。例えば、1～4の数字を2x2行列に変換したとき、Octaveでは縦→横の順で並びますが、

Octave

>> reshape([1:4],2,2)
ans =

   1   3
   2   4

Pythonの場合は、横→縦になります。（これに気づかなくて1時間ぐらいハマった）。

Python

>>> np.arange(1,5).reshape(2,2)
array([[1, 2],
       [3, 4]])

Octaveと同じreshapeをやりたければ、Pythonではorder='F'を指定しないといけません。こういう細かな言語仕様の違いが面倒くさい。上記のテストの出力は以下のようになります。

出力

Testing lrCostFunction() with regularization
Cost: 2.53481939611
Expected cost: 2.534819
Gradients:
[ 0.14656137 -0.54855841  0.72472227  1.39800296]
Expected gradients:
 0.146561
 -0.548558
 0.724722
 1.398003

##5.最急降下法の実装
Octaveのときはfmincgを使っていたのですが、Scipyのfmin_cgやfmin_ncgを使うとコスト値の浮動小数点の桁の問題で、反復が少ない段階で終わってしまうというちょっといただけない現象があります。
本来はこのモデルは精度95%近く出るのですが、この現象のせいで精度は90%までしか上がりません。なので、自分で最急降下法を実装したほうが精度は上がりました。もし他に良い組み込みの最適化関数があったら教えてください。

# 最急降下法（組み込みが遅いので自分で実装）
def gradient_descent(initial_theta, X, y, lambda_, eta, maxiter = 10000, tol=1e-3):
    theta_before = initial_theta
    for i in range(maxiter):
        J, grad = lr_cost_function(theta_before, X, y, lambda_)
        theta = theta_before - eta * grad
        norm = np.linalg.norm(theta - theta_before)
        if(i%100==0) : print("i =",i,", norm =", norm, "J =",J)
        if norm < tol:
            print("収束完了", i)
            break
        theta_before = theta
    return theta

maxiterが最大反復回数、その前にθのL2ノルムがtol以下になったら収束とみなして早期終了しています。etaは学習率ですが、この例では大きめの値（η＝1）を代入してOKです。収束まで相当時間かかるので。

#6.学習

# One-vs-allの訓練
def one_vs_all(X, y, num_labels, lambda_):
    m = X.shape[0]
    n = X.shape[1]
    all_theta = np.zeros((num_labels, n + 1))
    X = np.c_[np.ones(m), X]
    for i in range(num_labels):
        print("One vs all :", i+1, "/", num_labels)
        initial_theta = np.zeros(n+1)
        y_param = y == i+1
        theta = gradient_descent(initial_theta, X, y_param, lambda_, 1)
        all_theta[i, :] = theta
    return all_theta

num_labels = 10
lambda_ = 0.1
all_theta = one_vs_all(X_data, y, num_labels, lambda_)

自分の環境では10分ぐらいで計算が終わりました。

###7.精度

# 予測
def predict_one_vs_all(all_theta, X):
    m = X.shape[0]
    num_labels = all_theta.shape[0]
    XX = np.c_[np.ones(m), X]
    pred_array = sigmoid(np.dot(XX, all_theta.T))
    print(pred_array)
    p = np.argmax(pred_array, axis=1)+1 #行単位で集計
    return p

pred = predict_one_vs_all(all_theta, X_data)
print(np.bincount(pred))
print("Training Set Accuracy: ", np.mean(pred == y) * 100)

出力

[  0 512 489 485 503 483 501 497 510 511 509]
Training Set Accuracy:  94.78

精度は94.8%となりました。テキストによると94.9%出るとのことなので、ほぼ期待通りですね。

##・ニューラルネットワーク（Forward propagation）を自分で実装
これまでは隠れ層なしの多クラス分類をやっていましたが、入力と出力の間に隠れ層を1つ用意し、その層にはノードが26個（バイアスノードを抜くと25個）あるモデルを使います。そのノードが何を意味をするのかはよくわかりません。以下のような伝播をします。

入力層：$x=a^{(1)}=$5000x400次元、$\Theta^{(1)}=$25x401次元（バイアスノードが入るので列が1次元増えます）
隠れ層：$a^{(2)}=$5000x25次元、$\Theta^{(2)}=$10x26次元
出力層：$a^{(3)}=h_\theta(x)=$5000x10次元

この問題では、$\Theta^{(1)}$と$\Theta^{(2)}$が既に与えられているものとします。

###1.データの読み込み

import numpy as np
from scipy.io import loadmat

# データの読み込み
def load_data1():
    data = loadmat("ex3data1")
    return np.array(data['X']), np.ravel(np.array(data['y']))

X_data, y = load_data1()

# 計算済みの係数を読み込み
def load_weights():
    data = loadmat("ex3weights")
    return np.array(data['Theta1']), np.array(data['Theta2'])

Theta1, Theta2 = load_weights()

###2.Forward propagation

# Forward propagation
def predict(Theta1, Theta2, X):
    XX = np.c_[np.ones(X.shape[0]), X]
    Z2 = 1 / (1 + np.exp(-np.dot(XX, Theta1.T)))
    XX = np.c_[np.ones(Z2.shape[0]), Z2]
    Z3 = 1 / (1 + np.exp(-np.dot(XX, Theta2.T)))
    return np.argmax(Z3, axis=1)+1

多クラス分類を2回行っているだけです。行列の内積だらけで何やってるのか頭がこんがらがってきますが、「m×nの行列とn×pの行列の内積は、m×pの行列」という基本公式に戻ればどんなプログラムを書けばいいかおおよそ見当がつくはずです。

参考：行列の積　ＡＢの定義

Forward propagationというののいい訳語がわかりませんでしたが、前方誤差伝播ということでしょうか。

###3.精度

# 精度
pred = predict(Theta1, Theta2, X_data)
print("Training Set Accuracy: ", np.mean(pred == y) * 100)

出力

Training Set Accuracy:  97.52

隠れ層がなかった場合の精度が95前後だったのに対し、隠れ層を入れると97.5まで上昇しました。隠れ層を入れると精度が上昇することが確認できました。

以上です。一見難しそうなニューラルネットワークも組み込みの関数を使えば簡単にできますね。
これらの処理を自分で実装する必要はなくて、組み込みの関数がおおよそどんな処理をやっているかぐらいの理解ができれば十分だと思いますよ（逆に言うとそこを理解しないで精度の高い予測ができる魔法の道具として使うのは個人的にはよくないと思う）。

次回に進む
Coursera Machine LearningをPythonで実装 - [Week5]ニューラルネットワーク(2)