Looking up one additional array element increases runtime by three orders of magnitude

Question

The Problem

I have a list myList, which is a 487135 x 3 x 2 array of integers between -1000 and 1000. I want to be able to gather elements from a different list by their indices using elements from myList, like so:

GatherBy[Range[99], f[myList[[#]] ] &];
(* {0.000691, Null} *)

but even with a simple f, such as OddQ, calling any more than 99 elements caused the time to increase dramatically:

GatherBy[Range[100], f[myList[[#]]] &]; // AbsoluteTiming
(* {1.15235, Null} *)

At first, I thought this was an error with GatherBy or with my f, but after experimenting, I discovered the same issue happening with Map itself.

myList[[#]] & /@ Range[99]; //AbsoluteTiming
(* {0.000113, Null} *)

myList[[#]] & /@ Range[100]; // AbsoluteTiming
(* {1.1329, Null} *)

I figured out a way to work around the problem, but I have no idea what's causing this, or even how to make test code to replicate the problem fully. I uploaded the first thousand elements of myList on PasteBin here; the same qualitative behavior occurs on my computer and version of Mathematica (11.3 for Linux) for that partial list as for the full myList, though the slowdown is only a factor of ~40.

What I've tried:

Using a different set of indices.

myList[[#]] & /@ Range[1001, 1099]; // AbsoluteTiming
(* {0.000101, Null} *)

myList[[#]] & /@ Range[1001, 1100]; // AbsoluteTiming
(* {1.16167, Null} *)

I tried several more times, including taking numbers at random rather than using Range (in case I had a bizarre issue on every 100th element or something).

Replicating with RandomInteger.

myRandomList = RandomInteger[{-1000, 1000}, {487135, 3, 2}];
myRandomList[[#]] & /@ Range[100]; // AbsoluteTiming
(* {0.000149, Null} *)

This newly-generated array does not have the same problem.

Appending a RandomInteger to the end of myList.

myList2 = Join[myList, RandomInteger[{-1000, 1000}, {1, 3, 2}]];
myList2[[#]] & /@ Range[100]; // AbsoluteTiming
(* {0.000179, Null} *)

This was even more baffling - why would adding something to the end change behavior of the indices at the beginning?

Appending a RandomInteger to the end of myList, and then deleting it.

myList3 = Drop[Join[myList, RandomInteger[{-1000, 1000}, {1, 3, 2}] ], -1];
myList3 === myList
(* True *)

myList3[[#]] & /@ Range[100]; // AbsoluteTiming
(* {0.000161, Null} *)

myList3 is identical to myList, but does not have the problem.

Appending a non-random element to the end of myList, and then deleting it.

myList4 = Drop[Join[myList, {{{1, 1}, {1, 1}, {1, 1}}}], -1];
myList4 === myList
(* True *)

myList4[[#]] & /@ Range[100]; // AbsoluteTiming
(* {1.13683, Null} *)

myList4 is identical to myList, and does have the problem.

Running ByteCode on the three identical arrays.

ByteCount /@ {myList, myList3, myList4}
(* {163677440, 23382640, 163677440} *)

The arrays that have the problem are the same size, and 6.99996 times bigger than the array without.

Taking a slice of myList to replicate the problem on a smaller scale.

myList5 = myList[[1 ;; 1000]];
myList5[[#]] & /@ Range[99]; // AbsoluteTiming
(* {0.000131, Null} *)

myList5[[#]] & /@ Range[100]; // AbsoluteTiming
(* {0.005996, Null} *)

Even with only a thousand elements, it still takes many times longer to pull 100 elements than to pull 99.

My Questions

So, I have a workaround now for my actual program - add a random element to the end, and then delete it. But I very much want to figure out:

• Why is myList so much larger than the same array with one element added and then removed?
• What's so special about the number 100 when calling indices?
• Why is the slowdown more than three orders of magnitude?

Answer 1

The change you're seeing when increasing the list size to 100 is due to the autocompilation feature. There is a system options that controls when this switchover occurs:

SystemOptions["CompileOptions"->"MapCompileLength"]

{"CompileOptions" -> {"MapCompileLength" -> 100}}

So, at length 100, Mathematica takes time to compile the mapping function, and then it uses the compiled function instead of the original uncompiled function. The timing discrepancy is all due to the time to compile your function. The compilation process is faster when the original function is packed. Your original function was not packed, and creating the compiled function took a very long time (it is possible that Mathematica should automatically pack the function before compilation, but there may be good reasons not to do so). Note that this extra compilation time only happens once, although that is only a small consolation.

At any rate, the solution is to work with packed arrays when possible. You can use Developer`ToPackedArray to convert unpacked arrays to packed arrays. Here is a comparison:

m = RandomInteger[10^6, {10^6, 3, 2}];
um = Developer`FromPackedArray[m];

Developer`PackedArrayQ /@ {m, um}

{True, False}

(note that RandomInteger produces packed arrays be default, so I convert the packed array to an unpacked array instead)

Then:

GatherBy[Range[100], OddQ[m[[#]]]&]; //AbsoluteTiming
GatherBy[Range[100], OddQ[um[[#]]]&]; //AbsoluteTiming

{0.000755, Null}

{1.57199, Null}

I don't know what is actually being compiled, but here's an example of the possible difference:

With[{mat=m}, Compile[{{pos,_Integer}}, OddQ[mat[[pos]]]]]; //AbsoluteTiming
With[{mat=um}, Compile[{{pos,_Integer}}, OddQ[mat[[pos]]]]]; //AbsoluteTiming

{0.000452, Null}

{16.2089, Null}

If the above is representative to the actual compilation process, it was probably aborted to achieve the Gather timings shown above.

Answer 2

To complement Carl Woll's answer it may be useful to know that the Map operation in GatherBy, the source of the compilation, can actually be intercepted, as also shown by Carl Woll with his GatherByList, reproduced below:

GatherByList[list_, representatives_] := Module[{func},
    func /: Map[func, _] := representatives;
    GatherBy[list, func]
]

This knowledge and function allows us to perform the entire Part extraction at once, and if f can be written as a vectorized operation too (as OddQ is) the whole thing should be faster. Let's test it:

n = 487135;
m = RandomInteger[10^6, {n, 3, 2}];
r = RandomChoice[Range @ n, 50000];

GatherBy[r, OddQ @ m[[#]] &];    // RepeatedTiming
GatherByList[r, OddQ @ m[[r]] ]; // RepeatedTiming

{0.182, Null}
{0.125, Null}