Catastrophic cancellation

Formally, catastrophic cancellation happens because subtraction is ill-conditioned at nearby inputs: even if approximations ${\displaystyle {\tilde {x}}=x(1+\delta _{x})}$  and ${\displaystyle {\tilde {y}}=y(1+\delta _{y})}$  have small relative errors ${\displaystyle |\delta _{x}|=|x-{\tilde {x}}|/|x|}$  and ${\displaystyle |\delta _{y}|=|y-{\tilde {y}}|/|y|}$  from true values ${\displaystyle x}$  and ${\displaystyle y}$ , respectively, the relative error of the difference ${\displaystyle {\tilde {x}}-{\tilde {y}}}$  of the approximations from the difference ${\displaystyle x-y}$  of the true values is inversely proportional to the difference of the true values:

$${\displaystyle {\begin{aligned}{\tilde {x}}-{\tilde {y}}&=x(1+\delta _{x})-y(1+\delta _{y})=x-y+x\delta _{x}-y\delta _{y}\\&=x-y+(x-y){\frac {x\delta _{x}-y\delta _{y}}{x-y}}\\&=(x-y){\biggr (}1+{\frac {x\delta _{x}-y\delta _{y}}{x-y}}{\biggr )}.\end{aligned}}}$$ 

Thus, the relative error of the exact difference ${\displaystyle {\tilde {x}}-{\tilde {y}}}$  of the approximations from the difference ${\displaystyle x-y}$  of the true values is

$${\displaystyle \left|{\frac {x\delta _{x}-y\delta _{y}}{x-y}}\right|.}$$ 

which can be arbitrarily large if the true values ${\displaystyle x}$  and ${\displaystyle y}$  are close.

Subtracting nearby numbers in floating-point arithmetic does not always cause catastrophic cancellation, or even any error—by the Sterbenz lemma, if the numbers are close enough the floating-point difference is exact. But cancellation may amplify errors in the inputs that arose from rounding in other floating-point arithmetic.

Given numbers ${\displaystyle x}$  and ${\displaystyle y}$ , the naive attempt to compute the mathematical function ${\displaystyle x^{2}-y^{2}}$  by the floating-point arithmetic ${\displaystyle \operatorname {fl} (\operatorname {fl} (x^{2})-\operatorname {fl} (y^{2}))}$  is subject to catastrophic cancellation when ${\displaystyle x}$  and ${\displaystyle y}$  are close in magnitude, because the subtraction can expose the rounding errors in the squaring. The alternative factoring ${\displaystyle (x+y)(x-y)}$ , evaluated by the floating-point arithmetic ${\displaystyle \operatorname {fl} (\operatorname {fl} (x+y)\cdot \operatorname {fl} (x-y))}$ , avoids catastrophic cancellation because it avoids introducing rounding error leading into the subtraction.^[2]

For example, if ${\displaystyle x=1+2^{-29}\approx 1.0000000018626451}$  and ${\displaystyle y=1+2^{-30}\approx 1.0000000009313226}$ , then the true value of the difference ${\displaystyle x^{2}-y^{2}}$  is ${\displaystyle 2^{-29}\cdot (1+2^{-30}+2^{-31})\approx 1.8626451518330422\times 10^{-9}}$ . In IEEE 754 binary64 arithmetic, evaluating the alternative factoring ${\displaystyle (x+y)(x-y)}$  gives the correct result exactly (with no rounding), but evaluating the naive expression ${\displaystyle x^{2}-y^{2}}$  gives the floating-point number ${\displaystyle 2^{-29}=1.8626451{\underline {4923095703125}}\times 10^{-9}}$ , of which less than half the digits are correct and the other (underlined) digits reflect the missing terms ${\displaystyle 2^{-59}+2^{-60}}$ , lost due to rounding when calculating the intermediate squared values.

When computing the complex arcsine function, one may be tempted to use the logarithmic formula directly:

$${\displaystyle \arcsin(z)=i\log {\bigl (}{\sqrt {1-z^{2}}}-iz{\bigr )}.}$$ 

However, suppose ${\displaystyle z=iy}$  for ${\displaystyle y\ll 0}$ . Then ${\displaystyle {\sqrt {1-z^{2}}}\approx -y}$  and ${\displaystyle iz=-y}$ ; call the difference between them ${\displaystyle \varepsilon }$ —a very small difference, nearly zero. If ${\displaystyle {\sqrt {1-z^{2}}}}$  is evaluated in floating-point arithmetic giving

$${\displaystyle \operatorname {fl} {\Bigl (}{\sqrt {\operatorname {fl} (1-\operatorname {fl} (z^{2}))}}{\Bigr )}={\sqrt {1-z^{2}}}(1+\delta )}$$ 

with any error ${\displaystyle \delta \neq 0}$ , where ${\displaystyle \operatorname {fl} (\cdots )}$  denotes floating-point rounding, then computing the difference

$${\displaystyle {\sqrt {1-z^{2}}}(1+\delta )-iz}$$ 

of two nearby numbers, both very close to ${\displaystyle -y}$ , may amplify the error ${\displaystyle \delta }$  in one input by a factor of ${\displaystyle 1/\varepsilon }$ —a very large factor because ${\displaystyle \varepsilon }$  was nearly zero. For instance, if ${\displaystyle z=-1234567i}$ , the true value of ${\displaystyle \arcsin(z)}$  is approximately ${\displaystyle -14.71937803983977i}$ , but using the naive logarithmic formula in IEEE 754 binary64 arithmetic may give ${\displaystyle -14.719{\underline {644263563968}}i}$ , with only five out of sixteen digits correct and the remainder (underlined) all incorrect.

In the case of ${\displaystyle z=iy}$  for ${\displaystyle y<0}$ , using the identity ${\displaystyle \arcsin(z)=-\arcsin(-z)}$  avoids cancellation because ${\textstyle {\sqrt {1-(-z)^{2}}}={\sqrt {1-z^{2}}}\approx -y}$  but ${\displaystyle i(-z)=-iz=y}$ , so the subtraction is effectively addition with the same sign which does not cancel.

Numerical constants in software programs are often written in decimal, such as in the C fragment double x = 1.000000000000001; to declare and initialize an IEEE 754 binary64 variable named x. However, ${\displaystyle 1.000000000000001}$  is not a binary64 floating-point number; the nearest one, which x will be initialized to in this fragment, is ${\displaystyle 1.0000000000000011102230246251565404236316680908203125=1+5\cdot 2^{-52}}$ . Although the radix conversion from decimal floating-point to binary floating-point only incurs a small relative error, catastrophic cancellation may amplify it into a much larger one:

double x = 1.000000000000001;  // rounded to 1 + 5*2^{-52}
double y = 1.000000000000002;  // rounded to 1 + 9*2^{-52}
double z = y - x;              // difference is exactly 4*2^{-52}

The difference ${\displaystyle 1.000000000000002-1.000000000000001}$  is ${\displaystyle 0.000000000000001=1.0\times 10^{-15}}$ . The relative errors of x from ${\displaystyle 1.000000000000001}$  and of y from ${\displaystyle 1.000000000000002}$  are both below ${\displaystyle 10^{-15}=0.0000000000001\%}$ , and the floating-point subtraction y - x is computed exactly by the Sterbenz lemma.

But even though the inputs are good approximations, and even though the subtraction is computed exactly, the difference of the approximations ${\displaystyle {\tilde {y}}-{\tilde {x}}=(1+9\cdot 2^{-52})-(1+5\cdot 2^{-52})=4\cdot 2^{-52}\approx 8.88\times 10^{-16}}$  has a relative error of over ${\displaystyle 11\%}$  from the difference ${\displaystyle 1.0\times 10^{-15}}$  of the original values as written in decimal: catastrophic cancellation amplified a tiny error in radix conversion into a large error in the output.

Cancellation is sometimes useful and desirable in numerical algorithms. For example, the 2Sum and Fast2Sum algorithms both rely on such cancellation after a rounding error in order to exactly compute what the error was in a floating-point addition operation as a floating-point number itself.

The function ${\displaystyle \log(1+x)}$ , if evaluated naively at points ${\displaystyle 0<x\lll 1}$ , will lose most of the digits of ${\displaystyle x}$  in rounding ${\displaystyle \operatorname {fl} (1+x)}$ . However, the function ${\displaystyle \log(1+x)}$  itself is well-conditioned at inputs near ${\displaystyle 0}$ . Rewriting it as $${\displaystyle \log(1+x)=x{\frac {\log(1+x)}{(1+x)-1}}}$$  exploits cancellation in ${\displaystyle {\hat {x}}:=\operatorname {fl} (1+x)-1}$  to avoid the error from ${\displaystyle \log(1+x)}$  evaluated directly.^[2] This works because the cancellation in the numerator ${\displaystyle \log(\operatorname {fl} (1+x))={\hat {x}}+O({\hat {x}}^{2})}$  and the cancellation in the denominator ${\displaystyle {\hat {x}}=\operatorname {fl} (1+x)-1}$  counteract each other; the function ${\displaystyle \mu (\xi )=\log(1+\xi )/\xi }$  is well-enough conditioned near zero that ${\displaystyle \mu ({\hat {x}})}$  gives a good approximation to ${\displaystyle \mu (x)}$ , and thus ${\displaystyle x\cdot \mu ({\hat {x}})}$  gives a good approximation to ${\displaystyle x\cdot \mu (x)=\log(1+x)}$ .