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

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

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$ とおく．偏導関数は，変数を添字にして表すことにする．

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)$の定義の中で使っているのは和と積だけである．和と積の微分は高校で習ったと思うが，次の式で計算できる．

では， 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]

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

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