Live data from Hacker News

Calculating the mean of a list of numbers (2016)

hypothesis.works

71–80 of 96 posts

Re: Calculating the mean of a list of numbers (2016)

#71
post #65
post #53

Earlier quoted context omitted.

No, in fact it's slower than Kahan summation and less accurate. In particular, sorting doesn't help if you sum up a large number of similarly small values.

How is it less accurate?

>> In particular, sorting doesn't help if you sum up a large number of similarly small values.

Elaborating based entirely on my reading of the wikipedia article:

Kahan summation manually allots some extra bits to calculate with. You get extra precision because you computed with more bits of precision. Conceptually, you can do the same thing with integers:

Imagine I'm adding up a bunch of unsigned bytes. Ordinary addition will give me the lowest 8 bits of the sum, which is pretty much worthless. But in the analogue of Kahan summation, I can do better: I allocate 16 bits to hold my running result, and compute an accurate sum of, say, 0x02B7. Then I zero the low-order bits until the high-order bits fit in a byte: the sum after precision loss is 0x02B4. Then I report the sum: it's 0xAD (= 0x02B4 >> 2) with a left shift of 2. This is much better than the ordinary result of 0xB7.

(In that example, I needed to return an extra value, the left shift of 2. But floating-point numbers already include an indication of their order of magnitude, which integers don't, so in Kahan summation the extra return value is not necessary.)

So the quick answer to "why is Kahan summation more accurate" is that it's more accurate because you used more accuracy. Sort-and-add doesn't give you any more precision than usual, but it will lower your accumulated error.

The root issue that both these strategies attempt to mitigate is that when you add a large value to a small value, you end up losing precision from the small value. Sort-and-add will address that problem as to any two elements of the input -- if you have a list of values of which some are very large and some are very small, sort-and-add will sum up all of the small numbers into one large number before adding that sum to another large number. This reduces the number of times you would add a large number to a small number, compared to summing the input in a random order.

Kahan summation means you do all your addition at a higher level of precision. As such, it also addresses the (large+small) problem that occurs when an element of the input is small in comparison to the sum of all the input values, as opposed to being small in comparison to the maximum individual input value.

If your input consists of a medium amount of small numbers and a few large numbers, such that the sum of the input is well approximated by each of the large numbers within the input, sort-and-add should work fine.

If your input contains a large amount of small numbers, sort-and-add will fail: by the time your running sum has become large compared to each element of the input, every further addition will suffer from loss of precision -- you are adding a small number (the next input element) to a large one (the running sum).

If your input contains many large numbers, sort-and-add will fail for exactly the same reason it fails on many small numbers: each element individually is small compared to the total.

Re: Calculating the mean of a list of numbers (2016)

#73
post #65

Earlier quoted context omitted.

How is it less accurate?

>> In particular, sorting doesn't help if you sum up a large number of similarly small values. Elaborating based entirely on my reading of the wikipedia article: Kahan summation manually allots some extra bits to calculate with. You get extra precision because you computed with more bits of precision. Conceptually, you can do the same thing with integers: Imagine I'm adding up a bunch of unsigned bytes. Ordinary addi…

Thanks for thoughtful reply.

But the Kahan method aside, the idea would be to start with the smallest numbers, not the largest.

Re: Calculating the mean of a list of numbers (2016)

#74
post #73

Earlier quoted context omitted.

>> In particular, sorting doesn't help if you sum up a large number of similarly small values. Elaborating based entirely on my reading of the wikipedia article: Kahan summation manually allots some extra bits to calculate with. You get extra precision because you computed with more bits of precision. Conceptually, you can do the same thing with integers: Imagine I'm adding up a bunch of unsigned bytes. Ordinary addi…

Thanks for thoughtful reply. But the Kahan method aside, the idea would be to start with the smallest numbers, not the largest.

Yes, I understand that.

>> if you have a list of values of which some are very large and some are very small, sort-and-add will sum up all of the small numbers into one large number before adding that sum to another large number

I don't quite follow the point you're making -- can you elaborate?

Re: Calculating the mean of a list of numbers (2016)

#75
post #31

Earlier quoted context omitted.

Similar correct result from sbcl.

The pedantic point of the article is that an error in this situation isn't correct. The average is 8e307 and 8e307 is 8e307.

Agreed. My limited point was that detection of out-of-range is available out of the box in some languages.

Re: Calculating the mean of a list of numbers (2016)

#76
Here is a fascinating post by Stefan Karpinski, one of the creators of Julia: https://discourse.julialang.org/t/array-ordering-and-naive-s...

He shows off a function named "sumsto" which takes a single positive double precision floating point number "x" as an argument.

"sumsto" always returns the same vector of 2046 floating point numbers -- just in an order so that naive left-to-right summation returns "x". Just changing the order of that vector lets you sum to almost any positive double precision number.

If you want to run the code, I didn't see where "realmax" was defined, but this works:

realmax() = prevfloat(typemax(Float64))

Re: Calculating the mean of a list of numbers (2016)

#77
post #76

Here is a fascinating post by Stefan Karpinski, one of the creators of Julia: https://discourse.julialang.org/t/array-ordering-and-naive-s... He shows off a function named "sumsto" which takes a single positive double precision floating point number "x" as an argument. "sumsto" always returns the same vector of 2046 floating point numbers -- just in an order so that naive left-to-right summation returns "x". Just cha…

You think that's bad? Theoretical math has "conditionally convergent series'" [1] Sum in order and it can convergence to a given value. Rearrange the terms and any real number is possible.

I have a hunch there could be a bit of relation between these things.

[1] https://en.wikipedia.org/wiki/Conditional_convergence

Re: Calculating the mean of a list of numbers (2016)

#78
post #39
post #22

What about a "reduce" technique? Average the numbers in equal-sized chunks, then average those averages. You could even chunk the chunk averages and repeat the process as many levels down as you want to, and chunks could be as small as 2 each. I guess this still assumes that the largest number in the original list is less than or equal to the maximum floating point value, but otherwise you stay roughly in the same sp…

I think this works perfectly if your list has 2^n elements. Otherwise, you have to resort to multiplying by imprecise fractions.

You would only need one special case: if a list (top level or not) has 2n+1 numbers, weight the last one differently.

Re: Calculating the mean of a list of numbers (2016)

#80
post #49
post #35

Earlier quoted context omitted.

Most of the issues discussed in the article are not from limits in precision, they're from overflows that arise from limits in the exponential range. Finding the average of 8e+307 and 8e+307 is a "low precision" problem, and even naive methods to avoid overflow will not hit limits on precision in this case (e.g. 8e+307/2 + 8e+307/2). You're right that there are issues with precision and floating point math (the artic…

Sorry, I should’ve been more clear: when I said “precision” I meant any of several things related to floating point arithmetic including, but not limited to: actual numerical precision of operations, non-associativity of operations, overflows (which affects the latter case as well), etc. (Which is why I mention dynamic range in the GP.) The point is that there is an “unwritten contract” between the user and the libra…

> using arbitrary precision fractions for Gaussian elimination

Here's a way to get exact results from just standard precision, say, integers 32 bits long:

For Gauss elimination, that's for solving a system of linear equations, say, m equations in n unknowns. So, we are given, in floating point, a matrix m x n A, an unknown vector 1 x n x, and a right side constant 1 x m b. So we are looking for x so that

Ax = b

Now, multiply each equation by an integer, say, a power of 2, so that all the numbers in A and b are integers.

Next get a list of positive prime numbers, each, say, stored in an integer 32 bits long. Say we have p such prime numbers; so we have prime number for i = 1, 2, ..., p.

Next for each i, solve Ax = b by using Gauss elimination but in the field of integers modulo prime number i. For division, use the Euclidean greatest common divisor algorithm. Yes for this arithmetic have to be able to form the 64 bit sum or product of two 32 bit whole numbers and then divide by 32 bit integer and keep the 32 bit remainder -- commonly we write a little assembler routine for this.

After doing that work p times (could be done in parallel), we use the Chinese remainder theorem to put together the rational numbers that are quotients of multi-precision whole numbers of the solution. With those, we can get floating point approximations as accurate as we please.

But if want to work in floating point arithmetic anyway, then there are three ways to do better:

(1) To start, in Gauss elimination, look down column 1 of A, find the row with the number of largest absolute value, and swap rows to put that number in row 1 and column 1. More generally after p - 1 rows, will want a pivot element in row p and column p. By swapping rows, use the number largest in absolute value in column p and rows p, p + 1, ..., m.

(2) When form an inner product, say,

u(1)v(1) + u(2)v(2) + ... + u(q)v(q)

form each product in double precision and add the double precision numbers. If want to do a little better, then before adding, sort the numbers so that are adding the smaller numbers first. Or

"The main sin in numerical analysis is subtracting two large numbers whose difference is small.", etc.

(3) Once have x so that Ax ~ b, find difference d = b - Ax and then solve for e in Ae = d and replace x by x + e so that A(x + e) = Ax + Ae = (b - d) + d = b. So, take x + e as the solution. After Gauss elimination, solving Ae = d can go relatively quickly. Can do this step more than once.

Post reply on HN