16
$\begingroup$

I thought Kahan's summation method would make a nice example for students to use to think about round-off error [W. Kahan, Pracniques: Further Remarks on Reducing Truncation Errors, Commun. ACM 8  (1965), 40]. The method is available (I surmise) through the option Method -> "CompensatedSummation" of Total, but it's an easy program to write that is nicely laid out in the half-page article.

Below the vector x0 is a list whose elements are to be added. Note that s0 is the "exact" sum, that is, the sum of the exact numbers represented by the floating-point numbers in x0, rounded to machine precision. I consider it the target of the summation methods below.

SeedRandom[2];
x0 = 2 + RandomReal[1, 1000000];
s0 = N@Total[SetPrecision[x0, Infinity]];

Here is a norm (with sign) that measures the relative error in x - x0 in units of $MachineEpsilon.

meps[x_, x0_, eps_: $MachineEpsilon] := ((x - x0)/x0)/eps;

Then here are some ways to add up the vector. The timings are not particularly important, but they're somewhat interesting. I unpack x0 before apply Plus, since it would be unpacked anyway by Apply (@@), in order that the unpacking not be included in the timing; the sum is the same we would get with Plus @@ x0.

{With[{x0 = Developer`FromPackedArray[x0]}, 
   sP = Plus @@ x0; // AbsoluteTiming],
  sF = Fold[Plus, x0]; // AbsoluteTiming,
  sT = Total[x0]; // AbsoluteTiming,
  sC = Total[x0, Method -> "CompensatedSummation"]; // AbsoluteTiming
  }[[All, 1]]
(*  {0.131215, 0.120635, 0.000571, 0.01399}  *)

The example was chosen to show that Plus does something nearly like compensated summation but not exactly the same.

meps[{sP, sF, sT, sC}, s0]
(*  Plus       Fold     Total     Compensated  Version
   {-0.838972, 149.337, -56.2111, 0.}          10.4.1
   {-0.838972, 149.337, -2.51692, 0.}          11.1.1  *)

Aside from being surprised that Plus is such a slow method, I got to wondering what Plus is doing. With some effort, one can find that the order of the elements in x0 can affect the result of Plus. The only results I can find have errors of either 1 or 0 in the last bit. They are compared to the results of iterative adding using Fold.

meps[{Plus @@ Sort@x0, Fold[Plus, Sort@x0]}, s0]
(*  {-0.838972, -33.5589}  *)

meps[{Plus @@ -Sort[-x0], Fold[Plus, -Sort[-x0]]}, s0]
(*  {-0.838972, 182.057}  *)

SeedRandom[0];
With[{x0 = RandomSample[x0]}, meps[{Plus @@ x0, Fold[Plus, x0]}, s0]]
(*  {0., -158.566}  *)

SeedRandom[1];
With[{x0 = RandomSample[x0]}, meps[{Plus @@ x0, Fold[Plus, x0]}, s0]]
(*  {-0.838972, -93.9648}  *)

Does anyone know how Plus works? It would be nice to be able to explain it.

Improve this question
$\endgroup$
12
  • 3
    $\begingroup$ Might want to check on progressively larger examples to see if it scales as n log(n), that might indicate sorting under the hood. I remark that I do not know the answer offhand and probably will not have time to peek into the code to find out. $\endgroup$
    – Daniel Lichtblau
    Commented Jan 10, 2016 at 23:31
  • 1
    $\begingroup$ @george2079 Perhaps you mean extended-precision registers on the CPU/FPU are used to accumulate the result? That may be possible if the summand are chunked and the intermediate sums rounded to MP from time to time. If the whole computation is done with as little as 6 extra bits, the result rounds to the "exact" 53-bit result. $\endgroup$
    – Michael E2
    Commented Jan 11, 2016 at 17:19
  • 1
    $\begingroup$ If it is a matter of extended precision registers, this is likely to vary based on operating system, whether Compile was used somewhere, and probably a few other incidentals (like, phase of the moon). $\endgroup$
    – Daniel Lichtblau
    Commented Jan 12, 2016 at 15:56
  • 2
    $\begingroup$ @AlexeyPopkov I did check out the Orderless attribute: Plus can give different results (of at most one ulp) for different orderings of the arguments (via RandomSample). If the Orderless attribute affects evaluation, it is not by the arguments being sorted into canonical order before evaluation. It makes sense to me that orderless numerical functions like Plus would go directly to evaluation on numeric input. Simple example: Plus[1., -1., $MachineEpsilon/2] vs. Plus[1., $MachineEpsilon/2, -1.]. Orderless is more important for symbolic computations. $\endgroup$
    – Michael E2
    Commented May 6, 2017 at 20:07
  • 1
    $\begingroup$ Interesting question. Incidentally you forgot to include Tr, my favorite method to sum a list. Including it in v10.1 I find it twice as fast as (plain) Total with the same error. $\endgroup$
    – Mr.Wizard
    Commented Aug 1, 2017 at 9:14

1 Answer 1

Reset to default
2
$\begingroup$

Following up on your own comment about chunking I was able to replicate the result of Plus using a simple Partition with a block size of 100,000:

myPlus[x__] := Total /@ Partition[{x}, 1*^5, 1*^5, 1, 0] // Total

meps[myPlus @@ x0, s0]

-0.838972
Improve this answer
$\endgroup$
3
  • $\begingroup$ Interesting, but I think it's the unpacking, not the chunking, that leads to the Plus result: Any (evenly-dividing) block size works and so does plain unpacking. For instance meps[Total[Partition[x0, 1*^6, 1*^6, 1, 0], 2], s0] or meps[Total[DeveloperFromPackedArray@x0], s0]`. That the algorithm might depend on the data representation didn't occur to me until now. Coincidentally, I made this possibly related observation last night in chat. $\endgroup$
    – Michael E2
    Commented Aug 1, 2017 at 12:36
  • $\begingroup$ @MichaelE2 I may have fooled myself, but I also think I lost the point of this code by using Total[..., 2] at the last minute. Please see the restored definition and try other block sizes, which now for me give different results. $\endgroup$
    – Mr.Wizard
    Commented Aug 1, 2017 at 15:42
  • $\begingroup$ Okay, it's now giving me different results for different sizes. Testing Total[.., 2] more extensively shows it doesn't quite always give -0.83... but usually does. It's still interesting to me that Total on the unpacked x0 gives a different result than on the packed x0. $\endgroup$
    – Michael E2
    Commented Aug 1, 2017 at 16:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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