Numpy rounds in a different way than python

Question

The code

import numpy as np
a = 5.92270987499999979065
print(round(a, 8))
print(round(np.float64(a), 8))

gives

5.92270987
5.92270988

Any idea why?

Found nothing relevant in numpy sources.

Update:
I know that the proper way to deal with this problem is to construct programs in such a way that this difference is irrelevant. Which I do. I stumbled into it in regression testing.

Update2:
Regarding the @VikasDamodar comment. One shouldn't trust the repr() function:

>>> np.float64(5.92270987499999979065)
5.922709875
>>> '%.20f' % np.float64(5.92270987499999979065)
'5.92270987499999979065'

Update3:
Tested on python3.6.0 x32, numpy 1.14.0, win64. Also on python3.6.4 x64, numpy 1.14.0, debian.

Update4:
Just to be sure:

import numpy as np
a = 5.92270987499999979065
print('%.20f' % round(a, 8))
print('%.20f' % round(np.float64(a), 8))

5.92270987000000026512
5.92270988000000020435

Update5:
The following code demonstrates on which stage the difference takes place without using str:

>>> np.float64(a) - 5.922709874
1.000000082740371e-09
>>> a - 5.922709874
1.000000082740371e-09
>>> round(np.float64(a), 8) - 5.922709874
6.000000496442226e-09
>>> round(a, 8) - 5.922709874
-3.999999442783064e-09

Clearly, before applying 'round' they were the same number.

Update6:
In contrast to @user2357112's answer, np.round is roughly 4 times slower than round:

%%timeit a = 5.92270987499999979065
round(a, 8)

1.18 µs ± 26.5 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)  

%%timeit a = np.float64(5.92270987499999979065)
round(a, 8)

4.05 µs ± 43.9 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Also in my opinion np.round did a better job rounding to the nearest even than builtin round: originally I got this 5.92270987499999979065 number through dividing 11.84541975 by two.

Answer 1

tldr

Builtin round and numpy.round use different algorithms for rounding. Their results are identical for most numbers, but substantially different for some corner cases.

Both are good for certain uses.

round is faster for scalars, np.round is faster for arrays.

explanation

  •  Builtin round uses a straightforward approach of checking the decimal digit right after the requested one. It rounds 4 downwards no matter what happens to be after (even if it is ...499999), and rounds 5 to the nearest even unless it has a single 1 thereafter (eg ...500001), in which case it rounds up.

  •  np.round multiplies the number to the requested the power of 10, rounds to the nearest int by ordinary rules, then divides back by the same power of 10.

It gives more predictable results for cases like 0.33/2:

>>> 0.33/2
0.165
>>> round(0.33/2, 2)
0.17
>>> np.round(0.33/2, 2)
0.16

Here 0.165 should've been rounded to the nearest even which is 0.16.

Update:
Yet it suffers from round-off errors for cases like 1.09/2 (as noted by Mark Dickinson in the comment):

>>> 1.09/2
0.545
>>> round(1.09/2, 2)
0.55
>>> np.round(1.09/2, 2)
0.55

The only workaround I could think of is

>>> round(round(1.09*100)/2)/100
0.54

which works but is far from being universal.

Answer 2

float.__round__ takes special care to produce correctly-rounded results, using a correctly-rounded double-to-string algorithm.

NumPy does not. The NumPy docs mention that

Results may also be surprising due to the inexact representation of decimal fractions in the IEEE floating point standard [R9] and errors introduced when scaling by powers of ten.

This is faster, but produces more rounding error. It leads to errors like what you've observed, as well as errors where numbers even more unambiguously below the cutoff still get rounded up:

>>> x = 0.33499999999999996
>>> x
0.33499999999999996
>>> x < 0.335
True
>>> x < Decimal('0.335')
True
>>> x < 0.67/2
True
>>> round(x, 2)
0.33
>>> numpy.round(x, 2)
0.34000000000000002

---

You're getting a slower time for NumPy's rounding, but that doesn't have anything to do with which rounding algorithm is slower. Any time comparison between NumPy and regular Python math will boil down to the fact that NumPy is optimized for whole-array operations. Doing math on single NumPy scalars has a lot of overhead, but rounding an entire array with numpy.round easily beats rounding a list of floats with round:

In [6]: import numpy

In [7]: l = [i/7 for i in range(100)]

In [8]: a = numpy.array(l)

In [9]: %timeit [round(x, 1) for x in l]
59.6 µs ± 408 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

In [10]: %timeit numpy.round(a, 1)
5.27 µs ± 145 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

As for which one is more accurate, that's definitely float.__round__. Your number is closer to 5.92270987 than to 5.92270988, and it's round-ties-to-even, not round-everything-to-even. There's no tie here.

Answer 3

In my particular case a pretty straightforward way of working around the difference between those two functions to get consistent results is through multiplying and dividing.

For my application it seems to be working better than native round giving the same result as np.round:

'%.20f' % (round(a*1e8)/1e8)
'5.92270988000000020435'
'%.20f' % (round(np.float64(a)*1e8)/1e8)
'5.92270988000000020435'

Update
Thanks to @user2357112 I found out that it is exactly what happens in np.round internally (multiarray/calculation.c#L665), so unless you're cross-testing your results between numpy and native python it's safe just to use the numpy version of round without those extra divisions and multiplications on python level.

Update2
When dealing with scalars, this method of division&multiplication on python level is somewhat (~30%) slower than native round but significantly faster (~3 times) than np.round (giving the same result as np.round):

%%timeit c = 11.84541975
round(c/2)
349 ns ± 10.2 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

%%timeit c = 11.84541975
round(c*1e8/2)/1e8
519 ns ± 13 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

%%timeit c = np.float64(11.84541975)
round(c/2)
1.67 µs ± 20.2 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

%%timeit c = np.float64(11.84541975)
round(c*1e8/2)/1e8
2.01 µs ± 37.9 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

Update3
Python's builtin round uses a straightforward approach which engages 'round-to-even' rule only for exactly representable dyadic rationals, like 0.375 (which is an integer divided by an exact power of two), thus effectively substituting this rule for all other numbers with 'round the tie number which happens to have 49999 representation down and happens to end with 50001 up'. I'm not really sure if this algorithm is better or worse, but definitely it is less suitable for checking manually.

Answer 4

And yes, yet another way to deal with such things is to use Decimal which is not that slow in python3:

%%timeit d = D('11.84541975'); q = D('0.00000001')
(d/2).quantize(q)

485 ns ± 10.2 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)