物理を習ったことがある人向けの誤差逆伝播法

機械学習のニューラルネットワークにおける誤差逆伝播について、物理を習ったことがある人にわかりやすいようにまとめてみました。
なお、物理を習ったことがある人、というのは未定義で、「なんとなく物理をやったことがあるならわかりやすいかな」という程度です。
物理をやったことがあるなら、偏微分とか連鎖律とかはよく知っているわけです。また、あまりに煩雑な計算はうまいこと避けたいなとか思っていると思います。そして、変に図を描くよりも式変形した方がわかりやすいなとか思ったりすると思います。

参考文献
http://nnadl-ja.github.io/nnadl_site_ja/chap2.html

追記：全結合型以外の一般的な形の時の誤差逆伝播法についても追記しました。
追記：$\delta$の定義をより自然にした方法を追記しました。

# ニューラルネットワーク
まず、ニューラルネットワークについて考えます。
ある長さnのベクトルxがあるとします。そして、出力として、スカラーyがあるとします。この時、両者の関係が知りたいとします。つまり、

```math
y = f(x)
```
という関係式のfが知りたい関数です。
yがスカラー、xがベクトルなので、一番簡単な関数は

```math
y = w x + b
```
でしょう。wは1xnの行列です。
関数が知りたい、ということは、wとbを決めればよい、ということになります。
しかし、この関数はあまりにも単純です。
ですので、

```math
y = w_2 \sigma (w_1 x + b_1) + b_2
```
という関数を考えます。
今度は、w1はmxnの行列、b1はm次元ベクトル、w2は1xmの行列となります。mは自分で選びます。mが大きいと行列のサイズが大きくなるのでパラメータの数が増えます。ということで関数の表現力が上がることになります。
ここでσ(x)はベクトルxのそれぞれの要素に作用させる関数で、x = (2,3)などであれば、σ(x) = (σ(2),σ(3))となる感じです。σ(x)は非線形関数を使います。例えばtanh(x)などを使います。これを活性化関数と呼んでいます。

上の式をもう少し複雑にすることを考えます。
つまり、

```math
y = w_3 \sigma (w_2 \sigma (w_1 x + b_1) + b_2) + b_3
```
とします。
以下同文でどんどん複雑にできます。
最後に非線形関数をかけてもいいので、例えば、

```math
y = \sigma(w_3 \sigma (w_2 \sigma (w_1 x + b_1) + b_2) + b_3)
```
でも良いでしょう。

このままどんどん書いていっても大変なので、漸化式を使って書いておきます。
つまり、

```math
\begin{align}
a^l &= \sigma(w^l a^{l-1} + b^l) \\
y &= a^L \\
a^0 &= x
\end{align}
```
とします。
これで、L回同じことを繰り返したことを表現できていますね。

さらに、

```math
z^l = w^l a^{l-1} + b^l
```
というものを定義しておくと、

```math
\begin{align}
y(x) &= \sigma(z^L) \\
a^0 &= x
\end{align}
```
と簡単にかけます。この形式をL層のニューラルネットワークと呼びます。
それぞれ途中に出てくるベクトルの各要素をユニットと呼びます。$a^l_j$であれば、l層のj番目のユニットです。

# 最適化
インプットをベクトルx、アウトプットをスカラーyとする関数の表現方法を手に入れました。
次は、どんな関数を作るか、について考えます。
あるベクトルxがあった時、出力としてyoutを得たとします。この時、

```math
y = \sigma(z^l) \sim y_{out}
```
となるような関数が欲しいとします。
一番単純そうなのは、

```math
C = \sum_x (y_{out}(x) - y(x))^2
```
が最小となるような関数を作ることでしょう。もし完全に一致していればCはゼロになりますので。
xは色々なベクトルにして、たくさんのインプットベクトルxに対して対応するアウトプットyoutを集めてきて、それらが対応するように関数y(x)を作る、というわけです。

関数y(x)のパラメータはwやbです。ですので、Cを最小とするためには、パラメータ微分：

```math
\begin{align}
\frac{\partial C}{\partial w_{ij}^l} \\
\frac{\partial C}{\partial b_{j}^l} 
\end{align}
```
がゼロとなるようなwやbを持って来ればいいですね。
この微分を計算することが目的です。

# 誤差逆伝播法
Cはwやbに依存していますので、頑張れば微分できそうです。
しかし、むやみやたらに頑張りたくはないので、少し工夫します。

まず、

```math
\delta_j^l \equiv \frac{\partial C}{\partial z_j^l}
```
というl層でのj番目のユニットからくる誤差、というものを定義しておきます。

```math
z^l = w^l a^{l-1} + b^l
```
ということを思い返すと、パラメータbに関する微分は

```math
\begin{align}
\frac{\partial C}{\partial b_{j}^l} = \sum_k \frac{\partial z_k^l}{\partial b_{j}^l} \frac{\partial C}{\partial z_k^l} \\
= \sum_k \delta_{kj} \frac{\partial C}{\partial z_k^l} = \delta_j^l
\end{align}
```
となります。ここで$\delta_{kl}$はクロネッカーのデルタです。

パラメータwに関する微分も同様に、

```math
\begin{align}
\frac{\partial C}{\partial w_{ij}^l} = \sum_k \frac{\partial z_k^l}{\partial w_{ij}^l} \frac{\partial C}{\partial z_k^l} \\
= \sum_k \frac{\partial C}{\partial z_k^l} \delta_{ki} a^{l-1}_j = \delta_i^l a_j^{l-1}
\end{align}
```
となりますので、$\delta_j^l$さえなんとか計算できれば計算できそうです。

## 誤差逆伝播
まず、最後の層であるL番目の層に着目し、$\delta_j^L$を計算してみましょう。

```math
\delta_j^L = \sum_k \frac{\partial a_k^L}{\partial z_j^L} \frac{\partial C}{\partial a_k^L}
```
ここで、

```math
a^L = \sigma(z^L)
```
ですから、ベクトル$a^L$のj番目の要素はzjにしか依存しません。よって、

```math
\begin{align}
\delta_j^L = \sum_k \delta_{kj} \frac{\partial a_k^L}{\partial z_j^L} \frac{\partial C}{\partial a_k^L} = \frac{\partial \sigma(z_j^L)}{\partial z_j^L}\frac{\partial C}{\partial a_j^L} \\
=  \frac{\partial C}{\partial a_j^L}  \sigma'(z_j^L)
\end{align}
```
と書けます。つまり、活性化関数σの微分の値です。活性化関数は解析的な形を与えていますので、その微分は簡単に計算できます。また、最小化したい関数Cは、最終的な出力の関数ですから、$a^L$の関数として解析的な形でかけており、微分は簡単に計算できます。

残りは、l層での$\delta^l_j$の値です。最終層での値が分かっているので、漸化式的な形で求められれば、計算できることになります。よって、

```math
\begin{align}
\delta^l_j = \sum_k \frac{\partial z_k^{l+1}}{\partial z_j^l} \frac{\partial C}{\partial z_k^{l+1}} \\
=  \sum_k \frac{\partial z_k^{l+1}}{\partial z_j^l} \delta^{l+1}_k
\end{align}
```
となり、さらに、

```math
z^l = w^l a^{l-1} + b^l = w^l \sigma(z^{l-1}) + b^l
```
を用いれば、

```math
\frac{\partial z_k^{l+1}}{\partial z_j^l} = w^l_{kj} \sigma'(z_j^l)
```
となり、

```math
\delta^l_j = \sigma'(z_j^l) \sum_k w^{l+1}_{kj}  \delta^{l+1}_k
```
が得られます。
これで、$\delta_j^L$から遡る形で$\delta_j^j$が計算できるので、Cをwやbで微分できることになります。これを、誤差逆伝播法と呼びます。

# 一般化
物理を習ったことがあるなら、やはり一般化はしてみたいところです。
というわけで、

```math
\begin{align}
y(x) &= \sigma(z^L) \\
a^0 &= x \\
a^l &= \sigma(z^l) \\
z^l &= w^l a^{l-1} + b^l
\end{align}
```
を一般化してみます。
これはl層の値はl-1層の値から作られることを意味していて、パラメータがwとbです。
ここで、

```math
z^l = g_{W^l} (a^{l-1}) 
```
としてしまいましょう。パラメータは$W^l$としておきます。このWは複数の要素を持っているとします。ベクトルとか行列とかですね。今は簡単のため、$W^l_n$とします。l層でのn番目のパラメータです。
つまり、関数gはWをパラメータとしているわけです。
よって、何らかの関数を最小化したいのであれば、パラメータWで微分すれば良いことになります。これは、

```math
\begin{align}
\frac{\partial C}{\partial W^l_n} &= \sum_k \frac{\partial z_k^l}{\partial W_n^l} \frac{\partial C}{\partial z_k^l} \\
&= \sum_k \frac{\partial [g_{W^l} (a^{l-1})]_k}{\partial W_n^l} \delta^l_k
\end{align}
```
と書けます。gは自分で設定した関数で微分が計算できますし、$\delta_k^l$は前と同じ定義です。

よって、L層での値は前と同じで

```math
\begin{align}
\delta_j^L = \sum_k \frac{\partial a_k^L}{\partial z_j^L} \frac{\partial C}{\partial a_k^L} = \frac{\partial C}{\partial a_j^L}  \sigma'(z_j^L)
\end{align}
```
です。

l層での値は、途中まで同じで、

```math
\begin{align}
\delta^l_j 
=  \sum_k \frac{\partial z_k^{l+1}}{\partial z_j^l} \delta^{l+1}_k \\
= \sum_k \frac{\partial  [g_{W^l} (a^{l})]_k }{\partial z_j^l} \delta^{l+1}_k
\end{align}
```
となります。さらに、$a^l = \sigma(z^l) $を使えば、

```math
\begin{align}
\frac{\partial  [g_{W^l} (a^{l})]_k }{\partial z_j^l} = 
\frac{\partial  \sigma(z_j^l)}{\partial z_j^l} \frac{\partial  [g_{W^l} (\sigma(z^l))]_k }{\partial \sigma(z_j^l)} \\
= \frac{\partial  [g_{W^l} (\sigma(z^l))]_k }{\partial \sigma(z_j^l) } \sigma'(z_j^l) 
\end{align}
```
なので、

```math
\begin{align}
\delta^l_j 
= \sigma'(z_j^l) \sum_k \frac{\partial  [g_{W^l} (\sigma(z^l))]_k }{\partial \sigma(z_j^l)} \delta^{l+1}_k \\
= 
\sigma'(z_j^l) \sum_k \frac{\partial  [g_{W^l} (a^l)]_k }{\partial a_j^l} \delta^{l+1}_k 
\end{align}
```
となります。
このように書くと、l層とl+1層の関係がどんなものであっても、誤差逆伝播法が使えることになります。
ニューラルネットワークの言葉で言えば、全結合型ではないニューラルネットワークもこれで微分が計算できます。

## 確認
確認してみましょう。
一般化する前の式であれば、

```math
[g_{W^l} (a^l))]_k = \sum_j w^{l+1}_{kj} a^{l}_j + b^{l+1}_k
```
なので、

```math
\frac{\partial  [g_{W^l} (a^l)]_k }{\partial a_j^l} = w^{l+1}_{kj}
```
となりますので、ちゃんと同じ結果になりますね。

# 2023年追記：少しモヤる人のために
この記事を書いてから数年が経ちましたが、上の$\delta$が唐突に登場することにモヤモヤする人がいるかもしれません。というのは、少し天下りに定義しているようにみえるからです。ですので、少し追記することにします。
まず、計算したい量が何かをはっきりしておきます。まず、インプットを$\vec{x}$、アウトプットを$\vec{y}$とした時、$\vec{y}$は
```math
\begin{align}
\vec{a}^0(\vec{x}) &= \vec{x} \\
\vec{z}^l(\vec{a}^{l-1}) &= g_{\vec{W}^1}(\vec{a}^{l-1}) \\
\vec{a}^l(\vec{z}^l) &= \sigma^l(\vec{z}^l) \\
\vec{y} &= \vec{a}^L(\vec{z}^L)
\end{align}
```
という$L$層のニューラルネットワークで書けているとします。ここで、$z^l$と$z^{l-1}$の関係は上の式を用いると、
```math
\vec{z}^l = g_{\vec{W}^1}(\sigma^{l-1}(\vec{z}^{l-1}))
```
となっています。
そして、最小化したい量$C$は
```math
C(\vec{x}) = C(\vec{y})
```
です。見ての通り、関数$C$はインプットベクトル$\vec{x}$の関数です。そこで、インプットベクトル$\vec{x}$の要素での微分を考えます。
```math
\frac{\partial C}{\partial x_n} 
```
ここで、$x_n$で解析的に微分ができるのは$\vec{a}^0(\vec{x}) = \vec{x}$ですので、
```math
\frac{\partial C}{\partial x_n} = \frac{\partial C}{\partial a^0_n}\frac{\partial a^0_n}{\partial x_n}
```
とすれば、二つの積の一つ$\frac{\partial a^0_m}{\partial x_n}$は解析的に微分ができます。次に、$\frac{\partial C}{\partial a^0_n}$の$a^0_n$微分が解析的にできるのは$\vec{z}^1(\vec{a}^{0})$ですので、
```math
\frac{\partial C}{\partial a^0_n} = \sum_{m^1} \frac{\partial C}{\partial z^1_{m^1}} \frac{\partial z^1_{m^1}}{\partial a^0_n}
```
とすると$\frac{\partial z^1_{m^1}}{\partial a^0_m}$は解析的に微分できます。そして、
```math
\frac{\partial C}{\partial z^1_{m^1}} = \sum_{m^2} \frac{\partial C}{\partial z^2_{m^2}} \frac{\partial z^2_{m^2}}{\partial z^1_{m^1}} 
```
とすると、
```math
\begin{align}
\frac{\partial z^2_{m^2}}{\partial z^1_{m^1}}  = \frac{\partial g_{\vec{W}^1}(\sigma^{1}(\vec{z}^{1}))　}{\partial z^1_{m^1}} =  \frac{\partial g_{\vec{W}^1}(\sigma^{1}(\vec{z}^{1}))　}{\partial \sigma^{1}} \frac{\partial \sigma^{1}(\vec{z}^{1})}{\partial z^1_{m^1}} \\
= \frac{\partial g_{\vec{W}^1}(a^{1})　}{\partial a^{1}} \frac{\partial \sigma^{1}(\vec{z}^{1})}{\partial z^1_{m^1}}
\end{align}
```
となるので、解析的に微分を実行できます。そして、
```math
\frac{\partial C}{\partial z^2_{m^2}} = \sum_{m^3} \frac{\partial C}{\partial z^3_{m^3}} \frac{\partial z^3_{m^3}}{\partial z^2_{m^2}} 
```
となりますから、これを
```math
\frac{\partial C}{\partial z^l_{m^l}} = \sum_{m^{l+1}} \frac{\partial C}{\partial z^{l+1}_{m^{l+1}}} \frac{\partial z^{l+1}_{m^{l+1}}}{\partial z^l_{m^l}} 
```
とどんどん繰り返すことができます。そして、最後には、
```math
\frac{\partial C}{\partial z^L_{m^L}} =  \sum_{m} \frac{\partial C}{\partial y_{m}} \frac{\partial y_{m}}{\partial z^L_{m^L}} 
```
となります。ここで、$\frac{\partial C}{\partial y_{m}} $も$\frac{\partial y_{m}}{\partial z^L_{m^L}} $の両方とも解析的に微分を実行できるようになりましたので、全ての量を解析的に微分できるようになったことになります。つまり、$\frac{\partial C}{\partial x_n} $を解析に微分を実行できています。計算の流れとしては、
```math
\frac{\partial C}{\partial y_{m}},\frac{\partial y_{n}}{\partial z^L_{m^L}}  \rightarrow \frac{\partial C}{\partial z^L_{m^L}} \rightarrow \frac{\partial C}{\partial z^L_{m^{L-1}}} \rightarrow \cdots   \rightarrow \frac{\partial C}{\partial z^1_{m^1}} \rightarrow \frac{\partial C}{\partial a^0_m} \rightarrow \frac{\partial C}{\partial x_n} 
```
です。$L$番目の値から逆に進んでいくので、逆伝播法と呼ばれています。
以上のような計算では、重要な量は
```math
\delta^l_m \equiv \frac{\partial C}{\partial z^l_{m}} 
```
であることがわかります。そして、上で述べているように、パラメータ微分も$\delta^l_m$を使って計算できます。

## 自動微分との関連
自動微分に関する説明は例えば[Juliaでの例](https://qiita.com/cometscome_phys/items/4940eddaaa9086169ef1)の記事があります。この記事との関係を述べておきます。
まず、
```math
\frac{\partial C}{\partial z^l_{m^l}} = \sum_{m^{l+1}} \frac{\partial C}{\partial z^{l+1}_{m^{l+1}}} \frac{\partial z^{l+1}_{m^{l+1}}}{\partial z^l_{m^l}} 
```
の式は$\delta$を用いて書くと、
```math
\delta^l_{m^l} = \sum_{m^{l+1}} \delta^{l+1}_{m^{l+1}} \frac{\partial z^{l+1}_{m^{l+1}}}{\partial z^l_{m^l}} 
```
ですが、これは、$\delta^l_{m^l}$が
```math
\delta^l_{m^l} = K^{l+1}(\vec{\delta}^{l+1})
```
という関数$K$で書けていることを意味しています。これは記事でのpullbackそのものです。層ごとに関数$K$が実装されていれば、どんな層でも自動微分が可能になります。例えばJuliaでは[この記事](https://qiita.com/cometscome_phys/items/683bfb73160c54366636)に任意のpullbackの実装方法が書いてありますし、PyTorchの場合には、[公式ページのPYTORCH: DEFINING NEW AUTOGRAD FUNCTIONS](https://pytorch.org/tutorials/beginner/examples_autograd/two_layer_net_custom_function.html)にあります。