8
$\begingroup$

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.

$\endgroup$
4
  • 2
    $\begingroup$ 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. $\endgroup$
    – kirma
    Oct 12, 2014 at 8:20
  • $\begingroup$ @kirma Is there a way to set globally the machine precision at higher values? $\endgroup$ Oct 12, 2014 at 8:30
  • 2
    $\begingroup$ 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.5`20), 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. $\endgroup$
    – kirma
    Oct 12, 2014 at 8:57
  • $\begingroup$ @kirma Many Thanks! ;) $\endgroup$ Oct 12, 2014 at 9:42

1 Answer 1

21
$\begingroup$

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

enter image description here

$\endgroup$
5
  • 1
    $\begingroup$ Didn't know about Method -> "CompensatedSummation" - thanks, +1! $\endgroup$ Oct 12, 2014 at 14:05
  • 5
    $\begingroup$ 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. $\endgroup$ Oct 12, 2014 at 17:19
  • 4
    $\begingroup$ But exactly what is "compensated summation"?? Is it what's described in hal.archives-ouvertes.fr/file/index/docid/517618/filename/… ? $\endgroup$
    – murray
    Oct 15, 2014 at 20:55
  • 1
    $\begingroup$ @murray, it's a nice floating point trick due to Velvel Kahan. Search around; this is in fact a very well-studied procedure. $\endgroup$ May 6, 2015 at 14:33
  • $\begingroup$ @J.M. related (79174) $\endgroup$ Sep 3, 2015 at 8:20

Your Answer

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

Not the answer you're looking for? Browse other questions tagged or ask your own question.