Composable total derivative operator

D(f) returns the derivative of f as a new function. D(f, x) evaluates the derivative at x. D(D(f)) composes for second-order derivatives, and so on.

D(f, x = NULL, order = 1L)

Arguments

f: A function taking a parameter vector (via x[i] indexing) and returning a scalar, list, or dual_vector.
x: Optional numeric vector. If provided, evaluates D(f)(x).
order: Derivative order (default 1). order = k applies D k times.

Value

If x is NULL, a function. Otherwise, a numeric vector, matrix, or array of the appropriate tensor shape.

Details

Each application of D appends one n-dimension to the output shape, where n = length(x):

For f: R^n -> R: D gives (n) gradient, D^2 gives (n,n) Hessian, D^3 gives (n,n,n), etc.
For f: R^n -> R^m: D gives (m,n) Jacobian, D^2 gives (m,n,n), etc.

The composability works because the S4 dispatch for dualr arithmetic handles nested duals recursively. When D(f) is called with dual inputs (from an outer D), derivative propagation is automatic.

gradient(), hessian(), and jacobian() are convenience wrappers: gradient(f, x) is D(f, x), hessian(f, x) is D(f, x, order = 2), and jacobian(f, x) is D(f, x).

Examples

f <- function(x) x[1]^2 * x[2]
D(f, c(3, 4))                     # gradient: c(24, 9)
#> [1] 24  9
D(f, c(3, 4), order = 2)          # Hessian matrix
#>      [,1] [,2]
#> [1,]    8    6
#> [2,]    6    0

g <- function(x) list(x[1] * x[2], x[1]^2)
D(g, c(2, 3))                     # 2x2 Jacobian
#>      [,1] [,2]
#> [1,]    3    2
#> [2,]    4    0

# Composition: D(D(f)) == D(f, order = 2)
Df <- D(f)
DDf <- D(Df)
DDf(c(3, 4))                      # same as D(f, c(3,4), order=2)
#>      [,1] [,2]
#> [1,]    8    6
#> [2,]    6    0