Forward-mode automatic differentiation via dual numbers

Dual numbers compute a function’s value and its exact derivative together, with no finite-difference step-size tradeoff between truncation error (large h) and floating-point cancellation (small h).

import math

import numpy as np

from mathematicskit.calculus.systems.dual_numbers import derivative
from mathematicskit.calculus.systems.finite_differences import central_difference

Exact derivative, compared against central differences at various h

f_dual = lambda x: (x * x).sin()  # d/dx sin(x^2) = 2x cos(x^2)
f_plain = lambda x: math.sin(x * x)
x0 = 1.3
exact = 2.0 * x0 * math.cos(x0**2)

dual_result = derivative(f_dual, x0)
print(f"dual-number derivative: {dual_result:.15f} (exact: {exact:.15f}, error {abs(dual_result - exact):.2e})")

for h in (1e-1, 1e-3, 1e-5, 1e-7, 1e-9, 1e-11):
    cd = central_difference(f_plain, x0, h=h)
    print(f"central diff h={h:.0e}: {cd:.15f}, error {abs(cd - exact):.2e}")

dual-number derivative: -0.309196080171192 (exact: -0.309196080171192, error 0.00e+00)
central diff h=1e-01: -0.331234137017592, error 2.20e-02
central diff h=1e-03: -0.309198313356795, error 2.23e-06
central diff h=1e-05: -0.309196080394702, error 2.24e-10
central diff h=1e-07: -0.309196080405805, error 2.35e-10
central diff h=1e-09: -0.309196113157384, error 3.30e-08
central diff h=1e-11: -0.309197112358106, error 1.03e-06

Composite functions and the chain rule fall out automatically

g = lambda x: (x.exp() / (x + 1.0)).log()
x1 = 2.0
grad = derivative(g, x1)
print(f"d/dx log(exp(x)/(x+1)) at x={x1}: {grad:.10f}")

d/dx log(exp(x)/(x+1)) at x=2.0: 0.6666666667

A NumPy sanity check against a nearby closed-form derivative

numpy_gradient_estimate = np.gradient([f_plain(x0 - 1e-6), f_plain(x0), f_plain(x0 + 1e-6)], 1e-6)[1]
print("close to numpy finite-difference estimate:", bool(np.isclose(dual_result, numpy_gradient_estimate, atol=1e-3)))

close to numpy finite-difference estimate: True