fastdual

Fast forward-mode automatic differentiation via dual numbers, implemented as a CPython C extension.

Computes exact gradients, Jacobians, and Hessians with minimal overhead — no taping, no graph construction, just numbers that carry their derivatives.

Install

pip install fastdual

Drop-in Gradients

Any function that works with floats works with Dual — no rewriting, no framework, no JIT warmup:

from fastdual import Dual, der

def my_function(x):
    return x**3 - 2*x + 1

x = Dual(3.0)
y = my_function(x)
dy_dx = der(y, x)  # 25.0 (exact derivative)

Quick Start

from fastdual import Dual, der
import numpy as np

x = Dual(3.0)
y = Dual(5.0)

z = x * y + np.sin(x)
print(z.val)        # 15.1411...
print(der(z, x))    # 5.99 (dz/dx = y + cos(x))
print(der(z, y))    # 3.0  (dz/dy = x)

Array Operations

from fastdual import DualArray, val, jac
import numpy as np

xs = DualArray([1.0, 2.0, 3.0])  # independent seeds

result = np.sin(xs) + xs ** 2
print(val(result))     # [sin(1)+1, sin(2)+4, sin(3)+9]
print(jac(result, xs)) # diagonal Jacobian

Automatic Jacobians

from fastdual import autojac
import numpy as np

@autojac
def f(x, y):
    return np.array([x**2 + y, x * y**2])

result, J = f(2.0, 3.0)
# result = [7.0, 18.0]
# J = [[4.0,  1.0],
#      [9.0, 12.0]]

Automatic Hessians

Second-order derivatives via hyper-dual numbers (also a C extension):

from fastdual import autohess

@autohess
def rosenbrock(x, y):
    return (1.0 - x)**2 + 100.0 * (y - x**2)**2

result, H = rosenbrock(1.0, 1.0)
# result = 0.0
# H = [[802, -400],
#      [-400, 200]]

Optimization

Automatic gradient computation for scipy.optimize.minimize:

from fastdual import minimize

def objective(x):
    return (x[0] - 1)**2 + 100 * (x[1] - x[0]**2)**2

result = minimize(objective, [0.0, 0.0])
print(result.x)  # [1.0, 1.0]

Pass hess=True for exact Hessians via hyper-dual numbers, enabling second-order methods (trust-ncg, Newton-CG, etc.) that converge in fewer iterations:

result = minimize(objective, [0.0, 0.0], hess=True)

Requires pip install fastdual[optimize].

Sparse Jacobians

For large systems with known sparsity, avoid redundant computation via graph coloring:

from fastdual import sparse_jac
import numpy as np

def f(x):
    return np.array([x[i-1] + x[i] + x[i+1] for i in range(1, len(x)-1)])

sparsity = ...  # boolean (m, n) array of known nonzero entries
J = sparse_jac(f, x0, sparsity)  # only n_colors << n forward passes

NumPy Integration

DualArray supports __array_ufunc__ and __array_function__ protocols:

from fastdual import DualArray
import numpy as np

x = DualArray([1.0, 2.0, 3.0])

np.sin(x)              # ufuncs
np.dot(x, x)           # dot product
np.sum(x)              # reduction
np.linalg.norm(x)      # L2 norm
np.linalg.solve(A, b)  # linear solve with gradient propagation

Supported Operations

Both Dual (first-order) and HyperDual (second-order) types are C extensions. All operations work as methods and (for Dual) as NumPy ufuncs.

Arithmetic

Transcendental Functions

Available as methods (.sin()) on both types and via NumPy ufuncs (np.sin()) on Dual/DualArray.

Binary Functions (Dual only)

These are available via NumPy ufuncs on Dual/DualArray.

Utility Functions

Predicates (Dual only)

np.isfinite(), np.isinf(), np.isnan() — check the primal value, return bool.

Comparisons

<, <=, ==, !=, >, >= — compare on primal value only.

API Reference

Performance

All hot paths are in C — both Dual and HyperDual types are C extensions with zero Python object allocation in the inner loop.

Dual: overhead vs plain floats

HyperDual: overhead vs plain floats

HyperDual carries 4 fixed doubles — no sparse gradient bookkeeping. Per-element arithmetic is nearly free compared to floats.

Jacobian: fastdual vs finite differences

Jacobians use the C extension for forward-mode AD — one pass computes all partials simultaneously, vs n+1 function evaluations for finite differences.

Hessian: fastdual vs finite differences

Hessians require n(n+1)/2 function evaluations (each with HyperDual arithmetic). For small n, finite differences with simple functions can be competitive. The hyper-dual approach shines when derivatives must be exact (no step-size tuning) or when the function involves transcendentals where finite-difference errors grow.

Gradient vs Hessian

How much more does a Hessian cost compared to a gradient for the same function?

Dual computes the full gradient in a single forward pass but carries a sparse gradient vector that grows with the number of variables. HyperDual uses 4 fixed doubles per element (no per-variable scaling), but needs n(n+1)/2 passes for the full Hessian. The ratio reflects this: Hessians are roughly O(n²) more expensive than gradients.

Test

pytest tests/ -v