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

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.

• Without digging into details, this is almost certainly related to machine precision limitations, and rounding in machine precision. \$MachinePrecision is about 16, which means rounding occurs at that resolution (in practice, 53 bits). You can work in higher precision by using WorkingPrecision option in RandomReal. – kirma Oct 12 '14 at 8:20
• @kirma Is there a way to set globally the machine precision at higher values? – An old man in the sea. Oct 12 '14 at 8:30
• By checking that all reals in your computations are arbitrary-precision numbers (which WorkingPrecision does in this case, and for explicit numbers backtick like 1.520), numerical computations are performed in arbitrary-precision arithmetic. This is a bit different from machine precision, which is directly based on hardwired implementation on the CPU. This is a slightly convoluted topic to understand correctly; I suggest looking at the documentation on (arbitrary) precision and questions on precision and arbitrary-precision. – kirma Oct 12 '14 at 8:57
• @kirma Many Thanks! ;) – An old man in the sea. Oct 12 '14 at 9:42

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"}]


• Didn't know about Method -> "CompensatedSummation"` - thanks, +1! – dr.blochwave Oct 12 '14 at 14:05
• This answer (which I upvoted, of course) is a nice illustration of what I was going to put in a comment (and will still do so now): finite precision addition is not associative. – Daniel Lichtblau Oct 12 '14 at 17:19
• But exactly what is "compensated summation"?? Is it what's described in hal.archives-ouvertes.fr/file/index/docid/517618/filename/… ? – murray Oct 15 '14 at 20:55
• @murray, it's a nice floating point trick due to Velvel Kahan. Search around; this is in fact a very well-studied procedure. – J. M.'s ennui May 6 '15 at 14:33
• @J.M. related (79174) – dr.blochwave Sep 3 '15 at 8:20