Tutorials › Pharmacometric Workflows and Reproducible Modeling › Finite-Difference Step-Size Selection

Initial Estimates and Estimation Settings · Tutorial 113 of 399

Finite-Difference Step-Size Selection

Learn how finite-difference accuracy depends on step size and how to select and check perturbations that respect each parameter’s scale.

Advanced 12 min read

What You'll Learn

  • Distinguish truncation error from roundoff and objective-evaluation noise.
  • Derive the forward- and central-difference step-size trade-offs.
  • Choose perturbations in dimensionless, parameter-scaled coordinates.
  • Work through a two-parameter objective with widely differing parameter scales.
  • Diagnose step sizes using error trends, repeatability, and derivative checks.

Why the Step Size Matters

A finite-difference gradient is an approximation: it estimates how an objective function changes when a parameter is perturbed. The perturbation must be large enough for that change to stand out from numerical error, but small enough that the function is locally close to its linear approximation. Those requirements compete. A step that is too large causes truncation error; a step that is too small magnifies roundoff error and any noise in the objective evaluation.

In estimation, a gradient can guide parameter updates, so an inaccurate component may distort the direction or size of an update. The earlier material on analytic and finite-difference gradients distinguishes ways to obtain gradients. Here the focus is how to choose the finite-difference step once a numerical approximation is being used. Step size is a numerical accuracy choice, not a model parameter.

Difference Formulas and Their Error

Let \(f(\theta)\) be a scalar objective, and suppose it is evaluated near the point \(\theta\). A forward difference uses a positive step \(h\); a central difference uses equal steps on either side. The derivative is assumed to exist in a neighborhood of \(\theta\), and the evaluation points must be available.

Definition: The forward-difference estimate is \(D_F(h)=[f(\theta+h)-f(\theta)]/h\). The central-difference estimate is \(D_C(h)=[f(\theta+h)-f(\theta-h)]/(2h)\). The truncation error is the difference between the exact derivative and the estimate when function values are exact. Roundoff error arises from finite-precision arithmetic; objective-evaluation noise can add another error source.

Taylor’s theorem explains the different truncation orders. If \(f''\) is continuous near \(\theta\), the first-order expansion gives \(f(\theta+h)=f(\theta)+f'(\theta)h+f''(\theta)h^2/2+o(h^2)\). Dividing the difference by \(h\) yields \(D_F(h)=f'(\theta)+f''(\theta)h/2+o(h)\). If the third derivative is bounded near \(\theta\), the expansion through the quadratic term has remainder \(O(h^3)\); consequently the forward difference has an \(O(h^2)\) remainder after its leading \(f''(\theta)h/2\) term. A bounded third derivative does not, by itself, justify a Taylor remainder of \(O(h^4)\) after a cubic expansion.

For the central formula, if \(f'''\) is bounded near \(\theta\), Taylor expansion through the quadratic term at each side has remainder \(O(h^3)\). Subtracting the two expansions cancels the constant and quadratic terms; after dividing by \(2h\), the central-difference truncation error is \(O(h^2)\). If \(f'''\) is continuous, the more specific expansion is \(D_C(h)=f'(\theta)+f'''(\theta)h^2/6+o(h^2)\). The cancellation of even powers is why the central formula typically has smaller truncation error than the forward formula at the same sufficiently small step.

Neither truncation result says that decreasing \(h\) indefinitely improves the computed derivative. If each computed function value has an absolute error no greater than \(\delta\), the error introduced into the forward quotient is bounded by \(2\delta/h\), and that in the central quotient by \(\delta/h\). These bounds follow by applying the triangle inequality to the errors in the function values and dividing by the denominator. In floating-point calculations, \(\delta\) is often roughly proportional to machine precision times the magnitude of the function evaluation, but the actual error also depends on the implementation. If the objective is itself computed approximately, its numerical integration or stochastic evaluation error can exceed ordinary floating-point roundoff.

Deriving the Truncation–Roundoff Balance

A simple local error model makes the compromise explicit. For a forward difference, suppose the magnitude of the leading truncation error is approximately \(A h\), with \(A>0\), while roundoff contributes approximately \(B\delta/h\), with \(B>0\). The total error is modeled by \(E_F(h)=Ah+B\delta/h\). Differentiating with respect to \(h\) and setting the derivative to zero gives \(A-B\delta/h^2=0\), so the minimizing step in this model is \(h_F=\sqrt{B\delta/A}\).

For a central difference, suppose the corresponding truncation term is \(C h^2\), with \(C>0\), and retain the same inverse-step roundoff term. Then \(E_C(h)=C h^2+B\delta/h\), and its derivative is \(2Ch-B\delta/h^2\). Setting this equal to zero gives \(h_C=(B\delta/(2C))^{1/3}\). These are approximate balance rules, not universal constants: \(A\), \(C\), and \(\delta\) depend on the function, the evaluation procedure, and the point.

Worked Example: Balancing Errors in a Scaled Coordinate

Suppose a parameter has first been expressed in a dimensionless coordinate \(z\). For a hypothetical forward-difference calculation, let the local error model be \(E_F(h)=4h+\epsilon/h\), where \(\epsilon=2.22\times10^{-16}\) is the machine epsilon used in this illustration. The predicted minimizing step is

$$ h_F=\sqrt{\frac{\epsilon}{4}} =\sqrt{5.55\times10^{-17}} \approx 7.45\times10^{-9}. $$

For a hypothetical central-difference calculation with \(E_C(h)=2h^2+\epsilon/h\), the model instead gives

$$ h_C=\left(\frac{\epsilon}{4}\right)^{1/3} =(5.55\times10^{-17})^{1/3} \approx 3.81\times10^{-6}. $$

These values are not prescriptions for every objective. They demonstrate that the optimal scale depends on the truncation coefficient as well as the precision, and that the central and forward formulas balance different powers of \(h\). In practice, local coefficients and evaluation noise are rarely known precisely, so the model motivates a step-size sweep rather than replacing one.

Choose Steps in Each Parameter’s Scale

A single numerical step for every parameter is usually inappropriate when parameter magnitudes or units differ. For a vector of parameters \(\theta=(\theta_1,\ldots,\theta_p)\), assign each parameter a positive characteristic scale \(s_j\) in the same units as \(\theta_j\). Define the dimensionless coordinate \(z_j=\theta_j/s_j\). A common relative perturbation \(\eta\) then corresponds to \(h_j=\eta s_j\) in the original parameter units. The partial derivative with respect to \(\theta_j\) is approximated by perturbing only that coordinate and dividing by \(h_j\); in the dimensionless coordinate, it is the derivative with respect to \(z_j\), divided by \(s_j\).

A useful \(s_j\) reflects the expected magnitude or a meaningful change in parameter \(j\), and must be chosen in that parameter’s units. It might be based on a typical estimate, a scientifically meaningful increment, or a scale used by the estimation software. A rule such as \(s_j=\max(|\theta_j|,1)\) is not automatically meaningful: the number 1 has different implications for a parameter measured in litres, hours, or a dimensionless ratio. When a parameter value is near zero, use a meaningful nonzero scale rather than letting the perturbation collapse with the current value.

Worked Example: Different Scales in a Two-Parameter Objective

Consider the hypothetical objective \(f(x,y)=\exp(x/1000)+0.5\exp(y/0.01)\), evaluated at \(x=0\), \(y=0\). The parameters have different units and characteristic scales: take \(s_x=1000\) units of \(x\) and \(s_y=0.01\) units of \(y\). The exact derivatives at this point are \(\partial f/\partial x=1/1000=0.001\) per unit of \(x\), and \(\partial f/\partial y=0.5/0.01=50\) per unit of \(y\).

Choose dimensionless central step \(\eta=0.001\). Thus \(h_x=\eta s_x=1\) unit of \(x\), and \(h_y=\eta s_y=0.00001\) unit of \(y\). Perturb one coordinate at a time, holding the other fixed. For \(x\), the central estimate is

$$ D_{C,x} =\frac{\exp(1/1000)-\exp(-1/1000)}{2(1)} =\frac{\sinh(0.001)}{1000(0.001)} \approx 0.001000000167 $$

per unit of \(x\). For \(y\), the factor \(0.5\) in the objective must be retained:

$$ D_{C,y} =\frac{0.5\exp(0.001)-0.5\exp(-0.001)}{2(0.00001)} =50\frac{\sinh(0.001)}{0.001} \approx 50.00000833 $$

per unit of \(y\). In both coordinates the relative truncation error is approximately \(\sinh(0.001)/0.001-1\approx1.667\times10^{-7}\). The same dimensionless step has produced very different perturbations in the original units while keeping the relative perturbation comparable. The calculation also shows that the two gradient components have different units; comparing their raw magnitudes alone is not a check of their numerical accuracy.

Step-Size Diagnostics and Failure Modes

Because the ideal step is not generally known in advance, evaluate each gradient component at several nearby step sizes. A practical diagnostic is whether estimates remain reasonably stable over a range of steps before smaller steps begin to fluctuate. If the step is too large, truncation error often appears as a systematic change in the estimate as \(h\) decreases. If it is too small, cancellation in the subtraction or noise in the objective can cause erratic changes. For central differences in a smooth region, a truncation-dominated error often decreases approximately fourfold when \(h\) is halved, since its leading term is proportional to \(h^2\). Once roundoff or evaluation noise dominates, that pattern breaks down.

Worked Example: Reading a Step-Size Sweep

Let the hypothetical scalar objective be \(g(z)=\exp(z)\), evaluated at \(z=0\), where the exact derivative is 1. The forward estimate is \((\exp(h)-1)/h\), while the central estimate is \((\exp(h)-\exp(-h))/(2h)=\sinh(h)/h\). At \(h=0.1\), the estimates are approximately 1.051709 for forward and 1.001668 for central. At \(h=0.01\), they are approximately 1.005017 and 1.00001667. At \(h=0.00001\), they are approximately 1.000005 and 1.000000000017, respectively, before accounting for roundoff.

The forward errors decrease roughly in proportion to \(h\), while the central truncation errors decrease roughly in proportion to \(h^2\). But at the smallest step, a function-value error on the order of \(10^{-16}\), divided by a step of \(10^{-5}\), can contribute quotient error on the order of \(10^{-11}\). The central truncation error there is only about \(1.67\times10^{-11}\), so roundoff can already be comparable. The displayed exact-arithmetic estimates therefore do not guarantee that a computer will return the most accurate derivative at the smallest step. The appropriate conclusion is to inspect neighboring steps and the actual computed values, not to declare the smallest step best.

Other useful checks include repeating objective evaluations, comparing forward and central estimates when evaluations on both sides are available, and comparing with an analytic derivative or an independent numerical implementation when one exists. For stochastic objectives, a gradient may vary because of sampling noise rather than finite precision; a smaller step can amplify that variability. If evaluation noise is material, assess the repeatability of the function values and gradient estimates, and use an evaluation procedure whose precision is adequate for the intended step.

A one-sided formula may be necessary when the objective cannot be evaluated on both sides of a point. It has a different truncation behavior from the central formula, so central-difference step guidance should not be transferred to it without checking the formula’s error order. This tutorial does not prescribe how parameter constraints are handled; the essential requirement here is that every perturbed evaluation used by the selected formula is valid.

Choosing Among Practical Approaches

ApproachStrengthLimitation
Fixed absolute stepSimple and reproducible when parameters share units and scales.Can be too large for one parameter and too small for another.
Scaled step \(h_j=\eta s_j\)Adapts perturbations to parameter units and characteristic magnitudes.Requires defensible scales \(s_j\); a poor scale can still mislead.
Step-size sweepReveals truncation trends, instability, and a region of estimate stability.Costs additional objective evaluations and does not remove objective noise.
Analytic derivativeAvoids finite-difference step selection and subtraction error.Requires a correct derivative implementation and may not be available for every objective.

A reasonable workflow combines these approaches: use parameter-specific scales, choose a moderate dimensionless starting step, and test a small range around it. Prefer a stable range of derivative estimates over a single apparently precise value. Record the step convention and scales so that the gradient calculation is reproducible.

Common Mistakes and Limitations

  • Making every step smaller. This reduces truncation error only while truncation dominates. Once roundoff or evaluation noise is important, further reduction can worsen the derivative.
  • Using one absolute step for differently scaled parameters. A step appropriate for a parameter near 1000 may be enormous for a parameter near 0.01. Scale the steps in parameter-specific units.
  • Assuming machine precision is the only noise source. Approximate integrations, iterative calculations, and stochastic likelihood evaluations can introduce errors larger than floating-point roundoff.
  • Comparing gradient components without units or scaling. Derivatives with respect to parameters measured in different units have different units and raw magnitudes. Check their accuracy relative to a meaningful parameter scale.
  • Treating the balance formula as exact. The local coefficients and effective evaluation error are not usually known precisely. Use the formulas to understand the trade-off, then inspect a step-size sweep.

Finite-difference step selection is limited by local smoothness and objective reliability. If the objective changes irregularly over the perturbation scale, the smooth Taylor error model may not describe the observed behavior. In that case, a stable and reproducible derivative may require improving objective evaluation or using an alternative derivative calculation, rather than further tuning \(h\) alone.

Check Your Understanding

Answer each question, then compare with the solutions below.

  1. For a forward difference with error model \(E(h)=3h+12\delta/h\), derive the minimizing step.
  2. A parameter is measured in milligrams and has characteristic scale \(s=40\) mg. What perturbation corresponds to dimensionless step \(\eta=0.002\)?
  3. Why can a central-difference estimate become less accurate when its step is made very small?
  4. For a smooth central difference, what approximate change in truncation error is expected when the step is halved?

Solutions. For question 1, differentiate: \(E'(h)=3-12\delta/h^2\). Setting this to zero gives \(h^2=4\delta\), so \(h=2\sqrt{\delta}\). For question 2, \(h=\eta s=0.002(40\text{ mg})=0.08\text{ mg}\). For question 3, the subtraction of nearly equal objective values magnifies function-value roundoff or evaluation noise after division by \(h\); the error can grow approximately as \(\delta/h\). For question 4, the leading central truncation error is proportional to \(h^2\), so halving \(h\) reduces that component by approximately a factor of four, provided truncation error still dominates.

Clinical Biostats