深層学習／ニューラルネットワークをスクラッチで組んでみる

#1.はじめに
　今回は、順伝播、誤差逆伝播を理解する為に、実際に動くニューラルネットワークをスクラッチで実装してみます。データセットはMNISTを使います。

#2.ニューラルネットワークの仕様
　入力層28×28=784個、中間層10個、出力層2個のニューラルネットワークとし、活性化関数・誤差関数は sigmoid 、最適化手法は勾配降下法とします。

　データセットは、MNISTから「１」と「７」を抽出したものを使い、2値分類を行います。
![スクリーンショット 2020-03-16 09.00.03.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/6b6dd35c-7e77-4708-b3a8-48a843dd259d.png)

#3.順伝播
　まず、$a^0_0$を計算すると、$x^0_0$から$x^0_{783}$まで784個の入力が入っていて、それぞれ重み$w^0_{00}$から$w^0_{783,0}$が掛かっているので、
![スクリーンショット 2020-03-16 09.12.07.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/651e8c09-e543-6eb0-cdd5-439f47f76730.png)

　行列で表すと、$a^0_0$から$a^0_9$の全ての計算を簡単に表せます。
![スクリーンショット 2020-03-16 09.19.28.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/306373e9-f4b2-aaba-3357-a611391f4d59.png)

　$x^1_0$から$x^1_9$は、$a^0_0$から$a^0_9$を活性化関数sigmoidに通したものになるので、
![スクリーンショット 2020-03-16 09.21.38.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/7faef362-2768-d568-0c14-202ba3b212a4.png)

　次に、$a^1_0$を計算すると、$x^1_0$から$x^1_9$まで10個の入力が入っていて、それぞれ重み$w^1_{00}$から$w^1_{90}$が掛かってるので、
![スクリーンショット 2020-03-16 09.25.48.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/bf41d84c-8df5-0951-088e-6ca227bee6d8.png)

　先程同様、$a^1_0$, $a^1_1$を行列で表すと、
![スクリーンショット 2020-03-16 18.33.38.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/5fe0cdf4-3423-ba48-0fec-91d9844c29b5.png)

　最後に$y^0$と$y^1$は、
![スクリーンショット 2020-03-16 19.29.15.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/c5c5ccce-c9ba-d66b-3b88-aa46c896ae23.png)

この様に、順伝播は行列の内積や足し算等で簡単に行えます。

#4.誤差逆伝播（中間層から出力層)
　まず、中間層から出力層の重みとバイアスの更新です。

　重みwの更新式は、$w=w-\eta*\frac{\partial E}{\partial w}$で表せます。ここで、$\eta$は学習率、$\frac{\partial E}{\partial w}$は誤差Eを重みwで微分したものです。

　$\frac{\partial E}{\partial w}$について具体的な例を上げ、それを元に実装するための一般式で表してみます。まず、中間層から出力層の重みです。

![スクリーンショット 2020-03-15 19.10.21.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/68edbb1a-e583-0601-6006-dc91232c70d0.png)
　重み$w^1_{00}$を更新するために、$\frac{\partial E^0}{\partial w^1_{00}}$を求めます。微分の連鎖率から、
![スクリーンショット 2020-03-15 19.14.30.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/9509e30a-9ed9-cf73-d548-18a0fb78181f.png)
　一般式で表すと、k=0〜9, j=0〜1 で、
![スクリーンショット 2020-03-16 08.44.27.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/935ba320-c38c-d9d6-5d30-ddc8c409a7d2.png)
　となり、これによって重み$w^1_{kj}$の更新が可能になります。バイアスについては、$b^1$が1なので、上記式の$x^1_k$が1に代わるだけで、
![スクリーンショット 2020-03-16 10.02.58.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/729db79c-5b72-93d8-e0e4-9ce3041c2779.png)
　これによってバイアス$b^1_j$の更新が可能になります。


#5.誤差逆伝播（入力層から中間層)
　次に、入力層から中間層の重みとバイアスの更新です。

![スクリーンショット 2020-03-16 10.15.30.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/9eff1238-41f4-51f8-5b94-b6de5236c561.png)

　重み$w^0_{00}$を更新するためには、$\frac{\partial E^0}{\partial w^0_{00}}$と$\frac{\partial E^1}{\partial w^0_{00}}$ を求める必要があります。微分の連鎖率から、
![スクリーンショット 2020-03-16 08.32.17.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/51070a34-0388-b8dc-ae9e-79bbb51afd53.png)
　一般式で表すと、k=0〜783, j=0〜9 で、
![スクリーンショット 2020-03-16 08.33.56.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/428c1795-f4d6-d426-d746-bdcbd47cbdcd.png)
　となり、これによって重み$w^0_{kj}$の更新が可能になります。バイアスについては、$b^0$が1なので、上記式の$x^0_k$が1に代わるだけなので、
![スクリーンショット 2020-03-16 10.11.34.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/b82691c8-3006-f0b9-8425-e74b313a22e5.png)
　これによってバイアス$b^0_j$の更新が可能になります。

#6.順伝播、誤差逆伝播部分の実装
　先程求めた一般式を元に、順伝播、誤差逆伝播の部分を実装します。

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

# シグモイド関数の微分
def sigmoid_d(a):
    return (1 - sigmoid(a)) * sigmoid(a)

# 誤差逆伝播
def back(l, j):
    if l == max_layer - 1:
        return (y[j] - t[j]) * sigmoid_d(A[l][j])
    else:
        output = 0
        m = A[l+1].shape[0]   
        for i in range(m):
            output += back(l+1, i) * W[l+1][i,j] * sigmoid_d(A[l][j])
        return output
```
　`def back(l, j):`の具体的な動きは、

`l=1` の時は、
　`(y[j]-t[j])＊sigmoid_d(A[1][j])`が返されます。

`l=0` の時は、
　`(y[0]-t[0])＊sigmoid_d(A[1][0])*W[1][0,j]*sigmoid_d(A[0][j])
　+(y[1]-t[1])＊sigmoid_d(A[1][1])*W[1][1,j]*sigmoid_d(A[0][j])`
が返されます。


```python
# 重みWの設定
np.random.seed(seed=7)
w0 = np.random.normal(0.0, 1.0, (10, 784))
w1 = np.random.normal(0.0, 1.0, (2, 10))
W = [w0, w1]

# バイアスbの設定
b0 = np.ones((10, 1))
b1 = np.ones((2, 1))
B = [b0, b1]

# その他の設定
max_layer = 2 # 層数の設定
n = 0.5  # 学習率の設定
```
　重みWとバイアスb及びその他の設定をします。

![スクリーンショット 2020-03-16 19.48.21.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/613ead43-a443-c9d7-a55d-d7436785ee5c.png)
　学習がスムーズに開始できる様に、重み行列w0, w1の各項は、0〜1の正規分布に従う乱数にしています。ちなみに、`np.random.seed(seed=7)`のseed=の番号を変えると、学習のスタート具合（スムーズにスタートするか、少しもたつくか）が変化します。バイアス行列b0, b1は各項が1です。


```python
#学習ループ
count = 0 
acc = []

for x, t in zip(xs, ts):
    
    # 順伝播
    x0 = x.flatten().reshape(784, 1)
    a0 = W[0].dot(x0) + B[0]
    x1 = sigmoid(a0)
    a1 = W[1].dot(x1) + B[1]
    y = sigmoid(a1)

    # パラメータ更新のためにx,aをリスト化
    X = [x0, x1]
    A = [a0, a1]

    # パラメータ更新
    for l in range(len(X)):
        for j in range(W[l].shape[0]):
            for k in range(W[l].shape[1]):
                W[l][j, k] = W[l][j, k] - n * back(l, j) * X[l][k]  
            B[l][j] = B[l][j] - n * back(l, j) 
```
　学習ループです。順伝播は、行列の内積、足し算等で簡単に行えます。
パラメータ更新では、

`l = 0` の時は、`j = 0〜9, k = 0〜783`の範囲で、
　`W[0][j,k] = W[0][j,k] - n ＊ back(0,j) ＊ X[0][k]`
　`B[0][j] = B[0][j] - n ＊ back(0,j)`
　と更新されます。

`l = 1` の時は、`j = 0〜1, k = 0〜9`の範囲で、
　`W[1][j,k] = W[1][j,k] - n ＊ back(0,j) ＊ X[0][k]`
　`B[1][j] = B[1][j] - n ＊ back(0,j)`
　と更新されます。


#7.データセットの準備
　Keras で MNISTデータセットを読み込み、「１」, 「7」のみを抽出します。

```python
import numpy as np
from keras.datasets import mnist
from keras.utils import np_utils
import matplotlib.pyplot as plt

# 数字表示
def show_mnist(x):
    fig = plt.figure(figsize=(7, 7))   
    for i in range(100):
        ax = fig.add_subplot(10, 10, i+1, xticks=[], yticks=[])
        ax.imshow(x[i].reshape((28, 28)), cmap='gray')
    plt.show()

# データセット読み込み        
(x_train, y_train), (x_test, y_test) = mnist.load_data()
show_mnist(x_train)

# 1, 7　を抽出
x_data, y_data = [], []
for i in range(len(x_train)):  
    if y_train[i] == 1 or y_train[i] == 7:
       x_data.append(x_train[i])
       if y_train[i] == 1:
          y_data.append(0)
       if y_train[i] == 7:
          y_data.append(1)

show_mnist(x_data)

# list形式からnumpy形式に変換
x_data = np.array(x_data)
y_data = np.array(y_data)

# x_dataの正規化、y_dataのワンホット表現化
x_data = x_data.astype('float32')/255
y_data = np_utils.to_categorical(y_data)

# 学習、テストデータを取得
xs = x_data[0:200]
ts = y_data[0:200]  
xt = x_data[2000:3000]  
tt = y_data[2000:3000] 
```
![スクリーンショット 2020-03-16 14.42.10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/dcd323f8-9c58-1008-009f-31ee681a4e84.png)
![スクリーンショット 2020-03-16 14.42.25.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/3c49d9f6-e099-39d3-9754-335fcde0adeb.png)
　0〜9のデータと、そこから1,7を抽出したデータを、それぞれ先頭から100個表示しています。

　学習用データ xs, ts はそれぞれ200個、テスト用データ xt, tt はそれぞれ1,000個用意します。

#8.実装全体
　学習毎にテストデータによる精度確認と精度推移グラフ表示をプラスした、実装の全体です。

```python
import numpy as np
from keras.datasets import mnist
from keras.utils import np_utils

# データセット読み込み
(x_train, y_train), (x_test, y_test) = mnist.load_data()

# 1, 7 の数字のみ抽出
x_data, y_data = [], []
for i in range(len(x_train)):  
    if y_train[i] == 1 or y_train[i] == 7:
       x_data.append(x_train[i])
       if y_train[i] == 1:
          y_data.append(0)
       if y_train[i] == 7:
          y_data.append(1)

# list形式からnumpy形式に変換
x_data = np.array(x_data)
y_data = np.array(y_data)

# x_dataの正規化, y_dataのワンホット化
x_data = x_data.astype('float32')/255
y_data = np_utils.to_categorical(y_data)

# 学習データ、テストデータの取得
xs = x_data[0:200]  
ts = y_data[0:200]  
xt = x_data[2000:3000]  
tt = y_data[2000:3000]  


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

# シグモイド関数の微分
def sigmoid_d(a):
    return (1 - sigmoid(a)) * sigmoid(a)

# 誤差逆伝播
def back(l, j):
    if l == max_layer - 1:
        return (y[j] - t[j]) * sigmoid_d(A[l][j])
    else:
        output = 0
        m = A[l+1].shape[0]   
        for i in range(m):
            output += back(l + 1, i) * W[l + 1][i, j] * sigmoid_d(A[l][j])
        return output

# 重みWの設定
np.random.seed(seed=7)
w0 = np.random.normal(0.0, 1.0, (10, 784))
w1 = np.random.normal(0.0, 1.0, (2, 10))
W = [w0, w1]

# バイアスbの設定
b0 = np.ones((10, 1))
b1 = np.ones((2, 1))
B = [b0, b1]

# その他の設定
max_layer = 2 # 層数の設定
n = 0.5  # 学習率の設定

#学習ループ
count = 0 
acc = []

for x, t in zip(xs, ts):
    
    # 順伝播
    x0 = x.flatten().reshape(784, 1)
    a0 = W[0].dot(x0) + B[0]
    x1 = sigmoid(a0)
    a1 = W[1].dot(x1) + B[1]
    y = sigmoid(a1)

    # パラメータ更新のためにx,aをリスト化
    X = [x0, x1]
    A = [a0, a1]

    # パラメータ更新
    for l in range(len(X)):
        for j in range(W[l].shape[0]):
            for k in range(W[l].shape[1]):
                W[l][j, k] = W[l][j, k] - n * back(l, j) * X[l][k]  
            B[l][j] = B[l][j] - n * back(l, j) 
            
    # テストデータによる精度チェック       
    correct, error = 0, 0

    for i in range(1000):

        # 学習したパラメータで推論
        x0 = xt[i].flatten().reshape(784, 1)
        a0 = W[0].dot(x0) + B[0]
        x1 = sigmoid(a0)
        a1 = W[1].dot(x1) + B[1]
        y = sigmoid(a1)
    
        if np.argmax(y) == np.argmax(tt[i]):
           correct += 1
        else:
           error += 1
    calc = correct/(correct+error)
    acc.append(calc)
    count +=1
    print("\r[%s] acc: %s"%(count, calc))
   
# 精度推移グラフ表示    
import matplotlib.pyplot as plt
plt.plot(acc, label='acc')
plt.legend()
plt.show()   
```
![スクリーンショット 2020-03-16 17.21.46.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/209705/36c41189-8397-6229-acbf-69e01c9524a4.png)
　200ステップで、分類精度97.8％となりました。自分でスクラッチで実装したニューラルネットワークがちゃんと動くと嬉しいものですね。