誤差逆伝播法をPythonで実装しながら理解する

「誤差伝播法ってなぁに？」  
「ニューラルネットワークのパラメータの勾配を求める為の手法だよ」

## はじめに

誤差逆伝播法を、Pythonで実装しながら理解していくよ。最終的には簡単なニューラルネットワークを作る。

実装は割とシンプルにしたつもり。Pythonがある程度わかる人は、コードを見ることで理解が深まるかも！

## 合成関数の微分

高校でやるやつ。これが分かれば誤差逆伝播法なんてほとんど理解したようなもん

<br>

試しに、以下の関数を微分してみよう。

```math
y = (x + 1)^2
```

これは普通に展開しても解けるけど、合成関数の微分を使っても解けるね

```math
\begin{align}

u &= x + 1 \\
y &= u^2 \\
\frac{dy}{dx}
&= \frac{dy}{du} \frac{du}{dx} \\
&= 2u \cdot 1 \\
&= 2(x + 1)

\end{align}
```

こんな感じで。**関数の関数**(合成関数)を微分するときは、関数ごとに微分をしたものをかけ合わせればよかった。<br><br>


じゃあここまでの流れをPythonで実装してみよう

以下の二つの関数を`class`として実装する。
- $f(x) = x + 1$
- $g(x) = x^2$

`class`にする必要ある？関数でよくね？と思うかもしれんけど、まあ読んでみてよ。  
まずは$f(x) = x + 1$から


```Python
class Plus1:
    def __call__(self, x):
        return x + 1
```

出来た。入力した値に1を足して出力するだけの関数。`__call__()`というのは特殊メソッドで、関数のように`()`をつけて呼び出したときに実行されるヤツ。  
こんな感じ

```Python
plus1 = Plus1() # インスタンス生成
y = plus1(3) # __call__メソッドの呼び出し
print(y)
>> 4
```

入力した3に1を足した4が出力された。

では、今度は微分を行うメソッドを書いてみよう。「微分を行う」というのを、「ある入力された値に、微分した値をかけて出力する」と捉えるとこうなる

```Python
class Plus1:
    def __call__(self, x):
        return x + 1

    def backward(self, d):
        return d * 1
```

`backward()`というメソッドを追加した。  
$x + 1$を$x$で微分すると1になるので、入力値(d)に1をかけて出力させる。

こんなノリで、$g(x) = x^2$の方も書いちゃおう

```Python
class Square:
    def __call__(self, x):
        self.x = x
        return x ** 2

    def backward(self, d):
        return d * self.x*2
```

できた。$g(x)$を$x$で微分すると$2x$になるので、`backward()`では、それを入力値(d)にかけて出力する。微分するときに使うので、初めに入力された値(x)は変数に保存しておく。  
こんな感じに、`backward()`を追加したり、値を一時的に保存したりしたかったから、`class`として実装したのだ

<br>

ちなみに、`__call__()`で行なっている演算は**順伝播**  
`backward()`で行なっている演算は**逆伝播**という

あと、順伝播(`__call__()`)への入力と逆伝播(`backward()`)への入力が混ざる気がするので、今後以下のように区別する。  
- 順伝播への入力: **入力(x)**
- 逆伝播への入力: **入力(d)**

<br><br>

ではこれらを使って、実際に計算してみよう。

そしてここから、

```math
h(x) = (x + 1)^2
```

としよう

まずは順伝播から

```Python
# インスタンス生成
plus1 = Plus1()
square = Square()

# 計算
x = 3
z = plus1(x)
y = square(z)
print(y)
>> 16
```

ほい。$h(3) = 16$ということで、正解！

じゃあ今度は逆伝播で、$h$の$x=3$での傾きを求めてみよう。$h'(3)$のことだね。  
そしてこれはこんな感じで求められる

```Python
d = square.backward(1)
d = plus1.backward(d)
print(d)
>> 8
```

8とでてきた。計算してみると

```math
\begin{align}
h'(x) &= 2(x + 1) \\
h'(3) &= 8
\end{align}
```

となるので正解！  
これは合成関数の微分に基づいていて、正に「関数ごとの微分を掛け合わせる」という部分に当たるね！

> 入力された値に微分した値をかけて出力する

さっきこう捉えた意味が分かったかな...?  
微分した結果を逆方向に伝えていく感じだね。一番初めは特に何もないので1を入力しておく。

そして今後この「傾き」を「勾配」と呼ぶことにするよ。


## ニューラルネットワークの構築

見てもらった通り、複数の関数を経て出力された値を微分した値は、関数ごとに微分をすれば簡単に求まる。  
んで、これはニューラルネットワークが持つパラメータの勾配を求めるときにも使える。

損失をパラメータで微分した値(勾配)を求めるとき、損失を出す際に通った**層**や**損失関数**を一つ一つ微分すればいいよねという話。

じゃあ、誤差逆伝播法で学習を行うニューラルネットワークを実際に作ってみようじゃないか

### 層の定義

NNは複数の層から構成されるので、まずは層を作る。

#### ReLU

重みとかバイアスみたいなパラメータを持つ層は工夫が必要なので、一旦パラメータを持たない層を作ってみよう。  
ReLUはこういう関数

```math
f(x) = \left\{
\begin{array}{ll}
x & (x > 0) \\
0 & (x \leq 0)
\end{array}
\right.
```

さっきの関数と同じように作ればOK

```Python
import numpy as np

class ReLU:
    def __call__(self, x):
        self.x = x # 逆伝播で使うので変数に保存
        return np.maximum(0, x) # 0より大きければx、小さければ0

    def backward(self, d):
        return d * (self.x > 0) # xが0より大きければ傾き1、小さけれ0にする

    def update(self, lr):
        pass
```

できた。bool型は`+`とか`*`みたいな演算子と一緒に使うと0,1として扱ってくれるので、`backward()`はこういう書き方でOK。

さっき実装したような一般的な関数の入出力は「数値」だけど、NNの層の入出力はベクトルで行うのでnumpyを使うよ  
あと、後々のことを考えて、パラメータを更新するメソッド`update()`を書いている。ReLUはパラメータを持たないので何もしないけど

#### 全結合層

全結合層。kerasでいうDense。PyTorchでいうLinear。Affineと呼ぶこともあるし、Fully Connected(全結合)からfcと指すこともある。ここではPyTorchに倣ってLinearにしよー

```Python
class Linear:
    def __init__(self, n_input, n_output):
        self.w = np.random.normal(size=(n_input, n_output), scale=(2/n_output)**0.5)
        self.b = np.random.randn(n_output)

    def __call__(self, x):
        self.x = x
        return np.dot(x, self.w) + self.b

    def backward(self, d):
        self.grad_w = np.multiply(*np.meshgrid(d, self.x))
        self.grad_b = d
        return np.dot(d, self.w.T)

    def update(self, lr):
        self.w -= lr * self.grad_w
        self.b -= lr * self.grad_b
```

できた。**ミニバッチは非対応**にした。そのせいで`grad_w`が複雑になっているけど、「逆伝播」についてはミニバッチじゃない方が分かりやすい気がしたので。

パラメータを持つ層なので、`__init__()`でパラメータを初期化している。  
重みもバイアスも正規分布に従って初期化してるけど、重みの初期化は`randn()`ではなく`normal()`を使っている。  
`randn()`だと標準偏差が1になるんだけど、ReLUを使う場合、それだと重みが大きくなりすぎて上手くいかないときがあるので、小さくしてる。詳しくは**Xaivierの初期化**とかでググってみて。

`backward()`では、入力(d)を各パラメータで微分して変数に入れておく(`grad_なんちゃら`)。  
それらはパラメータの勾配となるので、`update()`(パラメータ更新)の時は、それらに学習率をかけてパラメータから引く感じ。<br><br><br><br>


**`backward()`(逆伝播)** について、もう少し詳しく見ていこう

まず、あるニューロンに対する順伝播時の入力(x)を$x$、出力を$y$、重みを$w$、バイアスを$b$として、これらの関係を表す。

$$
x \cdot w + b = y
$$

内積にバイアスを足しているね

これを踏まえて、ここ(逆伝播)で求める値を説明する。求める値は3つ。

- 重みの勾配(`grad_w`)

$l$番目の入力($x_l$)から$m$番目の出力($y_m$)への重みを$w_{lm}$とすると、求めたいのは$y_m$を$w_{lm}$で微分したもの。  
で、それは入力(d)に$x_l$をかけたものになる。$l, m$ではない部分やバイアスは微分をする上で関係ないので無視してよくて、残るのは$x_l \times w_{lm}$の部分だけという感じ。これを$w_{lm}$で微分すると$x_l$になる。

関係のある部分だけ図にするとこうなる。

![img1.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2235373/e8ee4ce7-409f-806b-9b08-15c6c83a2f85.png)


これを全ての$l, m$で求めるのがこのコード: `np.multiply(*np.meshgrid(self.x, d))`  
これは、全部の組み合わせで掛け算をしているだけ。↓の例を見て理解してくれ。

```Python
a = [1, 2, 3]
b = [4, 5]
y = np.multiply(*np.meshgrid(a, b))
print(y)
>> [[ 4  8 12]
    [ 5 10 15]]
```

(これ一発でできる関数とか知っている人いたら教えてほしい...)<br><br>


- バイアスの勾配(`grad_b`)

バイアスで微分する上で内積の部分は関係ないので無視してよい。上の式を$b$で微分すると1なので、バイアス(b)の勾配は、入力(d)に1をかけたd<br><br>

- 入力(x)の勾配(`return`)

ある一つの入力(x)に着目する。この入力(x)は$x_i$としよう。  
ここで、層のニューロンの数が一つだと仮定する。そうすると話は簡単。重みの勾配を求めたときと逆のことをすればいいので、$x_i$の勾配はそこに対応する重み$w_i$となる。  
でも実際は層のニューロンが複数あるので、勾配が複数出てくる。そんなときにどうすればよいか。

これは、実はとってもシンプルで、出てきた勾配を全部足せばよい。で、それは行列の積で上手くまとめられて、それがこれ: `np.dot(d, self.w.T)`

### ニューラルネットワーク

ではこれらを組み合わせてNNをつくろー

```Python
class NeuralNetwork:
    def __init__(self, *layers):
        self.layers = layers

    def __call__(self, x):
        for layer in self.layers:
            x = layer(x)
        return x

    def backward(self, d):
        for layer in self.layers[::-1]:
            d = layer.backward(d)

    def update(self, lr):
        for layer in self.layers:
            layer.update(lr)
```

できた。これだけ。かんたん。  
使い方は、インスタンス生成時にレイヤーを入れていくだけ。kerasとかPyTorchでいうSequentialみたいな感じ。以下メソッドの説明

- `__init__()`
与えられたレイヤーを変数に入れる。

- `__call__()`  
順伝播。与えた層で順番に演算しているだけ。

- `backward()`  
逆伝播。リストに`[::-1]`をつけると逆順になる。

- `update()`  
パラメータ更新。逆伝播で求めた勾配に学習率をかけてパラメータから引く。ここをシンプルにするために、パラメータを持たない層にもこのメソッドを書いておいた。

<br>

試しに使ってみよう

```Python
nn = NeuralNetwork(
    Linear(5, 32),
    ReLU(),
    Linear(32, 10,)
)

x = np.random.randn(5)
y = nn(x)
print(y)
>> [ 2.70924181 -1.15491763  2.64228401  1.00399676 -1.36650978 -2.33312378
     2.08624605 -0.2996883   0.48605541  2.26515702]
```

5個の値を入ると10個の値が返ってくるモデルを作った。でそれに適当な5個の乱数を入れると適当な10個の値が返ってきた。

### 損失関数

損失関数もクラスとして書いてみよう

#### 二乗和誤差

差の二乗の和...を、2で割ったもの。

$$
\frac{1}{2} \sum_{i} (y_i - t_i)^2
$$


なんで2で割るかって？ 微分した時に綺麗になるからだよ〜

```Python
class RSS:
    def __call__(self, y, t):
        self.y = y
        self.t = t
        return np.sum((y - t) ** 2) / 2

    def backward(self):
        return self.y - self.t
```

`backward()`の中身、綺麗でしょ〜。微分すると二乗の部分が前に出てくるので、$\frac{1}{2}$と打ち消し合っていい感じになる  
あと、ここの逆伝播では何も入力しない。逆伝播の一番最初の部分だから入力はないね。

#### 交差エントロピー

分類タスクで使うやつ。

```math
\sum_i - t_i \, log \, y_i 
```

これを使うときは、出力層の活性化関数はsoftmaxを使う。ので、それも一緒に書いちゃおう！  
softmaxはこれ

```math
y_k = \frac{e^{y_k}}{\sum_i e^{y_i}}
```

```Python
class CrossEntropy:
    def __call__(self, y, t):
        y = self._softmax(y)
        self.y = y
        self.t = t
        loss = -np.sum(t * np.log(y))
        return loss

    def backward(self):
        return self.y - self.t

    def _softmax(self, y):
        return np.exp(y) / np.sum(np.exp(y))
```

できた。これも、逆伝播がとってもきれいな形。普通活性化関数はNN側に書くけど、交差エントロピーと一緒に使うとこの部分の微分が綺麗になるので、こっちに書いた。  
微分がきれいなのはたまたまではなく、頭のいい人が綺麗になるように考えたからだと思う。<br><br><br>

二つの損失関数を書いたけど、どちらも逆伝播が「予測値-正解」になっているね。  
勾配を求めるときは、これをNNに入れて逆伝播を行うんだけど、「予測値-正解」って、正解と予測値の差だから、「誤差」と表せるよね。

「誤差逆伝播法」の由来が分かった気がするね


## 学習

では実際に学習させてみよう。定番のMNIST。

```Python
from keras.datasets import mnist
(x_train, y_train), (x_test, y_test) = mnist.load_data()

# 平坦化・正規化
x_train = x_train.reshape(-1, 784) / 255
x_test = x_test.reshape(-1, 784) / 255

# one-hotベクトルに
y_train = np.eye(10)[y_train]
y_test = np.eye(10)[y_test]
```

確率的勾配降下法で学習させる。ミニバッチに対応していないので、バッチサイズは1

```Python
def train(model, x, y, criterion, lr, n_epochs):
    for epoch in range(1, n_epochs+1):
        loss = 0
        for _ in range(len(x)):
            # ランダムに一つ抽出
            idx = np.random.randint(0, len(x))
            sample_x = x[idx]
            sample_y = y[idx]

            out = model(sample_x) # 順伝播
            loss += criterion(out, sample_y) # 損失の計算
            d = criterion.backward() # 損失関数の逆伝播
            model.backward(d) # ニューラルネットワークの逆伝播
            model.update(lr) # パラメータの更新
        print(f'{epoch}epoch loss:{loss / len(x)}')
```

ネットワーク構成は適当！

```Python
nn = NeuralNetwork(
    Linear(784, 512),
    ReLU(),
    Linear(512, 128),
    ReLU(),
    Linear(128, 10),
)
```

レッツゴー

```Python
train(nn, x_train, y_train, CrossEntropy(), 0.001, 5)
>> 1epoch loss: 0.2695821597613696
   2epoch loss: 0.10846042071623153
   3epoch loss: 0.07909254203673934
   4epoch loss: 0.05543957445079766
   5epoch loss: 0.04660230423813134
```

私の環境では、1エポックに2分ぐらいだった。cpuなので、こんなもん。

いい感じにlossが減っているので上手くいってそう。確認してみよう。

```Python
import matplotlib.pyplot as plt

sample_x = x_test[0]
sample_y = y_test[0]
t = np.argmax(sample_y)
pred = np.argmax(nn(sample_x))
plt.title(f'true: {t}, pred: {pred}')
plt.imshow(sample_x.reshape(28, 28), cmap='gray')
```
![output.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/2235373/43a0911f-6885-753d-2591-e900278eb392.png)

ちゃんと7と予測できているね。

精度も確かめてみよう

```Python
predicts = np.array([nn(x).argmax() for x in x_test])
y_test = y_test.argmax(axis=1)
print('accuracy:', (predicts == y_test).mean())
>> accuracy: 0.9728
```

バッチリ

## おわり

おわりです

**参考**
- [O'Reilly Japan - ゼロから作るDeep Learning](https://www.oreilly.co.jp/books/9784873117584/)