はじめに

自分で微分ってしたいですか？
値さえわかれば別に関数形いらなくないですか？

そんなあなたに最適、そうそれが我々が微分するのでは無く機会に微分させる自動微分です！
まあ正直自動微分は最近機械学習方面で流行ってるから探せばいくらでも解説出てくるし、導入書くの面倒臭くなってきたのでこんなところでいいですよね。
この自動微分をdual number (二重数または双対数)を使って実装したいわけですが、場合によっては多変数の高階微分まで必要だという事もあるのではないでしょうか。
今回はその多変数高階の自動微分をhyper-dual numberを使って実装できないかなあと思うわけです。

dual number

まず、dual numberとはなにかという話ですが、これは

d = a + b \varepsilon

という形を取る数の事です。
ただし、ここで$a$と$b$は実数ですが、$\varepsilon$は二乗を取ると$0$になる、すなわち$\varepsilon^2 = 0$となる特殊な何かです。
ここで、この$\varepsilon$をこの文章では勝手に「二重単位」と呼び、$a$を実部、$b$を「二重部」と呼んでおきます。

dual numberは複素数と同じように、結合法則や交換法則を満たします。
四則演算のルールですが、これも複素数と同じようにかけます。
例えば、2つのdual number

\begin{align}
d_1 & = a_1 + b_1 \varepsilon\\
d_2 & = a_2 + b_2 \varepsilon
\end{align}

に対して、

\begin{align}
d_1 + d_2 & = a_1 + b_1 \varepsilon + a_2 + b_2 \varepsilon\\
          & = \left(a_1 + a_2\right) + \left(b_1 + b_2\right) \varepsilon\\
d_1 - d_2 & = a_1 + b_1 \varepsilon - (a_2 + b_2 \varepsilon)\\
          & = \left(a_1 - a_2\right) + \left(b_1 - b_2\right) \varepsilon\\
d_1 d_2 & = \left(a_1 + b_1 \varepsilon \right) \left(a_2 + b_2 \varepsilon\right)\\
        & = a_1 a_2 + \left(a_1 b_2 + a_2 b_1\right) \varepsilon\\
\frac{d_2}{d_1} & = \frac{a_2 + b_2 \varepsilon}{a_1 + b_1 \varepsilon}\\
                & = \frac{\left( a_2 + b_2 \varepsilon \right) \left( a_1 - b_1 \varepsilon \right)}{\left( a_1 + b_1 \varepsilon \right) \left( a_1 - b_1 \varepsilon \right)}\\
                & = \frac{a_2}{a_1} + \frac{a_1 b_2 - a_2 b_1}{a_1^2} \varepsilon
\end{align}

となります。

さて、ではなぜこのdual numberが自動微分に関わるのでしょう。
話を進める前に、dual numberでは無く普通の実数のテイラー級数展開をまずは思い出してみましょう。
ある関数$f(x)$があったとして、その$x = x_0 + h$での級数展開を考えます。

\begin{align}
f(x_0 + h) & = f(x_0) + f'(x_0) h + \frac{1}{2!} f''(x_0) h^2 + \cdots \\
\end{align}

となります。
ここで、$h$に$b \varepsilon$を、$x_0$に$a$を代入することで、dual number $d$に対する$f(d)$の値を求めてみます。

\begin{align}
f(d) & = f(a + b\varepsilon) \\
     & = f(a) + f'(a) b \varepsilon + \frac{1}{2!} f''(a) (b \varepsilon)^2 + \cdots \\
     & = f(a) + f'(a) b \varepsilon
\end{align}

さらに、ここで$b = 1$を代入してしまいましょう。
すると、結果は

\begin{align}
f(a + \varepsilon) = f(a) + f'(a) \varepsilon
\end{align}

となります。
これはつまり、(当たり前ですが)「任意の関数$f(x)$に対し、二重部が1なdual number $a + \varepsilon$を入力すると、その結果の実部はdual numberの実部をその関数に入力した物に、二重部はその導関数になる。」という事です。

なんだかちょっと狐につつまれたような気分にならなくもないので、少し念の為確認しておきましょう。
関数

\begin{align}
f(x) = x^2 + 2x + 3
\end{align}

を考えます。
これの導関数は、

\begin{align}
f'(x) = 2x + 2
\end{align}

となります。
さて、一方、

\begin{align}
f(a + \varepsilon) & = (a + \varepsilon)^2 + 2(a + \varepsilon) + 3 \\
                   & = a^2 + 2 a \varepsilon + \varepsilon^2 + 2 a + 2 \varepsilon + 3\\
                   & = a^2 + 2 a + 3 + (2 a + 2) \varepsilon
\end{align}

となり、確かに$f(a) + f'(a)\varepsilon$になっています。

実装C++

大雑把にこんな感じ。
これをdualnumber.hppとかにしておいて、

#include <cmath>

class dualnumber{
	//real part and dual part
	double real, dual;
	public:
	dualnumber(){
	};
	dualnumber(const double real_, const double dual_): real(real_), dual(dual_){
	}
	dualnumber operator+ () const{
		return{
			real,
			dual,
		};
	}
	dualnumber operator- () const{
		return{
			- real,
			- dual,
		};
	}
	friend dualnumber operator+ (const dualnumber lhs, const dualnumber rhs){
		return{
			lhs.real + rhs.real,
			lhs.dual + rhs.dual,
		};
	}
	friend dualnumber operator+ (const dualnumber lhs, const double rhs){
		return{
			lhs.real + rhs,
			lhs.dual,
		};
	}
	friend dualnumber operator+ (const double lhs, const dualnumber rhs){
		return{
			lhs + rhs.real,
			rhs.dual,
		};
	}
	friend dualnumber operator- (const dualnumber lhs, const dualnumber rhs){
		return{
			lhs.real + rhs.real,
			lhs.dual + rhs.dual,
		};
	}
	friend dualnumber operator- (const dualnumber lhs, const double rhs){
		return{
			lhs.real - rhs,
			lhs.dual,
		};
	}
	friend dualnumber operator- (const double lhs, const dualnumber rhs){
		return{
			lhs - rhs.real,
			rhs.dual,
		};
	}
	friend dualnumber operator* (const dualnumber lhs, const dualnumber rhs){
		return{
			lhs.real * rhs.real,
			lhs.dual * rhs.real + lhs.real * rhs.dual,
		};
	}
	friend dualnumber operator* (const dualnumber lhs, const double rhs){
		return{
			lhs.real * rhs,
			lhs.dual * rhs,
		};
	}
	friend dualnumber operator* (const double lhs, const dualnumber rhs){
		return{
			lhs * rhs.real,
			lhs * rhs.dual,
		};
	}
	friend dualnumber operator/ (const dualnumber lhs, const dualnumber rhs){
		return{
			lhs.real / rhs.real,
			(lhs.dual * rhs.real - lhs.real * rhs.dual) / (rhs.real * rhs.real),
		};
	}
	friend dualnumber operator/ (const dualnumber lhs, const double rhs){
		return{
			lhs.real / rhs,
			lhs.dual / rhs,
		};
	}
	friend dualnumber operator/ (const double lhs, const dualnumber rhs){
		return{
			lhs / rhs.real,
			lhs * rhs.dual / (rhs.real * rhs.real),
		};
	}
	double getRealPart() const{
		return real;
	}
	double getDualPart() const{
		return dual;
	}
};

あとはmain.cppなどで

#include "dualnumber.hpp"
#include <iostream>

int main(){
	dualnumber x(2.0, 1.0);
	dualnumber f = x * x * x + 2.0 * x + 1.0;
	std::cout << f.getRealPart() << ", " << f.getDualPart() << std::endl;;
	return 0;
}

2変数関数

では、関数$f$が$f(x, y)$のような2変数関数だったらどうなるでしょうか。
まず、(書くのが面倒臭いので)$\frac{\partial f}{\partial x}$を$f_x$のように書く事にします。
また、dual numberとして、$d_1 = a_1 + \varepsilon_1$、$d_2 = a_2 + \varepsilon_2$を用意します。
ここで、$\varepsilon_1^2 = 0$、$\varepsilon_2^2 = 0$ですが、$\varepsilon_1 \varepsilon_2 \neq 0 $に注意してください。
また、$f_{xy} = f_{yx}$を一旦仮定しておきます。
2変数関数のテイラー級数展開は、

\begin{align}
f(d_1, d_2) & = f(a_1 + \varepsilon_1, a_2 + \varepsilon_2)\\
            & = f(a_1, a_2) + f_x(a_1, a_2) \varepsilon_1 + f_y(a_1, a_2) \varepsilon_2 + \frac{1}{2} \left\{f_{xx}(a_1, a_2) \varepsilon_1^2 + f_{xy}(a_1, a_2) \varepsilon_1 \varepsilon_2 + f_{yx}(a_1, a_2) \varepsilon_2 \varepsilon_1 + f_{yy}(a_1, a_2) \varepsilon_2^2 \right\}\\
            & = f(a_1, a_2) + f_x(a_1, a_2) \varepsilon_1 + f_y(a_1, a_2) \varepsilon_2 + f_{xy}(a_1, a_2) \varepsilon_1 \varepsilon_2
\end{align}

となります。
2変数関数の場合、個々の変数での微分項$f_x, f_y$に加え、交差微分$f_{xy}$もどうやらこれで求める事ができそうです。

二階導関数への拡張

では、二階の微分がほしい場合はどうでしょうか。
この場合、hyper dual numberと呼ばれる物を定義します。

d = a_0 + a_1 \varepsilon_1 + a_2 \varepsilon_2 + a_{12} \varepsilon_1 \varepsilon_2

ここで、前述の通り$\varepsilon_1^2 = 0$、$\varepsilon_2^2 = 0$ですが、$\varepsilon_1 \varepsilon_2 \neq 0 $に注意しましょう。

\begin{align}
f(d) & = f(a_0 + a_1 \varepsilon_1 + a_2 \varepsilon_2 + a_{12}\varepsilon_1 \varepsilon_2)\\
     & = f(a_0) + f'(a_0) (a_1 \varepsilon_1 + a_2 \varepsilon_2 + a_{12}\varepsilon_1 \varepsilon_2) + \frac{1}{2} f''(a_0) (a_1 \varepsilon_1 + a_2 \varepsilon_2 + a_{12}\varepsilon_1 \varepsilon_2)^2\\
     & = f(a_0) + f'(a_0) (a_1 \varepsilon_1 + a_2 \varepsilon_2 + a_{12}\varepsilon_1 \varepsilon_2) + f''(a_0) a_1 a_2 \varepsilon_1 \varepsilon_2
\end{align}

となります。
さて、二階微分がほしいので、$a_1 = a_2 = 1$として選ぶと$f''(a_0)$を引っ張り出すのに都合が良さそうです。
$a_{12}$は一階の導関数部分に出てきますが、これを$0$にしておくとこれまた綺麗そうです。

したがって、まとめますと二階の導関数が欲しい場合、

\begin{align}
d & = a + \varepsilon_1 + \varepsilon_2\\
f(d) & = f(a) + f'(a) \varepsilon_1 + f'(a)\varepsilon_2 + f''(a) \varepsilon_1 \varepsilon_2
\end{align}

とまあ、それなりに綺麗に対応が取れました。
二階の微分の場合、$f'(a)$が二回求まる(二階だけにね！(激ウマ爆笑ギャグ))という無駄が発生しますが、これで無事$f''(a)$も求めることができそうです。

なお、このhyper dual numberの四則演算に関しては、以下のようになります。

\begin{align}
d_1 & = a_0 + a_1 \varepsilon_1 + a_2 \varepsilon_2 + a_{12} \varepsilon_1 \varepsilon_2\\
d_2 & = b_0 + b_1 \varepsilon_1 + b_2 \varepsilon_2 + b_{12} \varepsilon_1 \varepsilon_2
\end{align}

として、

\begin{align}
d_1 + d_2 & = a_0 + b_0 + (a_1 + b_1) \varepsilon_1 + (a_2 + b_2) \varepsilon_2 + (a_{12} + b_{12}) \varepsilon_1 \varepsilon_2\\
d_1 - d_2 & = a_0 - b_0 + (a_1 - b_1) \varepsilon_1 + (a_2 - b_2) \varepsilon_2 + (a_{12} - b_{12}) \varepsilon_1 \varepsilon_2\\
d_1 d_2 & = a_0 b_0 + (a_0 b_1 + a_1 b_0) \varepsilon_1 + (a_0 b_2 + a_2 b_0) \varepsilon_2 + (a_0 b_{12} + a_1 b_2 + a_2 b_1 + a_{12} b_0) \varepsilon_1 \varepsilon_2\\
\frac{d_2}{d_1} & = \left( \frac{1}{a_0} - \frac{a_1}{a_0^2} \varepsilon_1 - \frac{a_2}{a_0^2} \varepsilon_2 + \left( - \frac{a_{12}}{a_0^2} + \frac{2 a_1 a_2}{a_0^3} \right) \varepsilon_1 \varepsilon_2 \right) (b_0 + b_1 \varepsilon_1 + b_2 \varepsilon_2 + b_{12} \varepsilon_1 \varepsilon_2)\\
                & = \frac{b_0}{a_0} + \frac{a_0 b_1 - a_1 b_0}{a_0^2} \varepsilon_1 + \frac{a_0 b_2 - a_2 b_0}{a_0^2} \varepsilon_2 + \left( \frac{b_{12}}{a_0} - \frac{a_1 b_2 + a_2 b_1}{a_0^2} - \frac{a_{12} b_0}{a_0^2} + \frac{2 a_1 a_2 b_0}{a_0^3} \right) \varepsilon_1 \varepsilon_2
\end{align}

とまあこんな感じです多分(最後の除算はなんか微妙に自信がない…)。

2変数関数の二階導関数

もうだんだんわけがわからなくなってきましたがやってみましょう。
文字がややこしくなってきたので$\varepsilon$に加えて$\epsilon$も導入します。
$d_1 = a_1 + \varepsilon_1 + \epsilon_1$、$d_2 = a_2 + \varepsilon_2 + \epsilon_2$
です。
くどいようですが$\varepsilon$などなどは自分との2乗の時のみ$0$になり、$\varepsilon$同士の掛け算では0になりません。
よし気合入れていくぞ！

\begin{align}
f(d_1, d_2) & = f(a_1 + \varepsilon_1 + \epsilon_1, a_2 + \varepsilon_2 + \epsilon_2)\\
            & = f(a_1, a_2) + f_x(a_1, a_2) (\varepsilon_1 + \epsilon_1) + f_y(a_1, a_2) (\varepsilon_2 + \epsilon_2) + \frac{1}{2} \left\{f_{xx}(a_1, a_2) (\varepsilon_1 + \epsilon_1)^2 + f_{xy}(a_1, a_2) (\varepsilon_1 + \epsilon_1) (\varepsilon_2 + \epsilon_2) + f_{yx}(a_1, a_2) (\varepsilon_2 + \epsilon_2) (\varepsilon_1 + \epsilon_2) + f_{yy}(a_1, a_2) (\varepsilon_2 + \epsilon_2)^2 \right\}\\
            & = f(a_1, a_2) + f_x(a_1, a_2) (\varepsilon_1 + \epsilon_1) + f_y(a_1, a_2) (\varepsilon_2 + \epsilon_2) + f_{xx}(a_1, a_2) \varepsilon_1 \epsilon_1 + f_{xy} (\varepsilon_1 + \epsilon_1) (\varepsilon_2 + \epsilon_2) + f_{yy}(a_1, a_2) \varepsilon_2 \epsilon_2
\end{align}

とまあひどく面倒くさいのですがなんとか求まりそうな感じです。

結論

hyper dual numberを使って、複数変数の高階導関数を自動微分で求められるかどうか、試してみました。
まずは2変数二階から始めましたがまあどうやらできそうです。
流石にそこまでは面倒くさいので実装はしませんが(というのは実は嘘で、したんですがコードが汚くてとても公にできるものではないというのが本音です)。