dual の設計

このページでは、Dual[T] の数学と、この型が dual API に記載された演算とインスタンスをちょうど持つ理由を説明します。dual チュートリアル は、この背景知識なしで型を使います。

設計目標

Luna Flow のトレイトに対して書かれた通常のプログラムを、別の数値型の上で実行することで、その正確な 1 階導関数を計算します。この型は、インスタンスが表明する代数法則を満たし、定義域の失敗を共有の arithmetic エラー値で報告し、コンテナや多項式のライブラリから独立していなければなりません。

数学的背景

双対数の代数

$R$ を可換環とします（浮動小数点型の場合は、それが近似する実数）。$R$ 上の双対数とは、多項式環 $R[x]$ を $x^2$ で生成されるイデアルで割った剰余環です：

$$R[\varepsilon] = R[x]/(x^2), \qquad \varepsilon = x + (x^2), \qquad \varepsilon^2 = 0 .$$

各剰余類には次数 1 以下の代表元がちょうど 1 つあるので、すべての元は一意な $a, b \in R$ を用いて $a + b\varepsilon$ と表されます。これが組 (value, tangent) です。$R[\varepsilon]$ は可換環をイデアルで割った剰余環なので、それ自体が可換環であり、その演算は多項式の演算を行った後に $\varepsilon^2$ を落としたものです：

$$\begin{aligned} (a + b\varepsilon) + (c + d\varepsilon) &= (a + c) + (b + d)\varepsilon, \\ -(a + b\varepsilon) &= -a - b\varepsilon, \\ (a + b\varepsilon)(c + d\varepsilon) &= ac + (ad + bc)\varepsilon + bd\,\varepsilon^2 = ac + (ad + bc)\varepsilon . \end{aligned}$$

単元と除算

$a + b\varepsilon$ は、$a$ が可逆であるときにちょうど可逆です。$a$ が単元ならば

$$(a + b\varepsilon)(a^{-1} - a^{-2} b\,\varepsilon) = 1 - a^{-1} b\,\varepsilon + a^{-1} b\,\varepsilon - a^{-2} b^2 \varepsilon^2 = 1 .$$

逆に、$(a + b\varepsilon)(c + d\varepsilon) = 1$ ならば $ac = 1$ なので、$a$ は単元です。特に $\varepsilon$ は逆元を持たない 0 でない元であり、$\varepsilon \cdot \varepsilon = 0$ より零因子です。したがって $R$ が体であっても $R[\varepsilon]$ は決して体になりません。単元 $c$ に対して、商は

$$\frac{a + b\varepsilon}{c + d\varepsilon} = (a + b\varepsilon)(c^{-1} - c^{-2} d\,\varepsilon) = \frac{a}{c} + \frac{bc - ad}{c^2}\,\varepsilon ,$$

となり、これが Dual::div と Dual::div_checked の実装する式です。

多項式：導関数はどこから来るのか

$n \ge 1$ に対し、二項定理と $\varepsilon^2 = 0$ から

$$(a + b\varepsilon)^n = \sum_{k=0}^{n} \binom{n}{k} a^{n-k} (b\varepsilon)^k = a^n + n a^{n-1} b\,\varepsilon .$$

線形性により、$R$ に係数を持つすべての多項式 $p(x) = \sum_k c_k x^k$ は次を満たします

$$p(a + b\varepsilon) = \sum_k c_k a^k + \Big(\sum_k k\,c_k a^{k-1}\Big) b\,\varepsilon = p(a) + p'(a)\,b\,\varepsilon ,$$

ここで $p'$ は形式的導関数です。極限はまったく必要ありません。この恒等式は整数を含むすべての可換環で成り立ちます。

滑らかな関数と連鎖律

多項式でない微分可能な関数 $f$ について、パッケージは同じ恒等式によって双対数への拡張を定義します。

$$f(a + b\varepsilon) := f(a) + f'(a)\,b\,\varepsilon ,$$

これは無限小 $h = b\varepsilon$ による 1 次のテイラー展開 $f(a + h) = f(a) + f'(a)h + O(h^2)$ であり、$\varepsilon^2 = 0$ なので剰余項は消えます。すると合成によって、追加のコードなしに連鎖律が得られます。微分可能な $g$ と $f$ について：

$$\begin{aligned} f\big(g(a + b\varepsilon)\big) &= f\big(g(a) + g'(a)\,b\,\varepsilon\big) \\ &= f(g(a)) + f'(g(a))\,g'(a)\,b\,\varepsilon \\ &= (f \circ g)(a) + (f \circ g)'(a)\,b\,\varepsilon . \end{aligned}$$

和・積・商の規則は、上の環演算を $b = u'(x)$、$d = v'(x)$ として読んだものです：

$$(u + v)' = u' + v', \qquad (uv)' = u'v + uv', \qquad \Big(\frac{u}{v}\Big)' = \frac{u'v - uv'}{v^2} .$$

パッケージの初等関数の規則は、この定義を既知の導関数に適用したものです：

メソッド	$f(a)$	$f'(a)$	実際に計算される接成分
sqrt	$\sqrt a$	$\frac{1}{2\sqrt a}$	$r = \sqrt a$ として $b / (2r)$
exp	$e^a$	$e^a$	$v = e^a$ として $b \cdot v$
exp2	$2^a$	$2^a \ln 2$	$v = 2^a$ として $b \cdot v \cdot \ln 2$
ln	$\ln a$	$1/a$	$b / a$
log2	$\log_2 a$	$\frac{1}{a \ln 2}$	$b / (a \ln 2)$
log10	$\log_{10} a$	$\frac{1}{a \ln 10}$	$b / (a \ln 10)$
sin	$\sin a$	$\cos a$	$b \cos a$
cos	$\cos a$	$-\sin a$	$-(b \sin a)$
tan	$\tan a$	$\sec^2 a$	$b / (\cos a \cdot \cos a)$

値への射影は準同型である

写像 $\pi : R[\varepsilon] \to R$、$a + b\varepsilon \mapsto a$ は $0$、$1$、$+$、$-$、$\times$ を保ちます。$\pi\big((a + b\varepsilon)(c + d\varepsilon)\big) = ac = \pi(a + b\varepsilon)\,\pi(c + d\varepsilon)$ であり、他の演算も同様です。初等関数の規則もすべて、定義により $\pi(f(z)) = f(\pi(z))$ を満たします。したがって、双対数上の任意の計算の value は、値に対して同じ計算を行った結果とちょうど一致します。接成分を加えてもプライマルの結果は変わらず、その丸めさえ変わりません。

設計上の判断

スカラーに対する 1 つのジェネリック型

問題。 導関数は Double、Float、および他の Luna Flow パッケージで定義されたスカラー型について必要です。

選択肢。 Double 専用の双対数型、または各演算が使う最小のトレイト集合を要求するジェネリックな Dual[T]。

選択。 Dual[T] はジェネリックであり、各メソッドが独自の境界を持ちます（Dual::mul は Add + Mul だけを、Dual::exp2 は Exponential + Logarithmic + IntegralHomomorphism + Mul を必要とします）。これは最小のトレイト合成に依存するという Luna Flow の規則に従っています。また T 自体を双対数にすることもでき、入れ子にすることで高階導関数が得られます（forward の設計 を参照）。

環レベルのインスタンスのみ

問題。 Dual[T] は luna-generic のどの構造トレイトを実装してよいでしょうか。

選択。 Zero、One、AddMonoid、AddGroup、MulMonoid、Semiring、Ring を、それぞれ T が同じ構造を持つときに実装します。これらは等式クラスです。その法則は項の間の恒等式であり、恒等式は剰余環 $T[x]/(x^2)$ に継承されます。確認として、積の規則の結合法則を示します：

$$\begin{aligned} \big((a + b\varepsilon)(c + d\varepsilon)\big)(e + f\varepsilon) &= (ac + (ad + bc)\varepsilon)(e + f\varepsilon) = ace + (acf + ade + bce)\varepsilon, \\ (a + b\varepsilon)\big((c + d\varepsilon)(e + f\varepsilon)\big) &= (a + b\varepsilon)(ce + (cf + de)\varepsilon) = ace + (acf + ade + bce)\varepsilon . \end{aligned}$$

Field、MulGroup、Inverse は却下しました。単元の計算 が示すとおり $\varepsilon \ne 0$ は逆元を持たないので、Inverse インスタンスはそれに対して誤った値を返さざるを得ないからです。

順序を持たない

問題。 ジェネリックなコードはしばしば < で分岐します。

選択肢。 value だけで順序付ける、(value, tangent) で辞書式に順序付ける、順序を持たない。

選択。 Compare インスタンスは持ちません。値による順序は、導出された Eq に関して反対称ではありません（$a + b\varepsilon$ と $a + c\varepsilon$ は等しくないのに互いに $\le$ になります）。辞書式順序は全順序ですが環の順序ではありません。$\varepsilon > 0$ となり、順序環では $x > 0, y > 0 \Rightarrow xy > 0$ が要求されますが、$\varepsilon \cdot \varepsilon = 0$ だからです。分岐するコードは x.value() を明示的に比較します。これにより、導関数が選ばれた分岐の導関数であることも見てわかるようになります。

体ではないがチェックなしの Div を持つ

問題。 通常の数式は / を使い、MoonBit の / 演算子には Div トレイトが必要です。しかし $R[\varepsilon]$ 上の除算は部分的にしか定義されません。

選択。 Dual[T] は商の公式を用いてチェックなしで Div を実装し、除数が 0 の場合の振る舞いは T のものを引き継ぎます（Double では無限大または NaN）。これは arithmetic のチェックなしの層を反映したもので、そこでは Sqrt などがスカラーの IEEE の意味論に従います。チェック付きの層は DivChecked で、T 自身の DivChecked の上に実装されています。Div だけでは構造の主張にはなりません。この型は依然として Field を実装しません。

チェック付きの形式は arithmetic を再利用する

問題。 ゼロ除算や負の数の平方根を、データとして報告できなければなりません。

選択肢。 autodiff 専用のエラー型、または arithmetic のトレイトとエラー値。

選択。 Dual[T] は Luna-Flow/arithmetic の DivChecked と SqrtChecked を実装し、その ArithmeticError を返します。これにより呼び出し側は、双対数とスカラーの失敗を 1 つの語彙で扱えます。値と接成分の計算はどちらもチェックされ、最初の失敗が返されます。チェック付きの形式を持つのは除算と平方根だけです。これは arithmetic が定義するチェック付きトレイトがその 2 つだからです。対数関数と三角関数は T のチェックなしの意味論に従います。

整数を経由した定数

問題。 sqrt、exp2、log2、log10 の規則には、T における定数 $2$ と $10$ が必要です。

選択肢。 Double からの変換を要求する、one + one として構築する、整数からの標準的な写像を使う。

選択。 IntegralHomomorphism::from_integral(2) と from_integral(10) です。写像 $\mathbb Z \to T$ は唯一の環準同型なので、どの環でも正しい定数を指し、$2$ や $10$ のような小さな整数はどの数値インスタンスでも厳密に表せます。Double からの変換は厳密な型には存在せず、繰り返しの加算は $10$ に対してはコストが大きくなります。

計算した値を再利用する

exp と exp2 は $v = f(a)$ を 1 回計算し、$f' = f$（係数 $\ln 2$ を除いて）なので $v$ から接成分を作ります。sqrt も同様に計算した平方根で割ります。これにより 1 回の呼び出しにつき初等関数の評価が 1 回節約され、接成分が返される値と整合します。

tan の導関数

$\tan' a$ を与える式は 2 つあります。値を再利用する $1 + \tan^2 a$ と、余弦をもう 1 回必要とする $\sec^2 a = 1/\cos^2 a$ です。パッケージは $b / (\cos a \cdot \cos a)$ を使います。両者は数学的には等しく、選んだ形式は $\tan$ 自体とちょうど同じ点で未定義になります。

明示的なメソッド昇格

MoonBit 0.10 では、トレイトインスタンスはもはや暗黙的にメソッドを生成しません。src/dual/extends.mbt は算術演算子、equal、zero、one、div_checked を昇格し、既存の呼び出し側のために not_equal、to_repr、from_nat、from_integral、pi、e、tau を隠された非推奨の形式として残しています（非推奨 を参照）。

正しさと不変条件

プログラム全体の導関数

プログラムが有限個のステップで $y = f(x)$ を計算するとし、各ステップは環演算、除算、または表中の初等関数のいずれかであるとします。これを $x + 1\varepsilon$（Dual::variable(x)）上で、他の入力はすべて定数として実行します。すると、すべての中間値 $v_k$ は $v_k(x) + v_k'(x)\,\varepsilon$ として表されます。

ステップに関する帰納法による証明。 入力：変数は $x + 1\varepsilon$ で $x' = 1$、定数は $c + 0\varepsilon$ で $c' = 0$ です。ステップ：オペランドが $u + u'\varepsilon$ と $w + w'\varepsilon$ ならば、和・積・商の公式は上で導いた規則により $(u \circ w) + (u \circ w)' \varepsilon$ を与え、初等関数は連鎖律により $g(u) + g'(u)\,u'\,\varepsilon = g(u) + (g \circ u)'\varepsilon$ を与えます。最後の中間値は $y$ なので、tangent は $f'(x)$ です。$\square$

入力の接成分を $b$ とした同じ議論により $f'(x)\,b$ が得られ、複数の入力をベクトル $v$ でシードすると方向微分 $\nabla f(x) \cdot v$ が得られます（linalg の設計 を参照）。得られる導関数はプログラムの導関数であって、それが近似する数学的関数の導関数ではありません。value() に基づく分岐では、選ばれた分岐が微分されます。

丸め誤差

各浮動小数点演算について $|\delta| \le u$ として $\mathrm{fl}(x \circ y) = (x \circ y)(1 + \delta)$ と書き（Double では $u = 2^{-53}$）、$k$ 個の因子 $(1 + \delta_i)^{\pm 1}$ の積は $|\theta_k| \le \gamma_k = ku/(1 - ku)$ を満たす $1 + \theta_k$ に等しいという標準的な補題を用います。¹1 N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002, Lemma 3.1.

積の接成分。 Dual::mul は $\mathrm{fl}(\mathrm{fl}(ad) + \mathrm{fl}(cb))$ を計算します：

$$\begin{aligned} \hat t &= \big(ad(1 + \delta_1) + cb(1 + \delta_2)\big)(1 + \delta_3) = ad(1 + \theta_2) + cb(1 + \theta_2'), \\ |\hat t - (ad + bc)| &\le \gamma_2 \big(|ad| + |bc|\big). \end{aligned}$$

商の接成分。 Dual::div は $\mathrm{fl}\big(\mathrm{fl}(bc - ad) / \mathrm{fl}(c \cdot c)\big)$ を計算します。丸めは 5 回で、そのうち 1 回は分母にあります：

$$\begin{aligned} \hat t &= \frac{\big(bc(1 + \delta_1) - ad(1 + \delta_2)\big)(1 + \delta_3)}{c^2 (1 + \delta_4)}\,(1 + \delta_5) = \frac{bc(1 + \theta_4) - ad(1 + \theta_4')}{c^2}, \\ \Big|\hat t - \frac{bc - ad}{c^2}\Big| &\le \gamma_4\,\frac{|bc| + |ad|}{c^2} . \end{aligned}$$

どちらの上界も、導関数の式を直接評価する場合の誤差と同じ形をしています。相対誤差が大きくなるのは $ad$ と $bc$ が打ち消し合う場合だけであり、これは導関数そのものの条件数の問題です。これらはどの中間値もオーバーフローやアンダーフローを起こさないことを仮定しています。特に $c^2$ は $c$ よりずっと早くアンダーフローします（Dual::div_checked の警告を参照）。初等関数の規則では、接成分の誤差は T による $f'(a)$ の実装の誤差（たとえば sin における cos の誤差）に高々 2 回の丸めを加えたものです。

有限差分との比較

前進差分 $D_h f(x) = (f(x + h) - f(x))/h$ には 2 つの誤差要因があります。テイラーの定理から打ち切り誤差が生じ、さらに浮動小数点で誤差高々 $u|f|$ で $f$ を評価することで桁落ち誤差が加わります：

$$\big|D_h \hat f(x) - f'(x)\big| \le \frac{h}{2}\,|f''(\xi)| + \frac{2u\,|f(x)|}{h} .$$

右辺の $h$ に関する導関数を 0 とおくと、

$$\frac{|f''|}{2} - \frac{2u|f|}{h^2} = 0 \quad\Longrightarrow\quad h^\ast = 2\sqrt{\frac{u|f|}{|f''|}}, \qquad \text{error}(h^\ast) = 2\sqrt{u\,|f|\,|f''|} ,$$

そのため最良の刻み幅でも有効数字の約半分が失われます(Double では $\sqrt u \approx 10^{-8}$)。中心差分 $(f(x + h) - f(x - h))/(2h)$ の打ち切り誤差は $\frac{h^2}{6}|f'''|$ で、約 $u^{2/3} \approx 10^{-11}$ に達します。双対数には打ち切り項がまったくありません。上の帰納法により、双対数はプログラムの正確な導関数を計算し、唯一の誤差は前節の上界に従う丸めだけです。²2 A. Griewank and A. Walther, Evaluating Derivatives, 2nd ed., SIAM, 2008 の第 2–3 章では、フォワードモードとその誤差解析を詳しく扱っています。 テストスイートと dual チュートリアル で両方の手法を比較しています。

コスト

加算、減算、符号反転は T の演算 1 回ではなく 2 回かかります。乗算は乗算 3 回と加算 1 回、除算は乗算 3 回、減算 1 回、除算 2 回かかります。初等関数はその関数と導関数の因子の分だけかかります。したがってプログラムは Dual[T] 上では T 上より高々小さな定数倍（約 3〜4 倍）遅く実行され、定数倍多くのメモリを必要とします。

守るべき法則

• 任意の結果の x.value() は、値に対して同じ計算を行った結果に等しい（射影は準同型である）。
• 定数、zero()、one()、from_integral、from_nat、および Constants インスタンスの接成分は 0 である。
• インスタンスは Ring までで、Field、MulGroup、Inverse、Compare はない。
• 積の規則は T の乗算が可換であることを仮定しており、Luna Flow のすべての数値インスタンスはこれを満たす。

却下した代替案

• 有限差分。 単純ですが、上で導いたように桁の半分以上を失い、問題ごとに刻み幅が必要です。
• 記号微分。 導関数の式を構築するには項の言語と簡約が必要で、式は急速に膨らむことがあります。これは CAS のレイヤーに属するもので、スカラー型の役割ではありません。
• リバースモード。 多数の入力に対する勾配には適していますが、計算を記録するテープやクロージャが必要です。実装されていません。
• 接成分のベクトルを格納する。 1 回のパスで勾配全体が得られますが、ベクトル値の型と、演算ごとのメモリ確保というコストがかかります。スカラーの接成分にすることで Dual[T] は単純な 2 フィールドの値に保たれます。ベクトルの接成分は今後の課題です。

境界

• 1 階導関数のみです。高階導関数は Dual[Dual[T]] を入れ子にすることで得られ、1 階ごとに $2^k$ のコストがかかります。打ち切りテイラー（ジェット）型はありません。
• リバースモードも記号微分もありません。
• Field、MulGroup、Inverse、順序のインスタンスはなく、Show もありません。
• 定義域をチェックするのは div_checked と sqrt_checked だけです。対数関数、三角関数、チェックなしの演算子は T に従います。
• 検証付きの包含はありません。上記の丸めの上界は事前評価であり、実行時に計算されるものではありません。T としての区間やボールのスカラーは、このリポジトリではテストも文書化もされていません。

Footnotes

1. N. J. Higham, Accuracy and Stability of Numerical Algorithms, 2nd ed., SIAM, 2002, Lemma 3.1. ↩

2. A. Griewank and A. Walther, Evaluating Derivatives, 2nd ed., SIAM, 2008 の第 2–3 章では、フォワードモードとその誤差解析を詳しく扱っています。 ↩