Feature suggestion: sum() ought to use a compensated summation algorithm

[Szabolcs Horvát]

I did the following calculation: Generated a list of a million random
numbers between 0 and 1, constructed a new list by subtracting the mean
value from each number, and then calculated the mean again.

The result should be 0, but of course it will differ from 0 slightly
because of rounding errors.

However, I noticed that the simple Python program below gives a result
of ~ 10^-14, while an equivalent Mathematica program (also using double
precision) gives a result of ~ 10^-17, i.e. three orders of magnitude
more precise.

Here's the program (pardon my style, I'm a newbie/occasional user):

from random import random

data = [random() for x in xrange(1000000)]

mean = sum(data)/len(data)
print sum(x - mean for x in data)/len(data)

A little research shows that Mathematica uses a "compensated summation"
algorithm. Indeed, using the algorithm described at
http://en.wikipedia.org/wiki/Kahan_summation_algorithm
gives us a result around ~ 10^-17:

def compSum(arr):
s = 0.0
c = 0.0
for x in arr:
y = x-c
t = s+y
c = (t-s) - y
s = t
return s

mean = compSum(data)/len(data)
print compSum(x - mean for x in data)/len(data)


I thought that it would be very nice if the built-in sum() function used
this algorithm by default. Has this been brought up before? Would this
have any disadvantages (apart from a slight performance impact, but
Python is a high-level language anyway ...)?

Szabolcs Horvát

 

[Arnaud Delobelle]

[...]

A little research shows that Mathematica uses a "compensated
summation" algorithm. Indeed, using the algorithm described at
http://en.wikipedia.org/wiki/Kahan_summation_algorithm
gives us a result around ~ 10^-17:

def compSum(arr):
s = 0.0
c = 0.0
for x in arr:
y = x-c
t = s+y
c = (t-s) - y
s = t
return s

mean = compSum(data)/len(data)
print compSum(x - mean for x in data)/len(data)


I thought that it would be very nice if the built-in sum() function
used this algorithm by default. Has this been brought up before?
Would this have any disadvantages (apart from a slight performance
impact, but Python is a high-level language anyway ...)?

Szabolcs Horvát

sum() works for any sequence of objects with an __add__ method, not
just floats! Your algorithm is specific to floats.

 

[Ivan Illarionov]

I did the following calculation: Generated a list of a million random
numbers between 0 and 1, constructed a new list by subtracting the mean
value from each number, and then calculated the mean again.

The result should be 0, but of course it will differ from 0 slightly
because of rounding errors.

However, I noticed that the simple Python program below gives a result
of ~ 10^-14, while an equivalent Mathematica program (also using double
precision) gives a result of ~ 10^-17, i.e. three orders of magnitude
more precise.

Here's the program (pardon my style, I'm a newbie/occasional user):

from random import random

data = [random() for x in xrange(1000000)]

mean = sum(data)/len(data)
print sum(x - mean for x in data)/len(data)

A little research shows that Mathematica uses a "compensated summation"
algorithm. Indeed, using the algorithm described at
http://en.wikipedia.org/wiki/Kahan_summation_algorithm gives us a result
around ~ 10^-17:

def compSum(arr):
s = 0.0
c = 0.0
for x in arr:
y = x-c
t = s+y
c = (t-s) - y
s = t
return s

mean = compSum(data)/len(data)
print compSum(x - mean for x in data)/len(data)


I thought that it would be very nice if the built-in sum() function used
this algorithm by default. Has this been brought up before? Would this
have any disadvantages (apart from a slight performance impact, but
Python is a high-level language anyway ...)?

Szabolcs Horvát

Built-in sum should work with everything, not just floats. I think it
would be useful addition to standard math module.

If you know C you could write a patch to mathmodule.c and submit it to
Python devs.

 

[Szabolcs Horvát]

sum() works for any sequence of objects with an __add__ method, not
just floats! Your algorithm is specific to floats.

This occurred to me also, but then I tried

sum(['abc', 'efg'], '')

and it did not work. Or is this just a special exception to prevent the
misuse of sum to join strings? (As I said, I'm only an occasional user.)

Generally, sum() seems to be most useful for numeric types (i.e. those
that form a group with respect to __add__ and __neg__/__sub__), which
may be either exact (e.g. integers) or inexact (e.g. floating point
types). For exact types it does not make sense to use compensated
summation (though it wouldn't give an incorrect answer, either), and
sum() cannot decide whether a user-defined type is exact or inexact.

But I guess that it would still be possible to make sum() use
compensated summation for built-in floating point types (float/complex).

Or, to go further, provide a switch to choose between the two methods,
and make use compensated summation for float/complex by default. (But
perhaps people would consider this last alternative a bit too messy.)

(Just some thoughts ...)

 

[hdante]

sum() works for any sequence of objects with an __add__ method, not
just floats!  Your algorithm is specific to floats.

This occurred to me also, but then I tried

sum(['abc', 'efg'], '')

and it did not work.  Or is this just a special exception to prevent the
  misuse of sum to join strings?  (As I said, I'm only an occasional user.)

Generally, sum() seems to be most useful for numeric types (i.e. those
that form a group with respect to __add__ and __neg__/__sub__), which
may be either exact (e.g. integers) or inexact (e.g. floating point
types).  For exact types it does not make sense to use compensated
summation (though it wouldn't give an incorrect answer, either), and
sum() cannot decide whether a user-defined type is exact or inexact.

But I guess that it would still be possible to make sum() use
compensated summation for built-in floating point types (float/complex).

Or, to go further, provide a switch to choose between the two methods,
and make use compensated summation for float/complex by default.  (But
perhaps people would consider this last alternative a bit too messy.)

(Just some thoughts ...)

The benefits should be weighted considering the common case. For
example, do you find an error of 10^-14 unbearable ? How many times
someone will add 1.000.000 numbers in a typical scientific application
written in python ?

I believe that moving this to third party could be better. What about
numpy ? Doesn't it already have something similar ?

 

[Arnaud Delobelle]

sum() works for any sequence of objects with an __add__ method, not
just floats! Your algorithm is specific to floats.

This occurred to me also, but then I tried

sum(['abc', 'efg'], '')

and it did not work. Or is this just a special exception to prevent
the misuse of sum to join strings? (As I said, I'm only an occasional
user.)

I think that's right: anything with an __add__ method, apart from
string, can be sum()ed.

 

[Ivan Illarionov]

sum() works for any sequence of objects with an __add__ method, not
just floats! Your algorithm is specific to floats.

This occurred to me also, but then I tried

sum(['abc', 'efg'], '')

Interesting, I always thought that sum is like shortcut of
reduce(operator.add, ...), but I was mistaken.

reduce() is more forgiving:
reduce(operator.add, ['abc', 'efg'], '' ) # it works
'abcefg'

 

[sturlamolden]

I believe that moving this to third party could be better. What about
numpy ? Doesn't it already have something similar ?

Yes, Kahan summation makes sence for numpy arrays. But the problem
with this algorithm is optimizing compilers. The programmer will be
forced to use tricks like inline assembly to get around the optimizer.
If not, the optimizer would remove the desired features of the
algorithm. But then we have a serious portability problem.

 

[Thomas Dybdahl Ahle]

sum() works for any sequence of objects with an __add__ method, not
just floats! Your algorithm is specific to floats.

This occurred to me also, but then I tried

sum(['abc', 'efg'], '')

Interesting, I always thought that sum is like shortcut of
reduce(operator.add, ...), but I was mistaken.

reduce() is more forgiving:
reduce(operator.add, ['abc', 'efg'], '' ) # it works
'abcefg'

Hm, it works for lists:
sum(([1], [2]), [])
[1, 2]

However I find the seccond argument hack ugly.
Does the sum way have any performance advantages over the reduce way?

 

[Ivan Illarionov]

Arnaud Delobelle wrote:

sum() works for any sequence of objects with an __add__ method, not
just floats! Your algorithm is specific to floats.

This occurred to me also, but then I tried

sum(['abc', 'efg'], '')

Interesting, I always thought that sum is like shortcut of
reduce(operator.add, ...), but I was mistaken.

reduce() is more forgiving:
reduce(operator.add, ['abc', 'efg'], '' ) # it works 'abcefg'

Hm, it works for lists:
sum(([1], [2]), [])
[1, 2]

However I find the seccond argument hack ugly. Does the sum way have any
performance advantages over the reduce way?

Yes, sum() is faster:

$ python -m timeit "" "sum([[1], [2], [3, 4]], [])"
100000 loops, best of 3: 6.16 usec per loop

$ python -m timeit "import operator" \
"reduce(operator.add, [[1], [2], [3, 4]])"
100000 loops, best of 3: 11.9 usec per loop

 

[George Sakkis]

sum() works for any sequence of objects with an __add__ method, not
just floats!  Your algorithm is specific to floats.
This occurred to me also, but then I tried
sum(['abc', 'efg'], '')
Interesting, I always thought that sum is like shortcut of
reduce(operator.add, ...), but I was mistaken.
reduce() is more forgiving:
reduce(operator.add, ['abc', 'efg'], '' ) # it works 'abcefg'

Hm, it works for lists:
sum(([1], [2]), [])
[1, 2]

So it's not reduce() that is more forgiving, it's sum() that makes an
exception for strings only.

However I find the seccond argument hack ugly. Does the sum way have any
performance advantages over the reduce way?

Yes, sum() is faster:

$ python -m timeit "" "sum([[1], [2], [3, 4]], [])"
100000 loops, best of 3: 6.16 usec per loop

$ python -m timeit "import operator" \
"reduce(operator.add, [[1], [2], [3, 4]])"
100000 loops, best of 3: 11.9 usec per loop

Not really; you measure the import and the attribute access time in
the second case. sum() is 2x faster for adding numbers only:

# Adding integers
python -mtimeit --setup="x=[1]*100" "sum(x)"
100000 loops, best of 3: 4.87 usec per loop
python -mtimeit --setup="from operator import add; x=[1]*100"
"reduce(add,x)"
100000 loops, best of 3: 10.1 usec per loop

# Adding floats
python -mtimeit --setup="x=[1.0]*100" "sum(x)"
100000 loops, best of 3: 5.14 usec per loop
python -mtimeit --setup="from operator import add; x=[1.0]*100"
"reduce(add,x)"
100000 loops, best of 3: 10.1 usec per loop

# Adding tuples
python -mtimeit --setup="x=[(1,)]*100" "sum(x,())"
10000 loops, best of 3: 61.6 usec per loop
python -mtimeit --setup="from operator import add; x=[(1,)]*100"
"reduce(add,x,())"
10000 loops, best of 3: 66.3 usec per loop

# Adding lists
python -mtimeit --setup="x=[[1]]*100" "sum(x,[])"
10000 loops, best of 3: 79.9 usec per loop
python -mtimeit --setup="from operator import add; x=[[1]]*100"
"reduce(add,x,[])"
10000 loops, best of 3: 84.3 usec per loop

George

 

[Ivan Illarionov]

44:19 +0200, Szabolcs Horvát wrote:

Arnaud Delobelle wrote:

sum() works for any sequence of objects with an __add__ method,
not just floats!  Your algorithm is specific to floats.

This occurred to me also, but then I tried

sum(['abc', 'efg'], '')

Interesting, I always thought that sum is like shortcut of
reduce(operator.add, ...), but I was mistaken.

reduce() is more forgiving:
reduce(operator.add, ['abc', 'efg'], '' ) # it works 'abcefg'
Hm, it works for lists:
sum(([1], [2]), [])
[1, 2]

So it's not reduce() that is more forgiving, it's sum() that makes an
exception for strings only.

However I find the seccond argument hack ugly. Does the sum way have
any performance advantages over the reduce way?

Yes, sum() is faster:

$ python -m timeit "" "sum([[1], [2], [3, 4]], [])" 100000 loops, best
of 3: 6.16 usec per loop

$ python -m timeit "import operator" \ "reduce(operator.add, [[1], [2],
[3, 4]])" 100000 loops, best of 3: 11.9 usec per loop

Not really; you measure the import and the attribute access time in the
second case. sum() is 2x faster for adding numbers only:

# Adding integers
python -mtimeit --setup="x=[1]*100" "sum(x)" 100000 loops, best of 3:
4.87 usec per loop python -mtimeit --setup="from operator import add;
x=[1]*100" "reduce(add,x)"
100000 loops, best of 3: 10.1 usec per loop

# Adding floats
python -mtimeit --setup="x=[1.0]*100" "sum(x)" 100000 loops, best of 3:
5.14 usec per loop python -mtimeit --setup="from operator import add;
x=[1.0]*100" "reduce(add,x)"
100000 loops, best of 3: 10.1 usec per loop

# Adding tuples
python -mtimeit --setup="x=[(1,)]*100" "sum(x,())" 10000 loops, best of
3: 61.6 usec per loop python -mtimeit --setup="from operator import add;
x=[(1,)]*100" "reduce(add,x,())"
10000 loops, best of 3: 66.3 usec per loop

# Adding lists
python -mtimeit --setup="x=[[1]]*100" "sum(x,[])" 10000 loops, best of
3: 79.9 usec per loop python -mtimeit --setup="from operator import add;
x=[[1]]*100" "reduce(add,x,[])"
10000 loops, best of 3: 84.3 usec per loop

George

Thanks for correction. Forget about --setup.

 

[hdante]

Yes, Kahan summation makes sence for numpy arrays. But the problem
with this algorithm is optimizing compilers. The programmer will be

No, optimizing compilers shouldn't discard floating point operations
by default, since it changes program behavior. If they do, at least
there should be a flag to turn it off.

 

[Gabriel Genellina]

There's a thread you might find interesting:

http://groups.google.com/group/comp...read/thread/9eda29faf92f532e/027cef7d4429aa3a

In that thread I suggested that Python ought to implement sum by adding
together each pair of values, then each pair of results and so on. This
means that for input values of a similar magnitude the intermediate results
stay balanced. The implementation doesn't actually do all the first level

Python doesn't require __add__ to be associative, so this should not be used as a general sum replacement. But if you know that you're adding floating point numbers you can use whatever algorithm best fits you. Or use numpy arrays; I think they implement Kahan summation or a similar algorithm.

 

[Szabolcs Horvát]

There's a thread you might find interesting:

http://groups.google.com/group/comp...read/thread/9eda29faf92f532e/027cef7d4429aa3a

In that thread I suggested that Python ought to implement sum by adding
together each pair of values, then each pair of results and so on. This
means that for input values of a similar magnitude the intermediate results
stay balanced. The implementation doesn't actually do all the first level
sums first, it builds up a stack containing the sum of 2^k, 2^(k-1)...2
values. Also it works for other types as well, and for strings or lists has
the advantage of minimising the number of large string or list concatenations.

If you follow that thread far enough you'll find a post by Jussi Salmela
which shows that summing adjacent pairs performs as well as Kahan summation
(he says 'almost as good' but in fact the only time they get different
results the 'correct' answer is 500000.000005 and Kahan and sumpairs get
the two nearest representable results either side: 500000.00000499998 and
500000.00000500004 respectively).

I never did get round to re-coding my sumpairs() function in C, I would be
interested to see just how much overhead it introduces (and I suspect it
would be faster than the existing sum in some cases such as for lists).

Thanks for this link! Though some of the thread is missing, it was an
interesting read.

This sounds like a very good idea: it is more generally applicable than
the Kahan summation because it does not require subtraction/negation, so
it would work with user-defined types, too.

While now I am convinced that using Kahan summation by default is not
such a good idea, I cannot see any serious disadvantages/problems with
pairwise summation!

 

[Szabolcs Horvát]

Python doesn't require __add__ to be associative, so this should not be used as a general sum replacement.

It does not _require_ this, but using an __add__ that is not commutative
and associative, or has side effects, would qualify as a serious misuse,
anyway. So I think that this isn't a real disadvantage (it can always
be documented that sum() expects __add__ to have these properties).

But you are right that summing floats with a requirement of high
precision is a quite specific use case, so there are no serious
arguments to do this, either.

 

[Raymond Hettinger]

I did the following calculation:  Generated a list of a million random
numbers between 0 and 1, constructed a new list by subtracting the mean
value from each number, and then calculated the mean again.

The result should be 0, but of course it will differ from 0 slightly
because of rounding errors.
See:
http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/393090



Here's the program (pardon my style, I'm a newbie/occasional user): . . .
mean = sum(data)/len(data)
print sum(x - mean for x in data)/len(data)

See:
http://svn.python.org/view/sandbox/trunk/statistics/statistics.py?rev=40981&view=markup


Raymond

 

[Raymond Hettinger]

However I find the seccond argument hack ugly. Does the sum way have any

....

Not really; you measure the import and the attribute access time in
the second case. sum() is 2x faster for adding numbers only:

Try it with Py2.6 or Py3.0. The sum() function has been optimized and
should be at least 6x faster for summing ints and floats. It should
even beat psyco optimized sums!


Raymond

 

[Szabolcs]

It does not _require_ this, but using an __add__ that is not commutative
and associative, or has side effects, would qualify as a serious misuse,
anyway.

Sorry, I meant associative only, not commutative.

 

[Szabolcs]

Unfortunately an __add__ which is not associative is actually perfectly
reasonable.

Well, 'reasonable' is a subjective notion in this case I have
never seen a non-associative use of the plus sign in maths, so I would
tend to take associativity for granted when seeing a + b + c in code,
therefore I would avoid such uses. (Multiplication signs are a
different story.) Of course others may feel differently about this,
and technically there's nothing in the way of using a non-associative
__add__!