二重数を用いた自動微分

はじめに

先日、面白い記事を見つけた。二重数(下記ページでは「双対数」と呼ばれているが、Wikipediaでは「二重数」)を用いると、初等的な関数($e^x,\sin x,\cos x,\log x,\dots$)の加減乗除と合成関数で作られた関数の微分係数を機械的に計算するプログラムを比較的簡単に実装できるという話。

双対数を利用した自動微分についてのメモ

数学的な妥当性の説明が簡単に述べられているが、私の知識不足故か、素直に納得できないところがあったので、自分が気になった箇所について考えてみた。以下それについて述べる。

紹介されている数式の導出

$f(a+b\epsilon )=f(a)+b\epsilon f^{\prime }(a)$の導出

関数$f:x\in G\subseteq \mathbb{C}\mapsto f(x)\in \mathbb{C}$は$a_0\in G$を中心に冪級数展開可能であるとする。$a,b\in \mathbb{C}$とし、$a$は収束領域に入っているものとする。

$$\begin{array}{rl}f(a+b\epsilon )& =\sum _{k=0}^{\mathrm{\infty }}\frac{1}{k!}{f}^{(k)}(a_0)(a+b\epsilon –a_0{)}^k\\ & =f(a_0)+\sum _{k=1}^{\mathrm{\infty }}\frac{1}{k!}{f}^{(k)}(a_0)(a+b\epsilon –a_0{)}^k\\ & =f(a_0)+\sum _{k=1}^{\mathrm{\infty }}\frac{1}{k!}{f}^{(k)}(a_0)\sum _{l=0}^k(\genfrac{}{}{0ex}{}{k}{l})(a-a_0{)}^{k-l}(b\epsilon {)}^l\\ & =f(a_0)+\sum _{k=1}^{\mathrm{\infty }}\frac{1}{k!}{f}^{(k)}(a_0)[(a-a_0{)}^k+b\epsilon k(a-a_0{)}^{k-1}]\\ & =\sum _{k=0}^{\mathrm{\infty }}\frac{1}{k!}{f}^{(k)}(a_0)(a-a_0{)}^k+b\epsilon \sum _{k=1}^{\mathrm{\infty }}\frac{1}{k!}{f}^{(k)}(a_0)k(a-a_0{)}^{k-1}\\ & =f(a)+b\epsilon {(\frac{\mathrm{d}}{\mathrm{d}x}\sum _{k=0}^{\mathrm{\infty }}\frac{1}{k!}{f}^{(k)}(a_0)(x-a_0{)}^k)|}_{x=a}\\ & =f(a)+b\epsilon f^{\prime }(a)\end{array}$$

$f(\mathit{a}+\epsilon \mathit{p})=f(\mathit{a})+{\mathit{p}}^{\mathrm{\top }}\mathrm{\nabla }f(a)$ の導出

$f:\mathit{x}\in G\subseteq {\mathbb{R}}^n\mapsto f(\mathit{x})\in \mathbb{C}$は${\mathit{a}}_0$を中心として冪級数展開可能であるとする。このとき$f$は無限回微分可能であり、Youngの定理より$\partial /\partial x_i$と$\partial /\partial x_j$が可換であることに注意する。$\mathit{a},\mathit{p}\in \mathbb{R}$とし、$\mathit{a}$は収束領域に入っているものとする。

$$\begin{array}{rl}f(\mathit{a}+\epsilon \mathit{p})& =f({\mathit{a}}_0+(\mathit{a}–{\mathit{a}}_0+\epsilon \mathit{p}))\\ & =f({\mathit{a}}_0)+\sum _{k=1}^{\mathrm{\infty }}\frac{1}{k!}{[(\mathit{a}–{\mathit{a}}_0+\epsilon \mathit{p}{)}^{\mathrm{\top }}\mathrm{\nabla }]}^k(f,{\mathit{a}}_0)\\ & =f({\mathit{a}}_0)+\sum _{k=1}^{\mathrm{\infty }}\frac{1}{k!}{[(\mathit{a}–{\mathit{a}}_0{)}^{\mathrm{\top }}\mathrm{\nabla }+\epsilon {\mathit{p}}^{\mathrm{\top }}\mathrm{\nabla }]}^k(f,{\mathit{a}}_0)\\ & =f({\mathit{a}}_0)+\sum _{k=1}^{\mathrm{\infty }}\frac{1}{k!}\{{[(\mathit{a}–{\mathit{a}}_0{)}^{\mathrm{\top }}\mathrm{\nabla }]}^k+k{[(\mathit{a}–{\mathit{a}}_0{)}^{\mathrm{\top }}\mathrm{\nabla }]}^{k-1}\epsilon {\mathit{p}}^{\mathrm{\top }}\mathrm{\nabla }\}(f,{\mathit{a}}_0)\\ & =f({\mathit{a}}_0)+\sum _{k=1}^{\mathrm{\infty }}\frac{1}{k!}{[(\mathit{a}–{\mathit{a}}_0{)}^{\mathrm{\top }}\mathrm{\nabla }]}^k(f,{\mathit{a}}_0)\\ & \phantom{=}+\{(\epsilon {\mathit{p}}^{\mathrm{\top }}\mathrm{\nabla })\sum _{k=1}^{\mathrm{\infty }}\frac{1}{(k-1)!}{[(\mathit{a}–{\mathit{a}}_0{)}^{\mathrm{\top }}\mathrm{\nabla }]}^{k-1}\}(f,{\mathit{a}}_0)\\ & =f(\mathit{a})+(\epsilon {\mathit{p}}^{\mathrm{\top }}\mathrm{\nabla })(\sum _{k=0}^{\mathrm{\infty }}\frac{1}{k!}{[(\mathit{a}–{\mathit{a}}_0{)}^{\mathrm{\top }}\mathrm{\nabla }]}^{k}f,{\mathit{a}}_0)\\ & =f(\mathit{a})+(\epsilon {\mathit{p}}^{\mathrm{\top }}\mathrm{\nabla })(g,{\mathit{a}}_0)\text{where}g:\mathit{x}\mapsto f(\mathit{x}+(\mathit{a}-{\mathit{a}}_0))\\ & =f(\mathit{a})+\epsilon {\mathit{p}}^{\mathrm{\top }}\mathrm{\nabla }(f,\mathit{a})\end{array}$$

計算機での実装を意識

計算機で先述の冪級数展開を計算するわけにはいかない。実装上は、まず二重数クラスを定義し、$f(a+b\epsilon )=f(a)+b\epsilon f^{\prime }(a)$を初等的な関数($\sqrt{x},e^x,\sin x,\cos x,\tan x,\log x,\dots$)について実装する。

双対数を利用した自動微分についてのメモの解説にあるように、$f+f^{\prime }\epsilon ,g+g^{\prime }\epsilon$型の二重数同士の四則演算が好都合な結果をもたらすので、$f(a+b\epsilon )=f(a)+b\epsilon f^{\prime }(a)$を実装した関数同士の四則演算で構成されたユーザー定義関数$h$については、引数に$a+\epsilon$を渡すだけで自動的に$h(a)+h^{\prime }(a)$が計算される。

Juliaでの実装

Juliaで実装したコードをGistで公開している: Automatic differentiation using Dual Number