README.md
In nabla: Exact Derivatives via Automatic Differentiation

nabla

Arbitrary-order exact derivatives at machine precision

nabla provides a single composable operator D that differentiates any R function to any order — exactly, at machine precision, through loops, branches, and all control flow:

library(nabla)
f <- function(x) x[1]^2 * exp(x[2])

D(f, c(1, 0))                # gradient
D(f, c(1, 0), order = 2)     # Hessian
D(f, c(1, 0), order = 3)     # 2×2×2 third-order tensor
D(f, c(1, 0), order = 4)     # 2×2×2×2 fourth-order tensor

Each application of D adds one dimension to the output. D(D(f)) gives the Hessian, D(D(D(f))) gives the third-order tensor, and so on — no limit on order, no loss of precision, no symbolic algebra.

| | Finite Differences | Symbolic Diff | AD (nabla) | |---|---|---|---| | Accuracy | O(h) or O(h²) truncation error | Exact | Exact (machine precision) | | Higher-order | Error compounds rapidly | Expression swell | Composes cleanly to any order | | Control flow | Works | Breaks on if/for/while | Works through any code |

Finite differences lose precision at higher orders (each order multiplies the error). Symbolic differentiation suffers from expression swell. nabla composes D via nested dual numbers — each order is as precise as the first.

# Install from CRAN
install.packages("nabla")

# Or install development version from GitHub
remotes::install_github("queelius/nabla")

The `D` operator

D is the core of nabla. It differentiates any function f and returns a new function — which can itself be differentiated:

f <- function(x) x[1]^2 * x[2] + sin(x[2])

Df   <- D(f)          # first derivative (function)
DDf  <- D(Df)         # second derivative (function)
DDDf <- D(DDf)        # third derivative (function)

Df(c(3, 4))           # gradient vector
DDf(c(3, 4))          # Hessian matrix
DDDf(c(3, 4))         # 2×2×2 tensor

Equivalently, evaluate directly at a point:

D(f, c(3, 4))              # gradient
D(f, c(3, 4), order = 2)   # Hessian
D(f, c(3, 4), order = 3)   # third-order tensor

gradient(), hessian(), and jacobian() are convenience wrappers:

gradient(f, c(3, 4))       # == D(f, c(3, 4))
hessian(f, c(3, 4))        # == D(f, c(3, 4), order = 2)

A dual number extends the reals with an infinitesimal ε where ε² = 0:

$$f(x + \varepsilon) = f(x) + f'(x)\,\varepsilon$$

For higher orders, nabla nests dual numbers: a dual whose components are themselves duals. Each level of nesting extracts one additional order of derivative — so D(D(D(f))) propagates through triply-nested duals to produce exact third derivatives. This works through lgamma, psigamma, trig functions, and all of R's math — no special cases needed.

Optimization — supply exact gradients to optim() and nlminb()
Maximum likelihood — Hessians for standard errors, third-order tensors for asymptotic skewness of MLEs
Sensitivity analysis — how outputs change with respect to inputs
Taylor approximation — exact coefficients to any order
Curvature analysis — second-order geometric properties

Introduction to nabla — dual numbers, arithmetic, composition, comparison with finite differences
MLE Workflow — gradient, Hessian, Newton-Raphson on statistical models
Higher-Order Derivatives — the D operator, curvature, Taylor expansion
Higher-Order MLE Analysis — third-order derivative tensors and asymptotic skewness of MLEs
Optimizer Integration — using gradient() and hessian() with optim() and nlminb()