Note: this question is not strictly related to Matlab or any other environment (even though I would prefer that you refer to my Matlab code).

I've n = 10000 equally spaced elements between -pi and pi in an array produced like this (in Matlab):

n = 10000;    
x = linspace(-pi, pi, n);

Then I create a second array x_r, which is basically contains the elements of x rounded to the third decimal place:

x_r = chop(x, 3);

Then, I calculate the sin of all elements of x and I call the resulting array y. I do the same for x_r.

y = sin(x);
y_r = sin(x_r);

Now, I try to sum all elements of y and then all elements of y_r, that is:

s = sum(y);
s_r = sum(y_r);

If I try to print s and s_r with 100 decimal digits after the dot, I obtain these huge numbers (and probably the first one is not finished):

s = -0.0000000000003170828901848210175171876477281318270923932134408573801920283585786819458007812500000000
s_r = -0.0000000000001331495677603378169351344695314764976501464843750000000000000000000000000000000000000000

Now, I want to find the mean of y and y_r, which I save respectively in m and m_r:

m = mean(y);
m_r = mean(y_r);

If I try to print m and m_r with 100 decimal digits, I obtain other two huge numbers:

m = -0.0000000000000000317082890184820990794524483366356918770088064486473850012160369260527659207582473755
m_r = -0.0000000000000000133149567760337822432508902790402516613129923689416538035157344666004064492881298065

My question is (as the title of this question indicates) which of the two operations is more robust to digit cancellation, the sum or the mean operation? And why?

My guess is that the mean operation is less robust to digit cancellation because it does more calculations, but I'm not sure, and I would not know if this would even be a right explanation. From the results, I'm not able to deduce anything useful honestly.

  • 1
    $\begingroup$ Mathematically, the only difference between sum and mean is that the latter involves a single additional division by a (in your case) moderate-sized number, so the stability properties can't be very different. Without knowing how Matlab implements the mean, it's hard to say much more -- but there are more stable approaches than the simple "sum-and-divide" approach, see, e.g., stackoverflow.com/questions/21918641/… $\endgroup$ Mar 20 '16 at 22:45

The first thing to realize is that reality doesn't work in decimal. And sin(x) in particular has no respect for decimals - it works in radians. So when you truncate x_r by up to 0.001, that really is 6.36619772368E-4'th of a right angle.

When I then see that y_r has so many trailing zeroes, I assume there is a massive loss of precision. This probably is because MATLAB has to revert to direct numerical approximation with finite precision, while the non-rounded version allows a more direct formula.

This means digit cancellation is just straight-out irrelevant in the second case.

  • $\begingroup$ @Kirill: Of course - with one bit that would be one digit (0.5), with 2 bits 2 digits (0.25), and in general with N bits you need N digits. But note that these results are the outcome of sin(x_r), in contrast to sin(x) so there's a definite difference in precision used. The non-rounded example appears to use more than double precision, or possibly another closed-form formula. $\endgroup$
    – MSalters
    Mar 21 '16 at 14:49
  • $\begingroup$ Note that s_r is exactly $153511\times 2^{-60}$, and s is exactly $1570118422128477\times 2^{-92}$, so we're just looking at the finite decimal expansions of some floating-point numbers. I don't think the trailing zeros say anything about lost precision. Since the sum would be $0$ in exact arithmetic (the summand is antisymmetric), we're just looking at accumulation of different roundoff errors. $\endgroup$
    – Kirill
    Mar 21 '16 at 16:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy