Proof that (x-y)(x+y) is more accurate than x²-y²

Question

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.

Answer 1

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

Answer 2

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.