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

「値と微分係数を同時に計算する型」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\)の偏導関数

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