# 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?
• Many of the differences you found would be caused by packed arrays (see e.g. this Q/A) -- a lot would be explained if myList is not packed (you can check with DeveloperPackedArrayQ), in which case most (if not all) of the things you tried ended up packing it. You should also see a difference between mapping over Range[99] and Range[100] if you do On["Packing"] first. Commented Oct 23, 2018 at 21:34
• You're correct - myList is packed and myList3 is unpacked according to DeveloperPackedArrayQ. Doing On["Packing"] first on the myList[[#]] & /@ Range[99] code gave 'DeveloperFromPackedArray::unpack: Unpacking array in call to HoldForm.', while for Range[100] it gave "DeveloperFromPackedArray::unpack: Unpacking array in call to List". Which seems odd, since myList wasn't packed to begin with. Commented Oct 23, 2018 at 21:51
• Are you trying to find all numbers in a certain range or just a random selection of numbers? I have a method that carries out OddQ on the entire array in {0.683298, Null} EDIT* this timing is baring in mind I have a slower system Commented Oct 23, 2018 at 22:56
• @AwkwardPanda I'm trying to find all numbers in a certain range. My use case is somewhat more convoluted than the example - I'm trying to do something like this: GatherBy[otherList[[i]], Union[Flatten[myList[[{i, #}]],1]] &], and then select the groups of a specific size. The problem doesn't lend itself (as far as I can tell) to straight mapping, and might not be the most efficient way to do it, but when using packed arrays (as I just learned), or when Length[otherList[[i]] ] is less than 100, it works fast enough for my purposes. Commented Oct 23, 2018 at 23:04

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 DeveloperToPackedArray to convert unpacked arrays to packed arrays. Here is a comparison:

m = RandomInteger[10^6, {10^6, 3, 2}];
um = DeveloperFromPackedArray[m];

DeveloperPackedArrayQ /@ {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.

• I believe this explains everything. My function f is a pure function which changes on every different Gather call. If Mathematica is compiling a different pure function every time, that explains why I didn't just see a single bit of overhead followed by fast responses. Thank you! Commented Oct 24, 2018 at 0:33

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}
`