Numpy Subtraction Precision

Question

I have a numpy array, x, that is calculated from a complicated equation:

x = (QT-m_μ_Q*M_T)/(m*σ_Q*Σ_T)
print(x)
print(x[0], x[1], x[2])
print(1.0-x)

This prints:

[ 1.  1.  1.]
1.0 1.0 1.0
[ -2.22044605e-16   3.33066907e-16  -4.44089210e-16]

Notice that the last line prints out something small and very close to zero but is non-zero. My next step is to take the square root of each value, so it should not contain negative values.

Explicitly starting with an array containing ones:

y = np.array([1., 1., 1.])
print(y)
print(y[0], y[1], y[2])
print(1.0-y)

produces the correct result so I'm not sure what the difference is:

[ 1.  1.  1.]
1.0 1.0 1.0
[ 0.  0.  0.]

Possible duplicate of [Is floating point math broken?](http://stackoverflow.com/questions/588004/is-floating-point-math-broken) — Tadhg McDonald-Jensen, Mar 02 '17 at 19:26
Try `print((x[0], x[1], x[2]))` - by putting them in a tuple, you force `repr` instead of `str`- `repr` of scalars always tries to print the full precision in numpy. Don't think is a duplicate - the problem here is that `str(f)` is hiding information — Eric, Mar 02 '17 at 20:57

score 3 · Answer 1 · edited May 23 '17 at 12:17

In short the ones in x are not truely 1. That's because you operate (the formula you used) on inexact values (floats are inexact, see also Is floating point math broken?).

However the default precision doesn't show the difference. But when you subtract one you get what is generally known as catastrophic cancellation (see for example Wikipedia or Why is 'catastrophic cancellation' called so?).

For example:

>>> import numpy as np
>>> np.array([1, 1.00000000000002, 1.000000000000002, 1.0000000000000002])
array([ 1.,  1.,  1.,  1.])

>>> 1 - np.array([1, 1.00000000000002, 1.000000000000002, 1.0000000000000002])
array([  0.00000000e+00,  -1.99840144e-14,  -1.99840144e-15, -2.22044605e-16])

Generally there are two approaches to fix this:

Use a formula that avoids catastrophic cancellation (may not be possible). However in your case you could try sqrt((y-x)/y) instead of sqrt(1-x/y) - just insert the correct expressions for x (QT-m_μ_Q*M_T) and y (m*σ_Q*Σ_T). Mathematically both are equivalent but due to the floating point math these may give different results.
Round your results and/or set negative values to zero.

Numpy Subtraction Precision

1 Answers1