step.K {pnd}R Documentation

Kink-based step selection (experimental!)

Description

Optimal step-size search using the full range of practical error estimates and numerical optimisation to find the spot where the theoretical total-error shape is best described by the data, and finds the step size where the ratio of rounding-to-truncation error is optimal.

Usage

step.K(
  FUN,
  x,
  h0 = NULL,
  deriv.order = 1,
  acc.order = 2,
  range = NULL,
  shrink.factor = 0.5,
  max.rel.error = .Machine$double.eps^(7/8),
  plot = FALSE,
  cores = 1,
  preschedule = getOption("pnd.preschedule", TRUE),
  cl = NULL,
  ...
)

Arguments

FUN

Function for which the optimal numerical derivative step size is needed.

x

Numeric scalar: the point at which the derivative is computed and the optimal step size is estimated.

h0

Numeric scalar: initial step size, defaulting to a relative step of slightly greater than .Machine$double.eps^(1/3) (or absolute step if x == 0).

deriv.order

Integer or vector of integers indicating the desired derivative order, \mathrm{d}^m / \mathrm{d}x^m, for each element of x.

acc.order

Integer or vector of integers specifying the desired accuracy order for each element of x. The final error will be of the order O(h^{\mathrm{acc.order}}).

range

Numeric vector of length 2 defining the valid search range for the step size.

shrink.factor

A scalar less than 1 that is used to create a sequence of step sizes. The recommended value is 0.5. Change to 0.25 for a faster search. This number should be a negative power of 2 for the most accurate representation.

max.rel.error

Error bound for the relative function-evaluation error (\frac{\hat f(\hat x) - f(x)}{f(x)}). Measures how noisy a function is. If the function is relying on numerical optimisation routines, consider setting to sqrt(.Machine$double.eps). If the function has full precision to the last bit, set to .Machine$double.eps/2.

plot

Logical: if TRUE, plots the estimated truncation and round-off errors.

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 TRUE, disables pre-scheduling for mclapply() or enables load balancing with parLapplyLB(). Recommended for functions that take less than 0.1 s per evaluation.

cl

An optional user-supplied cluster object (created by makeCluster or similar functions). If not NULL, the code uses parLapply() (if preschedule is TRUE) or parLapplyLB() on that cluster on Windows, and mclapply (fork cluster) on everything else.

...

Passed to FUN.

Details

This function computes the optimal step size for central differences using the statistical kink-search approach. The optimal step size is determined as the minimiser of the total error, which for central differences is V-shaped with the left-branch slope equal to the negative derivation order and the right-branch slope equal to the accuracy order. For standard simple central differences, the slopes are -1 and 2, respectively. The algorithm uses the least-median-of-squares (LMS) penalty and searches for the optimal position of the check that fits the data the best on a bounded 2D rectangle using derivative-free (Nelder–Mead) optimisation. Unlike other algorithms, if the estimated third derivative is too small, the function shape will be different, and two checks are made for the existence of two branches.

Value

A list similar to the one returned by optim():

References

There are no references for Rd macro ⁠\insertAllCites⁠ on this help page.

Examples

step.K(sin, 1, plot = TRUE)
step.K(exp, 1, range = c(1e-12, 1e-0), plot = TRUE)
step.K(atan, 1, plot = TRUE)

# Edge case 1: function symmetric around x0, zero truncation error
step.K(sin, pi/2, plot = TRUE)
step.K(sin, pi/2, shrink.factor = 0.8, plot = TRUE)

# Edge case 1: the truncation error is always zero and f(x0) = 0
suppressWarnings(step.K(function(x) x^2, 0, plot = TRUE))
# Edge case 2: the truncation error is always zero
step.K(function(x) x^2, 1, plot = TRUE)
step.K(function(x) x^4, 0, plot = TRUE)
step.K(function(x) x^4, 0.1, plot = TRUE)
step.K(function(x) x^6 - x^4, 0.1, plot = TRUE)
step.K(atan, 3/4, plot = TRUE)
step.K(exp, 2, plot = TRUE, range = c(1e-16, 1e+1))

[Package pnd version 0.1.0 Index]