# Understanding truncation and rounding error in IEEE floating point system?

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$?

• 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). – gnasher729 May 13 '18 at 16:17

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.

• 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$? – winnie99 May 12 '18 at 19:08
• @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. – D.W. May 12 '18 at 21:11
• 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. – gnasher729 May 12 '18 at 22:11
• @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). – D.W. May 13 '18 at 1:17