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

For example

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

Out[342]= {0.0198735, Null}


returns PackedArray (examined by DeveloperPackedArrayQ[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!!

• Er… personally I think this question is a little scattered, maybe you can consider divide it into 2 or 3 questions? – xzczd Jan 14 '16 at 12:53
• @xzczd Well, on the contraray, I would call "dividing it into 2 or 3" as scattered : ), and I suffered a lot from scattered questions when searching solutions on SE site. I organized them together, because all of them generate the same list, and they are what one can think of when dealing this kind of problem (I mean calculating data on a mesh) – matheorem Jan 14 '16 at 14:05
• Related: (33624) – Mr.Wizard Jun 11 '17 at 8:26

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:

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:

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

And a more recent question on this site:

• Thank you Mr.Wizard. What do you mean by "free" performance? – matheorem Jan 14 '16 at 15:44
• @matheorem What I meant is that where the developers saw fit to provide low-level optimizations you get a very nice performance boost without any additional effort like compiling to C. However the fast performance in such a case cannot be taken to mean that the function (like Outer) itself is fast. This has tripped me up a number of times as I tend to think, subconsciously if nothing else, things like "ListCorrelate is fast" then I am surprised when I see that my program is slow, only to realize that I used ListCorrelate in a non-default way, which is not optimized like the default. – Mr.Wizard Jan 14 '16 at 16:14
• I understand. Another question, Why Table is much faster than Map in this case ? Is it general? – matheorem Jan 14 '16 at 16:22
• @matheorem I added a few comments on Table that I hope you find illuminating. – Mr.Wizard Jan 14 '16 at 16:36
• Wow, I am shocked too ! Thank you so much for providing such valuable links, I learned a lot. And I find Array[#1 + #2 &, {1000, 1000}] is 10 times faster than Outer[#1 + #2 &, Range[1000], Range[1000]], and this is not easily read from the TraceInternal`, since the trace output are almost the same. It seems that they both works pairwise, right? – matheorem Jan 15 '16 at 9:02