過去の記事：

• 自動微分(0) — 導入
• 自動微分(1) — 1変数の場合
• 自動微分(2) — 偏微分係数への一般化
• 自動微分(3) — 高階導関数の計算

これまで数回に渡って自動微分を取り上げてきたが、それでは、自動微分は，数学的に見るとどのような概念と対応しているのだろうか？

簡単のため，ここでは実数 \(\mathbf{R}\) についての関数の、1変数による1階微分の自動微分のみを考えるものとする．1変数関数の自動微分を表す型を \(\mathrm{AutoDiff}(\mathbf{R})\) とする．\(\mathrm{AutoDiff}(\mathbf{R})\) には，\(\mathbf{R}\) から誘導された加法や乗法などの演算が入っている．また，\(\mathbf{R}^n\to\mathbf{R}\) の微分可能な関数から，\(\mathrm{AutoDiff}(\mathbf{R}^n)\to\mathrm{AutoDiff}(\mathbf{R})\) の関数が誘導される．加法や乗法などの演算は，\(\mathbf{R}^2\to\mathbf{R}\) の微分可能な関数と見ることができるので，関数の誘導だけを考えれば良い．

つまり，自動微分を考えるということは，実数 \(\mathbf{R}\) あるいはユークリッド空間 \(\mathbf{R}^n\) について，それぞれ \(\mathrm{AutoDiff}(\mathbf{R})\), \(\mathrm{AutoDiff}(\mathbf{R}^n)\) という空間を考えるということである．そして，(微分可能な)関数 \(\mathbf{R}^n\to\mathbf{R}^m\) に対しては，誘導された関数 \(\mathrm{AutoDiff}(\mathbf{R}^n)\to\mathrm{AutoDiff}(\mathbf{R}^m)\) が存在する．
これは圏論の言葉で言う「関手」(functor)によく似ている．しかし，\(\mathrm{AutoDiff}\) を圏の関手として見る場合，対象とする空間をユークリッド空間 \(\mathbf{R}^n\) に制限しているのは気持ち悪い．1階微分が定義できるようなもっと広い空間に対して，自動微分を一般化して考えられないだろうか？

数学では，ユークリッド空間を一般化したような空間を「多様体」(manifold)と呼ぶ．特に，「滑らか」で「関数の微分が定義できる」ものを「滑らかな多様体」とか「可微分多様体」と呼ぶ．多様体についての詳細はここでは省くので、数学の授業を受けるか、専門書を参照するかして欲しい。何が言いたいかというと，自動微分のようなものを，可微分多様体の圏から可微分多様体の圏への関手として定義できるはずである．

実は，この概念はすでに多様体論ではよく知られていて，「接バンドル」(tangent bundle)と呼ばれている．\(\mathrm{AutoDiff}(\mathbf{R}^n)\) に対応する空間は \(T\mathbf{R}^n\) と書かれる．一般に，多様体 \(M\) の接バンドルは \(TM\) と表される．写像の方は，\(f\colon M\to N\) から誘導される写像は \(df\colon TM\to TN\) と書かれる．

\(\mathbf{R}\) の自動微分、\(\mathrm{AutoDiff}(\mathbf{R})\) は \(\mathbf{R}\times\mathbf{R}\) のような直積で表された。このような場合、 \(T\mathbf{R}\) は自明(trivial)であるというような言い方をする。また、\(\mathrm{AutoDiff}(\mathbf{R}^n)\cong\mathrm{AutoDiff}(\mathbf{R})^n\)（自然な同型）となるが、これは接バンドルについて \(T(M\times N)\cong TM\times TN\)（自然な同型）となることに対応する。

筆者の知識が足らないので今書けるのはここまでだが、要するにプログラミングのテクニックの一つである自動微分は、幾何学の概念と対応していそうだということである。あとそれから、幾何ガチ勢からまさかりが飛んできそうで怖い。

前回まででは、1階の(偏)微分係数を考察した。今度は、高階の導関数の値(微分係数)も計算するように拡張しよう。簡単のため、1変数関数について考えることにする。

高階微分と演算、関数の合成との関係を見る。変数は \(t\) とする。

もとの関数	1階導関数、2階導関数
\(h(t)=\mathrm{const.}\)	\(h'(t)=0,\quad h^{\prime\prime}(t)=0\)
\(h(t)=t\)	\(h'(t)=1,\quad h^{\prime\prime}(t)=0\)
\(h(t)=f_1(t)\pm f_2(t)\)	\(h'(t)=f_1′(t)\pm f_2′(t),\quad h^{\prime\prime}(t)=f_1^{\prime\prime}(t)\pm f_2^{\prime\prime}(t)\)
\(h(t)=f_1(t)f_2(t)\)	\(\begin{array}{c} h'(t)=f_1′(t)f_2(t)+f_1(t)f_2′(t), \\ h^{\prime\prime}(t)=f_1^{\prime\prime}(t)f_2(t)+2f_1′(t)f_2′(t)+f_1(t)f_2^{\prime\prime}(t) \end{array}\)
\(h(t)=\frac{f_1(t)}{f_2(t)}\)	\(\begin{array}{c} h'(t)=\frac{f_1′(t)}{f_2(t)}-\frac{f_1(t)f_2′(t)}{f_2(t)^2}, \\ h^{\prime\prime}(t)=\frac{f_1^{\prime\prime}(t)}{f_2(t)}-\frac{2f_1′(t)f_2′(t)}{f_2(t)^2}-\frac{f_1′(t)f_2^{\prime\prime}(t)}{f_2(t)^2}+\frac{2f_1(t)f_2′(t)^2}{f_2(t)^3} \end{array}\)
\(h(t)=g(f(t))\)	\(\begin{array}{c} h'(t)=g'(f(t))f'(t), \\ h^{\prime\prime}(t)=g^{\prime\prime}(f(t))f'(t)^2+g'(f(t))f^{\prime\prime}(t) \end{array}\)
\(h(t)=g(f_1(t),f_2(t))\)	\(\begin{array}{c} h'(t)=g_x(f_1(t),f_2(t))f_1′(t)+g_y(f_1(t),f_2(t))f_2′(t), \\ h^{\prime\prime}(t)=g_{xx}(f_1(t),f_2(t))f_1′(t)^2+2g_{xy}(f_1(t),f_2(t))f_1′(t)f_2′(t)+g_{yy}(f_1(t),f_2(t))f_2′(t)^2 \\ \quad +g_x(f_1(t),f_2(t))f_1^{\prime\prime}(t)+g_y(f_1(t),f_2(t))f_2^{\prime\prime}(t) \end{array}\)

2階微分の自動微分では、「値」「1階の微分係数」「2階の微分係数」の3つの値を保持する必要があるので、実装するには \(\mathbf{R}\times\mathbf{R}\times\mathbf{R}\) に適切な演算を入れれば良い。上の表の \(f(t),f'(t),f^{\prime\prime}(t)\) を \(\mathbf{R}\times\mathbf{R}\times\mathbf{R}\) の成分で置き換えると、\(\mathbf{R}\times\mathbf{R}\times\mathbf{R}\) に入れるべき演算が分かる。すなわち、\begin{align*}
(x,x’,x^{\prime\prime})\pm(y,y’,y^{\prime\prime})&=(x\pm y,x’\pm y’,x^{\prime\prime}\pm y^{\prime\prime}) \\
(x,x’,x^{\prime\prime})\cdot(y,y’,y^{\prime\prime})&=(xy,x’y+xy’,x^{\prime\prime}y+2x’y’+xy^{\prime\prime}) \\
\frac{(x,x’,x^{\prime\prime})}{(y,y’,y^{\prime\prime})}&=\left(\frac{x}{y},\frac{x’}{y}-\frac{xy’}{y^2},\frac{x^{\prime\prime}}{y}-\frac{2x’y’}{y^2}-\frac{x’y^{\prime\prime}}{y^2}+\frac{2xy’^2}{y^3}\right),
\end{align*}\(g\colon\mathbf{R}\to\mathbf{R}\)(微分可能) に対し、\[g_*((x,x’,x^{\prime\prime}))=(g(x),g'(x)x’,g^{\prime\prime}(x)x’^2+g'(x)x^{\prime\prime}),\]\(g\colon\mathbf{R}\times\mathbf{R}\to\mathbf{R}\)(微分可能) に対し、\[g_*((x,x’,x^{\prime\prime}),(y,y’,y^{\prime\prime}))=\left(\begin{array}{c}g(x,y), \qquad
g_x(x,y)x’+g_y(x,y)y’, \\
g_{xx}(x,y)x’^2+2g_{xy}(x,y)x’y’+g_{yy}(x,y)y’^2+g_x(x,y)x^{\prime\prime}+g_y(x,y)y^{\prime\prime}\end{array}\right),\]
とすればよい。ただし、変数についている \({}^\prime\) は微分ではなく変数名の一部と考える。関数についている \({}^\prime\) や添字は微分の意である。

これを実際のプログラミング言語で実装するのは容易だろう。

前回まで考察したのは，1変数関数とその導関数の値を計算する場合だった．今回は，2変数関数とその偏導関数を計算する場合について考える．

自動微分を考える前に、2変数関数の微分(偏微分)と演算、関数の合成との関係を見ておこう。2つの変数を \(u\), \(v\) とおく．偏導関数は，変数を添字にして表すことにする．

もとの関数	偏導関数
\(h(u,v)=\mathrm{const.}\)	\(h_u(u,v)=0,\quad h_v(u,v)=0\)
\(h(u,v)=u\)	\(h_u(u,v)=1,\quad h_v(u,v)=0\)
\(h(u,v)=v\)	\(h_u(u,v)=0,\quad h_v(u,v)=1\)
\(h(u,v)=f_1(u,v)\pm f_2(u,v)\)	\(\begin{aligned} h_u(u,v)&=f_{1u}(u,v)\pm f_{2u}(u,v), \\ h_v(u,v)&=f_{1v}(u,v)\pm f_{2v}(u,v) \end{aligned}\)
\(h(u,v)=f_1(u,v)f_2(u,v)\)	\(\begin{aligned} h_u(u,v)&=f_{1u}(u,v)f_2(u,v)+f_1(u,v)f_{2u}(u,v), \\ h_v(u,v)&=f_{1v}(u,v)f_2(u,v)+f_1(u,v)f_{2v}(u,v) \end{aligned}\)
\(h(u,v)=f_1(u,v)/f_2(u,v)\)	\(\begin{aligned} h_u(u,v)&=(f_{1u}(u,v)f_2(u,v)-f_1(u,v)f_{2u}(u,v))/f_2(u,v)^2, \\ h_v(u,v)&=(f_{1v}(u,v)f_2(u,v)-f_1(u,v)f_{2v}(u,v))/f_2(u,v)^2 \end{aligned}\)
\(h(u,v)=g(f(u,v))\)	\(\begin{aligned} h_u(u,v)&=g'(f(u,v))f_u(u,v), \\ h_v(u,v)&=g'(f(u,v))f_v(u,v) \end{aligned}\)
\(h(u,v)=g(f_1(u,v),f_2(u,v))\)	\(\begin{aligned} h_u(u,v)&=g_x(f_1(u,v),f_2(u,v))f_{1u}(u,v)+g_y(f_1(u,v),f_2(u,v))f_{2u}(u,v), \\ h_v(u,v)&=g_x(f_1(u,v),f_2(u,v))f_{1v}(u,v)+g_y(f_1(u,v),f_2(u,v))f_{2v}(u,v) \end{aligned}\)

2変数の自動微分を表すときに必要な情報は

• 関数の値 \(f(u,v)\)
• 関数の \(u\) による偏微分係数 \(f_u(u,v)\)
• 関数の \(v\) による偏微分係数 \(f_v(u,v)\)

の3つである。そこで、\(\mathbf{R}\) の自動微分を表す型を \(\mathrm{AutoDiff2}\mathbf{R}=\mathbf{R}\times\mathbf{R}\times\mathbf{R}\) とおき、演算を定義する。上の関係式の \(f_1\),\(f_2\) とその微分係数を \((x,x_u,x_v),(y,y_u,y_v)\in\mathrm{AutoDiff2}\mathbf{R}\) で置き換えると、\(\mathrm{AutoDiff2}\mathbf{R}\) に入れるべき演算が分かる：\begin{align*}
(x,x_u,x_v)\pm(y,y_u,y_v)&=(x\pm y,x_u\pm y_u,x_v\pm y_v) \\
(x,x_u,x_v)\cdot(y,y_u,y_v)&=(x y,x_u y+x y_u,x_v y+x y_v) \\
(x,x_u,x_v)/(y,y_u,y_v)&=(x/y,(x_u y-x y_u)/y^2,(x_v y-x y_v)/y^2)
\end{align*}\(g\colon\mathbf{R}\to\mathbf{R}\) に対し，\[
g_*((x,x_u,x_v))=(g(x),g'(x)x_u,g'(x)x_v),\]
\(g\colon\mathbf{R}\times\mathbf{R}\to\mathbf{R}\) に対し，\[
g_*((x,x_u,x_v),(y,y_u,y_v))=(g(x,y),g_x(x,y)x_u+g_y(x,y)y_u,g_x(x,y)x_v+g_y(x,y)y_v)
\]

Haskellで実装すると，たとえば次のようになる:
[sourcecode lang=”plain”]
data AutoDiff2 a = AutoDiff2 a a a deriving(Eq)
instance Num a => Num (AutoDiff2 a) where
(AutoDiff2 f f_u f_v) + (AutoDiff2 g g_u g_v) = AutoDiff2 (f+g) (f_u+g_u) (f_v+g_v)
(AutoDiff2 f f_u f_v) – (AutoDiff2 g g_u g_v) = AutoDiff2 (f-g) (f_u-g_u) (f_v-g_v)
(AutoDiff2 f f_u f_v) * (AutoDiff2 g g_u g_v) = AutoDiff2 (f*g) (f_u*g+f*g_u) (f_v*g+f*g_v)
abs _ = undefined
signum _ = undefined
fromInteger n = AutoDiff2 (fromInteger n) 0 0
[/sourcecode]

さて、1変数の場合の自動微分(\(\mathbf{R}\times\mathbf{R}\) に演算を入れる)についてもう一度見てみよう。\begin{align*}
(x_0,x_1)\pm(y_0,y_1)&=(x_0\pm y_0,x_1\pm y_1) \\
(x_0,x_1)\cdot(y_0,y_1)&=(x_0 y_0,x_1 y_0+x_0 y_1) \\
(x_0,x_1)/(y_0,y_1)&=(x_0/y_0,(x_1 y_0-x_0 y_1)/y_0^2),
\end{align*}\(g\colon\mathbf{R}\to\mathbf{R}\)(微分可能) に対し、\[g_*((x_0,x_1))=(g(x_0),g'(x_0)x_1),\]\(g\colon\mathbf{R}\times\mathbf{R}\to\mathbf{R}\)(微分可能) に対し、\[g_*((x_0,x_1),(y_0,y_1))=(g(x_0,y_0),g_x(x_0,y_0)x_1+g_y(x_0,y_0)y_1),\]
であった。
よく見てみると、第2成分の \(x_1\), \(y_1\) は \(\mathbf{R}\) でなくても，\(\mathbf{R}\)-ベクトル空間なら何でもよいことが分かる。2変数の偏微分の自動微分は、\(x_1\), \(y_1\) として \(\mathbf{R}^2\) の元をとったものだと考えることができる。一般に、\(n\) 変数の自動微分は \(x_1\), \(y_1\) として \(\mathbf{R}^n\) の元をとることにより得られる。

前回の記事が適当すぎたので仕切り直しを。

「値と微分係数を同時に計算する型」AutoDiff型は，\(\mathbf{R}\times\mathbf{R}\)に適当な演算を入れたものと考えられる。微分を考える変数を \(x\)，微分を計算したい点を \(t\) とするとき，第一成分を「\(x=t\) における値」，第二成分を「\(x=t\) における微分係数」と考える。AutoDiff型の値を \(y=(y_0,y_1)\) とすると，\(y\) はある関数 \(f\) について \(y=(f(t),f'(t))\) となっている。\(h(x)=g(f(x))\) という関数があったとき，すでに得られている \(f(t),f'(t)\) の値を使って \(h(t),h'(t)\) を計算したいとする。合成関数の微分は\[h'(x)=g'(f(x))f'(x)\]で与えられるので、ここに \(x=t\) を代入すれば\[h'(t)=g'(f(t))f'(t)=g'(y_0)y_1\]となる。\(h(t)\) はもちろん \(h(t)=g(y_0)\) で与えられる。

今度は \(k(x)=g(l(x))\) という関数を考えてみよう。\(x=t\) における \(l\) の値と \(l’\) の値はすでに計算されているとして、これを \(z=(z_0,z_1)=(l(t),l'(t))\) とおく。このとき，先ほどと同じように \(k(t)=g(z_0)\)，\(k'(t)=g'(z_0)z_1\) がわかる。

このようにして，\(g\) という関数は、AutoDiff型からAutoDiff型への関数を定める。この関数を \(g_*\) とすると，AutoDiff型の値 \((y_0,y_1)\) を与えた時の \(g_*\) の値は\[g_*(y_0,y_1)=(g(y_0),g'(y_0)y_1)\]となる。AutoDiff型が「これまで計算した関数の値と微分係数」を表す型だとすれば，\(g_*\) という関数は「これまで計算した関数に \(g\) という関数を合成した時の、関数の値と微分係数」を表す型である。

上の説明では \(g\) と \(g_*\) を区別して書いたが、実際のプログラミングでは関数や演算子のオーバーロード、あるいは型クラスを使って、\(g_*\) に相当する関数も同じ演算子・関数名で書けるようにするのが普通である。

いくつかの演算について、AutoDiff型を \(\mathbf{R}\times\mathbf{R}\) とみなしたときの計算方法を与えておく。\begin{align*}
(x_0,x_1)\pm(y_0,y_1)&=(x_0\pm y_0,x_1\pm y_1) \\
(x_0,x_1)\cdot(y_0,y_1)&=(x_0 y_0,x_1 y_0+x_0 y_1) \\
(x_0,x_1)/(y_0,y_1)&=(x_0/y_0,(x_1 y_0-x_0 y_1)/y_0^2),
\end{align*}\(g\colon\mathbf{R}\to\mathbf{R}\)(微分可能) に対し、\[g_*((x_0,x_1))=(g(x_0),g'(x_0)x_1),\]\(g\colon\mathbf{R}\times\mathbf{R}\to\mathbf{R}\)(微分可能) に対し、\[g_*((x_0,x_1),(y_0,y_1))=(g(x_0,y_0),g_x(x_0,y_0)x_1+g_y(x_0,y_0)y_1),\]となる。

文章書くのだるい。気が向いたら続く。

関数の値 \(f(t)\) とその微分係数 \(f'(t)\) を同時に計算したい状況があるとする．ここでは，具体的な2次関数
\[f(x)=x^2+2x\]
で話を進める．この関数は，C++を使うと次のように実装できる．

[sourcecode lang=”cpp”]
template<typename T>
constexpr T f(T x) {
return x*x+2*x;
}
[/sourcecode]

あるいは，Haskellでは次のように書けるだろう．

[sourcecode lang=”plain”]
f :: Num a => a -> a
f x = x*x+2*x
[/sourcecode]

さて，賢明な読者は，\(f\)の導関数が\(2x+2\)であることに気づいたかもしれない．そうすると，\(f(x)\)と\(f'(x)\)を同時に計算する関数\(g\)は次のように実装できる．

[sourcecode lang=”cpp”]
// in C++
template<typename T>
std::pair<T,T> g(T x) {
return {x*x+2*x, 2*x+2};
}
[/sourcecode]

 

[sourcecode lang=”plain”]
— in Haskell
g :: Num a => a -> (a,a)
g x = (x*x+2*x, 2*x+2)
[/sourcecode]

しかし，せっかく\(f(x)\)の定義式をプログラム中に書いているのに，\(f'(x)\)の計算式を改めて書くのは無駄，と言わざるを得ない．また，プログラムを変更して\(f(x)\)の定義式を書き換えたときに\(f'(x)\)の式を直し忘れる，というミスが起こる可能性もある．できることなら，\(f(x)\)の定義式を使って，\(f'(x)\)を自動的に計算させたい．

そこで登場するのが，自動微分(automatic differentiation)と呼ばれるテクニックである。\(f(x)\)の定義の中で使っているのは和と積だけである．和と積の微分は高校で習ったと思うが，次の式で計算できる．

もとの関数	導関数
\(h(x)=f(x)\pm g(x)\)	\(h'(x)=f'(x)\pm g'(x)\)
\(h(x)=f(x)g(x)\)	\(h'(x)=f'(x)g(x)+f(x)g'(x)\)

では， f の中での和と積の挙動をいじくって，微分もついでに計算させることはできないだろうか？

幸い，C++には演算子オーバーロードという便利な機能がある．先のコードで， f はテンプレートを使って定義したので，引数の型は int や double でなくてもよい．演算子をオーバーロードしてやった自前の型を与えても良い．そこで， AutoDiff というクラスを次のように定義しよう：

[sourcecode lang=”cpp”]
template<typename T>
class AutoDiff
{
public:
T value; // 値
T diff; // 微分係数の値
constexpr AutoDiff() : value(), diff() {} // 定数0
constexpr AutoDiff(T value) : value(value), diff() {} // 定数で初期化（微分係数は0）
constexpr AutoDiff(T value, T diff) : value(value), diff(diff) {} // 値 value, 微分係数 diff で初期化

friend constexpr AutoDiff<T> operator+(AutoDiff<T> f, AutoDiff<T> g) {
return {f.value+g.value, f.diff+g.diff};
}
friend constexpr AutoDiff<T> operator*(AutoDiff<T> f, AutoDiff<T> g) {
return {f.value*g.value, f.diff*g.value+f.value*g.diff};
}
};
[/sourcecode]

このクラスを使って，先ほどの関数 g を実装しなおしてみると，次のようになる：

[sourcecode lang=”cpp”]
template<typename T>
std::pair<T,T> g(T x) {
AutoDiff<T> y = f(AutoDiff<T>(x, 1));
return {y.value, y.diff};
}
[/sourcecode]

実際にこの定義し直した g が関数の値と微分係数を正しく計算することは，各自確かめてほしい．

Haskellでは，型クラスを使うと同じようなことができる．

[sourcecode lang=”plain”]
data AutoDiff a = AutoDiff a a deriving(Eq)
instance Num a => Num (AutoDiff a) where
(AutoDiff f f’) + (AutoDiff g g’) = AutoDiff (f+g) (f’+g’)
(AutoDiff f f’) – (AutoDiff g g’) = AutoDiff (f-g) (f’-g’)
(AutoDiff f f’) * (AutoDiff g g’) = AutoDiff (f*g) (f’*g+f*g’)
abs _ = undefined
signum _ = undefined
fromInteger = flip AutoDiff 0 . fromInteger

g :: Num a => a -> (a,a)
g x = case f (AutoDiff x 1) of AutoDiff f f’ -> (f,f’)
[/sourcecode]

ここでのコード例は加法と乗法しか定義しなかったが，残りの演算は，次の表を参考にすれば容易に実装できるだろう．

もとの関数	導関数
\(h(x)=f(x)\pm g(x)\)	\(h'(x)=f'(x)\pm g'(x)\)
\(h(x)=f(x)g(x)\)	\(h'(x)=f'(x)g(x)+f(x)g'(x)\)
\(h(x)=f(x)/g(x)\)	\(h'(x)=(f'(x)g(x)-f(x)g'(x))/g(x)^2\)
\(h(x)=g(f(x))\)	\(h'(x)=g'(f(x))f'(x)\)
\(h(x)=g(f_1(x),f_2(x))\)	\(h'(x)=g_u(f_1(x),f_2(x))f_1′(x)+g_v(f_1(x),f_2(x))f_2′(x)\) ただし\(g_u,g_v\)は\(g\)の偏導関数

文章を書くのに疲れたので，今日はここまで．（続く）