関数の値 \(f(t)\) とその微分係数 \(f'(t)\) を同時に計算したい状況があるとする.ここでは,具体的な2次関数
\[f(x)=x^2+2x\]
で話を進める.この関数は,C++を使うと次のように実装できる.
template<typename T> constexpr T f(T x) { return x*x+2*x; }
あるいは,Haskellでは次のように書けるだろう.
f :: Num a => a -> a f x = x*x+2*x
さて,賢明な読者は,\(f\)の導関数が\(2x+2\)であることに気づいたかもしれない.そうすると,\(f(x)\)と\(f'(x)\)を同時に計算する関数\(g\)は次のように実装できる.
// in C++ template<typename T> std::pair<T,T> g(T x) { return {x*x+2*x, 2*x+2}; }
-- in Haskell g :: Num a => a -> (a,a) g x = (x*x+2*x, 2*x+2)
しかし,せっかく\(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
というクラスを次のように定義しよう:
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}; } };
このクラスを使って,先ほどの関数 g
を実装しなおしてみると,次のようになる:
template<typename T> std::pair<T,T> g(T x) { AutoDiff<T> y = f(AutoDiff<T>(x, 1)); return {y.value, y.diff}; }
実際にこの定義し直した g
が関数の値と微分係数を正しく計算することは,各自確かめてほしい.
Haskellでは,型クラスを使うと同じようなことができる.
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')
ここでのコード例は加法と乗法しか定義しなかったが,残りの演算は,次の表を参考にすれば容易に実装できるだろう.
もとの関数 | 導関数 |
---|---|
\(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\)の偏導関数 |
文章を書くのに疲れたので,今日はここまで.(続く)