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

Question

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

@5xum Because it gives more correct answer? The first one gives 0.74... the latter 0.99... The correct answer is $\approx 1$. — mavavilj, Jan 24 '18 at 14:12
Using what software? How did you calculate the bottom expression? Because if you just write log(...) / (1+3e-16-1), then the expression 1+3e-16-1 gets evaluated first, then the division happens, and since 1+3e-16-1 will get evaluated to approximatelly, but not exactly, 3e-16, you should, in general, expect a bigger mistake. — 5xum, Jan 24 '18 at 14:15
Please explain what results give rise to your description of stability vs. Instability. Since these are "one-off" computations, it is unclear how you meant that. — hardmath, Jan 24 '18 at 14:15
@mavavilj The code snippet, please. I'm trying to recreate the problem so I can help you. Please, be helpful and give as much information as possible. — 5xum, Jan 24 '18 at 14:16
just when you evaluate 1+x-1 in the console happens to error out! from 3e-16 to 2.22e-16 — Abr001am, Jan 24 '18 at 16:03

5xum · Accepted Answer · 2018-01-24T14:30:11.713

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

A very nice answer, including a warning to draw conclusions too fast! (+1) — Peter, Jan 24 '18 at 14:29
In fact, you can circumvent the issue with the log in this particular situation by using "log1p" (available in the math library of a lot of languages), which will show similar behavior to the sine example. — Ian, Jan 24 '18 at 14:34

score 1 · Answer 2 · answered Jan 24 '18 at 14:34

1

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))

answered Jan 24 '18 at 14:34

Does this mean that it's not numerical instability? Or is it? – mavavilj Jan 24 '18 at 14:36
@mavavilj Why would it not be numerical instability, if this answer is talking about rounding errors? – 5xum Jan 24 '18 at 14:38
@mavavilj Since you didn't give a definition of "numerical stabiliy", I don't know what you mean. I've told you what the reason of the error was, and how exaxctly it can be avoided. – Jan 24 '18 at 14:38
@5xum Okay it's. I didn't read the "rounded to". So this is a round-off error? – mavavilj Jan 24 '18 at 14:39
I'm referring to https://en.wikipedia.org/wiki/Numerical_stability – mavavilj Jan 24 '18 at 14:39
@mavavilj Such references belong into the question, not some comment down the thread. – Jan 24 '18 at 14:42

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

2 Answers2