Which algorithm is more accurate for computing the sum of a sorted array of numbers?

Question

Given is an increasing finite sequence of positive numbers $z_{1} ,z_{2},.....z_{n}$. Which of the following two algorithms is better for computing the sum of the numbers?

s=0; 
for \ i=1:n 
    s=s + z_{i} ; 
end

Or:

s=0; 
for \ i=1:n 
s=s + z_{n-i+1} ; 
end

In my opinion it would be better to start adding the numbers from the largest to the smallest number, because the error gets smaller and smaller. We also know that when we add a very large number to a very small number, the approximate result can be the large number.

Is this correct? What else can be said?

Answer 1

Adding arbitrary floating point numbers will usually give some rounding error, and the rounding error will be proportional to the size of the result. If you calculate a single sum and start by adding the largest numbers first, the average result will be larger. So you would start adding with the smallest numbers.

But you get better result (and it runs faster) if you produce four sums, for example: Start with sum1, sum2, sum3, sum4 and add four array elements in turn to sum1, sum2, sum3, sum4. Since each result is on average only 1/4th of the original sum, your error is four times smaller.

Better still: Add the numbers in pairs. Then add the results in pairs. Add those results in pairs again, and so on until you are left with two numbers to add.

Very simple: Use higher precision. Use long double to calculate a sum of doubles. Use double to calculate a sum of floats.

Close to perfect: Look up Kahan's algorithm, described before. Best still used by adding starting with the smallest number.

Answer 2

Are these integers or floating point numbers? Assuming it's floating point, I would go with the first option. It's better to add the smaller numbers to each other, then add the bigger numbers later. With the second option, you'll end up adding a small number to a big number as i increases, which can lead to problems. Here's a good resource on floating point arithmetic: What Every Computer Scientist Should Know About Floating-Point Arithmetic

Answer 3

animal_magic's answer is correct that you should add the numbers from smallest to largest, however I want to give an example to show why.

Assume we are working in a floating point format that gives us a staggering 3 digits of accuracy. Now we want to add ten numbers:

[1000, 1, 1, 1, 1, 1, 1, 1, 1, 1]

Of course the exact answer is 1009, but we can't get that in our 3 digit format. Rounding to 3 digits, the most accurate answer we get get is 1010. If we add smallest to largest, on each loop we get:

Loop Index        s
1                 1
2                 2
3                 3
4                 4
5                 5
6                 6
7                 7
8                 8
9                 9
10                1009 -> 1010

So we get the most accurate answer possible for our format. Now lets assume that we add from largest to smallest.

Loop Index        s
1                 1000
2                 1001 -> 1000
3                 1001 -> 1000
4                 1001 -> 1000
5                 1001 -> 1000
6                 1001 -> 1000
7                 1001 -> 1000
8                 1001 -> 1000
9                 1001 -> 1000
10                1001 -> 1000

Since the floating point numbers are rounded after each operation, all of the additions are rounded away, increasing our error from 1 to 9 from the exact. Now imagine if your set of numbers to add had a 1000, and then a hundred 1's, or a million. Note that to be truly accurate, you would want to sum the smallest two numbers, then resort the result into your set of numbers.

Answer 4

For the general case, I'd use compensated summation (or Kahan summation). Unless the numbers are already sorted, sorting them will be much more expensive than adding them. Compensated summation is also more accurate than sorted summation or naive summation (see the previous link).

As for references, What every programmer should know about floating-point arithmetic covers the basic points in enough detail that someone could read it in 20 (+/- 10) minutes and understand the basics. "What every computer scientist should know about floating-point arithmetic" by Goldberg is the classical reference, but most people I know who recommend that paper haven't read it in detail themselves, because it's around 50 pages (more than that, in some printings), and written in dense prose, so I have trouble recommending that as a first-line reference for people. It is good for a second look at the subject. An encyclopedic reference is Higham's Accuracy and Stability of Numerical Algorithms, which covers this material, as well as the accumulation of numerical errors in many other algorithms; it's also 680 pages, so I wouldn't look at this reference first either.

Answer 5

The previous answers already discuss the matter at large and give sound advice, but there is an additional quirk that I'd like to mention. On most modern architectures, the for loop that you have described would be performed anyway in 80-bit extended precision, which guarantees additional accuracy, since all temporary variables will be put in registers. So you already have some form of safeguard from numerical errors. However, in more complicated loops, the intermediate values will be stored in memory in between the operations, and hence truncated to 64 bits. I guess that

s=0; 
for \ i=1:n 
    printf("Hello World");
    s=s + z_{i} ; 
end

suffices to get lower precision in your summation (!!). So be very careful if you want to printf-debug your code while checking for accuracy.

For the interested, this paper describes a problem in a widely used numerical routine (Lapack's rank-revealing QR factorization) whose debugging and analysis was very tricky precisely because of this issue.

Answer 6

Of the 2 options, adding from smaller to larger will produce less numerical error then adding from larger to smaller.

However, >20 years ago in my "Numerical Methods" class the instructor stated this and it occurred to me that this was still introducing more error than necessary because of the relative difference in value between the accumulator and the values that were being added.

Logically, a preferable solution is to add the 2 smallest numbers in the list, then re-insert the summed value into the sorted list.

To demonstrate it, I worked out an algorithm that could do that efficiently (in space and time) by using the space freed up as elements were removed from the primary array to build a secondary array of the summed values which were inherently ordered since the additions were of the sums of values that were always increasing. On each iteration the "tips" of both arrays are then checked to find the 2 smallest values.

Answer 7

Since you didn't restrict the data-type to be used, to achieve a perfectly accurate result, simply use arbitrary length numbers... in which case the order won't matter. It will be much slower, but obtaining perfection does take time.

Answer 8

Use binary tree addition , ie , Choose the mean of the distribution(closest number) as the root of the binary tree, and create a sorted binary tree by adding lesser values to the left of the graph and larger ones to the right and so on. The add all child nodes of a single parent recursively in a bottom up approach. This will be efficient as the avg error increases with the number of summations and in a binary tree approach, the number of summations are in the order of log n in base 2. Hence the avg error would be lesser.

Answer 9

What Hristo Iliev said above about 64-bit compilers preferring the SSE and AVX instructions over the FPU (AKA NDP) is absolutely true, at least for Microsoft Visual Studio 2013. However, for the double precision floating-point operations I was using I found it actually faster, as well as in theory more accurate, to use the FPU. If it is important to you, I would suggest testing various solutions first, before choosing a final approach.

When working in Java, I very frequently use the arbitrary-precision BigDecimal data type. It is just too easy, and one usually does not notice the speed decrease. Computing the transcendental functions with infinite series and sqrt using Newton's method can take a millisecond or more, but it is doable and quite accurate.

Answer 10

I only left this here https://stackoverflow.com/a/58006104/860099 (when you go there, click to 'show code snippet' and run it by button

It is JavaScript example which clearly shows that sum starting from largest gives bigger error

arr=[9,.6,.1,.1,.1,.1];

sum     =             arr.reduce((a,c)=>a+c,0);  // =  9.999999999999998
sortSum = [...arr].sort().reduce((a,c)=>a+c,0);  // = 10

console.log('sum:     ',sum);
console.log('sortSum:',sortSum);