2
$\begingroup$

Why is $$h = \epsilon \cdot(1 + |x|)$$ a good step size in approximating numerically the derivative of a function $f:\mathbb{R} \to \mathbb{R}$ with $$\frac{f(x+h)-f(x)}{h}?$$

$\endgroup$
5
  • $\begingroup$ Well, you want $h$ to be small but you don't want it to be too small. Machine epsilon is the measure of too small in computer arithmetic. The $(1+|x|)$ is included because machine numbers spread out as their absolute value increases. $\endgroup$ – Mark McClure Apr 1 '17 at 16:22
  • $\begingroup$ @MarkMcClure Do you mind adding some details to explain where $(1+|x|)$ comes from? $\endgroup$ – user411609 Apr 1 '17 at 16:27
  • $\begingroup$ A numerical example: if we for example take $f(x) = \log(x)$ then $f'(x) = 1/x$ so we should have $f'(0.1) = 10$ and $f'(10) = 0.1$. Computing this on my system (in a simple c code) gives me $9.09$ and $0.0$ with this $h$ (I'm assuming as above that $\epsilon$ is machine epsilon). $\endgroup$ – Winther Apr 3 '17 at 0:13
  • $\begingroup$ Whereas if we let $h = \sqrt{\epsilon}$ in double precision, we have $\log'(10) = 0.0999999860$, recovering about half the correct digits possible for double precision. $\endgroup$ – user14717 Apr 3 '17 at 0:35
  • $\begingroup$ Numerical Recipes pages 229-231 has a good discussion on this. In particular, they show that the optimal $h$ to minimize $e_r + e_t$ is $\sqrt \epsilon\, x_c$ for one-sided, $\epsilon^{1/3} \, x_c$ for central differences, where "one often assumes $x_c = x$" . $\endgroup$ – denis Jun 9 '18 at 16:24
3
$\begingroup$

The premise of the question is incorrect: The choice of $h = \epsilon(1+|x|)$ is much too small.

There are two sources of error in approximating a derivative. The first comes from calculating it as a difference quotient, and it is easily seen that the error from this is \begin{align*} \left| \frac{f(x+h) - f(x)}{h} - f'(x) \right| \le \frac{h}{2}|f''(x)| \end{align*} However, we cannot calculate $f(x)$ and $f(x+h)$, we calculate some representable numbers close to them (call them $\hat{f}(x+h)$ and $\hat{f}(x)$). For floating point arithmetic with a guard digit we know that $\hat{f}(x+h) = f(x+h)(1+\epsilon_{1})$ and $\hat{f}(x) = f(x)(1+\epsilon_{2})$, where $|\epsilon_1|, |\epsilon_2| \le \epsilon$, where $\epsilon$ is the unit roundoff.1 So the total error is \begin{align*} \left| \frac{\hat{f}(x+h) - \hat{f}(x)}{h} - f'(x) \right| &= \left| \frac{f(x+h)(1+\epsilon_1) - f(x)(1+\epsilon_2)}{h} - f'(x) \right| \\ &\le\frac{h}{2}|f''(x)| + \frac{1}{h}(|f(x+h)|+|f(x)|)\epsilon =: g(h) \end{align*} We can see that taking $h$ as the smallest representable is likely to lead to castatrophe as the error from the roundoff of the function evaluation becomes unbounded. Instead we want to minimize the maximum error $g(h)$. You can differentiate $g$ and set it to zero to discover that the choice of $h$ that minimizes the maximum error is \begin{align*} \sqrt{\frac{2(|f(x+h)| + |f(x)|)\epsilon}{|f''(x)|}} \end{align*} For this choice of $h$, $g(h) = \sqrt{\epsilon}(|f''(x)| + (|f(x+h)| + |f(x)|)/2)$. Since we do not know $f''(x)$ (and cannot compute it for less than the cost of using a higher order method) then we can choose \begin{align*} h = 2\sqrt{\epsilon} \end{align*} This choice preserves the $O(\sqrt{\epsilon})$ error estimate irrespective of the ratio $(|f(x+h)| + |f(x)|)/|f''(x)|$.

There is a pathological case: If $x$ is so large that $x + 2\sqrt{\epsilon} == x$, then our choice of $h$ calculates the derivative to be precisely zero. Then the correct choice of $h$ is h = float_next(x) - x, which is roughly $|x|\epsilon$. However, in this regime $|x|\epsilon \gg 2\sqrt{\epsilon}$. The error is then dominated by the fact that the distance between floats increases as $|x|$ increases.

All the previous analysis depends on $x+h$ being representable. If $x+h$ is not representable, then you can show that the maximum error is \begin{align*} \frac{h}{2}|f''(x)| + \frac{1}{h}(|f(x+h)|+|f(x)| + |xf'(x)|)\epsilon \end{align*} but this does not change the conclusion that we should choose $h \sim O(\sqrt{\epsilon})$.

As to choosing $h = \epsilon(1+|x|)$, we can see that the error of this choice is \begin{align*} \epsilon(1+|x|)|f''(x)|/2 + \frac{|f(x+h)| + |f(x)|}{1+|x|} \approx \frac{2|f(x)|}{1+|x|} \end{align*} So the error is enormous.

1 Note that this assumes that the algorithm to compute our function produces the best representable estimate for the given precision. For C standard library functions there are sometimes guarantees that this is the case but most special function libraries do not compute to this accuracy as the compute time to extract the last few digits of precision represent the bulk of the computation time. Hence we should assume that $\hat{f}(x) \approx f(x)(1+k\epsilon)$ for some $k \in \mathbb{N}$.

$\endgroup$
3
  • $\begingroup$ And where does $(1+|x|)$ come from? $\endgroup$ – user411609 Apr 1 '17 at 16:24
  • 1
    $\begingroup$ This is an interesting explanation regarding the choice of an $h\ge\sqrt\epsilon$ where $\epsilon$ is the unit roundoff, but I think the question was more oriented toward chosing an $h=h(\varepsilon)$ given $\varepsilon$ fixed, such that $x+h\neq x$ especially when $x\gg1$ so the derivative is not just plain $0$. I could be wrong but I think there is possible divergence with OP intent. $\endgroup$ – zwim Apr 3 '17 at 0:09
  • $\begingroup$ @zwim: Yes, I agree (see edits). But nonetheless, the analysis which leads to OPs choice of $h$ is flawed, and in fact there is much literature which makes this mistake, so I used this opportunity to publicize the correct reasoning. $\endgroup$ – user14717 Apr 3 '17 at 0:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy