Why is Plus so much slower than Total?

Question

Why is Plus so much slower than Total for summing lists?

list = RandomReal[1, 10000000];
AbsoluteTiming[Plus @@ list]
AbsoluteTiming[Total[list]]

 {1.67408, 5.00073*10^6}
 {0.00570379, 5.00073*10^6}

It seems like Plus is a few hundred times slower than Total, but of course both give the same result.

Answer 1

Performance issues with large numeric arrays very often happen because of unpacking of PackedArrays. Numeric operations on Packed arrays are very well optimized, while unpacking and operations on huge unpacked lists are usually slow.

If you turn on packing warnings (On["Packing"]) before your tests, you will see that applying Plus to packed array leads to unpacking. On the contrary, Total can take PackedArrays as input and the calculation is very fast.

If you work with large numeric arrays you should avoid unpacking as much as possible.

More information on packed arrays can be found here:

What is a Mathematica packed array?

Guidelines for avoiding the unpacking of a packed array

Answer 2

In Mathematica, "everything is an expression". All data and code is represented in the same manner, as Mathematica expressions. Thus expressions must be very general. They need to be able to represent anything that appears in the language.

There is usually a tradeoff between generality and performance. Mathematica expressions, which can contain everything from symbols through strings to thousand-digit precision numbers, are not competitive on performance with MATLAB matrices, which can only contain machine numbers (i.e. low-precision numbers that computers can work very quickly with).

For this reason, Mathematica has an alternative, optimized way to store lists of machine-precision numbers internally. This storage method is called "packed array".

Only Lists can be packed, and only if they contain numbers which are all of the same type.

Thus

arr = {1,2,3} (* same as List[1,2,3] *)

can be packed, but

expr = something[1,2,3]

cannot.

When you Apply Plus to an array, like arr above, its head changes to Plus: Plus[1,2,3]. This means that it can no longer be packed. This severely impacts performance.

The following tests demonstrate that this unpacking accounts for the performance difference.

list = RandomReal[1, 10000000];

Developer`PackedArrayQ[list]
(* True *)

unpackedList = 
   Developer`FromPackedArray[list]; // AbsoluteTiming
(* {0.258629, Null} *)

Total[list] // AbsoluteTiming
(* {0.004191, 4.99851*10^6} *)

Total[unpackedList] // AbsoluteTiming
(* {1.24454, 4.99851*10^6} *)

Plus @@ list // AbsoluteTiming
(* {1.58597, 4.99851*10^6} *)

Total@Developer`FromPackedArray[list] // AbsoluteTiming
(* {1.52471, 4.99851*10^6} *)

Answer 3

Introduction

@Shadowray has already put a finger on the source of the time difference in the OP's example: unpacking. But it seems to me that the question points out just a symptom of a broader question of how best to add up lists and arrays. One might want to consider the problem of addition from a design perspective and consider the roles of the tools available. Total is my default go-to function for numerical summation of arrays, but there are some surprises such as the dot product with a vector of ones or folding Plus. It is interesting to consider why.

The computational environment

Some basic tools for adding arrays are the following five:

• Plus + Apply: Plus @@ list.
• Total[list, levelspec], which with both arguments is formally equivalent to a combination of Plus, Apply and Flatten, but not computationally equivalent.
• Others, such as, Sum, Fold, and Dot, which perhaps appear to have other purposes than simple summing, but with the right coding they will sum a finite list: Sum[x, {x, list}], Fold[Plus, list], and {1, 1,...,1}.list.

There are two types of summation, symbolic addition (both in mathematical formulas and in code writing/generation) and numerical addition (complex, real or integer). For symbolic addition, Plus is the primary tool. The syntax of Total or Sum may be more convenient in many cases, and they all seem of comparable speed (summing a million random monomials). Since the question concerns numerics, let's leave it at that and focus on numerical summation.

For numerical addition, let us first note that current CPUs are designed for scalar, vector and matrix arithmetic. One should expect standard libraries such as MKL, BLAS, and LAPACK to take advantage of the relevant features. Packed arrays are designed to also take advantage of these strengths. This, in fact, gives Dot advantages that one might not suspect.

Another factor are number of system limits such as these:

Cases["CompileOptions" /. SystemOptions["CompileOptions"], 
  HoldPattern[s_ /; StringMatchQ[s, ___ ~~ "CompileLength"] -> len_] :> {s, len}
  ] // TableForm

When these boundaries are crossed, there is usually a discontinuous jump in timing, which could be either up or down, depending on the time to compile and the speed-up, if any, of the compiled code. In the examples below, one sometimes perceives other possible boundaries at lengths of around 1000 and 10^6. I do not know the source of these jumps.

Finally there are packed and unpacked arrays. Arrays can be packed at different levels (from the bottom up). An unpacked array consists of a list of pointers to its elements. Unpacking an array involves making copies of its elements and constructing a list of pointers. "Unpacking to level 1" means copying array[[1]], array[[2],... as packed arrays and constructing a list of pointer to them. Copying entailed by unpacking can take a significant bit of time and should be avoided. On the other hand, if not packed or already unpacked, passing elements of the arrays is efficiently done by pointers without further copying. This is the reason for the surprising efficiency of Fold in an example below.

The five methods mentioned above

On unpacked arrays Plus, Total, and Sum are indistinguishable. Fold is a bit faster. On lists of scalars Dot is faster Fold, but on unpacked lists of packed vectors, Dot is equivalent to the first three. Lowering the system option "ApplyCompileLength" will speed up Plus @@ list, but it still will not be as fast as Total. Fold loses some speed when the length of the list reaches the system option "FoldCompileLength", at which point the compile time adds to the running time. Increase "FoldCompileLength" to Infinity and the performance in the example below actually improves, which may seem counterintuitive. Packing the array Developer`ToPackedArray can help; Dot seems to do this automatically.

On packed arrays, Plus has a distinct advantage and a distinct disadvantage. The advantage is that on vector arguments, Plus uses the strengths of the hardware and the standard software libraries. (It seems to do this better than Total, but I do not know why Total should be slower.) The disadvantage is that when there are many arguments (addends), the slots have to be filled. For a packed array, this requires unpacking the first level of the array (the issue in the OP's example).

Total pretty much cannot be beat on a packed list of scalars. On an unpacked list, it is virtually tied with Plus It is less successful on vectors and higher dimensional arrays. Its implementation has changed recently. By V10.3, Total "is now using several threads in parallel." (@ilian) The sum it computes is different in V10.4.1 and V11.1.1 and the average rounding error appears to be less (see here). Maybe it can be improved more (see comparisons below).

Fold and Sum are rarely the worst or the best and are included mainly for the sake of comparison. In some situations, Sum, Total, and Plus seem equivalent (see examples below). Interestingly, Fold seems faster than the other methods on unpacked lists of packed vectors.

Dot is frequently the fastest option, the overhead of constructing the array of ones notwithstanding. Total beats it on a packed list of scalars, but not on an unpacked one by far, unless the list is (re)packed first.

Examples

The methods of adding behave similarly on lists of vectors and lists of matrices but differently on lists of scalars; so we'll look at lists of scalars and lists of vectors only. The basic test computations are the following:

Plus @@ list
Total @ list
Dot[ConstantArray[1., Length@list], list]
Fold[Plus, list]
Sum[x, {x, list}]

For unpacked arrays, we also tested Total @ Developer`ToPackedArray @ list. These were timed with RepeatedTiming. A code dump is at the end.

For scalars, different length lists were used to test these methods. For lists of vectors, different dimension arrays were constructed so that 10^7 additions would be performed. This means that the number of addends (and memory) varies, since {a1, a2} + {b1, b2} is two additions just as does c1 + c2 + c3.

Packed Arrays

Unpacked Arrays

Plus @@ array or f @@ array will unpack a packed array down to the first level only. Below we'll look arrays that are already unpacked to level one ('unpacked") and that are completely unpacked ("UNPACKED").

Surprisingly preventing Fold from compiling with SetSystemOptions["CompileOptions" -> "FoldCompileLength" -> Infinity] speeds it up:

Code dump

Code for the examples (general functions follow):

(*** PACKED ***)
divisors = Flatten@Outer[Times, 10^Range[1, 6], {1, 3}];
data = Table[timing[div], {div, divisors}];
comp["scalars", "packed"] = llp[data, divisors, "scalars"]

divisors = Flatten@Outer[Times, 10^Range[0, 3], {1, 2, 5}];
data = Table[timing[div, 10^7/div], {div, divisors}];
comp["vectors", "packed"] = llp[data, divisors, "vectors"]

(*** UNPACKED ***)
divisors = Flatten@Outer[Times, 10^Range[1, 6], {1, 3}];
data = Table[timing[div, "unpacked"], {div, divisors}];
comp["scalars", "unpacked"] = llp[data, divisors, "scalars"]

divisors = Flatten@Outer[Times, 10^Range[0, 3], {1, 2, 5}];
data = Table[timing[div, 10^7/div, "unpacked"], {div, divisors}];
comp["vectors", "unpacked"] = llp[data, divisors, "vectors"]

divisors = Flatten@Outer[Times, 10^Range[0, 3], {1, 2, 5}];
data = Table[timing[div, 10^7/div, "UNPACKED"], {div, divisors}];
comp["vectors", "UNPACKED"] = llp[data, divisors, "vectors"]

(*** Uncompiled Fold ***)
With[{opts = SystemOptions["CompileOptions"]},
 Internal`WithLocalSettings[
  SetSystemOptions["CompileOptions" -> "FoldCompileLength" -> Infinity],
  divisors = Flatten@Outer[Times, 10^Range[0, 3], {1, 2, 5}];
  data = Table[timing[div, 10^7/div, "unpacked"], {div, divisors}];
  comp["vectors", "Fold"] = llp[data, divisors, "vectors"],
  SetSystemOptions[opts]
  ]]

General utility functions:

ClearAll[timing];
timing[list_List, "unpacked"] := timing[List @@ list, Automatic];
timing[list_List, "UNPACKED"] := timing[Developer`FromPackedArray@ list, Automatic];
timing[list_List, _] := {
   RepeatedTiming[Plus @@ list; "Plus @@", 0.1],
   RepeatedTiming[Total @ list; "Total", 0.1],
   If[! Developer`PackedArrayQ[list],
    RepeatedTiming[Total@ Developer`ToPackedArray@ list; "Total[Packed]", 0.1],
    Nothing],
   RepeatedTiming[Dot[ConstantArray[1., Length@list], list]; "Dot", 0.1],
   RepeatedTiming[Fold[Plus, list]; "Fold", 0.1],
   RepeatedTiming[Sum[x, {x, list}]; "Sum", 0.1]
   };
timing[n_Integer, type_String: "packed"] := 
  timing[RandomReal[1, {n + 1}], type];
timing[m_Integer, n_Integer, type_String: "packed"] := 
  timing[RandomReal[1, {m + 1, n}], type];
timing[m_Integer, n_Integer, p_Integer, type_String: "packed"] := 
  timing[RandomReal[1, {m + 1, n, p}], type];

llp[data_List, divisors_List, "scalars", opts___] :=
  llp[
   Transpose[data[[All, All, 1]]/divisors],
   data[[1, All, 2]],
   opts,
   FrameTicks -> {{Automatic, 
      Automatic}, {Transpose[{Range@Length@divisors, 
        10^(HoldForm /@ Round[Log10@divisors, 0.01])}], Automatic}},
   FrameLabel -> {"number of scalars", "time per addition"},
   PlotLabel -> "Addition of scalars"];
llp[data_List, divisors_List, t : "vectors" | "matrices", opts___] := 
  llp[Transpose[data[[All, All, 1]]], data[[1, All, 2]], 
   FrameTicks -> {{Automatic, 
      Automatic}, {Transpose[{Range@Length@divisors, 1 + divisors}], 
      Automatic}},
   FrameLabel -> {"number of " <> t, "time"},
   PlotLabel -> "Addition of " <> t, opts];
llp[data_List, labels_List, opts___] := ListLogPlot[data,
   opts, Frame -> True, PlotLegends -> labels,
   PlotStyle -> 
    With[{n = Length@labels}, 
     Table[ColorData["Rainbow"][i/(n + 1)], {i, n}]],
   Joined -> True, PlotMarkers -> Automatic, PlotRange -> All];