Backward Stability

563 Views Asked by At

For $x$ close to $0$ the computation

$y=\log(1-x)/x$

is numerically unstable but computing via

$y=\log(1-x)/(1-(1-x))$

is not. This question comes from http://www.cs.berkeley.edu/~demmel/cs267/lecture21/lecture21.html where they say that only the latter is backwards stable, but I can't figure out how to prove it.

1

There are 1 best solutions below

0
On

This is a neat example of the trickiness of the floating-point arithmetics. Here is a magnified plot of $\log(1-x)/x$ from the lecture:

weird plot

The horizontal scale is $x\approx 10^{-15}$. Suppose $x=1.23456789\cdot 10^{-15}$. The double-precision IEEE 754 arithmetics has no problem storing this value of $x$ by itself: the mantissa isn't too long, and the exponent is stored separately. But $1-x$ results in $0.9999999999999988$ because the mantissa only has $53$ bits, good for $16$ decimal digits at most. We lost a bit of precision here, rounding $1-x$ to the nearest representable float. But the subsequent division of $\log(1-x)$ by $x$ uses the original, non-truncated value of $x$. The result is $\approx -0.9892$, a noticeable deviation from $-1$.

Putting $1-(1-x)$ in the denominator truncates $x$ in the same way as it's truncated in the numerator. The result is $$\frac{\log(1-x)}{1-(1-x)} \approx -1.0000000000000007$$

The proof of instability could proceed as follows: the numerator $\log(1-x)$ is constant on each interval of the form $( j\cdot 2^{-54} , (j+1)\cdot 2^{-54})$, but $x$ changes by the factor of $(j+1)/j$ within this interval. Hence, the computed ratio will oscillate by a factor that grows toward $2$ as $j$ decreases. Eventually, $1-x$ underflows to $1$, and we get $0$.

The proof of stability: with the second approach, the same value of $ x$ is used in both places, so the mathematical fact about the limit of $\log(1-x)/x$ applies.