自動微分

**自動微分**（じどうびぶん、）や**アルゴリズム微分**（）とは、プログラムで定義された関数を解析し、関数の値と同時に偏導関数の値を計算するアルゴリズムである。

自動微分は複雑なプログラムであっても加減乗除などの基本的な算術演算や基本的な関数（指数関数・対数関数・三角関数など）のような基本的な演算の組み合わせで構成されていることを利用し、これらの演算に対して合成関数の偏微分の連鎖律を繰り返し適用することによって実現される。自動微分を用いることで偏導関数値を少ない計算量で自動的に求めることができる。

## 他の微分方式との違い

自動微分は以下のどちらとも異なる。
* 記号微分・数式微分（symbolic differentiation） - 原関数を表す数式から数式処理により導関数を導出する。数式処理システムで実装されている。
* 数値微分（numerical differentiation） - 原関数の値から数値的に微分係数を算出する

記号微分は効率が悪くなりやすく、プログラムで定義された関数から微分表現を導くのは困難であるという問題がある。一方、数値微分では離散化の際の丸め誤差や桁落ちによる精度の低下が問題である。さらに、どちらの手法も計算量や誤差の関係で高次の微分係数を求めることが難しい。また、勾配を用いた最適化で必要となる、多くの入力変数を持つ関数に対する偏微分値の計算を行うには速度が遅い。自動微分はこれらの古典的手法の問題を解決する。

また、自動微分は計算フローを追いかけることで計算できるので、分岐（if文）やループや再帰を含むようなアルゴリズムでも偏微分できる。

## 合成関数の偏微分の連鎖律
自動微分の基本原理は、合成関数の偏微分の連鎖律を用いた偏微分の分解である。

合成関数の偏微分の連鎖律とは y = f(w_1, w_2), w_1 = g(x_1, x_2), w_2 = h(x_1, x_2) の時、下記が成立することである。
:
\frac{\partial y}{\partial x_1} = \frac{\partial y}{\partial w_1} \frac{\partial w_1}{\partial x_1} + \frac{\partial y}{\partial w_2} \frac{\partial w_2}{\partial x_1}


## 2種類の自動微分
自動微分は2種類に分けられ、それぞれ
* **ボトムアップ型自動微分**（**フォーワード・モード**、**フォーワード・アキュムレーション**、**タンジェント・モード**、**狭義の自動微分**）
* **トップダウン型自動微分**（**リバース・モード**、**リバース・アキュムレーション**、**随伴モード**、**高速自動微分**）
と呼ばれる。

ボトムアップ型自動微分では連鎖律を内側から外側に計算し（を計算した後で  を計算する）、トップダウン型自動微分では外側から内側に計算する。

使い分けは、入力が n 次元、出力が m 次元とした場合、以下の違いがある。
* n  m ならばトップダウン型の方が計算量が少ない。トップダウン型の計算回数はm回。
機械学習において、損失関数の出力はほぼ常に m = 1 のスカラーなので、トップダウン型が使われる。機械学習で用いられる多層パーセプトロンのバックプロパゲーションはトップダウン型自動微分の特殊なケースである。

ボトムアップ型はR.E. Wengertが1964年に発表したが、2ページの論文で特に名前を付けていない。Andreas Griewank によると、トップダウン型を誰が最初に発明したのか判然としないが、計算量を減らす工夫として1960年代後半には提案されていた。

## ボトムアップ型自動微分
ボトムアップ型自動微分では最初に偏微分を行う入力変数を固定し、それぞれの部分式を再帰的に計算する。手計算においては連鎖律において内側の関数を繰り返し置き換えていくことに相当する。

:\frac{\partial y}{\partial x}
= \frac{\partial y}{\partial w_1} \frac{\partial w_1}{\partial x}
= \frac{\partial y}{\partial w_1} \left(\frac{\partial w_1}{\partial w_2} \frac{\partial w_2}{\partial x}\right)
= \frac{\partial y}{\partial w_1} \left(\frac{\partial w_1}{\partial w_2} \left(\frac{\partial w_2}{\partial w_3} \frac{\partial w_3}{\partial x}\right)\right)
= \cdots

多変数の場合はヤコビ行列の積として一般化できる。

トップダウン型自動微分と比較すると、ボトムアップ型自動微分は自然であり、偏微分に関する情報の流れが計算の順序と一致するため簡単に実行できる。それぞれの変数にその偏導関数値 \dot w_i

:\dot w_i = \frac{\partial w_i}{\partial x}

の計算値を保存する領域を付け加えるだけで、変数値の計算と同時に偏導関数値を計算することができる。

連鎖律より、w_i が計算グラフで先行ノードを持つ場合、
:\dot w_i = \sum_j \frac{\partial w_i}{\partial w_j} \dot w_j, j ∈ {i の先行ノード}

例として次の関数を考える。
:\begin{align}
y
&= f(x_1, x_2) \\
&= x_1 x_2 + \sin x_1 \\
&= w_1 w_2 + \sin w_1 \\
&= w_3 + w_4 \\
&= w_5
\end{align}

それぞれの部分式を中間変数 w_i としてラベル付けした。最初に w_1 = x_1, w_2 = x_2 としている。

どの入力変数で偏微分するかによって \dot w_1 や \dot w_2 の初期値が異なる。x_1 に関する偏微分を計算する場合の初期値は、

: \begin{align}
\dot w_1 = \frac{\partial w_1}{\partial x_1} = \frac{\partial x_1}{\partial x_1} = 1 \\
\dot w_2 = \frac{\partial w_2}{\partial x_1} = \frac{\partial x_2}{\partial x_1} = 0
\end{align}

となる。初期値が決まったら下の表に示す連鎖律に従って各偏導関数値を順に計算していく。図2はこの処理を計算グラフとして表している。\dot w_5 = \tfrac{\partial y}{\partial x_1} が求めたい値である。

:\begin{array}{l|l}
\text{Operations to compute value} &
\text{Operations to compute derivative}
\\
\hline
w_1 = x_1 &
\dot w_1 = 1 \text{ (seed)}
\\
w_2 = x_2 &
\dot w_2 = 0 \text{ (seed)}
\\
w_3 = w_1 \cdot w_2 &
\dot w_3 = w_2 \cdot \dot w_1 + w_1 \cdot \dot w_2
\\
w_4 = \sin w_1 &
\dot w_4 = \cos w_1 \cdot \dot w_1
\\
w_5 = w_3 + w_4 &
\dot w_5 = 1 \cdot \dot w_3 + 1 \cdot \dot w_4
\end{array}

この関数 f に対する勾配を求めるためには \tfrac{\partial y}{\partial x_1} だけではなく \tfrac{\partial y}{\partial x_2} も求める必要がある。そのため、初期値を \dot w_1 = 0; \dot w_2 = 1 とした同様の計算を追加で行わなければならない。

勾配を求める場合に必要なボトムアップ型自動微分の実行回数は入力変数の個数と等しく、トップダウン型自動微分では出力変数の個数に等しい。ボトムアップ型やトップダウン型の自動微分を1回行うのに必要な計算量は、元のプログラムの計算量に比例する。そのため、偏微分する関数 が  を満たす場合、ボトムアップ型自動微分はトップダウン型自動微分よりも効率的である。

## トップダウン型自動微分
トップダウン型自動微分では、はじめに偏微分される出力変数を固定して、それぞれの部分式に関する偏導関数値を再帰的に計算する。手計算においては部分式を連鎖律における外側の関数の偏微分で繰り返し置き換えていくことに相当する。

:\frac{\partial y}{\partial x}
= \frac{\partial y}{\partial w_1} \frac{\partial w_1}{\partial x}
= \left(\frac{\partial y}{\partial w_2} \frac{\partial w_2}{\partial w_1}\right) \frac{\partial w_1}{\partial x}
= \left(\left(\frac{\partial y}{\partial w_3} \frac{\partial w_3}{\partial w_2}\right) \frac{\partial w_2}{\partial w_1}\right) \frac{\partial w_1}{\partial x}
= \cdots

トップダウン型自動微分において、求める値はボトムアップ型の随伴であり、上線 \bar w_i で表す。これは中間変数 w_i に関する偏微分
:\bar w_i = \frac{\partial y}{\partial w_i}
である。連鎖律より、w_i が計算グラフで後続ノードを持つ場合、以下のように計算できる。\bar y は1で、y以外のノードで後続ノードがなければ0になる。
:\bar w_i = \sum_j \bar w_j \frac{\partial w_j}{\partial w_i}, j ∈ {i の後続ノード}

トップダウン型自動微分は連鎖律を外側から内側（図3の計算グラフでは上から下）にたどっていく。

偏微分する関数 が を満たすとき、トップダウン型自動微分はボトムアップ型自動微分よりも効率的である。例示した関数は m = 1 のスカラー値関数（1つだけ値を返す関数）なので、勾配を計算するには計算グラフを一度たどるだけでよい。これは n = 2回計算グラフをたどる必要があったボトムアップ型自動微分の計算量の半分である。

しかし、トップダウン型自動微分は、テープやWengert リスト（なお、Wengert はボトムアップ型を1964年に発表したが、トップダウン型を発表したわけではない）と呼ばれるサイズ変更が可能な配列に中間変数 w_i の計算結果を追記していく必要があるためメモリ使用量の点で不利であり、計算グラフが巨大になるとメモリ使用量が問題となる可能性がある。この問題は保存する中間変数を限定し、必要な中間変数を再計算することによって軽減される。

トップダウン型自動微分を用いて偏導関数値を計算するための演算は以下の通りである（関数値を求める時と順番が逆であることに注意）。\bar w_1 = \tfrac{\partial y}{\partial x_1} と \bar w_2 = \tfrac{\partial y}{\partial x_2} が求めたい値である。
:\begin{array}{l}
\text{Operations to compute derivative}
\\ \hline
\bar w_5 = 1 \text{ (seed)}
\\
\bar w_4 = \bar w_5 \cdot 1
\\
\bar w_3 = \bar w_5 \cdot 1
\\
\bar w_2 = \bar w_3 \cdot w_1
\\
\bar w_1 = \bar w_3 \cdot w_2 + \bar w_4 \cdot \cos w_1
\end{array}


最後の \bar{w_1} は以下の理由で成立する。図3で x_1 = w_1 から w_3 と w_4 に矢印が引かれていることより、w_3 と w_4 は w_1 の後続ノードである。合成関数の偏微分の連鎖律より 
\bar{w_1} = \bar{w_3} \frac{\partial w_3}{\partial w_1} + \bar{w_4} \frac{\partial w_4}{\partial w_1}
 であり、後は w_3 と w_4 の定義を代入すれば良い。どの中間変数がどの中間変数に影響を及ぼすか記録する必要がある。

### テープを使用しない場合
テープを使用しない場合は、以下のように Python で実装できる。

```python
import math

class Var:
    def __init__(self, value, children=None):
        self.value = value
        self.children = children or []
        self.grad = 0

    def __add__(self, other):
        return Var(self.value + other.value, [(1, self), (1, other)])

    def __mul__(self, other):
        return Var(self.value * other.value, [(other.value, self), (self.value, other)])

    def sin(self):
        return Var(math.sin(self.value), [(math.cos(self.value), self)])

    def calc_grad(self, grad=1):
        self.grad += grad
        for coef, child in self.children:
            child.calc_grad(grad * coef)

# 例: f(x, y) = x * y + sin(x)
x = Var(2)
y = Var(3)
f = x * y + x.sin()

# 偏微分の計算
f.calc_grad()

print("f =", f.value)
print("∂f/∂x =", x.grad)
print("∂f/∂y =", y.grad)
```


### テープを使用する場合
テープを使用する場合のアルゴリズムの実装の流れは以下の通りである。
# テープをサイズ変更が可能な動的配列として用意する。テープの要素は中間変数 w_i 毎に作り、以下の5項目を記録する。
## 計算内容（足し算や掛け算など）
## 計算グラフの入辺：計算に使用した引数のテープ上のインデックスのリスト。w_3 なら w_1 と w_2。
## 計算グラフの出辺：この中間変数がどの計算で使われたかのテープ上のインデックスのリスト。計算しながら追記していく。w_1 なら w_3 と w_4。
## w_i の値
## \bar w_i の値（最初は未定）
# w_i の値を計算していく。その際、上記の内容をテープに追記していく。また、テープには、どの中間変数の計算で、どの中間変数を使用したかも追記する（上記の3番）。同一の中間変数を2回以上計算しないように注意が必要。
# テープの最後の \bar w_i を 1 として、その1つ前から逆順に \bar w_i = \sum_j \bar w_j \frac{\partial w_j}{\partial w_i} を計算していく。

### 多次元配列（テンソル）の場合
スカラー値ではなく NumPy のような多次元配列（テンソル）を扱う場合も処理すべき内容は同じである。\bar w_j \tfrac{\partial w_j}{\partial w_i} を1つの処理単位として実装する。これを vector-Jacobian product (VJP) や L-operator (Lop) と言う。行列積の VJP は行列積で実装できる。ただし、ある軸での和（NumPyのsum）の VJP は、その軸で値を繰り返して \bar w_j の軸を復元する処理になるなど、直観的で無い物もある。ハーバード大学の HIPS が開発していた Autograd は、ほぼ全ての NumPy の関数の VJP を実装していて、簡潔に実装しているので参考になる。なお \tfrac{\partial w_i}{\partial w_j} \dot w_j は Jacobian-vector product (JVP) や R-operator (Rop) と言い、ボトムアップ型で使用する。

#### 行列積
行列積の vector-Jacobian product (VJP) は、以下のように求める。Z = XY を考え、\bar{X} を求める。行列積を行列の要素で表記する。
: z_{ik} = \sum_j x_{ij} y_{jk}
行列積はかけ算して足し算をするが、足し算のトップダウン型自動微分は \bar{z} を上から降ろし、かけ算の自動微分は相手側の項を \bar{z} とかけるということより、\bar{x}_{ij} が求まる。
: \bar{x}_{ij} = \sum_k \bar{z}_{ik} y_{jk}
これを要素表記から行列表記に直すと、転置行列との行列積になる。
: \bar{X} = \bar{Z} Y^T
同様に、\bar{y}_{jk} = \sum_i x_{ij} \bar{z}_{ik} より \bar{Y} = X^T \bar{Z} である。

#### 2次元畳み込み
畳み込みニューラルネットワークのバイアス項無しの2次元畳み込み（Conv2d）は、Hは縦、Wは横、Cはチャンネルとして、入力を \mathbf{X}(H_\text{in},W_\text{in},C_\text{in})、カーネルを \mathbf{K}(H_\text{kernel},W_\text{kernel},C_\text{in},C_\text{out})、出力を \mathbf{Y}(H_\text{out},W_\text{out},C_\text{out}) とした場合、要素表記で下記式で表される。
: \mathbf{Y}(i,j,n) = \sum_{k,l,m} \mathbf{K}(k,l,m,n) \mathbf{X}(i+k,j+l,m)
行列積と同じくこれも積和演算なので、\bar{\mathbf{K}} は下記となる。
: \bar{\mathbf{K}}(k,l,m,n) = \sum_{i,j} \bar{\mathbf{Y}}(i,j,n) \mathbf{X}(i+k,j+l,m)
\bar{\mathbf{X}} は擬似コードを使い、\bar{\mathbf{X}}=0 から始めて、下記の6重ループで計算できる。
 **for** i,j,k,l,m,n **in** それぞれの値域
     \bar{\mathbf{X}}(i+k,j+l,m) += \mathbf{K}(k,l,m,n) \bar{\mathbf{Y}}(i,j,n)
数式で書くと \bar{\mathbf{X}}(p,q,m) = \sum_{i,j,k,l,n}\mathbf{K}(k,l,m,n)\bar{\mathbf{Y}}(i,j,n), \text{ where } i+k=p \land j+l=q となる。

バイアス項 \mathbf{B} は \mathbf{Y}(i,j,k) = \mathbf{X}(i,j,k) + \mathbf{B}(k) であるが、\bar{\mathbf{X}}(i,j,k) = \bar{\mathbf{Y}}(i,j,k) と \bar{\mathbf{B}}(k) = \sum_{i,j}\bar{\mathbf{Y}}(i,j,k) となる。

#### 2次元最大値プーリング
畳み込みニューラルネットワークのカーネルの大きさが(p,q)の2次元最大値プーリング（MaxPool2d）は要素表記では下記となるが、
: \mathbf{Y}(i,j,m) = \max_{k,l} \mathbf{X}(pi+k,qj+l,m)
これの \bar{\mathbf{X}} は、最大値をとる所に \bar{\mathbf{Y}} を戻せば良いので下記となる。argmax が複数の k, l で最大となる場合は1つ選ぶ。
: \bar{\mathbf{X}}(pi+k,qj+l,m) = \begin{cases}
\bar{\mathbf{Y}}(i,j,m) & \text{if }(k,l) = \mathop{\rm argmax}\limits_{k,l} \mathbf{X}(pi+k,qj+l,m) \\
0 & \text{otherwise}
\end{cases}

#### ReLU
畳み込みニューラルネットワークの活性化関数の ReLU は要素表記では \mathbf{Y}(i,j,k) = \max \{ 0, \mathbf{X}(i,j,k) \} だが、これの \bar{\mathbf{X}} は、\mathbf{X} \ge 0 の所に \bar{\mathbf{Y}} を戻せば良いので下記となる。
: \bar{\mathbf{X}}(i,j,k) = \begin{cases}
\bar{\mathbf{Y}}(i,j,k) & \text{if }\mathbf{X}(i,j,k) \ge 0 \\
0 & \text{otherwise}
\end{cases}

## 二重数を用いた自動微分
ボトムアップ型自動微分は実数の代数に（元を）添加して新しい算術を導入することによって可能である。全ての数（通常の実数）に対して、その数における関数の微分を表現する追加の成分が足され、全ての算術演算がこの添加代数に拡張される。すなわち二重数の代数である。このアプローチはの理論（つまり双対空間のテンソル代数）によって一般化される（を見よ)。

各数 x を数 x + x'\varepsilon に置き換える。ここで x' は実数だが、\varepsilon は \varepsilon^2=0 を満たすである（無限小；滑らかな無限小解析も参照）。ちょうどこれだけを用いて通常の演算が得られる：
:\begin{align}
      (x + x'\varepsilon) + (y + y'\varepsilon) &= x + y + (x' + y')\varepsilon \\
  (x + x'\varepsilon) \cdot (y + y'\varepsilon) &= xy + xy'\varepsilon + yx'\varepsilon + x'y'\varepsilon^2 = xy + (x y' + yx')\varepsilon
\end{align}

引き算と割り算についても同様である。

いまやこの拡張算術のもとで多項式を計算できる。もし P(x) = p_0 + p_1 x + p_2x^2 + \cdots + p_n x^n ならば、

\begin{align}
  P(x + x'\varepsilon) &= p_0 + p_1(x + x'\varepsilon) + \cdots + p_n (x + x'\varepsilon)^n \\
                       &= p_0 + p_1 x + \cdots + p_n x^n + p_1x'\varepsilon + 2p_2xx'\varepsilon + \cdots + np_n x^{n-1} x'\varepsilon \\
                       &= P(x) + P^{(1)}(x)x'\varepsilon
\end{align}

ここで P^{(1)} は P の最初の引数に対する微分であり、x'（種と呼ばれる）は任意に取れる。

上に述べたように、この新しい算術は、順序対（\langle x, x' \rangle と書かれる元）と、最初の成分に対しては通常の算術を、第二の成分に対しては一階微分の算術を、それぞれ与えたものからなる。多項式に関する上の結果を解析関数に広げれば、新しい算術に対する、基本的な算術と幾つかの標準的な関数のリストが得られる：
:\begin{align}
  \left\langle u,u'\right\rangle + \left\langle v,v'\right\rangle &= \left\langle u + v, u' + v' \right\rangle \\
  \left\langle u,u'\right\rangle - \left\langle v,v'\right\rangle &= \left\langle u - v, u' - v' \right\rangle \\
  \left\langle u,u'\right\rangle * \left\langle v,v'\right\rangle &= \left\langle u v, u'v + uv' \right\rangle \\
  \left\langle u,u'\right\rangle / \left\langle v,v'\right\rangle &= \left\langle \frac{u}{v}, \frac{u'v - uv'}{v^2} \right\rangle \quad ( v\ne 0) \\
                               \sin\left\langle u,u'\right\rangle &= \left\langle \sin(u) , u' \cos(u) \right\rangle \\
                               \cos\left\langle u,u'\right\rangle &= \left\langle \cos(u) , -u' \sin(u) \right\rangle \\
                               \exp\left\langle u,u'\right\rangle &= \left\langle \exp u , u' \exp u \right\rangle \\
                               \log\left\langle u,u'\right\rangle &= \left\langle \log(u) , u'/u \right\rangle \quad (u>0) \\
                                 \left\langle u,u'\right\rangle^k &= \left\langle u^k , k u^{k - 1} u' \right\rangle \quad (u \ne 0) \\
                    \left| \left\langle u,u'\right\rangle \right| &= \left\langle \left| u \right| , u' \mbox{sign} u \right\rangle \quad (u \ne 0)
\end{align}
一般に、プリミティヴの関数 g に対して、
:g(\langle u,u' \rangle , \langle v,v' \rangle ) = \langle g(u,v) , g_u(u,v) u' + g_v(u,v) v' \rangle
ここで g_u と g_v はそれぞれ g の最初と2番目の引数に対する微分である。

基本的な二項算術演算を（実数と二重数の）混在した引数に対して、つまり順序対 \langle u, u' \rangle と実数 c に対して適用するとき、まずこの実数を \langle c, 0 \rangle に引き上げる（lifting）。関数 f : \mathbb{R}\rightarrow\mathbb{R} の x_0 に於ける微分はいまや f(\langle x_0, 1 \rangle) に上の算術を使って計算することによって得られる。これは \langle f ( x_0 ) , f' ( x_0 ) \rangle  を結果として与える。

### ベクトル引数と関数
多変数関数は、方向微分作用素を用いることで、一変数関数の場合と同様の効率と仕組みで取り扱える。つまり、f:\mathbb{R}^n\rightarrow\mathbb{R}^m の x \in \mathbb{R}^n における x' \in \mathbb{R}^n 方向微分 y' = \nabla f(x)\cdot x' を計算したい場合、それは (\langle y_1,y'_1\rangle, \ldots, \langle y_m,y'_m\rangle) = f(\langle x_1,x'_1\rangle, \ldots, \langle x_n,x'_n\rangle) を上と同様の算術を使って計算すればよい。もし \nabla f の全ての要素が望みならば、n 個の関数評価が要求される。ここで、多くの最適化アプリケーションでは、実際には方向微分があれば十分である。

### 高階・多変数
上の算術は多変数関数の二階やもっと高階の微分の計算の為に一般化出来る。しかし、その算術規則は直ちに極めて複雑なものとなる：複雑性は最高次の微分の次数に対して二次関数的となる。その代わりに、途中で打ち切った（truncated）テイラー多項式の代数を使用できる。結果得られる算術（一般化された二重数の上で定義された）は、関数を新しいデータ型であるかのように使って、効率的に計算することを可能にする。ひとたび関数のテイラー多項式が分かれば、その導関数たちは容易に抽出できる。
厳密で一般的な定式化はを通してを用いることにより達成される。

## 実装
自動微分の実装方法には大きく分けて、ソースコードの変換とオペレータオーバーローディングによる方法の2つがある。

### ソースコード変換

関数値を求める関数を記述した元のソースコードから、偏導関数値を計算する処理を含んだプログラムを自動的に生成する手法である。ソースコード変換はあらゆるプログラミング言語で実装でき、コンパイル時の最適化を行いやすいが、自動微分ツールの作成は難しい。

### オペレータオーバーローディング

この手法は演算子のオーバーロードがサポートされているプログラミング言語で記述されたソースコードに対してのみ適用可能である。元のソースコードの流れを大きく変更することなく実現できるが、基本データ型の変更などの小さな変更は必要である。

ボトムアップ型自動微分をオペレータオーバーロードで実現するのは容易である。トップダウン型自動微分についても可能であるが、現状のコンパイラではボトムアップ型自動微分と比べると最適化の面で不利である。

## ソフトウェア
自動微分を実装したライブラリなどのソフトウェアが多数存在する。2010年代の第3次人工知能ブームの際にディープラーニングに自動微分が必要なため、TensorFlowやPyTorchなどトップダウン型の自動微分を含むライブラリが多数作られた。

## 脚注

## 参考文献
* 
* 伊理正夫、「高速自動微分法(第2回年会特別講演)」『応用数理』 1993年 3巻 1号 p.58-66, , 日本応用数理学会
* 
* 
*  

## 外部リンク
* www.autodiff.org - Community Portal for Automatic Differentiation