matlab strange result for mean of single vs double

Hello,

I get vastly different result when avraging an array using single vs double precision. Here is a code to reproduce my issue

a=single([99 100 65 101]'+zeros(4,512,512,128)); 
display(' single mean : '); 
display([num2str(mean(a,[2 3 4]))]);
display( ' double mean : ');
display(num2str(mean(double(a),[2 3 4])));

The output on matlab R2020a and R2020b is :

 single mean : 
64
64
64
64
 double mean : 
 99
100
 65
101

I would expect small difference between single and double due to precision error but not that big.

I could also expect the value to become out of range if using the sum operator but this should not occure with a well implemented mean function .

This also only occures if the first dismention is not one; else i have no problem :

for exemple if using

display([num2str(mean(a(1,:,:,:),[2 3 4]))]);

I get the ouput :

99

which is the expected output

ps: i have already read that post : https://www.mathworks.com/matlabcentral/answers/5401-matlab-function-mean-returns-the-exact-same-value-for-uint16-and-double-values-not-for-single but the problem seems different (in the above post only a small difference can be seen which is explained by precision error while in my case results are not even similar)

Thank you

Herve

[Bruno Luong]

Hmm, interesting it looks like a serious MATLAB bug to me.

This is fine as workaround

>> mean(mean(mean(a,2),3),4)

ans =

  4×1 single column vector

    99
   100
    65
   101

[Bruno Luong]

SUM suffers the same issue

>> sum(a,[2 3 4])

ans =

  4×1 single column vector

   1.0e+09 *

   2.1474836
   2.1474836
   2.1474836
   2.1474836

>> sum(double(a),[2 3 4])

ans =

   1.0e+09 *

   3.321888768000000
   3.355443200000000
   2.181038080000000
   3.388997632000000

[Herve Hugonnet]

ok I will report the bug to matlab support

[Rik]

@Bruno, I suspect the reason for that is that mean uses sum under the hood.

[Bruno Luong]

"ok I will report the bug to matlab support"

Good idea. Thanks.

[Matt J]

Whatever it is, it is a highly array-shape dependent thing

a=ones(2,512^2,'single')*99; 

rowmean1=mean(a,2).', 

rowmean1 = 1×2

99.3535 99.3535

rowmean2=mean(a.',1),

rowmean2 = 1×2

99 99

[Quy]

This was posted in 2020, and I am using 2023b with the same problem. I don't expect my issue to be resolved anytime soon, although I suspect that my issue is probably related to precision?

K>> meanW = mean(wavelet_array3,2);
K>> meanS = zeros(height(wavelet_array3),1);
K>> for ii = 1:height(wavelet_array3); meanS(ii) = mean(wavelet_array3(ii,:)); end
K>> isequal(meanW,meanS)

ans =

  logical

   0

K>> d = meanW-meanS; max(abs(d))

ans =

   3.1572e-16

[Matt J]

@Quy It doesn't look like the same problem, nor any kind of problem. Both of your calculations are in doubles and differ by amounts that would be expected with that precision.

[Quy]

Maybe you're right and it is expected because of precision... But the output are still different answers when you would expect them to be the same.

[Matt J]

Feature, not bug.

[Paul]

Original issue still not fixed?

a=single([99 100 65 101]'+zeros(4,512,512,128)); 
display(' single mean : '); 

 single mean : 

display([num2str(mean(a,[2 3 4]))]);

64
64
64
64

display( ' double mean : ');

 double mean : 

display(num2str(mean(double(a),[2 3 4])));

 99
100
 65
101

[John D'Errico]

@Quy - the difference there is just the impact of floating point arithmetic. Welcome to the wonderful wacky world of floating point arithmetic. This is definitely not a bug though.

[Matt J]

Original issue still not fixed?

It's still a matter of debate whether there is anything to fix. Might not the user be responsible for choosing a data type with sufficient range and precision to accomodate the calculation? These kinds of hazards show up with int type operations all the time.

Alternatively, you could propose that mean(X,'double') be the default behavior of mean(X), but that might have speed implications.

[Paul]

I didn't realize it's still a matter of debate. Why is that? How did "int type operations" enter the discussion? The original problem illustrated by the OP was about single vs. double.

[Matt J]

Well, it's just a matter of the numeric type not having enough range or precision to accomodate the operation being done. It's the same kind of thing as what we see below:

mean(double(50:50:200).^2)

ans = 18750

mean(uint8(50:50:200).^2)

ans = 255

Here, it's the .^2 operation that causes an overflow in the case of uint8. Should power() be made to do overflow pre-checks and pre-convert its input to doubles?

[Bruno Luong]

I agree with Paul. integer overflow is not the same with floating round off.

The sum will be more accurate with the sum calculation is carried out hierarchically on floating points. This won't fix for fix-point array.

[Matt J]

I dont know what "carried out hierarchically" means. If it means "when converted to doubles", my example above shows explicitly that accuracy improves when integers are converted to doubles.

[Paul]

a=single([99 100 65 101]'+zeros(4,512,512,128));

mean is implemented in an m-file

which mean(a,[2 3 4])

/MATLAB/toolbox/matlab/datafun/mean.m

For this use case, I believe the observation arises from the call to sum on line 125 (which is line 127 in 2021b for whatever that's worth, though 2021b exhibits the same behavior) with flag being 'default'

dbtype /MATLAB/toolbox/matlab/datafun/mean.m 118:130

118       if omitnan     
119           % Compute sum and number of NaNs
120           m = sum(x, dim, flag, 'omitnan');
121           nr_nonnan = mysize(x, dim) - matlab.internal.math.countnan(x, dim);
122           % Divide by the number of non-NaNs.
123           y = m ./ nr_nonnan;
124       else
125           y = sum(x, dim, flag) ./ mysize(x,dim);
126       end
127   end
128       
129   end
130   

All of the elements of a are positive

all(a(:)>0)

ans = logical
   1

and realmax for a single is quite large

format long e
realmax('single')

ans = single
   3.4028235e+38

Callling sum over the specified dimensions yields different results

[sum(a,[2 3 4],'default') sum(double(a),[2 3 4],'default')]

ans = 4×2

1.0e+00 * 2.1474836e+09 3.3218888e+09 2.1474836e+09 3.3554432e+09 2.1474836e+09 2.1810381e+09 2.1474836e+09 3.3889976e+09

but summing over all dimensions yields the same result as @Bruno Luong noted three years ago

[sum(a,[1 2 3 4],'default') sum(double(a),[1 2 3 4],'default')]

ans = 1×2

1.0e+00 * 1.2247368e+10 1.2247368e+10

isequal(ans(1),ans(2))

ans = logical
   1

What could be "the operation being done" for which single does not have enough range or precision to sum the elements of a subset of 'a', all of which are postive, and where the sum of the whole 'a' matrix is less than realmax and gives the correct result for single or double?

"Should power() be made to do overflow pre-checks and pre-convert its input to doubles?"

Not in my opinion.

[Matt J]

All of the elements of a are positive...and realmax for a single is quite large

Yes, but I think @John D'Errico demonstrated with,

asum + 100 == asum

that this is not enough to prevent all precision issues. At some point, as the number you are accumulating becomes very large, a small increment in the running sum is no longer captured.

but summing over all dimensions yields the same result

My guess is that is a multithreading effect. Matlab's internal multi-threading decomposes the sum into sub-chunks and then combines the partial sums at the end. If the partial sums are not too different in magnitude, you will not have the precision issue that John describes. How the decomposition is done is a complicated internal decision by the software, and not something the user will be able to predict or rely on.

[Bruno Luong]

asum + 100 == asum

I'm not sure if such equality occurs systematically for large enough asum. It could be a LSB of matissa would incease by 1, depending on round off and internal CPU precision of plus operation.

I'm sure how the multi-thread sum engine decides the chunk size of serial sum. It could be that it would be more , suitable for double array and not for single array. Ideally advanced user should have a hand on such setting. I recorgnize that designing such robust sum/mean is not a trivial task.

[Stephen23]

"I'm not sure if such equality occurs systematically for large enough asum."

Here is a gigabyte of them:

asum = 1e19:1.0000000001e19;
numel(asum)

ans = 999999489

all((asum+100) == asum)

ans = logical
   1

"The original problem illustrated by the OP was about single vs. double."

For a sufficiently large array you would get the same effect using a DOUBLE array.

[Paul]

Is there any way to determine if the multi-threading, or chunking, decomposition of sum is different for double and single in general, or more specifically for this problem? Not to determine what the decomposition actually is, but just determine that it's different for the two data types.

[Matt J]

Is there any way to determine if the multi-threading, or chunking, decomposition of sum is different for double and single in general, or more specifically for this problem?

Not that I know of, but here's a little more insight. The problem only seems to occur when the elements being summed are discontiguous in memory. This is significant because it is easy to imagine that discontiguous data might not be partitioned in chunks.

Here's a simpler, 2D version of the problem to illustrate,

B=99*ones(2, 512*512*128,'single');

mean(B,2) %wrong

ans = 2×1

64 64

mean(B,2,'double') %right

ans = 2×1

99 99

However, if we perform the exact same computation, but with the data to be averaged placed along the 1st dimension (so that they are continguous), a conversion to double is not necessary,

mean(B',1)'

ans = 2×1

99 99

[Bruno Luong]

Interesting

[Paul]

it is easy to imagine that discontiguous data might not be partitioned in chunks.

I'd like to pull this thread further, but using sum because I think that's where the problem lies

B = single(99)*ones(2, 512*512*128,'single');

The "correct" answer for the sum of each row is

format long
dsum = 99*512*512*128

dsum = 
     3.321888768000000e+09

Summing both rows of discontiguous data

Bsum = sum(B,2)

Bsum = 2×1

1.0e+09 * 2.1474836 2.1474836

yields the exact same result as a naive, non-chunked loop

ssum = single(0); ninetynine = single(99);
for ii = 1:(512*512*128)
    ssum = ssum + ninetynine;
end
isequal(Bsum(1),ssum)

ans = logical
   1

But summing just one row of B, which is discontiguous, yields the correct result

Bsum = sum(B(1,:));
isequal(Bsum,dsum)

ans = logical
   1

Does this imply that B(1,:) is copied to contiguous memory and then summed?

My mental model of sum(B,2) is

Bsum = zeros(2,1,'single');
for ii = 1:2
    Bsum(ii) = sum(B(ii,:));
end
Bsum

Bsum = 2×1

1.0e+09 * 3.3218887 3.3218887

I guess my mental model needs to be adjusted.

[Matt J]

Does this imply that B(1,:) is copied to contiguous memory and then summed?

Yes, extracting a subset of a matrix through indexing, e.g., B(i,:) or B(:,j), always creates a copy (except possibly when the subset is just a scalar). And, because in this case B(1,:) is just a vector, it will be contiguous.

Answer 1

[John D'Errico]

You need to understand what you did, AND you need to understand floating point arithmetic.

a=single([99 100 65 101]'+zeros(4,512,512,128));
mean(a,[2 3 4])

ans = 4×1

64 64 64 64

So what happened there? What is the mean? A mean is just a sum of the list of numbers, divided by the number of elements in that list.

asum = sum(a,[2 3 4])

asum = 4×1

1.0e+09 * 2.1475 2.1475 2.1475 2.1475

You can see here that the sum is the SAME, for each of those cases. And that cannot be right, unless you understand what happened. When you compute those sums in single precision, at some point, we have a number where if you add a number on the order of 100 to it, nothing changes.

asum + 100 == asum

ans = 4×1 logical array
   1
   1
   1
   1

Do you see that the accumulation in the sum no longer works? We cannot increase those values, essentially an overflow condition. And that means we get your result.

asum/prod([512,512,128])

ans = 4×1

64 64 64 64

How do you fix this? First, I would ask why you are doing this in single precision. The simple solution is to do as @Matt J has suggested. Either compute the mean in multiple steps, or convert to double for the computation.

amean1 = mean(mean(mean(a,2),3),4)

amean1 = 4×1

99 100 65 101

amean2 = single(mean(double(a),[2 3 4]))

amean2 = 4×1

99 100 65 101

[Paul]

Why does summing all of the elements of a yield the expected result? Wouldn't that also suffer from the same problem?

a=single([99 100 65 101]'+zeros(4,512,512,128));
isequal(sum(a,'all'),sum(double(a),'all'))

ans = logical
   1

[Matt J]

Wouldn't that also suffer from the same problem?

That would depend on how the sum was multithreaded, see also my comment above. Depending on the multithreading, the elements a(i) would be summed in a different order, which would affect the occurence of these kinds of precision failures.

Here is a mini-example of the order-dependence:

e=eps('single')/2;
x=round([1,e,e]./e); 
y=flip(x);

isequal(sum(x), sum(double(x)))

ans = logical
   0

isequal(sum(y), sum(double(x)))

ans = logical
   1

[Herve Hugonnet]

I understand that all the elements are summed before beeing devided by the total number of elements. But my comment was that the mean should not be done this way and that this algorithm migh induce unexpected error in many consumer codes. Furthermore, as highlighted by the fact that

mean(mean(mean(a,2),3),4)

returns the correct result. A different and more precise implementation of the mean algorithm is possible using single precision. I would appreciate if mathlab updated the algorithm.

"First, I would ask why you are doing this in single precision." --> because single precision is faster and uses less memory

[Paul]

@Herve Hugonnet

Did you report this behavior to TMW Tech Support a few years ago (relevant comment)? If so, do you recall and would you mind sharing their response?

[Matt J]

because single precision is faster and uses less memory

Note that,

mean(a,[2,3,4],'double')

allows you to do the sum in double precision while keeping a in single precision form to conserve memory.

and that this algorithm migh induce unexpected error in many consumer codes.

True, but most people understand that when your array is large, you should not do such computations in single precision. That was the whole point of introducing the syntax mean(___'double') back in R2008 or so.

[Herve Hugonnet]

@Paul I am sorry but i did not report this to the tech support only uplaoded it here

[Herve Hugonnet]

@Matt J thanks for the tip regarding ,'double'

Answer 2

[Matt J]

This is fine as workaround...mean(mean(mean(a,2),3),4)

Another workaround is to force the mean calculation in double type,

a=single([99 100 65 101]'+zeros(4,512,512,128)); 
mean(a,[2,3,4],'double')

ans =

    99
   100
    65
   101

Answer 3

[Matt J]

I reported the alternative example below to Tech Support and their reply was that they do in fact consider it a bug. I'm not sure from the reply I got if the bug is the occurence of precision underflow itself, or if the bug is just the inconsistency in the result when taking means/sums across different dimensions.

B=99*ones(2, 512*512*128,'single');

mean(B,2) %wrong

ans = 2×1

64 64

mean(B,2,'double') %right

ans = 2×1

99 99

mean(B',1)' %right?

ans = 2×1

99 99