{"id":4675,"date":"2021-09-07T08:13:33","date_gmt":"2021-09-07T12:13:33","guid":{"rendered":"https:\/\/blogs.mathworks.com\/loren\/?p=4675"},"modified":"2021-09-07T13:27:45","modified_gmt":"2021-09-07T17:27:45","slug":"a-faster-and-more-accurate-sum","status":"publish","type":"post","link":"https:\/\/blogs.mathworks.com\/loren\/2021\/09\/07\/a-faster-and-more-accurate-sum\/","title":{"rendered":"A faster and more accurate sum"},"content":{"rendered":"
Today's guest blogger is Christine Tobler, who's a developer at MathWorks working on core numeric functions.<\/span><\/div>
Hi everyone! I'd like to tell you a story about round-off error, the algorithm used in <\/span>sum<\/span><\/a>, and compability issues. Over the last few years, my colleague Bobby Cheng has made changes to <\/span>sum<\/span> to make it both more accurate and faster, and I thought it would be an interesting story to tell here.<\/span><\/div>
Table of Contents<\/span><\/div>
Numerical issues in sum\r\n<\/span><\/a>Different ways of computing the sum\r\n<\/span><\/a>Making changes to sum\r\n<\/span><\/a>References about accurately computing the sum\r\n<\/span><\/a>Helper Functions<\/span><\/a><\/div><\/div>

Numerical issues in sum<\/span><\/h2>
Even for a simple function like <\/span>sum<\/span><\/a>, with floating-point numbers the order in which we sum them up matters:<\/span><\/div>
x = (1 + 1e-16) + 1e-16<\/span><\/span><\/div>
x = 1<\/div><\/div><\/div>
y = 1 + (1e-16 + 1e-16)<\/span><\/span><\/div>
y = 1.0000<\/div><\/div><\/div>
x - y<\/span><\/span><\/div>
ans = -2.2204e-16<\/div><\/div><\/div><\/div>
So the result of the <\/span>sum<\/span> command depends on the order in which the inputs are added up. The <\/span>sum<\/span> command in MATLAB started out with the most straight-forward implementation: Start at the first number and add them up one at a time. I have implemented that in <\/span>sumOriginal<\/span> at the bottom of this post.<\/span><\/div>
But we ran into a problem: For single precision, the results were sometimes quite wrong:<\/span><\/div>
sumOfOnes = sumOriginal(ones(1e8, 1, <\/span>'single'<\/span>))<\/span><\/span><\/div>
sumOfOnes = single\r\n 16777216<\/div><\/div><\/div><\/div>
What's going on here? The issue is that for this number, the round-off error is larger than 1, as we can see by calling <\/span>eps<\/span><\/a><\/div>
eps(sumOfOnes)<\/span><\/span><\/div>
ans = single\r\n 2<\/div><\/div><\/div><\/div>
Therefore, adding another 1 and then rounding to single precision returns the same number again, as the exact result is rounded down to fit into single precision:<\/span><\/div>
sumOfOnes + 1<\/span><\/span><\/div>
ans = single\r\n 16777216<\/div><\/div><\/div><\/div>
This was particularly noticeable when <\/span>sum<\/span> was used to compute the mean of a large array, since the mean computed could be smaller than all the elements of the input array.<\/span><\/div>

Different ways of computing the sum<\/span><\/h2>
Now we already had a threaded version of sum in MATLAB, where we would compute the sum of several chunks of an array in parallel, and then add those up in the end. It turned out that this version didn't have the same issues:<\/span><\/div>
numBlocks = 8;<\/span><\/span><\/div><\/div>
sumInBlocks(ones(1e8, 1, <\/span>'single'<\/span>), numBlocks)<\/span><\/span><\/div>
ans = single\r\n 100000000<\/div><\/div><\/div><\/div>
We have made this change in MATLAB's <\/span>sum<\/span> even for the non-threaded case, and it addresses this and other similar cases. But keep in mind that while this is an improvement for many cases, it's not a perfect algorithm and still won't always give the correct answer. We could modify the number of blocks that we split the input into, and not all work out great:<\/span><\/div>
sumInBlocks(ones(1e8, 1, <\/span>'single'<\/span>), 4)<\/span><\/span><\/div>
ans = single\r\n 67108864<\/div><\/div><\/div>
sumInBlocks(ones(1e8, 1, <\/span>'single'<\/span>), 128)<\/span><\/span><\/div>
ans = single\r\n 99999832<\/div><\/div><\/div><\/div>
Let's look at another case where we're not just summing up the number one (still summing up the same number everytime, because it makes it easier to compare with the exact value). We'll be looking at the relative error of the result here, which is more relevant than if we get the exact integer correct:<\/span><\/div>
x = repmat(single(3155), 54194, 1);<\/span><\/span><\/div><\/div>
exactSum = 3155*54194;<\/span><\/span><\/div><\/div>
numBlocks = 2.^(0:9)<\/span><\/span><\/div>
numBlocks = <\/span>1\u00d710<\/span><\/div><\/div>
1 2 4 8 16 32 64 128 256 512\r\n<\/div>
<\/div>
<\/div><\/div><\/div><\/div><\/div><\/div>
err = zeros(1, length(numBlocks));<\/span><\/span><\/div><\/div>
for <\/span>i=1:length(numBlocks)<\/span><\/span><\/div><\/div>
err(i) = abs(sumInBlocks(x, numBlocks(i)) - exactSum) \/ exactSum;<\/span><\/span><\/div><\/div>
end<\/span><\/span><\/div><\/div>
loglog(numBlocks, err)<\/span><\/span><\/div><\/div>
xlabel(<\/span>'Number of blocks'<\/span>)<\/span><\/span><\/div><\/div>
ylabel(<\/span>'Relative error in sumInBlocks function'<\/span>)<\/span><\/span><\/div>
<\/div><\/div><\/div><\/div><\/div>
So there's definitely a balancing act for choosing the exact number of blocks, with the plot above just representing one dataset.<\/span><\/div>
There are other possible algorithms for computing the sum, I'm including some links at the bottom. Some more complicated ones would result in a slowdown in computation time, which would put a burden on cases that aren't experiencing the issues with accuracy described above.<\/span><\/div>
We chose this \"computing by blocks\" algorithm because the change could be made in a such a way that sum became faster. This was done using loop-unrolling: Instead of adding one number at a time, several numbers are added to separate running totals at the same time. I've included an implementation in <\/span>sumLoopUnrolled<\/span> below.<\/span><\/div>
By the way: I'm not explicitly giving the choice of <\/span>numBlocks<\/span> we made, because we wouldn't want anyone to rely on this - it might change again in the future, after all.<\/span><\/div>

Making changes to <\/span>sum<\/span><\/h2>
In R2017b, we first shipped this new version of <\/span>sum<\/span>, for the case of single-precision numbers only. All practical issues with results being off by magnitudes for large arrays of single data had been addressed, and the function even got faster at the same time. However, we also got quite some feedback from people who were unhappy about the behavior change; while the new behavior gives as good or better results, their code was relying on the old behavior and adjusting to the new behavior was painful.<\/span><\/div>
While we don't want to get stuck in never being able to improve our functions for fear of breaking a dependency, we definitely want to make the process of updating as pain-free as possible. Generally speaking, we aim for run-to-run reproducibility: If a function is called twice with the same inputs, the outputs will be the same. However, if the machine, OS or MATLAB version (among other externals) change, this can also change the exact value that's returned by MATLAB, within round-off error. It's common for the output of matrix multiplication to change for performance' sake, for example. But with <\/span>sum<\/span>, no changes had been made for a long time, and so more people had come to rely on its exact behavior than we would have expected.<\/span><\/div>
So after making this change to <\/span>sum<\/span> for single values, we waited a few releases to evaluate customer feedback on the change. In R2020b, we added the same behavior for all other datatypes. This time, we added a <\/span>Release Note<\/span><\/a> that describes both the performance improvement and mentions the change as a \"Compatibility Consideration\". While we still had some feedback on this being an inconvenient change, it was less than in the first case.<\/span><\/div>
Have you had changes in round-off error with a new MATLAB release break some of your code? Do you find the compatibility considerations in the release notes useful? Please let us know <\/span>in the comments<\/span><\/a>.<\/span><\/div>

References about accurately computing the sum<\/span><\/h2>