3
$\begingroup$

I'm trying to understand the theory behind finding the optimal $h$ value for differentiation in this definition:

$$ \frac{f(x+h) - f(x)}{h}$$

as $h$ tends to 0.

Here is my understanding:

  1. Truncation error comes from say approximating the value of $f(x+h)$ using the taylor series. In this case, it could be that:

$$ f(x+h) \approx f(x) + hf'(x) + \frac{h^2}{2!}f''(x) + O(h^3)$$

  1. Letting $f'(x)$ to be the subject since we want to find gradient $f'(x)$, we have

$$f'(x) = \frac{f(x) - (hf'(x) + \frac{h^2}{2!}f''(x) + O(h^3))}{h}$$ $$\implies f'(x) = \frac{f(x+h)-f(x)}{h} - \frac{h}{2}f''(x) + O(h^2)$$

Intuitively, we can stop at here by letting $h \approx \sqrt\epsilon$ where $\epsilon$ is the machine epsilon, we can thus limit our truncation error. But suppose we continue, then our truncation error is $\frac{h}{2}f''(x)$.

  1. Then rounding error comes from the interaction of the sums in the differentiation equation:

$$ \frac{f(x+h) - f(x)}{h} = \frac{[f(x) + hf'(x) + O(h^2)] - f(x)}{h}$$

But for $f'(x) \approx 1$ and since we have $h \approx \sqrt \epsilon$, we have

$$ \frac{f(x+h) - f(x)}{h} \approx 1 \pm \frac{\epsilon}{h}$$

Which means the rounding error is approximately $\frac{\epsilon}{h}$, which concurs with the fact that the left hand side of the equation is actually $f'(x)$ as well, just that we have managed to reduce the rounding error to what is seen on the right hand side.

  1. Now, we see that our truncation error varies linearly with $h$ while our rounding error varies inversely with $h$. To find the optimal $h$ value, we equate both errors:

$$\frac{\epsilon}{h} = \frac{h}{2}f''(x) \approx \frac{h}{2}$$ since we have $f''(x) \approx 1$.

Letting $h$ be the subject, we have $$h = \sqrt{2\epsilon}$$ which means the optimal $h$ is around $h = \sqrt{\epsilon}$ as expected.

Is my understanding correct? Does the final factor of 2 in the end matter?

Big problem in the above argument: in point 3, one already use the fact that $h = \sqrt \epsilon$. Wouldn't that already defeat the purpose of finding the optimal $h$?

$\endgroup$
  • $\begingroup$ It doesn't really make sense to calculate f(x+h) using a Taylor series in order to find f'(x), because the Taylor series already requires f'(x). $\endgroup$ – gnasher729 May 13 '18 at 16:17
1
$\begingroup$

As listed your reasoning is indeed circular: you start by assuming $h \approx \sqrt{\epsilon}$, then use the results in point 3 to infer that we should take $h \approx \sqrt{\epsilon}$. So that's not convincing. So let's redo the analysis more carefully, then we'll look at the meaning of the constant.

A more careful analysis

As in your analysis, we'll assume $f'(x) \approx 1$ and $f''(x) \approx 1$.

Your analysis of truncation error is fine; the truncation error is indeed $\frac{h}{2} f''(x) + O(h^2)$, which given your assumptions is $\frac{h}{2} + O(h^2)$. In our analysis $h$ will be small enough that $O(h^2)$is negligible, so we'll call the truncation error $\frac{h}{2}$ for simplicity.

What about rounding error? The rounding error of representing $f(h)$ in floating point is $\epsilon$ in relative error, where $\epsilon$ is the machine epsilon; in absolute error, that is $\epsilon f(h)$. Since you assumed $f(h) \approx 1$, that is $\approx \epsilon$. Similarly, the rounding error of representing $f(x+h)$ is $\epsilon f(x+h)$ in absolute error; since $f(x+h) \approx f(x) \approx 1$, that is $\approx \epsilon$ as well. Now the absolute error of $f(x+h)-f(x)$ is at most about $2\epsilon + \epsilon (f(x+h)-f(x))$; the $2\epsilon$ comes from the error in each of the two subterms (and the triangle inequality), and the $\epsilon (f(x+h)-f(x))$ comes from roundoff error in representing the result. Since $f(x+h)-f(x) \approx h f'(x) \approx h$, that means the absolute error of $f(x+h)-f(x)$ is at most about $2\epsilon + \epsilon h$.

Then we divide by $h$. The absolute error becomes $2\epsilon/h + \epsilon$, plus a term due to representing the result in floating point; that extra term is $\epsilon ((f(x+h)-f(x)/h)$, and since $(f(x+h)-f(x)/h \approx f'(x) \approx 1$, we see that the extra term is $\approx \epsilon$. Therefore the absolute error in computing $(f(x+h)-f(x)/h$ is about $2\epsilon/h + \epsilon + \epsilon$. So, the total roundoff error is about $2\epsilon/h + 2\epsilon$, in absolute error.

The total error is at most the sum of the truncation error and the roundoff error, or about $$\frac{h}{2} + {2\epsilon \over h} + 2\epsilon.$$ We want to know what value of $h$ minimizes this expression. We can figure that out by letting $g(h) = \frac{h}{2} + 2\epsilon/h + 2\epsilon$ and minimizing $g(h)$; the minima must be at a point where the derivative is zero. So, differentiate: $$g'(h) = \frac{1}{2} - {2\epsilon \over h^2}.$$ Setting $g'(h)=0$ and solving for $h$ yields $h^2 =4 \epsilon$, i.e., the function $g(h)$ is minimized at $h=2\sqrt{\epsilon}$; at that point, we have a total (absolute) error of at most $g(2\sqrt{\epsilon}) \approx 2 \sqrt{\epsilon} + 2 \epsilon \approx 2\sqrt{\epsilon}$.

The summary: assuming $f(x)\approx 1$ and $f'(x) \approx 1$, the total error is minimized when you choose $h$ to be about $h\approx 2\sqrt{\epsilon}$. If you choose $h$ too much larger than that, the truncation error increases and starts dominating. If you choose $h$ too much smaller than that, the roundoff error increases and starts dominating.

So, if you know $f(x)\approx 1$ and $f'(x) \approx 1$, then it's reasonable to choose $h\approx 2\sqrt{\epsilon}$.

The constant

Now, you might be wondering a bunch of things. Why the constant 2? Where did that come from? How critical is it for the constant to be 2? What if I use $h \approx \sqrt{\epsilon}$ instead of $h \approx \sqrt{2\epsilon}$?

Well, if you redo the calculation with a different value of $h$, you'll see that the total error is not too sensitive to the constant 2. For instance, if you choose $h=\epsilon$, then the total error will be at most about $g(\epsilon) \approx 2.5 \sqrt{\epsilon} + 2 \epsilon \approx 2.5 \sqrt{\epsilon}$, which is not that different. So don't worry too much about the exact value of the constant.

Where did the constant 2 came from? Well, partly it is coming from the fact that we are doing multiple operations, and you can get a little bit of roundoff error from each operation, so the error can accumulate. And, partly it is coming from the assumptions that $f'(x) \approx 1$ and $f''(x) \approx 1$. Of course, those assumptions are pretty sketchy. In practice, we don't know what $f'(x)$ is before doing the computation (that's the whole point of this computation; we are trying to estimate/compute $f'(x)$), and we don't know what $f''(x)$ is. So, that assumption also looks circular.

But here's why in this case it might be OK. You could redo the calculation with different assumptions, e.g., $f'(x) \approx 2$ and $f''(x) \approx 3$. If you do that, you'll find a slightly different optimal value of $h$; instead of $2 \sqrt{\epsilon}$, you'll have a slightly different constant replacing the 2. So we can see that the 2 comes from the assumptions about the magnitude of $f'(x)$ and $f''(x)$. But, crucially, if $f'(x)$ and $f''(x)$ aren't too large or too small, then the resulting constant in the optimal choice of $h$ won't be too different. It might be $5 \sqrt{\epsilon}$ or $0.5 \sqrt{\epsilon}$, but that's not so bad.

And, in practice, we're often dealing with smooth functions where the derivatives don't get too large, so while we don't know exactly what $f'(x)$ or $f''(x)$ are, we might be able to assume they won't be too large. As a result, taking $h \approx c \sqrt{\epsilon}$ for some small value of $c$ will be optimal. We don't know what the optimal value of $c$ is exactly, but we can just heuristically pick a small value for $c$, like $c=1$ or $c=2$, and even if that's not exactly optimal, we'll probably do reasonably well -- it won't be too much worse.

So, the final answer is: no, the 2 is not essential. It's perfectly reasonable to use $h = \sqrt{\epsilon}$ instead of $h = 2 \sqrt{\epsilon}$ -- but you probably want to choose $h$ somewhere in that neighborhood.

$\endgroup$
  • $\begingroup$ A really great answer. Can you explain what's the triangle inequality in the statement " the 2ϵ comes from the error in each of the two subterms (and the triangle inequality)"? Why must we include this 2 $\epsilon$ value when the computations (the minus) already give an absolute error of 2$\epsilon$? $\endgroup$ – winnie99 May 12 '18 at 19:08
  • $\begingroup$ @winnie99, Glad it was helpful! If $|t-t^*| \le \epsilon$ and $|u-u^*| \le \epsilon$, then $|(t-u)-(t^*-u^*)| \le 2\epsilon$ (that's the triangle inequality). Now think of $t$ as the correct value (of $f(x+h)$) and $t^*$ as the machine representation in floating point, and think of $u$ as the correct value (of $f(x)$) and $u^*$ as its machine representation, and hopefully all will become clear. I think this is capturing that the prior computations already contribute $2\epsilon$ to the error. $\endgroup$ – D.W. May 12 '18 at 21:11
  • $\begingroup$ With IEEE floating point arithmetic, the difference of two nearby numbers will be calculated exactly. This is always the case for x - y where x/2 ≤ y ≤ 2x. This has two consequences: One, when calculating f (x+h), we already have a rounding error in the addition x + h. But (x + h) - x is the exact difference between x+h and x, so if we let h' = (x+h) - x, then x+h' will be calculated exactly, without rounding error. So we just replace h with (x+h) - x. The other effect is that f(x+h) - f(x) will be calculated without rounding errors. $\endgroup$ – gnasher729 May 12 '18 at 22:11
  • $\begingroup$ @gnasher729, cool, I didn't know that fact! Thank you. So my estimates of the error are overestimates (but only by a constant factor, so it only affects the constants in the final answer, not otherwise). $\endgroup$ – D.W. May 13 '18 at 1:17

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.