dual design

This page explains the mathematics of Dual[T] and why the type has exactly the operations and instances listed in the dual API. The dual tutorial uses the type without this background.

Design goal

Compute exact first derivatives of ordinary programs, written against the Luna Flow traits, by running the same program on a different number type. The type has to satisfy the algebraic laws its instances advertise, report domain failures through the shared arithmetic error values, and stay independent of any container or polynomial library.

Mathematical background

The algebra of dual numbers

Let $R$ be a commutative ring (for floating-point types, the real numbers they approximate). The dual numbers over $R$ are the quotient of the polynomial ring $R[x]$ by the ideal generated by $x^2$:

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

Every class has exactly one representative of degree at most one, so every element is $a + b\varepsilon$ with unique $a, b \in R$. This is the pair (value, tangent). Because $R[\varepsilon]$ is a quotient of a commutative ring by an ideal, it is itself a commutative ring, and its operations are those of polynomials followed by dropping $\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}$$

Units and division

$a + b\varepsilon$ is invertible exactly when $a$ is. If $a$ is a unit, then

$$(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 .$$

Conversely, if $(a + b\varepsilon)(c + d\varepsilon) = 1$ then $ac = 1$, so $a$ is a unit. In particular $\varepsilon$ is a non-zero element without an inverse, and $\varepsilon \cdot \varepsilon = 0$ makes it a zero divisor: $R[\varepsilon]$ is never a field, even when $R$ is. For a unit $c$ the quotient is

$$\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 ,$$

which is the formula Dual::div and Dual::div_checked implement.

Polynomials: where derivatives come from

For $n \ge 1$ the binomial theorem and $\varepsilon^2 = 0$ give

$$(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 .$$

By linearity, every polynomial $p(x) = \sum_k c_k x^k$ with coefficients in $R$ satisfies

$$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 ,$$

where $p'$ is the formal derivative. Nothing here needs limits: the identity holds in every commutative ring, including the integers.

Smooth functions and the chain rule

For a differentiable function $f$ that is not a polynomial, the package defines the extension to dual numbers by the same identity,

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

which is the first-order Taylor expansion $f(a + h) = f(a) + f'(a)h + O(h^2)$ with the infinitesimal $h = b\varepsilon$; the remainder vanishes because $\varepsilon^2 = 0$. Composition then produces the chain rule without any extra code. For differentiable $g$ and $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}$$

The sum, product and quotient rules are the ring operations above read with $b = u'(x)$ and $d = v'(x)$:

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

The elementary rules of the package are this definition applied to known derivatives:

Method	$f(a)$	$f'(a)$	Tangent as computed
sqrt	$\sqrt a$	$\frac{1}{2\sqrt a}$	$b / (2r)$ with $r = \sqrt a$
exp	$e^a$	$e^a$	$b \cdot v$ with $v = e^a$
exp2	$2^a$	$2^a \ln 2$	$b \cdot v \cdot \ln 2$ with $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)$

The value projection is a homomorphism

The map $\pi : R[\varepsilon] \to R$, $a + b\varepsilon \mapsto a$, preserves $0$, $1$, $+$, $-$ and $\times$: $\pi\big((a + b\varepsilon)(c + d\varepsilon)\big) = ac = \pi(a + b\varepsilon)\,\pi(c + d\varepsilon)$, and likewise for the others. Every elementary rule also has $\pi(f(z)) = f(\pi(z))$ by definition. Hence the value of any computation on dual numbers is exactly the result of the same computation on the values: adding tangents never changes the primal result, not even its rounding.

Design decisions

One generic type over the scalar

Problem. Derivatives are needed for Double, Float, and for scalar types defined in other Luna Flow packages.

Options. A Double-only dual type; a generic Dual[T] whose operations ask for the smallest trait set they use.

Choice. Dual[T] is generic, and every method carries its own bound (Dual::mul needs only Add + Mul, Dual::exp2 needs Exponential + Logarithmic + IntegralHomomorphism + Mul). This follows the Luna Flow rule of depending on the smallest trait composition. It also lets T itself be a dual number, which gives higher derivatives by nesting (see the forward design).

Ring-level instances only

Problem. Which structure traits from luna-generic may Dual[T] implement?

Choice. Zero, One, AddMonoid, AddGroup, MulMonoid, Semiring and Ring, each when T has the same structure. These are equational classes: their laws are identities between terms, and identities are inherited by the quotient $T[x]/(x^2)$. As a check, associativity of the product rule:

$$\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 and Inverse are rejected because the units computation shows that $\varepsilon \ne 0$ has no inverse: an Inverse instance would have to return a wrong value for it.

No ordering

Problem. Generic code often branches on <.

Options. Order by value only; order lexicographically by (value, tangent); no order.

Choice. No Compare instance. Ordering by value is not antisymmetric with respect to the derived Eq ($a + b\varepsilon$ and $a + c\varepsilon$ would each be $\le$ the other without being equal). The lexicographic order is total but is not a ring order: it makes $\varepsilon > 0$, and an ordered ring requires $x > 0, y > 0 \Rightarrow xy > 0$, while $\varepsilon \cdot \varepsilon = 0$. Code that branches compares x.value() explicitly, which also makes visible that the derivative is the derivative of the branch taken.

An unchecked Div despite not being a field

Problem. Ordinary formulas use /, and MoonBit’s / operator needs the Div trait; but division is partial on $R[\varepsilon]$.

Choice. Dual[T] implements Div with the quotient formula and no check, inheriting T’s behaviour on a zero divisor (for Double, infinities or NaN). This mirrors the unchecked tier of arithmetic, where Sqrt and friends follow the IEEE semantics of the scalar. The checked tier is DivChecked, implemented on top of T’s own DivChecked. Div alone is not a structure claim; the type still does not implement Field.

Checked forms reuse arithmetic

Problem. Division by zero and square roots of negative numbers must be reportable as data.

Options. A dedicated autodiff error type; the arithmetic traits and error values.

Choice. Dual[T] implements DivChecked and SqrtChecked from Luna-Flow/arithmetic and returns its ArithmeticError, so callers handle dual and scalar failures with one vocabulary. Both the value and the tangent computation are checked, and the first failure is returned. Only division and square root have checked forms, because those are the checked traits arithmetic defines; logarithms and trigonometric functions follow the unchecked semantics of T.

Constants through the integers

Problem. The rules for sqrt, exp2, log2 and log10 need the constants $2$ and $10$ in T.

Options. Require a conversion from Double; build them as one + one; use the canonical map from the integers.

Choice. IntegralHomomorphism::from_integral(2) and from_integral(10). The map $\mathbb Z \to T$ is the unique ring homomorphism, so it names the right constant in every ring, and small integers such as $2$ and $10$ are exact in every numeric instance. A Double conversion would not exist for exact types, and repeated addition costs more for $10$.

Reuse the computed value

exp and exp2 compute $v = f(a)$ once and form the tangent from $v$, because $f' = f$ (up to the factor $\ln 2$). sqrt likewise divides by the computed root. This saves one elementary evaluation per call and makes the tangent consistent with the returned value.

The derivative of tan

Two formulas give $\tan' a$: $1 + \tan^2 a$, which reuses the value, and $\sec^2 a = 1/\cos^2 a$, which needs one more cosine. The package uses $b / (\cos a \cdot \cos a)$. Both are mathematically equal; the chosen form is undefined at exactly the same points as $\tan$ itself.

Explicit method promotion

With MoonBit 0.10, trait instances no longer create methods implicitly. src/dual/extends.mbt promotes the arithmetic operators, equal, zero, one and div_checked, and keeps not_equal, to_repr, from_nat, from_integral, pi, e and tau as hidden deprecated forms for existing callers (see Deprecated).

Correctness and invariants

Derivatives of whole programs

Let a program compute $y = f(x)$ by a finite sequence of steps, each a ring operation, a division or one of the elementary functions of the table. Run it on $x + 1\varepsilon$ (Dual::variable(x)) with every other input a constant. Then every intermediate $v_k$ is represented as $v_k(x) + v_k'(x)\,\varepsilon$.

Proof by induction on the steps. Inputs: the variable is $x + 1\varepsilon$ and $x' = 1$; a constant is $c + 0\varepsilon$ and $c' = 0$. Step: if the operands are $u + u'\varepsilon$ and $w + w'\varepsilon$, the sum, product and quotient formulas give $(u \circ w) + (u \circ w)' \varepsilon$ by the rules derived above, and an elementary function gives $g(u) + g'(u)\,u'\,\varepsilon = g(u) + (g \circ u)'\varepsilon$ by the chain rule. The last intermediate is $y$, so tangent is $f'(x)$. $\square$

The same argument with input tangent $b$ gives $f'(x)\,b$, and with several inputs seeded by a vector $v$ it gives the directional derivative $\nabla f(x) \cdot v$ (see the linalg design). The derivative is that of the program, not of a mathematical function it approximates: a branch on value() differentiates the branch taken.

Rounding error

Write $\mathrm{fl}(x \circ y) = (x \circ y)(1 + \delta)$ with $|\delta| \le u$ for each floating-point operation ($u = 2^{-53}$ for Double), and use the standard lemma that a product of $k$ factors $(1 + \delta_i)^{\pm 1}$ equals $1 + \theta_k$ with $|\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.

Product tangent. Dual::mul computes $\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}$$

Quotient tangent. Dual::div computes $\mathrm{fl}\big(\mathrm{fl}(bc - ad) / \mathrm{fl}(c \cdot c)\big)$, five roundings with one of them in the denominator:

$$\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}$$

Both bounds have the shape of the error of evaluating the derivative formula directly; the relative error is only large when $ad$ and $bc$ cancel, which is the conditioning of the derivative itself. They assume that no intermediate overflows or underflows; in particular $c^2$ underflows long before $c$ does (see the warning on Dual::div_checked). For an elementary rule the tangent error is the error of T’s implementation of $f'(a)$ (for example of cos in sin) plus at most two more roundings.

Comparison with finite differences

The forward difference $D_h f(x) = (f(x + h) - f(x))/h$ has two error sources. Taylor’s theorem gives the truncation error, and evaluating $f$ in floating point with error at most $u|f|$ adds a cancellation error:

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

Setting the derivative of the right-hand side with respect to $h$ to zero,

$$\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''|} ,$$

so even the best step loses about half of the significant digits ($\sqrt u \approx 10^{-8}$ for Double). The central difference $(f(x + h) - f(x - h))/(2h)$ has truncation error $\frac{h^2}{6}|f'''|$ and reaches about $u^{2/3} \approx 10^{-11}$. Dual numbers have no truncation term at all: by the induction above they compute the exact derivative of the program, and the only error is rounding with the bounds of the previous section.²2 A. Griewank and A. Walther, Evaluating Derivatives, 2nd ed., SIAM, 2008, chapters 2–3, treat forward mode and its error analysis in full. The test suite and the dual tutorial compare both methods.

Cost

Addition, subtraction and negation cost two operations of T instead of one; multiplication costs three multiplications and one addition; division costs three multiplications, one subtraction and two divisions; an elementary function costs the function plus its derivative factor. A program therefore runs at most a small constant factor (about three to four) slower on Dual[T] than on T, and needs a constant factor more memory.

Laws to keep

• x.value() of any result equals the same computation on values (the projection is a homomorphism).
• Constants, zero(), one(), from_integral, from_nat and the Constants instance have tangent zero.
• Instances stop at Ring; no Field, MulGroup, Inverse or Compare.
• The product rule assumes commutative multiplication of T, which every numeric instance in Luna Flow has.

Alternatives rejected

• Finite differences. Simple, but they lose half the digits or more, as derived above, and need a step size per problem.
• Symbolic differentiation. Building derivative expressions needs a term language and simplification, and expressions can grow quickly; that belongs to a CAS layer, not to a scalar type.
• Reverse mode. Better for gradients of many inputs, but it needs a tape or closures that record the computation; it is not implemented.
• Storing a vector of tangents. One pass would give a whole gradient, at the cost of a vector-valued type and allocation on every operation. The scalar tangent keeps Dual[T] a plain two-field value; vector tangents are future work.

Boundaries

• First derivatives only. Higher derivatives come from nesting Dual[Dual[T]], which costs $2^k$ per order; there is no truncated Taylor (jet) type.
• No reverse mode and no symbolic differentiation.
• No Field, MulGroup, Inverse or order instance, and no Show.
• Only div_checked and sqrt_checked check their domain. Logarithms, trigonometric functions and the unchecked operators follow T.
• No validated enclosures: the rounding bounds above are a priori, not computed at run time. Interval or ball scalars as T are neither tested nor documented by this repository.

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, chapters 2–3, treat forward mode and its error analysis in full. ↩