Why does Outer sometimes return a packed array, and sometimes not?

Question

For example

In[342]:= 
test3 = Flatten[Outer[Plus, Range[1000], Range[1000]], 
    1]; // AbsoluteTiming

Out[342]= {0.0198735, Null}

returns PackedArray (examined by Developer`PackedArrayQ[test3] ), while

In[334]:= 
test3 = Flatten[Outer[#1 + #2 &, Range[1000], Range[1000]], 
    1]; // AbsoluteTiming

Out[334]= {0.822719, Null}

Unpacks. And we can see the huge timing difference.

Why different unpacking behaviour of Outer?

What is more, the first Outer example is the fastest possible code. Even faster than Compiled to C target code.

Oddly, Map pure function and Map built-in Total is slow

In[281]:= 
test2 = (#[[1]] + #[[2]] &) /@ 
    Tuples[{Range[1000], Range[1000]}]; // AbsoluteTiming

Out[281]= {0.214871, Null}

In[340]:= 
test22 = Total /@ Tuples[{Range[1000], Range[1000]}]; // AbsoluteTiming

Out[340]= {0.176941, Null}

Though both of the above returns PackedArray. And they are much slower than trivial Table

In[299]:= 
test4 = Flatten[Table[i + j, {i, 1, 1000}, {j, 1, 1000}], 
    1]; // AbsoluteTiming

Out[299]= {0.0530814, Null}

The performance is so inconsistent. I was trying to calculate something like 251001 (-1 + #1) + 501 (-1 + #2) + #3 & on Tuples[{{1}, Range[501], Range[501]}], found it is a bottleneck, Tried Map, Outer, ApplyAll(Apply unpacks), sophisticated Compile, and finally found trivial Table is already almost the same efficiency as Compile!!

Answer 1

I think this is pretty clearly a case where Outer is simply overloaded with an optimized definition for Plus. A look at the Trace shows that there are no intermediate steps:

Trace[
  Outer[Plus, Range[3], Range[3]],
  TraceInternal -> True
]

{{Range[3], {1, 2, 3}}, {Range[3], {1, 2, 3}}, 
 Outer[Plus, {1, 2, 3}, {1, 2, 3}], {{2, 3, 4}, {3, 4, 5}, {4, 5, 6}}}

By comparison observe the individual evaluation steps with a Function:

Trace[
  Outer[#1 + #2 &, Range[3], Range[3]],
  TraceInternal -> True
]

{{Range[3], {1, 2, 3}}, {Range[3], {1, 2, 3}}, 
 Outer[#1 + #2 &, {1, 2, 3}, {1, 2, 3}],
 {(#1 + #2 &)[1, 1], 1 + 1, 2},
 {(#1 + #2 &)[1, 2], 1 + 2, 3},
 {(#1 + #2 &)[1, 3], 1 + 3, 4},
 {(#1 + #2 &)[2, 1], 2 + 1, 3},
 {(#1 + #2 &)[2, 2], 2 + 2, 4},
 {(#1 + #2 &)[2, 3], 2 + 3, 5},
 {(#1 + #2 &)[3, 1], 3 + 1, 4},
 {(#1 + #2 &)[3, 2], 3 + 2, 5},
 {(#1 + #2 &)[3, 3], 3 + 3, 6},
 {{2, 3, 4}, {3, 4, 5}, {4, 5, 6}}}

Basically you get "free" performance for certain operators, but in general Outer is a pairwise operation that does not use vector optimizations. See Leonid's post on the topic:

• Fastest way to calculate matrix of pairwise distances

---

An examination of Table will also show that there are hidden optimizations taking place. Consider this potentially baffling example:

Table[i + j, {i, 1, 1000}, {j, 1, 1000}]; // RepeatedTiming

n = 1000;
Table[i + j, {i, 1, n}, {j, 1, n}]; // RepeatedTiming

{0.0127, Null}

{0.220, Null}

I first learned of this from a post on Stack Overflow and it really shocked me. But it has a reasonable explanation which is that when Table sees explicit numeric ranges it triggers an optimized routine. When it does not, as in the second case, it allows for the possibility of a changing value of n in the j iterator by using a slower but more literal algorithm. Here are two examples where the value of n does change; try to guess the output of each before you check it yourself:

n = 3;
Table[i + j, {i, 1, n}, {j, 1, n++}]

n = 3;
Table[(n++; i + j), {i, 1, n}, {j, 1, n}]

Another Stack Overflow post I actually answered before I understood what was going on with Table:

• Using Array and Table Functions in Mathematica. Which is best when

(As usual Leonid gave me a polite nudge in the right direction.)

And a more recent question on this site:

• Why does array packing in Table behave like this?