50行くらいで実装するニューラルネットワーク
これ何
classもdef文も書かずに, 隠れ層1つのニューラルネットワークを一から実装します.
実装するに当たって、順伝播と逆伝播の考え方は
Coureseraのdeep learning specializationが大変役立ちました*1。
サイトとしてはhttp://cs231n.github.io/neural-networks-case-study/を参考にしました。ただ、逆伝播について、より数学的に厳密に実装したつもりです。
実装したニューラルネットワークは下図のようになります。
上付きの[]数字は何番目の層か、ユニットの個数は特徴量の数とみなせます。
斜め文字W, bは行列です。
特徴は
- 出力層ユニット1の2値分類を扱う
- 最も単純な2層からなるニューラルネットワークを扱う
- 活性化関数はシグモイド関数
- 入力層ユニット数は2(つまり特徴量2), 隠れ層ユニット数は4
なお、ニューラルネットワークにおける層の数は, (隠れ層の数) + (出力層の数) としてカウントされ、入力層は含まれません。
2層以上で構成されるモデルは多層パーセプトロン(Multi Layer Perseptron: MLP)と呼ばれ、
深層学習(deep learning)とは、一般的に3層以上で構成されるMLPモデルを用いた機械学習のことを指します。
ニューラルネットワークは何をしているのか?
実は、そこまで難しくありません。
ニューラルネットワークは、ロジスティック回帰です。
流石に嘘です。でも、とても近いです。
ロジスティック回帰式は、線形回帰式を非線形活性化関数であるシグモイド関数の中に入れてあげることで、0以上1未満の2値分類を表すことができるのでした。
それと同じで、
ニューラルネットワーク = 線形結合 + 活性化関数による非線形変換
と言えます。
順伝播のセクションをみるともうちょっとイメージ湧くはずです。
さて、ニューラルネットワークの1回(1 epoch)あたりの学習の流れは以下の通りです。
- 重みを初期化する
- 順伝播により出力層に重みを伝え、確率を算出する
- 逆伝播により重みの差(勾配)を計算する
- 重みを更新する
2番目と3番目、一度は聞いたことがあるであろう「順伝播」と「逆伝播」を説明します
順伝播(forward propagation)とは
入力層から出力層に向けた伝播を表します。こちらはイメージしやすく、線型変換 + 非線形変換を経て、最後に確率を計算する流れです。
単層パーセプトロンで見てみます。
これは、ロジスティック回帰と全く同じ計算となっています。
2値分類では、出力層をsigmoid関数とすることで、0以上1未満の値を出力することができます(つまり活性化関数の出力aがそのまま予測値になる)。の場合、ラベル1となる確率が0.8ということです。
そして、予測値aと、真値yからコスト関数L(a, y)が求めます。
同じことが3層以上のニューラルネットワークでも言えます。
逆伝播(back propagation)とは
逆伝播はイメージつきにくいかもしれません。
ポイントは
重みを更新するために、誤差に基づく重みの差分(勾配:)を求めたい
ということです。差分がわかれば、
のようにして、初期値Wから重みを更新をすることができます。ここでのは学習率というやつです。
じゃあ、どうやって求めるのか。
「合成関数の連鎖積」を使います。
ここは数学的知識がある前提で進めます...
図で見るとこんな感じになります。
この例では、を求めたいとします。
右から見てみてください。
- コスト関数L(a, y)がわかっているので、Lのaに対する偏微分daが求まる
- 次に偏微分dZを求めたいが、を直接求めるのではなく、合成関数の連鎖積で(青とピンクに)まで分解し、既に既知のdaとを用いて計算する。
- 同様に、dw1も、合成関数の連鎖積で、既知のdZとを用いて計算する。
ここで用いている偏微分da およびは、
コスト関数として2値分類で用いられる交差エントロピー(1サンプルあたり)
dZまで実際に計算してみると、a - yと、こんなにスッキリした値になるんですね。
ぜひ、実際に手を動かして確認してみてください。
「逆伝播」と言っている通り、1→3、順方向とは逆の流れで重みが伝播しているのがわかります。
後はひたすら、
- 順伝播で確率aを計算
- コスト関数を計算
- 逆伝播により勾配を計算
- 重みを更新
これをループで繰り返します。
実装
# ライブラリのインポート import numpy as np import matplotlib.pyplot as plt %matplotlib inline import sklearn from sklearn.linear_model import LogisticRegressionCV from sklearn.metrics import accuracy_score # データの生成 from sklearn.datasets import make_moons X, y_ = make_moons(n_samples=100, random_state=0, noise=.3) y = y_.reshape((y_.shape[0], -1)) print(X.shape) # (100, 2) print(y.shape) # (100, 1) # 描写 plt.scatter(X[y_==0, 0], X[y_==0, 1], color='r', marker='^', alpha=0.5) plt.scatter(X[y_==1, 0], X[y_==1, 1],color='b', marker='o', alpha=0.5) plt.show()
データ数100, ラベルは0か1, 特徴量すなわち次元数は2のデータです。
ちなみにsklearnで色々なデータが入っててこの半月状のデータもあることは知ってましたが、引数にnoise
も設定できるんですね。便利。
まずは古典的なロジスティック回帰をみてみます
# Train the logistic regression classifier clf = LogisticRegressionCV() clf.fit(X, y_) # Print accuracy y_pred_lr = clf.predict(X) score = accuracy_score(y_, y_pred_lr) print ('Accuracy of logistic regression: %d ' % (score * 100) + '% ' + "(percentage of correctly labelled datapoints)")
色ぬりもしてみるとこんな感じになりました
必要な関数はシグモイド関数a = σ(Z)でこれは作っておきます
def sigmoid(x): return 1 / (1 + np.exp(-x))
ここからは、
- ニューラルネットワークの構造の定義
- 重みの初期化
- ループで重みを更新
です。
# back propagation
における勾配は全て、合成関数の連鎖積を計算しています。
コスト関数は、mサンプル全体だと次のようになります。
ここでのa[2]は実装におけるy_outと対応しています。
ここで、1サンプルあたりの、実装されている勾配の数式をまとめると次のようになります。
ちょっとわかりにくいですが、全ての値は行列だと思ってください。
(図2における2. a-yに対応)
(*は単純な(要素同士の)掛け算, da/dz = a_1 * (1-a_1)を利用)
注意点
- 図2では単層だったためZやWが一つしかないが、今回は2層なのでdW1やdb1も求める必要がある。ただ、同じ(線型結合+sigmoid)が加わっただけで逆伝播も同じように計算できる。
- サンプル全体のコストを求めているので、サンプル数(m, コードではn_samples)で割り一つの重みの変化に対応させる。勾配も。
- 行列演算上そのまま内積を求めるのではダメで、転置する必要がある
# 1. Define neural network stracture n_x = X.shape[1] # 特徴量の数がinput n_h = 4 # 適当, 今回は4ユニット n_y = 1 # 2クラス分類は1つの確率で決まるので、1ユニット num_iteration = 15000 # 何回計算するか n_samples = X.shape[0] # サンプル数 # 2. Initialize the model's parameters W1 = np.random.randn(n_x, n_h) * 0.01 W2 = np.random.randn(n_h, n_y) * 0.01 b1 = np.zeros((1, n_h)) b2 = np.zeros((1, n_y)) # 3. loop for i in range(num_iteration): # forward propagation Z1 = np.dot(X, W1) + b1 # 線型結合(隠れ層) A1 = sigmoid(Z1) # 活性化関数(隠れ層) Z2 = np.dot(A1, W2) + b2 # 線型結合(出力層) assert Z2.shape == (n_samples, 1) y_out = sigmoid(Z2) # probability(0~1), y_out = a2 # cost function cost = -1 / n_samples * np.sum((np.log(y_out)*y) + (np.log(1-y_out))*(1-y)) assert(isinstance(cost, float)) # back propagation dZ2 = y_out - y dW2 = np.dot(A1.T, dZ2) / n_samples assert dW2.shape == W2.shape db2 = np.sum(dZ2, axis=0, keepdims=True) / n_samples dZ1 = np.dot(dZ2, W2.T) * (A1 - np.power(A1, 2)) dW1 = np.dot(X.T, dZ1) / n_samples db1 = np.sum(dZ1, axis=0, keepdims=True) / n_samples # update parameters lr = 1e-0 W1 -= lr * dW1 W2 -= lr * dW2 b1 -= lr * db1 b2 -= b2 * db2 if i % 1000 == 0: print('iteration %d:loss %f' % (i, cost))
15000回目のイテレーションで出力された結果y_out
をみてみる
# 出力された確率が0.5以上だったらラベル=1, それ以下だったらラベル=0 y_pred = np.where(y_out < 0.5, 0, 1) print('training accuracy: %.2f' % (np.mean(y_pred == y)))
training accuracy: 0.97
こんな感じになりました。
最後に
簡単そうに見えますが、全然大変でした。
ひたすらassert
文で自分の期待する行列のshapeと一致するか、手探りで進んでいった感じです...
コツとしては、ひたすら「行列の形を意識すること」だと思います。
- ユニット数は次元数と同じとみなせるので、データの形は基本(n_samples, ユニット数)なはず
- コスト関数はfloatなはず
など、色々なポイントがあるので、逐一チェックするといいと思います。
1から理解しながら実装するのはかなり大変ですが、これをやると例えば
- 初期値や学習率をどう決めるのがいいのかが(解説読んですんなり)わかるようになる
- sigmoid関数以外のtanhやReLUが良いのはなぜかわかるようになる
- CNNや最新のディープラーニング技術も(雰囲気)わかるようになる
- Inputするデータの形をめちゃくちゃ意識するようになる
- 行列演算ちょっとできるようになる
などかなりメリットがあるなと思います。
ディープラーニングの闇の魔術に飲み込まれない程度には理解したいところです。
最後まで読んでくださってありがとうございました〜
*1:自分は1ヶ月かけて5つのコースを終えることができました.