How exactly does MATLAB calculate a sum of array elements by its sum() function? Does it use any compensated summation algorithm such as Kahan?

I'm currently working on a software project translating some MATLAB code to Java. During this process, I recognized, that the output of MATLAB's sum() function applied to some 'single' or 'double' array is not identical to a plain addition of consecutive array elements in a loop. I did not find any further information regarding the concrete implementation of sum() in any recent MATLAB version (e.g. for me R2019a), so I tried to replicate its functionality by implementing some common summation algorithms for floating-point addition such as the Kahan summation algorithm or some of its variants. However, the results did not match with the output of MATLAB's sum().

Does anyone know any details about the implementation of MATLAB's sum() function? Is any kind of compensated summation included? Are 'single' arrays cast to 'double' before the summation? Is the summation distributed across multiple threads per default? Are arrays sorted to reduce errors due to the addition of tiny and big floating-point numbers?

Answer 1

[Steven Lord]

You might find this post on Loren Shure's blog interesting and informative.

[Jonas K.]

Thanks a lot, very interesting and helpful information about the development of sum() in this blog post!

Answer 2

[dpb]

TMW does not document publicly algorithms beyond anything provided in the documentation Description section or, occasionally an Algorithms section may add some additional insight.

I've not poked around to investigate, the following thread has some Info https://www.mathworks.com/matlabcentral/answers/550-compensated-summation-in-sum?s_tid=answers_rc1-1_p1_MLT#answer_822 Of course, as noted there, while not likely to have changed, there's no guarantee TMW hasn't changed any heuristic rules.

[Jonas K.]

Thank you, unfortunately I'm already aware of the thread you noted, but it would be great to have any information about a current implementation. Especially, since TMW seems to have changed the multi-threading capabilities since then, as Bruno Luong pointed out, which is also in agreement with what I've tested so far.

Answer 3

[Edric Ellis]

The outtype parameter to sum controls the numeric type used for summation, described in the doc. The default is that sum on single values is performed in single, but other types are operated on in double. (The examples below do rely on the order of operations, which is not guaranteed)

singles = realmax('single') .* [1, 1, -1, -1]

singles = 1×4

1.0e+38 * 3.4028 3.4028 -3.4028 -3.4028

% Saturates
sum(singles)

ans = single Inf

% In 'double', doesn't saturate
sum(singles, 'double')

ans = 0

int8s = [intmax('int8'), intmax('int8'), intmin('int8')]

int8s = 1×3

127 127 -128

% Default summation in 'double' - doesn't saturate
sum(int8s)

ans = 126

% Saturates
sum(int8s, 'native')

ans = int8 -1

[Jonas K.]

Thank you, the information that 'single' arrays are summed up per default using their native type is definitely helpful, although any details about the concrete way / algorithm how the summation is performed, would be even more helpful.

Answer 4

[Bruno Luong]

According to my test it seems MATLAB does not sum by chunk when operating on vector, to ensure the result is consistent, i.e. not depending on number of threads.

>> a=rand(1,1e7);
>> maxNumCompThreads=1;
>> tic; s1=sum(a), toc

s1 =

   4.9999e+06

Elapsed time is 0.005212 seconds.
>> maxNumCompThreads=4;
>> tic; s4=sum(a), toc

s4 =

   4.9999e+06

Elapsed time is 0.005153 seconds.
>> s1-s4

ans =

     0
>> ss=sum(sort(a));
>> ss-s1

ans =

   3.7253e-09

IMO opinion the sum just carried out linearly from left to right with some internal internat result with fiw number of bits > 64.

IIRC, in some version (2015?) MATLAB implements a multi-thread on vector and that raises some questions and that has been discussed on the old newsgroup, then they switched back to single thread.

The multi-thread is used for sum on 2D or ND-array, where each thread is in charge a set of vectors.

All that is hypothetic as TMW does not document the algorirthm.

[Jonas K.]

Thank you, this change regarding multi-threading is in agreement with what I've read and tested so far. I also guess, that more than 64 bits are used internally. However, a concrete number of bits and any details about any used algorithm would be great, since I've already tested some implementation using a 128 bit representation and linear summation, which did not lead to equivalent results either.

[Bruno Luong]

Have you tried the x87 80-bit extended format which is historical the base of Intel FPU?

You migh also take a look at rounding mode _MCW_RC

[Bruno Luong]

I check with MATLAB coder and the function sum does not work linearly.

The algorithm depends on the length of the vector. The order seems to be deterministic and independent of values. The algorithm has few branches and quite complex. I cannot post for the reason that kind of reverse engineering.

If you have matlab coder, you can do your own reverse engineering if you are curious how it works.

[Paul]

Matlab Coder generates source code for built-in functions? I did not know that. I thought the built-in functions are called from the source code and linked in to the executable.

[Bruno Luong]

I can see the C code for sum. It might be just a replicate of the built-in function, but it supposes to return the same result.

[Jonas K.]

Nice to know, thanks for your efforts! I haven't used the MATLAB coder so far for this purpose, but I may try it out, if the information given by the blog post linked by Steven Lord still doesn't help.

[Bruno Luong]

I don't think Loren's blog include details about sum. And beside this post is aimed to alert the behavior change in R2021b, however you are using R2019a.