4
$\begingroup$

I was carrying on my reading of What Every Computer Scientist Should Know About Floating-Point Arithmetic but got stuck on the proof of Theorem 2 (page 34).

At some point it says: \begin{align} (x \otimes x) \ominus (y \otimes y) & = \left[x^2(1 + \delta_1) - y^2(1 + \delta_2)\right](1 + \delta_3) \\ & = \left[(x^2 - y^2)(1 + \delta_1) + (\delta_1 - \delta_2)y^2\right](1 + \delta_3) \\ \end{align} I'm ok with the rewriting, but I don't understand the argument that:

When $x$ and $y$ are nearby, the error term $(\delta_1 - \delta_2)y^2$ can be as large as the result $x^2 - y^2$.

This doesn't make much sens to me right now. I understand that if both quantities are close to each other, then the relative error is close to $1$. But not why they should be close to each other.

$\endgroup$
3
$\begingroup$

Suppose $y = x + \Delta$. Then $(x^2 - y^2) = x^2 - (x^2 + 2\Delta x + \Delta^2) = - (2\Delta x + \Delta^2)$ with leading term on the order of $2 \Delta x$. Multiply by $1 + \delta_1$ and that's still the leading term.

Compare to $(\delta_1 - \delta_2) y^2$ with leading term $(\delta_1 - \delta_2) x^2$.

$\endgroup$
3
$\begingroup$

I'm not 100% sure of everything, but here are some elements that are implicit in the proof. Note that I reuse the notations and common knowledge of the paper, therefore this answer isn't self-contained.

First, let's assume that $x \gt y \ge 0$. It's a safe assumption to make because the sign of the variables doesn't matter and swapping $x$ and $y$ only change the sign of the result.

Let's write the relative error for the formula $(x \otimes x) \ominus (y \otimes y)$.

\begin{align} relerr & = \frac{\left|(x \otimes x) \ominus (y \otimes y) - (x^2 - y^2)\right|}{\left|x^2 - y^2\right|} \\ & = \left|\frac{(x \otimes x) \ominus (y \otimes y) - (x^2 - y^2)}{x^2 - y^2}\right| \\ & = \left|\frac{(x \otimes x) \ominus (y \otimes y)}{x^2 - y^2} - 1\right| \\ & = \left|\frac{\left[(x^2 - y^2)(1 + \delta_1) + (\delta_1 - \delta_2)y^2\right](1 + \delta_3)}{x^2 - y^2} - 1\right| \\ & = \left|\frac{(\delta_1 - \delta_2)y^2}{x^2 - y^2}(1 + \delta_3) + (1 + \delta_1)(1 + \delta_3) - 1\right| \\ \end{align}

Since we want to bound this expression by some amount of the machine-$\epsilon$, we need to find when the first term is the largest. Which happen to be when $y$ is large, making the numerator large, and when $x$ is close to $y$, making the denominator small. We could probably stop there (as the paper does), but this doesn't provide an upper bound for the relative error.

Let's assume that $x = (1 + \delta_x) y$ with $0 < \delta_x \ll 1$.

\begin{align} relerr & = \left|\frac{(\delta_1 - \delta_2)y^2}{(1 + \delta_x)^2 y^2 - y^2}(1 + \delta_3) + (1 + \delta_1)(1 + \delta_3) - 1\right| \\ & = \left|\frac{\delta_1 - \delta_2}{(1 + \delta_x)^2 - 1}(1 + \delta_3) + (1 + \delta_1)(1 + \delta_3) - 1\right| \\ & = \left|\frac{\delta_1 - \delta_2}{2 \delta_x + \delta_x^2}(1 + \delta_3) + (1 + \delta_1)(1 + \delta_3) - 1\right| \\ \end{align}

In order to bound this expression, we will need a lower bound on $\delta_x$. It needs to be greater than $0$ so that $x \neq y$ and the equations make sens. The smallest non-zero value it can take is $\delta_x \ge \beta^{-p+1} = 2 \epsilon$.

Moreover we have: $-\epsilon \le \delta_1 \le \epsilon$ and $-\epsilon \le \delta_2 \le \epsilon$.

Therefore, we have: \begin{align} \delta_1 - \delta_2 & \le 2 \epsilon \\ 2 \delta_x + \delta_x^2 & \ge 4 \epsilon + 4 \epsilon^2 \\ \frac{\delta_1 - \delta_2}{2 \delta_x + \delta_x^2} & \le \frac{2 \epsilon}{4 \epsilon + 4 \epsilon^2} & = \frac{1}{2 + \epsilon} & \lt \frac{1}{2} \\ \end{align}

Hence: \begin{align} relerr & \le \left|\frac{\delta_1 - \delta_2}{2 \delta_x + \delta_x^2}(1 + \delta_3)\right| + \left|(1 + \delta_1)(1 + \delta_3) - 1\right| \\ & \lt \frac{1 + \delta_3}{2} + |\delta_1 + \delta_3 + \delta_1 \delta_3| \\ & \lt \frac{1 + 2 \epsilon}{2} + \epsilon + 2 \epsilon + 2 \epsilon^2 \\ & \lt \frac{1}{2} + 4 \epsilon + 2 \epsilon^2 \\ & \lt \frac{1}{2} + 5 \epsilon \end{align}

Meaning that we could get an answer between 50% and 150% of the true result. I was however unable to find an actual example by exhaustive search in the high numbers. The best I could find (or worse, depending how you look at it ^^) was: $$x = 1.340780783004664349 \times 10^{154}$$ $$y = 1.340780783004664051 \times 10^{154}$$ using 64 bit floats (53 bits mantissa). Which has a relative error of roughly $0.25$.

Given that the binary mantissa of these numbers looks pretty regular, there might be something special about them and my upper bound is not correct.

My guess is that when $x$ is the float right after $y$, then maybe the rounding error of $x^2$ is similar to that of $y^2$. Which would make $\delta_1 - \delta_2 \le \epsilon$ instead of $2 \epsilon$. If anyone can prove it, feel free to edit this answer.

$\endgroup$

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.