GenD {pnd} | R Documentation |
Numerical derivative matrices with parallel capabilities
Description
Computes numerical derivatives of a scalar or vector function using finite-difference methods.
This function serves as a backbone for Grad()
and Jacobian()
, allowing for detailed control
over the derivative computation process, including order of derivatives, accuracy, and step size.
GenD
is fully vectorised over different coordinates of the function argument,
allowing arbitrary accuracies, sides, and derivative orders for different coordinates.
Usage
GenD(
FUN,
x,
elementwise = NA,
vectorised = NA,
multivalued = NA,
deriv.order = 1L,
side = 0,
acc.order = 2L,
stencil = NULL,
h = NULL,
zero.tol = sqrt(.Machine$double.eps),
h0 = NULL,
control = list(),
f0 = NULL,
cores = 1,
preschedule = TRUE,
cl = NULL,
func = NULL,
method = NULL,
method.args = list(),
...
)
## S3 method for class 'GenD'
print(
x,
digits = 4,
shave.spaces = TRUE,
begin = "",
sep = " ",
end = "",
...
)
Arguments
FUN |
A function returning a numeric scalar or a vector whose derivatives are to be computed. If the function returns a vector, the output will be a Jacobian. |
x |
Numeric vector or scalar: the point(s) at which the derivative is estimated.
|
elementwise |
Logical: is the domain effectively 1D, i.e. is this a mapping
|
vectorised |
Logical: if |
multivalued |
Logical: if |
deriv.order |
Integer or vector of integers indicating the desired derivative order,
|
side |
Integer scalar or vector indicating the type of finite difference:
|
acc.order |
Integer or vector of integers specifying the desired accuracy order
for each element of |
stencil |
Optional custom vector of points for function evaluation.
Must include at least |
h |
Numeric or character specifying the step size(s) for the numerical
difference or a method of automatic step determination ( |
zero.tol |
Small positive integer: if |
h0 |
Numeric scalar of vector: initial step size for automatic search with
|
control |
A named list of tuning parameters passed to |
f0 |
Optional numeric: if provided, used to determine the vectorisation type to save time. If FUN(x) must be evaluated (e.g. second derivatives), saves one evaluation. |
cores |
Integer specifying the number of CPU cores used for parallel computation. Recommended to be set to the number of physical cores on the machine minus one. |
preschedule |
Logical: if |
cl |
An optional user-supplied |
func |
For compatibility with |
method |
For compatibility with |
method.args |
For compatibility with |
... |
Ignored. |
digits |
Positive integer: the number of digits after the decimal comma to round to (i.e. one less than the number of significant digits). |
shave.spaces |
Logical: if true, removes spaces to ensure compact output; if false, results in nearly fixed-width output (almost). |
begin |
A character to put at the beginning of each line, usually |
sep |
The column delimiter, usually |
end |
A character to put at the end of each line, usually |
Details
For computing gradients and Jacobians, use convenience wrappers Jacobian
and Grad
.
If the step size is too large, the slope of the secant poorly estimates the derivative; if it is too small, it leads to numerical instability due to the function value rounding.
The optimal step size for one-sided differences typically approaches Mach.eps^(1/2)
to balance the Taylor series truncation error with the rounding error due to storing
function values with limited precision. For two-sided differences, it is proportional
to Mach.eps^(1/3). However, selecting the best step size typically requires knowledge
of higher-order derivatives, which may not be readily available. Luckily,
using step = "SW"
invokes a reliable automatic data-driven step-size selection.
Other options include "DV"
, "CR"
, and "CRm"
.
The step size defaults to 1.5e-8
(the square root of machine epsilon)
for x
less than 1.5e-8
,
(2.2e-16)^(1/(deriv.order + acc.order)) * x
for x > 1
, and interpolates
exponentially between these values for 1.5e-8 < x < 1
.
The use of f0
can reduce computation time similar to the use of f.lower
and f.upper
in uniroot()
.
Unlike numDeriv::grad()
and numDeriv::jacobian()
, this function
fully preserves the names of x
and FUN(x)
.
Value
A vector or matrix containing the computed derivatives, structured according
to the dimensionality of x
and FUN
. If FUN
is scalar-valued,
returns a gradient vector. If FUN
is vector-valued, returns a Jacobian matrix.
See Also
gradstep()
for automatic step-size selection.
Examples
# Case 1: Vector argument, vector output
f1 <- sin
g1 <- GenD(FUN = f1, x = 1:100)
g1.true <- cos(1:100)
plot(1:100, g1 - g1.true, main = "Approximation error of d/dx sin(x)")
# Case 2: Vector argument, scalar result
f2 <- function(x) sum(sin(x))
g2 <- GenD(FUN = f2, x = 1:4)
g2.h2 <- Grad(FUN = f2, x = 1:4, h = 7e-6)
g2 - g2.h2 # Tiny differences due to different step sizes
g2.auto <- Grad(FUN = f2, x = 1:4, h = "SW")
print(attr(g2.auto, "step.search")$exitcode) # Success
# Case 3: vector input, vector argument of different length
f3 <- function(x) c(sum(sin(x)), prod(cos(x)))
x3 <- 1:3
j3 <- GenD(f3, x3, multivalued = TRUE)
print(j3)
# Gradients for vectorised functions -- e.g. leaky ReLU
LReLU <- function(x) ifelse(x > 0, x, 0.01*x)
system.time(replicate(10, suppressMessages(GenD(LReLU, runif(30, -1, 1)))))
system.time(replicate(10, suppressMessages(GenD(LReLU, runif(30, -1, 1)))))
# Saving time for slow functions by using pre-computed values
x <- 1:4
finner <- function(x) sin(x*log(abs(x)+1))
fouter <- function(x) integrate(finner, 0, x, rel.tol = 1e-12, abs.tol = 0)$value
# The outer function is non-vectorised
fslow <- function(x) {Sys.sleep(0.01); fouter(x)}
f0 <- sapply(x, fouter)
system.time(GenD(fslow, x, side = 1, acc.order = 2, f0 = f0))
# Now, with extra checks, it will be slower
system.time(GenD(fslow, x, side = 1, acc.order = 2))
# Printing whilst preserving names
x <- structure(1:3, names = c("Andy", "Bradley", "Ca"))
print(Grad(function(x) prod(sin(x)), 1)) # 1D derivative
print(Grad(function(x) prod(sin(x)), x))
print(Jacobian(function(x) c(prod(sin(x)), sum(exp(x))), x))