Why do I get different result in different versions of MATLAB (2016 vs 2021)?-CodePudding

Why do I get different results when using the same code running in different version of MATLAB (2016 vs 2021) for sum(b.*x1) where b is single and x1 is double. How to avoid such error between MATLAB version?

MATLAB v.2021:

sum(b.*x1)

ans =

  single

  -0.0013286

MATLAB 2016

sum(b.*x1)

ans =

  single

  -0.0013283

CodePudding user response：

In R2017b, they changed the behavior of sum for single-precision floats, and in R2020b they made the same changes for other data types too.

The change speeds up the computation, and improves accuracy by reducing the rounding errors. Simply put, previously the algorithm would just run through the array in sequence, adding up the values. The new behavior computes the sum over smaller portions of the array, and then adds up those results. This is more precise because the running total can become a very large number, and adding smaller numbers to it causes more rounding in those smaller numbers. The speed improvement comes from loop unrolling: the loop now steps over, say, 8 values at the time, and in the loop body, 8 running totals are computed (they don’t specify the number they use, the 8 here is an example).

Thus, your newer result is a better approximation to the sum of your array than the old one.

For more details (a better explanation of the new algorithm and the reason for the change), see this blog post.

Regarding how to avoid the difference: you could implement your own sum function, and use that instead of the builtin one. I would suggest writing it as a MEX-file for efficiency. However, do make sure you match the newer behavior of the builtin sum, as that is the better approximation.

Here is an example of the problem. Let's create an array with N 1 elements, where the first one has a value of N and the rest have a value of 1.

N = 1e8;
a = ones(N 1,1,'single');
a(1) = N;

The sum over this array is expected to be 2*N. If we set N large enough w.r.t. the data type, I see this in R2017a (before the change):

>> sum(a)
ans =
  single
   150331648

And I see this in R2018b (after the change for single-precision sum):

>> sum(a)
ans =
  single
   199998976

Both implementations make rounding errors here, but one is obviously much, much closer to the expected result (2e8, or 200000000).