(cache)双対数を利用した自動微分のしくみ

最近、SLAMの勉強をしてみたいなぁと思いceres-solverのコードを読んでいたのですが、ヤコビアンの導出のところで見慣れない操作をしていることに気付きました。調べてみると双対数というものを利用して自動微分なる計算をしているらしいです。

ということで、この記事では双対数（二重数, dual number）とそれを利用した自動微分についてまとめてみます。平たく言えば、プログラムで微分を計算する方法を紹介します。

ざっくり理解で書いている部分もあるので、間違いを見つけたらマサカリ投げてください。

1. 前置き

今回の最終目標は、与えられた関数のある点における1回微分を計算することです。

物理シミュレーションや最適化のアルゴリズムでは頻繁に微分値が要求されるのですが、偏導関数は変数の次元に応じて数が増えるので、全てを手計算してコーディングするには限界があります。

そこで「元の関数だけを定義しておいて、導関数はプログラムに自動的に計算させよう」という考え方をするのが自動微分です。

高校数学で習うように、基礎的な数学関数から構成される関数の微分は、ある程度のルールさえ覚えれば手計算できてしまいます。

たとえば、元の関数が

f (x, y) = \frac{\sin x y}{x}

なら、 $x$ についての偏導関数は

\begin{array}{rcl} \frac{\partial f}{\partial x} & = & {(\frac{1}{x} \cdot \sin x y)}^{'} \\ = & \frac{1}{x} \cdot y \cdot \cos x y - \frac{\sin x y}{x^{2}} \end{array}

ですね。計算の過程では、 $\sin x$ の導関数が $\cos x$ であるという知識と、関数の加減乗除の微分・合成関数の微分についてのルールを利用しています。

こうした基礎的な数学関数の導関数と四則演算のルールに基づいて機械的に計算するだけならプログラムでも可能なはずです。しかし、なんの指針もなく導関数の計算ルールをアルゴリズム化するのはなかなかツラいものがあります。

そこで利用するのが、今回のもうひとつのテーマである双対数（二重数・dual number）という概念です。

2. 導関数の演算則

まずは関数の加減乗除の微分についてルールを確認しておきます。

$f, g$ を微分可能な関数とします。

線形性

(α f + β g)^{'} = α f^{'} + β g^{'}

積の微分

(f g)^{'} = f g^{'} + g f^{'}

商の微分

{(\frac{f}{g})}^{'} = \frac{f^{'} g - f g^{'}}{g^{2}}

このあたりは高校で習った通りだと思います。

3. 双対数

次に双対数について見ていきます。

3.1. 双対数の定義

双対数は複素数と同じように二次元の値を持つ数です。

d = a + b ϵ

複素数では虚部を表すために $i$ を使いますが、双対数では非実部を表すために $ϵ$ を使います。

双対数は複素数と同じように交換法則や結合法則を満たします。違うのは複素数が $i^{2} = - 1$ としていたのに対し、双対数では $ϵ^{2} = 0$ と定義するという点です。つまり「2回掛けると0になる特殊な元」を導入して値を記述するのが双対数です。

3.2. 双対数の演算則

それでは、双対数の四則演算のルールを確認していきます。

線形和

α (a_{1} + b_{1} ϵ) + β (a_{2} + b_{2} ϵ) = (α a_{1} + β a_{2}) + (α b_{1} + β b_{2}) ϵ

積算

(a_{1} + b_{1} ϵ) \cdot (a_{2} + b_{2} ϵ) = a_{1} a_{2} + (a_{1} b_{2} + a_{2} b_{1}) ϵ

除算
少し迷うところですが、分母・分子に同じ数を掛け合わせることで計算できます。

\begin{array}{rcl} \frac{(a_{2} + b_{2} ϵ)}{(a_{1} + b_{1} ϵ)} & = & \frac{(a_{2} + b_{2} ϵ) (a_{1} - b_{1} ϵ)}{(a_{1} + b_{1} ϵ) (a_{1} - b_{1} ϵ)} \\ = & \frac{a_{1} a_{2} + (a_{1} b_{2} - a_{2} b_{1}) ϵ}{a_{1}^{2}} \\ = & \frac{a_{2}}{a_{1}} + \frac{a_{1} b_{2} - a_{2} b_{1}}{a_{1}^{2}} ϵ \end{array}

$ϵ^{2} = 0$ になるというルール以外は普通の多項式計算と同じですね。

3.3. 双対数と導関数の演算則の共通性

さて、ここで双対数と導関数の計算を比較すると性質が共通している部分があることに気付きます。

実部を関数 $f, g$ , 非実部を導関数 $f^{'}, g^{'}$ とする双対数同士の演算を考えてみると、

線形和

α (f + f^{'} ϵ) + β (g + g^{'} ϵ) = (α f + β g) + (α f^{'} + β g^{'}) ϵ

積算

(f + f^{'} ϵ) \cdot (g + g^{'} ϵ) = f g + (f g^{'} + g f^{'}) ϵ

除算

\frac{(f + f^{'} ϵ)}{(g + g^{'} ϵ)} = \frac{f}{g} + \frac{g f^{'} - f g^{'}}{g^{2}} ϵ

というように、四則演算を行う前後で「実部は関数、非実部はその導関数」という関係を保っていることがわかります。この事実から、任意の実係数多項式で与えられる関数とその合成関数についてはこの関係が保たれることが保証されます。

指数関数や三角関数もテイラー展開すると実係数多項式になるため、双対数に関する計算を矛盾なく定義できます。

例 : 指数関数・三角関数

\begin{array}{rcl} e^{a + b ϵ} & = & e^{a} + e^{a} b ϵ = e^{a} (1 + b ϵ) \\ \sin (a + b ϵ) & = & \sin a + (\cos a) b ϵ \\ \cos (a + b ϵ) & = & \cos a - (\sin a) b ϵ \end{array}

結果として、基礎的な数学関数のみで記述される関数は大抵この性質を保存することになります。

4. 双対数を利用した自動微分

4.1. 自動微分の原理

次に、この双対数の性質をどう利用すれば微分値が計算できるのか、という部分を見ていきます。

ここまでの話に従えば、実係数多項式の合成関数として記述される関数 $f (d)$ に双対数

\begin{array}{rcl} d & = & y + \dot{y} ϵ \\ = & y + \frac{\partial y}{\partial x} ϵ \end{array}

を代入すると、

\begin{array}{rcl} f (d) & = & f (y) + f^{'} (y) ϵ \\ = & f (y (x)) + \frac{\partial f}{\partial x} ϵ \end{array}

が出てきますね。

つまり $x_{0} + ϵ$ という双対数を目的の関数に代入して計算すると $f (x_{0} + ϵ) = f (x_{0}) + f^{'} (x_{0}) ϵ$ が出てくると言えます。

なんとも不思議な話ですが、 $x_{0} + ϵ$ という特定の双対数を代入して関数を計算するだけで関数値と微分値が同時に求められるわけです。同じく変数 $x$ が多次元の場合も、長さ1の方向ベクトル $p$ を用いて $f (x_{0} + p ϵ)$ を計算すれば、そのベクトル方向の微分値が計算できます。

4.2. 計算例

いまいちピンとこないと思いますので、適当な関数の計算例を見てみましょう。

目的の関数を $f (x, y) = \sin x + x y$ とします。

手計算で勾配ベクトルを求めると $[\cos x + y, x]^{T}$ ですね。

$x$ についての偏微分とは $[1, 0]^{T}$ 方向の勾配なので、 $[x_{0} + ϵ, y_{0}]^{T}$ を代入してみると、

\begin{array}{rcl} f (x_{0} + ϵ, y_{0}) & = & \sin (x_{0} + ϵ) + (x_{0} + ϵ) y_{0} \\ = & \sin x_{0} + x_{0} y_{0} + (\cos x_{0} + y_{0}) ϵ \end{array}

ちゃんと $x$ 方向の勾配が計算できています。

同じく $[x_{0}, y_{0} + ϵ]^{T}$ を代入してみると、

\begin{array}{rcl} f (x_{0}, y_{0} + ϵ) & = & \sin x_{0} + x_{0} (y_{0} + ϵ) \\ = & \sin x_{0} + x_{0} y_{0} + x_{0} ϵ \end{array}

ちゃんと $y$ 方向の勾配が計算できてますね。

4.3. 実装

以上の自動微分をプログラムとして実装する方法について考えてみます。

最も素直な方法は、双対数を表す型を定義しそれに対して四則演算子とcmath相当の数学関数を定義するというものです。この型を数値として扱いcos(x) + x * yのように元となる関数を記述するだけで自動微分のためのコードを実装したのと等価になります。

C++の線形代数ライブラリEigenやGoogleが開発している最適化ライブラリceres-solver内の自動微分はこの方式を採用しています。実装はそれぞれAutoDiffScalar.hとjet.hにまとまっているので読めば雰囲気をつかめると思います。

ただし自動微分を使ってヤコビアンを計算する場合、何も考えずにこの方法を使うと枝刈りが行われず同じ計算を何度も実行することになり効率が悪くなります。

たとえば $f (x, y, z) = g (x, y) h (z)$ という関数のヤコビアンを計算するとき、 $x, y$ の偏微分を計算する過程で $h (z)$ の値は変化しないので使い回すのが好ましいです。この例では変数が3次元なので計算コストの差はたかが知れていますが、変数の次元が増えると無視できないものになってきます。

また、そもそも双対数は実部で関数値を計算するのでヤコビアンの計算過程で同じ $f (x_{0})$ が何度も計算されることになります。同じ $x_{0}$ に対して計算する際は、途中で計算される実部をすべてキャッシュするのが理想です。

自動微分を実装するもうひとつの方法は、関数定義を木構造に分解し最適化をかけた上でコードを生成するようなツールを利用するというものです。関数の定義は、基本となる数学関数の四則演算と合成で記述されているので、木構造に分解して整理をかけることで効率的なコードを生成できます。このあたり筆者はあまり詳しくないのですが、そういう関数型言語の機能やDSLみたいなものってあったりするんでしょうかね…？

5. 双対数による自動微分と数値微分の比較

最後に、双対数による自動微分と数値微分を比較してみます。

5.1. 数値微分とは

数値微分というのは微分の定義に基づいて、

f^{'} (x_{0}, p) = lim_{h \to 0} \frac{f (x_{0} + h p) - f (x_{0})}{h}

f^{'} (x_{0}, p) = lim_{h \to 0} \frac{f (x_{0}) - f (x_{0} - h p)}{h}

f^{'} (x_{0}, p) = lim_{h \to 0} \frac{f (x_{0} + h p) - f (x_{0} - h p)}{2 h}

などのように微分値を計算する方法です。

上の式はそれぞれ前進差分/後退差分/中央差分と呼ばれており、前進差分と後退差分は計算コストが小さい一方、中央差分は計算精度が良いことが知られています。

5.2. 数値微分の利点

数値微分の最大の利点は、任意の関数の微分が計算できるという点です。先ほど説明したように、双対数を用いた自動微分ができるのは実係数多項式の合成関数だけであり、それ以外の関数には適用できません。その点、数値微分は $f (x)$ の定義さえあれば任意の関数の微分が可能です。

また実装しやすいという利点もあります。先ほど述べたように双対数の計算は比較的重いので、自動微分の計算コストは相応に高くなります。一方、数値微分では単に $f (x)$ を1 or 2回計算するだけなので単純ですし、計算コスト削減のための工夫もしやすいです。

5.3. 自動微分の利点

自動微分の最大の利点は数値計算上の誤差（丸め誤差や桁落ち）を抑えやすいという点です。数値微分の計算では、ほぼ同じ値を取る数同士の引き算 $f (x + h) - f (x - h)$ の部分で桁落ち誤差が出やすいです。また誤差をなるべく小さくするような $h$ の選定にも気を配らなければなりません。一方、自動微分では解析的に導出される手順に基づいて微分値を計算するので、純粋に導関数の計算上で生じる誤差のみが見込まれます。

以上、双対数を利用した自動微分についてでした。

こういう代数学の応用を考える人って本当にすごいと思います。

この記事を書くのに当たって以下の文献にお世話になりました。ありがとうございました。

かみのメモ

コンピュータビジョン・プログラムな話題中心の勉強メモ

双対数を利用した自動微分のしくみ

もくじ