step.K: Kink-based step selection

In pnd: Parallel Numerical Derivatives, Gradients, Jacobians, and Hessians of Arbitrary Accuracy Order

Kink-based step selection

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),
  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.
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():

• par – the optimal step size found.

• value – the estimated numerical first derivative (using central differences).

• counts – the number of iterations (each iteration includes two function evaluations).

• abs.error – an estimate of the truncation and rounding errors.

• exitcode – an integer code indicating the termination status:

  • 0 – Optimal termination; a minimum of the V-shaped function was found.

  • 1 – Third derivative is too small or noisy; a fail-safe value is returned.

  • 2 – Third derivative is nearly zero; a fail-safe value is returned.

  • 3 – There is no left branch of the V shape; a fail-safe value is returned.

  • 4 – Step trimmed to 0.1|x| when |x| is not tiny and within range.

• message – A summary message of the exit status.

• iterations – A list including the step and argument grids, function values on those grids, estimated derivative values, estimated error values, and predicted model-based errors.

References

\insertAllCited

Examples

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

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

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

pnd documentation built on Sept. 9, 2025, 5:44 p.m.