dual 设计

本页解释 Dual[T] 背后的数学，以及为什么这个类型恰好具有 dual API 中列出的运算和实例。dual 教程 在不涉及这些背景的情况下使用该类型。

设计目标

通过在另一种数值类型上运行同一个程序，计算针对 Luna Flow trait 编写的普通程序的精确一阶导数。该类型必须满足其实例所声明的代数定律，通过共享的 arithmetic 错误值报告定义域错误，并且不依赖任何容器库或多项式库。

数学背景

对偶数的代数

设 $R$ 为交换环（对于浮点类型，指它们所近似的实数）。$R$ 上的对偶数是多项式环 $R[x]$ 对由 $x^2$ 生成的理想所作的商：

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

每个等价类恰好有一个次数至多为一的代表元，因此每个元素都可唯一地写成 $a + b\varepsilon$，其中 $a, b \in R$。这就是二元组 (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$ 是一个没有逆元的非零元素，而 $\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$ 时的一阶泰勒展开 $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}$	$b / (2r)$，其中 $r = \sqrt a$
exp	$e^a$	$e^a$	$b \cdot v$，其中 $v = e^a$
exp2	$2^a$	$2^a \ln 2$	$b \cdot v \cdot \ln 2$，其中 $v = 2^a$
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 都恰好等于在值上进行同一计算的结果：附加切向分量绝不会改变原始结果，连其舍入也不会改变。

设计决策

对标量泛型的单一类型

问题。 Double、Float 以及其他 Luna Flow 包中定义的标量类型都需要求导。

选项。 只支持 Double 的对偶类型；或者泛型的 Dual[T]，其每个运算只要求它用到的最小 trait 集合。

选择。 Dual[T] 是泛型的，每个方法都带有自己的约束（Dual::mul 只需要 Add + Mul，Dual::exp2 需要 Exponential + Logarithmic + IntegralHomomorphism + Mul）。这遵循了 Luna Flow 依赖最小 trait 组合的规则。它还允许 T 本身就是对偶数，从而通过嵌套得到高阶导数（见 forward 设计）。

只提供环层面的实例

问题。 Dual[T] 可以实现 luna-generic 中的哪些结构 trait？

选择。 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 trait；但除法在 $R[\varepsilon]$ 上是部分定义的。

选择。 Dual[T] 用商公式实现 Div，不做检查，在除数为零时沿用 T 的行为（对于 Double，是无穷大或 NaN）。这与 arithmetic 中不带检查的那一层相对应，在那里 Sqrt 等 trait 遵循标量的 IEEE 语义。带检查的那一层是 DivChecked，它基于 T 自身的 DivChecked 实现。仅有 Div 并不构成结构上的声明；该类型仍然不实现 Field。

带检查的形式复用 arithmetic

问题。 除以零和负数的平方根必须能够作为数据报告。

选项。 专门的 autodiff 错误类型；arithmetic 的 trait 和错误值。

选择。 Dual[T] 实现 Luna-Flow/arithmetic 中的 DivChecked 和 SqrtChecked，并返回其 ArithmeticError，因此调用方可以用同一套词汇处理对偶数和标量的失败。值和切向分量的计算都会被检查，并返回第一个失败。只有除法和平方根有带检查的形式，因为 arithmetic 定义的带检查 trait 就是这两个；对数和三角函数沿用 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)$，并由 $v$ 构成切向分量，因为 $f' = f$（至多相差因子 $\ln 2$）。sqrt 同样除以已计算出的根。这样每次调用可省去一次初等函数求值，并使切向分量与返回的值保持一致。

tan 的导数

有两个公式可以给出 $\tan' a$：复用值的 $1 + \tan^2 a$，以及需要多算一次余弦的 $\sec^2 a = 1/\cos^2 a$。本包使用 $b / (\cos a \cdot \cos a)$。两者在数学上相等；所选形式恰好在与 $\tan$ 本身相同的点上无定义。

显式方法提升

从 MoonBit 0.10 起，trait 实例不再隐式生成方法。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() 的分支只会对所走的那个分支求导。

舍入误差

对每个浮点运算记 $\mathrm{fl}(x \circ y) = (x \circ y)(1 + \delta)$，其中 $|\delta| \le u$（对于 Double，$u = 2^{-53}$），并使用标准引理：$k$ 个因子 $(1 + \delta_i)^{\pm 1}$ 之积等于 $1 + \theta_k$，其中 $|\theta_k| \le \gamma_k = ku/(1 - ku)$。¹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)$，共五次舍入，其中一次在分母中：

$$\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 的误差），再加上至多两次舍入。

与有限差分的比较

前向差分 $D_h f(x) = (f(x + h) - f(x))/h$ 有两个误差来源。泰勒定理给出截断误差，而在浮点数中以至多 $u|f|$ 的误差计算 $f$ 会带来抵消误差：

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

令右端关于 $h$ 的导数为零，

$$\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 与 A. Walther 所著 Evaluating Derivatives（第 2 版，SIAM，2008）第 2–3 章完整论述了前向模式及其误差分析。 测试套件和 dual 教程 对两种方法进行了比较。

代价

加法、减法和取负需要两次 T 运算而非一次；乘法需要三次乘法和一次加法；除法需要三次乘法、一次减法和两次除法；初等函数需要计算函数本身及其导数因子。因此程序在 Dual[T] 上运行时至多比在 T 上慢一个较小的常数倍（约三到四倍），并需要常数倍的额外内存。

需要保持的定律

• 任何结果的 x.value() 都等于在值上进行同一计算的结果（投影是同态）。
• 常数、zero()、one()、from_integral、from_nat 以及 Constants 实例的切向分量都为零。
• 实例止于 Ring；没有 Field、MulGroup、Inverse 或 Compare。
• 乘积法则假定 T 的乘法满足交换律，Luna Flow 中的所有数值实例都满足这一点。

被否决的方案

• 有限差分。 简单，但如上所推导，它会损失一半甚至更多的有效数字，并且每个问题都需要选择步长。
• 符号求导。 构建导数表达式需要项语言和化简，而且表达式可能迅速膨胀；这属于 CAS 层，而不属于标量类型。
• 反向模式。 对于多输入的梯度更有优势，但它需要记录计算过程的磁带（tape）或闭包；目前尚未实现。
• 存储切向分量向量。 一遍计算就能得到整个梯度，代价是引入向量值类型，并且每次运算都要分配内存。标量切向分量使 Dual[T] 保持为简单的双字段值；向量切向分量留待将来实现。

边界

• 只支持一阶导数。高阶导数通过嵌套 Dual[Dual[T]] 得到，其代价按阶数呈 $2^k$ 增长；没有截断泰勒（jet）类型。
• 没有反向模式，也没有符号求导。
• 没有 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 与 A. Walther 所著 Evaluating Derivatives（第 2 版，SIAM，2008）第 2–3 章完整论述了前向模式及其误差分析。 ↩