1
$\begingroup$

I started learning about floating point by reading "What Every Computer Scientist Should know About Floating-Point Arithmetic" by David Goldberg. On page 4 he presents a proof for the maximum relative error suffered when rounding any real number to a floating-point number, as follows:

To compute the relative error that corresponds to $1/2\ \text{ulp}$, observe that when a real number is approximated by the closest possible floating-point number $\overbrace{d\ dd\ \cdots\ dd}^p \times \beta^e$, the absolute error can be as large as $\overbrace{0\ 00\ \cdots\ 00}^p\beta' \times \beta^e$ where $\beta'$ is the digit $\beta/2$. This error is $((\beta/2)\beta^{-p}) \times \beta^e$. Since numbers of the form $d.dd\ \cdots\ dd \times \beta^e$ all have this same absolute error but have values that range between $\beta^e$ and $\beta \times \beta^e$, the relative error ranges between $((\beta/2)\beta^{-p} \times \beta^e/\beta^e$ and $((\beta/2)\beta^{-p} \times \beta^e/\beta^{e+1}$. That is,

$$\frac 12 \beta^{-p} \le \frac 12 \text{ulp} \le \frac{\beta}{2} \beta^{-p}. \tag 2$$

[...] Setting $\epsilon = (\beta/2)\beta^{-p}$ to the largest of the bounds in $(2)$, we can say that when a real number is rounded to the closest floating-point number, the relative error is always bounded by $\epsilon$, which is referred to as machine epsilon.

In this text $\beta$ and $p$ are the base and precision of a given floating-point representation, and $e$ is the exponent used in the approximate representation of a real number. In $(2)$, I believe $\frac 12 \text{ulp}$ refers to the relative error associated with a $1/2\ \text{ulp}$ approximation error.

I understood this proof, however reading the Wikipedia article on the machine epsilon it suggests that this is a simplification without developing the unsimplified version. The above proof (equivalent to the one outlined on Wikipedia) assumes that the largest relative error occurs when representing the real number $1.0 \times \beta^e$. However, this number is actually exactly represented; in order to correct for this, a small term must be added to reach the real number half way to the next floating-point number, corresponding to $1/2\ \text{ulp} = ((\beta/2)\beta^{-p})$, which must be added to the denominator when calculating the relative error. This lead me to the following formula:

$$\begin{align} \frac{\frac{\beta}{2}\beta^{-p} \times \beta^e}{\beta^e + \frac{\beta}{2}\beta^{-p}} & = \frac{\frac 12 \beta^{1-p} \times \beta^e}{\frac 12 \beta^{1-p} + \beta^e} && \text{Relative error of $1.0 \times \beta^e + 1/2$ ulp} \tag 1 \newline \frac{\beta^{-e}}{\beta^{-e}} \left(\frac{\frac 12 \beta^{1-p} \times \beta^e}{\frac 12 \beta^{1-p} + \beta^e}\right) & = \frac{\frac 12 \beta^{1-p}}{\frac 12 \beta^{1-p-e} + 1} \tag 2 \newline \frac{2 \beta^{p-1}}{2 \beta^{p-1}} \left(\frac{\frac 12 \beta^{1-p}}{\frac 12 \beta^{1-p-e} + 1}\right) & = \frac{1}{\beta^{-e} + 2 \beta^{p-1}} \tag 3 \end{align}$$

Since we were no longer able to eliminate the term $e$, we know that the actual relative error when rounding a real number to a floating-point number depends on the exponent used to represent it. Taking the limit, we can proceed as follows:

$$\begin{align} \lim_{e \to +\infty} \frac{1}{\beta^{-e} + 2 \beta^{p-1}} & = \frac 12 \beta^{-(p-1)} = \frac{\beta}{2} \beta^{-p} \newline \lim_{e \to -\infty} \frac{1}{\beta^{-e} + 2 \beta^{p-1}} & = 0 \end{align}$$

The first limit proves that the simplified definition of the machine epsilon is true in the limit when the exponent $e$ used to represent a real number approaches infinity, meaning that it is an over-estimation of the relative error. This makes sense because adding the $1/2\ \text{ulp}$ term to the denominator can only decrease the relative error.

In practice, for $e \gt -p$ the difference between the relative error and this over-estimation is extremely small. However, the IEEE 754 single precision floating-point format has an $e_{min} = -126$ and $p = 24$, with double precision $e_{min} = -1022$ and $p = 53$, meaning that a significant amount of the representable floating-point numbers will have a nearly zero relative rounding error, according to this formula.

If this is true, then supposing that a computation doesn't require the entire dynamic range of the floating-point type being used, would it be possible to scale the values down to decrease rounding errors and potentially improve numerical stability? If the result had to be scaled back up at the end of the computation then it would be subject to higher $e$ rounding error, but intermediary computations would not. I strongly suspect that everything I've done is completely wrong as I am new to floating-point. I appreciate any help!

$\endgroup$
3
  • $\begingroup$ You might want to read again what a relative error is. $\endgroup$ Commented Feb 21, 2023 at 12:13
  • $\begingroup$ The calculations presented in my post are for the relative error. I showed that it is incorrect to take the relative error of the floating-point number $1.0 \times \beta^e$ (which is normally used to establish the maximum relative error bound) because this number is exactly represented, so you must take the relative error of the real number $1.0 \times \beta^e + \frac{\beta}{2} \beta^{-p}$, which is approximated by the floating point number $1.0 \times \beta^e$. $\endgroup$ Commented Feb 21, 2023 at 23:19
  • 2
    $\begingroup$ Ah, well I realized the error in my formula now, fortunately. I forgot to multiply the $1/2\ \text{ulp} = \frac{\beta}{2} \beta^{-p}$ figure by $\beta^e$, because of course the absolute error introduced by $1/2\ \text{ulp}$ is relative to the exponent used. I'll try to answer my own question shortly. $\endgroup$ Commented Feb 21, 2023 at 23:23

1 Answer 1

1
$\begingroup$

In step $(1)$ of my formula I defined an equation for the relative error of $1.0 \times \beta^e + 1/2\ \text{ulp}$. This equation was incorrect because I did not multiply the $1/2\ \text{ulp}$ term by $\beta^e$, which of course is necessary because the absolute value of $1/2\ \text{ulp}$ changes with $e$. Doing this allows $\beta^e$ to be factored out as before, such that the relative error is no longer a function of the exponent $e$ used to represent the floating-point number approximating a real number.

\begin{align} \frac{\frac{\beta}{2}\beta^{-p} \times \beta^e}{\beta^e} & = \frac{\beta}{2}\beta^{-p} && \text{Maximum relative error as defined by references} \newline \frac{\frac{\beta}{2}\beta^{-p} \times \beta^e}{\beta^e + \frac{\beta}{2}\beta^{-p}} & = \frac{1}{\beta^{-e} + 2\beta^{p-1}} && \text{Incorrect definition of relative error from question} \newline \frac{\frac{\beta}{2}\beta^{-p} \times \beta^e}{\beta^e + \frac{\beta}{2}\beta^{-p}\beta^e} & = \frac{\frac 12 \beta^{1-p} \times \beta^e}{\left(\frac 12 \beta^{1-p} + 1\right) \beta^e} && \text{Correct relative error of $1.0 \times \beta^e + 1/2$ ulp} \newline & = \frac{\frac 12 \beta^{1-p}}{\frac 12 \beta^{1-p} + 1} \newline & = \frac{1}{1 + 2\beta^{p-1}} \newline & = \frac{1}{1 + \frac{2}{\beta} \beta^{p}} \newline & \approx \frac{\beta}{2} \beta^{-p} && \text{for $\beta^p \gg 1$.} \end{align}

$\endgroup$

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .