Catastrophic cancellation

In numerical analysis, catastrophic cancellation[1][2] is the phenomenon that subtracting good approximations to two nearby numbers may yield a very bad approximation to the difference of the original numbers.

For example, suppose you have two wood studs, one ${\displaystyle L_{1}=254.5\,{\text{cm}}}$ long and the other ${\displaystyle L_{2}=253.5\,{\text{cm}}}$ long. If you measure them with a ruler that is good only to the centimeter, you might get approximations ${\displaystyle {\tilde {L}}_{1}=255\,{\text{cm}}}$ and ${\displaystyle {\tilde {L}}_{2}=253\,{\text{cm}}}$. Depending on your needs, these may be good approximations, in relative error, to the true lengths: the approximations are in error by less than 2% of the true lengths, ${\displaystyle |L_{1}-{\tilde {L}}_{1}|/|L_{1}|<2\%}$ .

However, if you subtract the approximate lengths, you will get ${\displaystyle {\tilde {L}}_{1}-{\tilde {L}}_{2}=255\,{\text{cm}}-253\,{\text{cm}}=2\,{\text{cm}}}$, even though the true difference between the lengths is ${\displaystyle L_{1}-L_{2}=254.5\,{\text{cm}}-253.5\,{\text{cm}}=1\,{\text{cm}}}$. The difference of the approximations, ${\displaystyle 2\,{\text{cm}}}$, is in error by 100% of the magnitude of the difference of the true values, ${\displaystyle 1\,{\text{cm}}}$.

Catastrophic cancellation may happen even if the difference is computed exactly, as in the example above—it is not a property of any particular kind of arithmetic like floating-point arithmetic; rather, it is inherent to subtraction, when the inputs are approximations themselves. Indeed, in floating-point arithmetic, when the inputs are close enough, the floating-point difference is computed exactly, by the Sterbenz lemma—there is no rounding error introduced by the floating-point subtraction operation.

Formal analysis

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 approximate difference ${\displaystyle {\tilde {x}}-{\tilde {y}}}$ from the true difference ${\displaystyle x-y}$ is inversely proportional to the true difference:

$${\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 numbers is

$${\displaystyle {\biggl |}{\frac {x\delta _{x}-y\delta _{y}}{x-y}}{\biggr |}.}$$

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

In numerical algorithms

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.

Example: Difference of squares

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 will amplify 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 nearest to ${\displaystyle 1.8626451{\underline {49230957}}\times 10^{-9}}$, of which only half the digits are correct and the other half (underlined) are garbage.

Example: Complex arcsine

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

In the case of ${\displaystyle z=iy}$ for ${\displaystyle y<0}$, using the identity ${\displaystyle \arcsin(z)=-\arcsin(-z)}$ avoids cancellation because ${\displaystyle {\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.

Example: Radix conversion

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.

Benign cancellation

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 \operatorname {fl} (1+x)-1}$ to undo the error.[2] This works when ${\displaystyle x}$ is near zero (though not exactly zero) because ${\displaystyle \operatorname {fl} (1+x)}$ is near ${\displaystyle 1}$, so ${\displaystyle \log(\operatorname {fl} (1+x))\approx \operatorname {fl} (1+x)-1}$, and thus the numerator and denominator of the quotient cancel, leaving only ${\displaystyle x\approx \log(1+x)}$ as desired.

References

1. Muller, Jean-Michel; Brunie, Nicolas; de Dinechin, Florent; Jeannerod, Claude-Pierre; Joldes, Mioara; Lefèvre, Vincent; Melquiond, Guillaume; Revol, Nathalie; Torres, Serge (2018). Handbook of Floating-Point Arithmetic (2nd ed.). Gewerbestrasse 11, 6330 Cham, Switzerland: Birkhäuser. p. 102. doi:10.1007/978-3-319-76526-6. ISBN 978-3-319-76525-9.CS1 maint: location (link)
2. Goldberg, David (March 1991). "What every computer scientist should know about floating-point arithmetic". ACM Computing Surveys. New York, NY, United States: Association for Computing Machinery. 23 (1): 5–48. doi:10.1145/103162.103163. ISSN 0360-0300. Retrieved 2020-09-17.