step.CR: Curtis-Reid automatic step selection

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

Curtis–Reid automatic step selection

Description

Curtis–Reid automatic step selection

Usage

step.CR(
  FUN,
  x,
  deriv.order = 1,
  acc.order = 1,
  h0 = NULL,
  max.rel.error = .Machine$double.eps^(7/8),
  aim = NULL,
  tol = NULL,
  range = NULL,
  maxit = 20L,
  seq.tol = 1e-04,
  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.
deriv.order	Order of the derivative (m in \frac{d^m f}{dx^m}) for which a numerical approximation is needed.
acc.order	Order of accuracy: defines how the approximation error scales with the step size h, specifically O(h^{a+1}), where a is the accuracy order and depends on the higher-order derivatives of the function.
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).
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.
aim	Positive real scalar: desired ratio of truncation-to-rounding error. The original version over-estimates the truncation error, hence a higher aim is recommended. For the modernised version, aim should be equal to deriv.order / acc.order.
tol	Numeric scalar greater than 1: tolerance multiplier for determining when to stop the algorithm based on the current estimate being between aim/tol and aim*tol.
range	Numeric vector of length 2 defining the valid search range for the step size.
maxit	Integer: maximum number of algorithm iterations to prevent infinite loops in degenerate cases.
seq.tol	Numeric scalar: maximum relative difference between old and new step sizes for declaring convergence.
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 \insertCitecurtis1974choicepnd algorithm. If the estimated third derivative is exactly zero, then, the initial step size is multiplied by 4 and returned.

If 4th-order accuracy is requested, then, two things happen. Firstly, since 4th-order-accurate differences requires a larger step size and the truncation error for the 2nd-order-accurate differences grows if the step size is larger than the optimal one, a higher ratio of truncation-to-rounding errors should be targeted. Secondly, a 4th-order-accurate numerical derivative is returned, but the truncation and rounding errors are still estimated for the 2nd-order-accurate differences. Therefore, the estimated truncation error is higher and the real truncation error of 4OA differences is lower.

TODO: mention that f must be one-dimensional

The arguments passed to ... must not partially match those of step.CR(). For example, if cl exists, then, attempting to avoid cluster export by using step.CR(f, x, h = 1e-4, cl = cl, a = a) will result in an error: a matches aim and acc.order. Redefine the function for this argument to have a name that is not equal to the beginning of one of the arguments of step.CR().

The original version of the article uses deriv.order = 1, acc.order = 1, aim = 100.

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; especially useful for computationally expensive functions).

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

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

• exitcode – an integer code indicating the termination status:

  • 0 – Optimal termination within tolerance.

  • 1 – Third derivative is zero; large step size preferred.

  • 2 – No change in step size within tolerance.

  • 3 – Solution lies at the boundary of the allowed value range.

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

  • 5 – Maximum number of iterations reached.

• message – A summary message of the exit status.

• iterations – A list including the full step size search path, argument grids, function values on those grids, estimated error ratios, and estimated derivative values.

References

\insertAllCited

Examples

f <- function(x) x^4
step.CR(x = 2, f)  # Wrap plot(...) around each of these functions
step.CR(x = 2, f, h0 = 1e-3)
step.CR(x = 2, f, acc.order = 2)
step.CR(x = 2, f, deriv.order = 2)
step.CR(x = 2, f, deriv.order = 3, acc.order = 4)

# A bad start: too far away
step.CR(x = 2, f, h0 = 1000)  # Bad exit code + a suggestion to extend the range
step.CR(x = 2, f, h0 = 1000, range = c(1e-10, 1e5))  # Problem solved

library(parallel)
cl <- makePSOCKcluster(names = 2, outfile = "")
abc <- 2
f <- function(x, abc) {Sys.sleep(0.02); abc*sin(x)}
x <- pi/4
system.time(step.CR(f, x, h = 1e-4, cores = 1, abc = abc))  # To remove speed-ups
system.time(step.CR(f, x, h = 1e-4, cores = 2, abc = abc))  # Faster
f2 <- function(x) f(x, abc)
clusterExport(cl, c("f2", "f", "abc"))
system.time(step.CR(f2, x, h = 1e-4, cl = cl))  # Also fast
stopCluster(cl)

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