Is adding and subtracting 1 some generic method of removing numerical instability?

Question

Is adding and subtracting 1 some generic method of removing numerical instability?

Like in:

Unstable: $$\frac{\log(1+3 \cdot 10^{-16})}{3 \cdot 10^{-16}}$$

Stable:

$$\frac{\log(1+3 \cdot 10^{-16})}{1+3 \cdot 10^{-16}-1}$$

The correct answer (of the first expression) is $\approx 1$.

Computing these in Python 3:

import math

math.log(1+3*pow(10,-16))/(3*pow(10,-16))
# > 0.7401486830834376

math.log(1+3*pow(10,-16))/(1+3*pow(10,-16)-1)
# > 0.9999999999999999

---

For a definition of numerical stability, I refer to:

https://en.wikipedia.org/wiki/Numerical_stability

Answer 1

What you are experiencing is pure coincidence, and a nice example of how one example cannot be used to draw general conclusion.

When you write

 math.log(1+3*pow(10,-16))/(1+3*pow(10,-16)-1)

the program first calculates the expression on the right and gets a value of (1+3*pow(10,-16)-1) = 2.220446049250313e-16, which is quite a ways off from the actual value of 3e-16.

Then it calculates the expression on the left and gets math.log(1+3*pow(10,-16)) = 2.2204460492503128e-16 which again, is a long way off from the actual value which should be close to 3e-16.

The reason the two errors are the same is that calculating the logarithm is probably done by using its power series, i.e.

$$\ln(1+x)=x+\frac12 x^2+\cdots$$

which is reworked into

$$\ln(x) = \ln(1+(x-1)) = (x-1) + \frac12 (x-1)^2 + \cdots$$

so you are, when calculating log(1+3e-16), forcing the software to calculate (1+3e-16) - 1, the exact same computation (resulting in the exact same error) it did when calculating the denominator.

It's by pure coincidence that you discovered a function (i.e., log), that makes the same numerical error as the addition function, so the two errors cancel out.

---

In general, you should expect a bigger error when calculating $$\frac{f(x)}{1+x-1}$$ when compared to $\frac{f(x)}{x}$ for small values of $x$, because calculating $1+x-1$ will already result in sometimes big errors.

For example, taking $f(x) = \sin x$ should show you that your modified calculation is not, in general, more stable.

>>> from math import sin
>>> x=3e-16
>>> sin(x)/x
1.0
>>> sin(x)/(1+x-1)
1.3510798882111488

Answer 2

The reason for the inaccuracy is that $1+3\cdot10^{-16}$ is rounded to the next machine number, which happens to be not very exactly $3\cdot10^{-16}$ away from $1$. That's why there is a special function in math calculating $\log(1+x)$ directly from $x$, it's called log1p(x). So just try calculating

 math.log1p(3*pow(10,-16))/(3*pow(10,-16))