Why does taking advantage of Listable change the results of a numerical computation slightly?

Question

I have two variables: t0, and teta0. The first is computed using several nested sums, the second is computed taking advantage to some listable properties of some functions. Both should return the same, when using the same data/sample.

For a small random sample size (e.g. 10) the difference between them is of the order of 10^-15 or 10^-16 or 0, and when I do t0==teta0 I get True. When I increase the sample size (e.g. 999), the difference between them is sometimes of the same magnitude. But when I do t0==teta0 I get False.

Here is the code:

n = 999;

Tauijk = RandomReal[{1, 10}, {2, n}];

y = RandomReal[{1, 10}, {1, n + 1}][[1]];
one = ConstantArray[1, n];

teta0 = (Plus @@ 
     Flatten[(Tauijk.y[[2 ;; n + 1]])*(Tauijk*{y[[1 ;; n]], 
           y[[1 ;; n]]})*(Tauijk.one)^(-1) - (Tauijk*{y[[
           2 ;; n + 1]], y[[2 ;; n + 1]]}*{y[[1 ;; n]], 
          y[[1 ;; n]]})])/(Plus @@ 
     Flatten[(Tauijk.y[[1 ;; n]])*(Tauijk*{y[[1 ;; n]], 
           y[[1 ;; n]]})*(Tauijk.one)^(-1) - (Tauijk*{y[[
           1 ;; n]], y[[1 ;; n]]}*{y[[1 ;; n]], y[[1 ;; n]]})]);


t0 = (Sum[
     Sum[(Sum[y[[j + 1]]*Tauijk[[i, j]], {j, 1, n}])*Tauijk[[i, 
         k]]*y[[k]]/(Sum[Tauijk[[i, j]], {j, 1, n}]) - 
       y[[k + 1]]*y[[k]]*Tauijk[[i, k]], {i, 1, 2}], {k, 1, 
      n}])/(Sum[
     Sum[(Sum[
           y[[j]]*Tauijk[[i, j]], {j, 1, 
            n}]/(Sum[Tauijk[[i, j]], {j, 1, n}]) - 
         y[[k]])*Tauijk[[i, k]]*y[[k]], {i, 1, 2}], {k, 1, n}]);

This is strange because when I try something similar with some two other simpler variables like

 teta1 = (Plus @@ 
     Flatten[Tauijk[[1, 
        1 ;; n]]*(y[[2 ;; n + 1]] - teta0*y[[1 ;; n]])])/(Plus @@ 
     Flatten[Tauijk[[1, 1 ;; n]]]);

t1 = Sum[Tauijk[[1, j]]*(y[[j + 1]] - teta0*y[[j]]), {j, 1, n}]/
   Sum[Tauijk[[1, j]], {j, 1, n}];

I get the same order of magnitude for the difference as I would get for t0-teta0 and when doing t1==teta1, I get True. I simply cannot understand why this happens.

Any help would be appreciated.

Answer 1

When you work with machine numbers the precision depends on the order of summation. The built-in functions Plus and Total are more intelligent then straightforward summation. It can be shown by the following example

x = RandomInteger[{-100, 100}, {10000}]/RandomInteger[{1, 100}, {10000}];
res = Accumulate[x];
nx = N[x];
s[n_] := Sum[nx[[k]], {k, 1, n}];
p[n_] := Plus @@ nx[[;; n]];
tc[n_] := Total[nx[[;; n]], Method -> "CompensatedSummation"];

r = Range[1, 10000, 10];

ListLinePlot[{s /@ r - res[[r]], p /@ r - res[[r]], tc /@ r - res[[r]]}, PlotRange -> All, 
 PlotLegends -> {"Sum/Accumulate", "Plus/Total", "Total with compensation"}]