はじめに

今までで一番今更感がある記事ですが、誤差逆伝播法について理解が弱かったので改めて勉強してみました。

注意事項
・誤差逆伝播法の行列ですが、私は高速化手法という位置づけだと思っているので、本質の理解を優先し本記事では除外しています。
・本記事はプログラムで実装できるまでの理解を目安にしており、その先の数学的な厳密性までは踏み込んでいません。
更にちゃんと理解したい場合は参考サイトや別の記事を参照してください。

人工ニューラルネットワーク

人工ニューラルネットワークは人の脳の神経網（生物のニューラルネットワーク）を模した数理モデルで、各ユニット（ニューロン）が重み（シナプス）を介して繋がっており、重みの強さを調整することで出力を変化させるモデルとなります。

ユニットは各入力に対して重みを持っており、入力×重みの合計に活性化関数を通して値を出力します。
活性化関数はユニットの性質を決める任意の関数です。

また各ユニットは複数連結させることを想定しており、一般的には層と呼ばれる単位で以下のように表現されます。

参考：https://ja.wikipedia.org/wiki/%E3%83%8B%E3%83%A5%E3%83%BC%E3%83%A9%E3%83%AB%E3%83%8D%E3%83%83%E3%83%88%E3%83%AF%E3%83%BC%E3%82%AF

普遍性定理/万能近似定理

ニューラルネットワークをここまで有名にした土台となる定理です。
この定理を簡単に言うと、どんな複雑な関数でもユニット（及び層）を適切に用意すればその関数を限りなく正確に表現できるというものです。

定理名ですがネット上では2種類引っ掛かり、普遍性定理（Universal Approximation Theorem）は主に1層を想定してユニット数を増やしたもの、万能近似定理（Universal Representation Theorem）は層も増やしたものに使われている印象を受けました。
（前者が3層ニューロンネットワークで後者がディープニューラルネットワークに関係しているのかな）
（後日本語訳が結構バラバラで英語名の方が正確かも？）

数学的な証明はかなり難しいらしいので省略し、直感的な理解にのみ説明します。

まず、ニューラルネットワークは活性化関数を通して出力されるので、以下のように切り取ると関数の足し算（合成関数）と見ることができます。

この合成関数ですが、無限の関数を合成した場合どんな関数でも表現できることが保証されているというのがこの定理の主張です。

どういうことかを活性化関数がReLUの場合で実際に見てみます。
ReLUはxが0より大きい値ならx、0以下なら0を返す関数で、pythonですと以下です。

def relu(x):
    if x > 0:
        return x
    else:
        return 0

ReLUは重みで傾きが変わり、biasで横の位置が変わります。

合成のイメージは以下です。

実際に２つのReLU関数を合成してみます。

左の図が $y_1=relu(x)$ と $y_2=-relu(x-1)$ の合成です。（緑の関数）
この緑の関数ですが、傾きを急にしていくとステップ関数に限りなく近づくことができます。（右図）
ステップ関数は0未満で0、0以上で1を返す関数です。

このステップ関数を組み合わせると曲線を描くことができます。

ディスプレイを細かく見ると赤青緑の点から色が表現されているみたいな感じで、ギザギザを細かくすれば直線でも曲線に見えるイメージですね。
厳密な証明ではないですが、感覚的にこれを無限に細かくすればどんな形でも表現できることが分かると思います。

出力コード

import matplotlib.pyplot as plt
import numpy as np


def relu(x):
    return np.maximum(0, x)


x = np.linspace(-1, 3, 1000)

fig, axes = plt.subplots(1, 2, figsize=(14, 6))

y1 = relu(x)
y2 = -relu(x - 1)

axes[0].plot(x, y1, label="y1 = relu(x)")
axes[0].plot(x, y2, label="y2 = -relu(x-1)")
axes[0].plot(x, y1 + y2, label="y1 + y2")
axes[0].set_title("Graph of y = relu(x)")
axes[0].set_xlabel("x")
axes[0].set_ylabel("y")
axes[0].set_ylim(-2, 3)
axes[0].axhline(0, color="black", linewidth=0.5)
axes[0].axvline(0, color="black", linewidth=0.5)
axes[0].grid(color="gray", linestyle="--", linewidth=0.5)
axes[0].legend()

y1 = 100 * relu(x)
y2 = -100 * relu(x - 1 / 100)

axes[1].plot(x, y1, label="y1 = 100*relu(x)")
axes[1].plot(x, y2, label="y2 = -100*relu(x-0.01)")
axes[1].plot(x, y1 + y2, label="y1 + y2")
axes[1].set_title("Graph of y = -relu(x) + 1")
axes[1].set_xlabel("x")
axes[1].set_ylabel("y")
axes[1].set_ylim(-2, 3)
axes[1].axhline(0, color="black", linewidth=0.5)
axes[1].axvline(0, color="black", linewidth=0.5)
axes[1].grid(color="gray", linestyle="--", linewidth=0.5)
axes[1].legend()

plt.tight_layout()
plt.show()

参考
・ニューラルネットワークにおけるUniversal Approximation Theorem(普遍性定理)について(Qiita)
・ニューラルネットワークの普遍性定理(Qiita)
・ニューラルネットワークが任意の関数を表現できることの視覚的証明

誤差逆伝播法（Backpropagation;BP）

ニューラルネットの表現力は分かったと思いますが、次にこれをどう学習するかという話です。
一般的には誤差逆伝播法と呼ばれる手法で学習されます。

誤差逆伝播法を理解するには「勾配降下法」と「連鎖率」を理解する必要があるのでそれを説明していきます。

勾配降下法

勾配降下法は関数の最小値（または最大値）を見つける最適化手法です。
ある関数を微分すると、特定の地点xの接線の傾き（勾配）を得ることができます。
この傾きはxに対する出力の方向と強さを表しています。

出力を増やしたい場合は勾配方向に値を増やし、出力を減らしたい場合は勾配とは逆方向に値を増やして出力をコントロールし、出力が最小となる点を探す手法が勾配降下法です。

微分の表記

私は数学に普段触れるような生活はしていないので軽く復習しておきます。

$f(x)=x^2+1$ を$x$に対して偏微分した表現を、$\frac{\partial \boldsymbol{f}}{\partial \boldsymbol{x}}f(x)$ または $\frac{\partial \boldsymbol{f}}{\partial \boldsymbol{x}} \left\{ x^2+1 \right\}$ と表現します。
また、$f(x)$ に対して偏微分を意識する必要がなくなり、ただの微分と同じになった場合は $f'(x)$ と表記します。
※分かりやすいように微分の範囲を $\left\{ \right\}$ で表します。
※また意味のある掛け算をあえて明示して書きます。

参考：https://ja.wikipedia.org/wiki/%E5%BE%AE%E5%88%86%E3%81%AE%E8%A8%98%E6%B3%95

連鎖率（chain rule）

合成関数を偏微分する場合以下の式が成り立ちます。
例えば $(x)$ から $(u,v)$ が定まり、$(u,v)$ から $f$ が定まるときの $f$ に対する $x$ の偏微分は以下です。

$$
\frac{\partial \boldsymbol{f}}{\partial \boldsymbol{x}}
=\frac{\partial \boldsymbol{f}}{\partial \boldsymbol{{\color{red}u}}}
\frac{\partial \boldsymbol{{\color{red}u}}}{\partial \boldsymbol{x}}
+\frac{\partial \boldsymbol{f}}{\partial \boldsymbol{{\color{red}v}}}
\frac{\partial \boldsymbol{{\color{red}v}}}{\partial \boldsymbol{x}}
$$

具体的に $f(x) = (x^2 + 1)^2$ で、$x$の偏微分で見てみます。
$u=x^2+1$ と置くと連鎖率を用いて以下のように計算できます。

\begin{align}
\frac{\partial \boldsymbol{f}}{\partial \boldsymbol{x}}
f(x) &= \frac{\partial \boldsymbol{f}}{\partial \boldsymbol{x}} \left\{u^2\right\} \\
&= \frac{\partial \boldsymbol{f}}{\partial \boldsymbol{{\color{red}u}}}\left\{u^2\right\} \times \frac{\partial \boldsymbol{{\color{red}u}}}{\partial \boldsymbol{x}} \left\{x^2 + 1 \right\} \tag{1} \\
&= 2u \times 2x \\
&= 2(x^2 + 1) \times 2x \\
&= 4x^3 + 4x \\
\end{align}

$(1)$が連鎖率を適用した形で「全体の微分 × 括弧内の微分」に分けることができます。

参考
・連鎖律（多変数関数の合成関数の微分）(高校数学の美しい物語)
・合成関数の偏微分における連鎖律(チェインルール)とその証明(数学の景色)

誤差逆伝播

モデルの出力を二乗和誤差（SSE：Sum of Squared Error）で最小化することを考えます。
（平均二乗誤差だと思い込んでいました…orz、基礎に戻るのは大事ですね）

モデルの出力を$o$、正解データを$t$と置くと以下です。

$$
E(o,t) = \frac{1}{2}\sum_i(o_i-t_i)^2 \
$$

※$\frac{1}{2}$ は微分後の値をすっきりさせるための係数です
このSSEを最小化するために各重み$w$を勾配方向の逆に変化させることを考えます。
例として以下の図で考えます。（図の記号と数式は一致させています）

この図における $w_3$ と $w_1$ の重みの勾配を計算してみます。
まずは$w_3$の勾配として$w_3$の偏微分を計算します。

\begin{align}
E(o,t) &= \frac{1}{2} \bigg( (o_1-t_1)^2 + (o_2-t_2)^2 \bigg) \\
o_1 &= g(y_2) \\
y_2 &= x_3w_3 + x_4w_4 + bias2 \\
\end{align}

\begin{align}
\frac{\partial \boldsymbol{E}}{\partial \boldsymbol{w_3}} E(o,t) &= \frac{\partial \boldsymbol{E}}{\partial \boldsymbol{w_3}} \left\{ \frac{1}{2} \bigg( (o_1-t_1)^2 + (o_2-t_2)^2 \bigg) \right\} \\
&= \frac{\partial \boldsymbol{E}}{\partial \boldsymbol{{\color{red}o_1}}} \left\{ \frac{1}{2} \bigg( (o_1-t_1)^2 + (o_2-t_2)^2 \bigg) \right\} 
\times \frac{\partial \boldsymbol{{\color{red}o_1}}}{\partial \boldsymbol{{\color{red}y_2}}} g(y_2)
\times \frac{\partial \boldsymbol{{\color{red}y_2}}}{\partial \boldsymbol{w_3}} (x_3w_3 + x_4w_4+bias2) \\
&= (o_1-t) \times  g'(y_2) \times x_3 \\
&= \Delta y_2 \times x_3
\end{align}

次の層への計算用に $(o_1-t) \times g'(y_2)$ を $\Delta y_2$ と表現しています。

次に$w_1$の偏微分です。
$w_1$を計算する場合、導線が$x_3$側と$x_5$側の2方向ありますが、この場合連鎖率は足し算をします。

\begin{align}
o_2 &= h(y_3) \\
y_3 &= x_5w_5 + bias3 \\
x_5 &= f(y_1) \\
x_3 &= f(y_1) \\
y_1 &= x_1w_1 + x_2w_2 + bias1
\end{align}

\begin{align}

\frac{\partial \boldsymbol{E}}{\partial \boldsymbol{w_1}} E(o,t) &= \frac{\partial \boldsymbol{E}}{\partial \boldsymbol{w_1}} \left\{ \frac{1}{2} \bigg( (o_1-t_1)^2 + (o_2-t_2)^2 \bigg) \right\} \\

&= \frac{\partial \boldsymbol{E}}{\partial \boldsymbol{{\color{red}o_1}}} \left\{ \frac{1}{2} \bigg( (o_1-t_1)^2 + (o_2-t_2)^2 \bigg) \right\} 
\times \frac{\partial \boldsymbol{{\color{red}o_1}}}{\partial \boldsymbol{{\color{red}y_2}}} g(y_2)
\times \frac{\partial \boldsymbol{{\color{red}y_2}}}{\partial \boldsymbol{\color{red}x_3}} \left\{ x_3w_3 + x_4w_4+bias2 \right\} \\
& \quad \times \frac{\partial \boldsymbol{{\color{red}x_3}}}{\partial \boldsymbol{\color{red}y_1}} f(y_1) \times \frac{\partial \boldsymbol{{\color{red}y_1}}}{\partial \boldsymbol{w_1}} \left\{ x_1w_1 + x_2w_2+bias1 \right\} \\
&\quad + \frac{\partial \boldsymbol{E}}{\partial \boldsymbol{{\color{red}o_2}}} \left\{ \frac{1}{2} \bigg( (o_1-t_1)^2 + (o_2-t_2)^2 \bigg) \right\} 
\times \frac{\partial \boldsymbol{{\color{red}o_2}}}{\partial \boldsymbol{{\color{red}y_3}}} h(y_3)
\times \frac{\partial \boldsymbol{{\color{red}y_3}}}{\partial \boldsymbol{\color{red}x_5}} \left\{ x_5w_5 + bias3 \right\} \\
& \quad \times \frac{\partial \boldsymbol{{\color{red}x_5}}}{\partial \boldsymbol{\color{red}y_1}} f(y_1) \times \frac{\partial \boldsymbol{{\color{red}y_1}}}{\partial \boldsymbol{w_1}} \left\{ x_1w_1 + x_2w_2+bias1 \right\} \\

&=\quad (o_1-t) \times g'(y_2) \times w_3 \times f'(y_1) \times x_1\\
& \quad + (o_2-t) \times h'(y_3) \times w_5 \times f'(y_1) \times x_1\\

=& (\Delta y2 \times w_3 + \Delta y3 \times w_5 ) \times f'(y_1) \times x_1 \\
=& \Delta y1 \times x_1 \\

\end{align}

計算式まとめ

各記号を以下と置いた場合のニューロンに関する計算式は以下です。（ニューロン視点で書いています）

ニューロンの接続元：$i$
ニューロンの入力 : $x_i$
ニューロンの重み : $w_i$
ニューロンのバイアス : $b$
ニューロンの活性化関数適用前の値 : $y=\sum_i x_i w_i + b$
ニューロンの活性化関数 : $f$
ニューロンの接続先：$j$
学習率 : $\alpha$

順伝搬（ニューロンの出力）

$$
f(y)
$$

勾配（重み更新時の前のニューロンへの入力）

$$
d_i = f'(y) w_i \sum_j( d_{j} )
$$

重みの更新

$$
w_i \leftarrow w_i - \alpha f'(y) x_i \sum_j(d_j)
$$

biasの更新

$$
b \leftarrow b - \alpha f'(y) \sum_j(d_j)
$$

参考
・誤差逆伝播法のノート(Qiita)
・誤差逆伝播法をはじめからていねいに(Qiita)
・ニューラルネットワークの数学（逆伝播）(キカガク)

実装例

pythonのlistのみで実装してみました。

import等

import math
import random

import numpy as np
from tqdm import tqdm

np.random.seed(30)
random.seed(40)

活性化関数

class Sigmoid:
    def __call__(self, x: float) -> float:
        return 1 / (1 + math.exp(-x))

    def derivative(self, x: float) -> float:
        return self(x) * (1 - self(x))


class ReLU:
    def __call__(self, x: float) -> float:
        return x if x > 0 else 0

    def derivative(self, x: float) -> float:
        return 1 if x > 0 else 0


class Liner:
    def __call__(self, x: float) -> float:
        return x

    def derivative(self, x: float) -> float:
        return 1

損失関数


class SSE:
    def __call__(self, y: list[float], target: list[float]) -> float:
        self.y = y
        self.target = target
        return sum([(y[i] - target[i]) ** 2 for i in range(len(y))]) / 2

    def backward(self) -> list[float]:
        return [self.y[i] - self.target[i] for i in range(len(self.y))]


class CrossEntropy:
    def __call__(self, y: list[float], target: list[float]) -> float:
        y = self._softmax(y)
        self.y = y
        self.target = target
        cross_entropy = sum([target[i] * math.log(y[i]) for i in range(len(y))])
        return -cross_entropy

    def backward(self) -> list[float]:
        return [self.y[i] - self.target[i] for i in range(len(self.y))]

    def _softmax(self, y: list[float]) -> list[float]:
        total = sum([math.exp(y2) for y2 in y])
        return [math.exp(y2) / total for y2 in y]

ニューロンとモデル


class Neuron:
    def __init__(self, in_num: int, activate) -> None:
        self.in_num = in_num
        self.act_w = activate
        self.weights = [random.random() / in_num for _ in range(in_num)]
        self.bias = 0.0

    def forward(self, x: list[float]) -> float:
        assert len(self.weights) == len(x)
        y = sum([self.weights[i] * x[i] for i in range(self.in_num)])
        y += self.bias
        self.x = x
        self.y = y
        return self.act_w(y)

    def backward(self, error: float) -> list[float]:
        dy = self.act_w.derivative(self.y)
        # 勾配は error*dy*x
        self.grad_w = [error * dy * self.x[i] for i in range(self.in_num)]
        self.grad_b = error * dy
        # 次に渡す勾配は error*dy*w
        return [error * dy * self.weights[i] for i in range(self.in_num)]

    def update(self, lr: float):
        for i in range(self.in_num):
            self.weights[i] -= lr * self.grad_w[i]
        self.bias -= lr * self.grad_b


class BPModel:
    def __init__(self, in_num: int, layers) -> None:
        # (layers, layer(neurons), neuron)
        self.layers: list[list[Neuron]] = []
        for size, activate in layers:
            layer = [Neuron(in_num, activate) for _ in range(size)]
            self.layers.append(layer)
            in_num = size

    def forward(self, x: list[float]):
        for layer in self.layers:
            x = [neuron.forward(x) for neuron in layer]
        return x

    def backward(self, errors: list[float]):
        for layer in reversed(self.layers):
            # errorは順伝搬で次につながっているニューロンの合計
            next_errors = [0.0 for _ in range(layer[0].in_num)]
            for i, neuron in enumerate(layer):
                e = neuron.backward(errors[i])
                assert len(e) == len(next_errors)
                for j in range(len(e)):
                    next_errors[j] += e[j]
            errors = next_errors

    def update(self, lr: float):
        for layer in self.layers:
            [neuron.update(lr) for neuron in layer]

実行


def _simple_run(dataset, model: BPModel, criterion, epochs, lr):
    for i in range(epochs):
        x, target = dataset[i % len(dataset)]

        y = model.forward(x)
        loss = criterion(y, target)
        error = criterion.backward()
        model.backward(error)
        model.update(lr)

        y2 = model.forward(x)
        print(f"{i} {x} {y} -> {y2}, target {target}, loss {loss}")

    print("--- weights ---")
    for layer in model.layers:
        for neuron in layer:
            print(neuron.weights, neuron.bias)

    print("--- result ---")
    for x, target in dataset:
        y = model.forward(x)
        print(f"{x} -> {y}, target {target}")


def main_not():
    model = BPModel(
        in_num=1,
        layers=[
            (4, Sigmoid()),
            (1, Sigmoid()),
        ],
    )
    dataset = [
        [[0], [1]],
        [[1], [0]],
    ]
    _simple_run(dataset, model, SSE(), 1000, lr=0.8)


def main_xor():
    model = BPModel(
        in_num=2,
        layers=[
            (4, ReLU()),
            (4, ReLU()),
            (1, Liner()),
        ],
    )
    dataset = [
        [[0, 0], [0]],
        [[1, 0], [1]],
        [[0, 1], [1]],
        [[1, 1], [0]],
    ]
    _simple_run(dataset, model, SSE(), 1000, lr=0.1)


def main_mnist():
    from keras.datasets import mnist

    (x_train, y_train), (x_test, y_test) = mnist.load_data()
    train_indices = np.where((y_train == 0) | (y_train == 1))[0]
    test_indices = np.where((y_test == 0) | (y_test == 1))[0]
    x_train = x_train[train_indices]
    y_train = y_train[train_indices]
    x_test = x_test[test_indices]
    y_test = y_test[test_indices]

    x_train = x_train.reshape(-1, 784) / 255
    x_test = x_test.reshape(-1, 784) / 255
    y_train = np.eye(10)[y_train]

    x_train = x_train[:1000]
    y_train = y_train[:1000]
    x_test = x_test[:10000]
    y_test = y_test[:10000]

    model = BPModel(
        in_num=784,
        layers=[
            (256, ReLU()),
            (64, ReLU()),
            (10, Liner()),
        ],
    )
    criterion = CrossEntropy()
    lr = 0.01

    for epoch in range(5):
        loss = 0
        for j in tqdm(range(len(x_train))):
            x = x_train[j].tolist()
            target = y_train[j].tolist()

            y = model.forward(x)
            loss += criterion(y, target)
            error = criterion.backward()
            model.backward(error)
            model.update(lr=lr)
        print(f"{epoch}epoch loss:{loss / len(x_train)}")

    # ---acc
    acc = 0
    for i in tqdm(range(len(x_test))):
        x = x_test[i].tolist()
        y = model.forward(x)
        predicts = np.argmax(y)
        if predicts == y_test[i]:
            acc += 1
    print("accuracy:", acc / len(x_test))

参考
・誤差逆伝播法をPythonで実装しながら理解する(Qiita)