# Issue with very large lists in Mathematica

I have a need to work with very large arrays of data in Mathematica. The length of the list is not known upfront. The typical size is a few hundred thousands elements and could easily go beyond a million elements in some cases. There are two major types of lists I need to work with: 1. list of real numbers, 2. list of some (different) symbolic expressions.

While the list of real numbers is internally converted into packed array and does not seem to cause much performance issues, it is the second type of list (very long list of some different symbolic expressions), which is causing a major problem. It turned out that access time (get or set) to a single element of such lists grows linearly (if not faster) with the length of the list.

I wrote a simple test to evaluate average access speed per element for such lists along with the workaround to deal with the issue: https://www.dropbox.com/s/7j9qlppmrsrmjug/CLM_DynArray_Test_204__008.nb?dl=0 .

I ran the test on Mathematica 7, 9, and 10. Version 7 does not have get issue but still has set issue and both versions 9 and 10 have both get and set issues. This definitely looks like an error in internal Mathematica implementation unless I missed something. I would appreciate any help with that problem.

Here is the excerpt from the code above to illustrate the problem:

ClearAll["Global*"];

nnnList = {1000, 1500, 2000, 3000, 4000, 5000, 7000, 10000, 15000, 20000, 30000,
40000 , 50000, 70000, 100000};

MinReperition = 5;
MaxLength = MinReperition*Max[nnnList];

CreationTiming[nnn_?IntegerQ, useVariables_?BooleanQ] := CreationTiming[nnn,useVariables, True];
CreationTiming[nnn_?IntegerQ, useVariables_?BooleanQ, useRepetitions_?BooleanQ] :=
Module[
{timing, noOfRepetitions, repetitionCount},
noOfRepetitions = If[useRepetitions, Ceiling[MaxLength/nnn], 1];
timing = 10^6*Mean[
Table[
If[useVariables,
Timing[vars = Table[{ToExpression["var" <> ToString[ii]]}, {ii, 1, nnn}];][[1]],
Timing[vars = Table[RandomReal[], {ii, 1, nnn}];][[1]]
],
{repetitionCount, 1, noOfRepetitions}
]
];
Return[timing];
];

AccessTiming[nnn_?IntegerQ, useVariables_?BooleanQ] :=
Module[
{timing, x, noOfRepetitions, repetitionCount},
noOfRepetitions = Ceiling[MaxLength/nnn];
CreationTiming[nnn, useVariables, False];
timing = 10^6*Mean[
Table[Timing[Do[(x = vars[[ii]]), {ii, 1, nnn}];][[1]],
{repetitionCount, 1, noOfRepetitions}]
];
Return[timing];
];

plotOptions = {PlotRange -> All, Frame -> True, GridLines -> Automatic};

Print["Access time per Real Table element in microseconds."];
Print[DiscretePlot[AccessTiming[nnn, False]/nnn, {nnn, nnnList}, Evaluate[plotOptions]]];

Print["Access time per variable Table element in microseconds."];
Print[DiscretePlot[AccessTiming[nnn, True]/nnn, {nnn, nnnList}, Evaluate[plotOptions]]];


nnnList is the list of lengths to be tested for access speed.

For small list lengths it is necessary to repeat the test more times so that to obtain a better average result. For the largest list length the test will be performed MinReperition times.

CreationTiming measures time in microseconds, which is necessary to create a Table of length nnn. If useVariables [Equal] False, then the function populates the Table with random real numbers. If useVariables [Equal] True, then the function populates the Table with generated symbolic variables (var1, var2, var3, etc...). useRepetitions [Equal] True (default) will result in multiple runs. useRepetitions [Equal] False will result in a single run.

AccessTiming measures time in microseconds, which is necessary to consecutively get all the elements from a Table of length nnn.

DiscretePlot plots the time per element to access all elements. As the time per elements grows at least linearly the total time to access all elements grows at least quadratically.

• I tried to test access times like this: arr = Table[RandomInteger[10, 5].x^Range[5], {1000000}]; Table[AbsoluteTiming@Do[arr[[i]], {i, Length[arr]}], {10}] I see a linear increase as arr grows. However, the very first time to iterate through the array is much slower than subsequent times. – Szabolcs Sep 9 '14 at 3:35
• Hi, Szabolcs. Thanks for comments. I updated to original post to include the piece of code with the issue of the problem. I thought GPL would be a good choice in case someone find use to the workaround until Wolfram fixes the issue. As the issue has been around since version 7 and actually got worse, I don't expect that anytime soon. I send them an example but it is probably a good idea to send them a shorter version. – Konstantin Konstantinov Sep 9 '14 at 4:17
• @Szabolcs In your code you forgot to divide the AbsoluteTiming by the number of calls. If I do so I get equal timings independently of the length of the list: Table[arr=Table[RandomInteger[10,5].x^Range[5],{nnn}];Mean[Table[AbsoluteTiming[Do[arr[[i]],{i,Length[arr]}]][[1]]/Length[arr],{10}]],{nnn,{10000,100000,1000000}}]. – Alexey Popkov Sep 9 '14 at 8:22
• @AlexeyPopkov I didn't forget, but it was late at night and I wasn't clear enough. My point was that I see a linear increase which is exactly what one would expect. – Szabolcs Sep 9 '14 at 14:22
• Thanks for the example(s). I'll get some of it into a bug report. – Daniel Lichtblau Sep 9 '14 at 16:29

## Update

This should have been obvious to me from the start but of course Part does not hold its first argument. This means that the entire list of Symbols must be checked against the changing value of the Do iterator, which I think easily explains the slow-down observed. Consider:

x = Table[Unique["var"], {100000}];
Do[x, {i, 10000}] // AbsoluteTiming // First

6.605378


That's the time for 100,000 Symbols to be resolved 10,000 separate times. For 20,000 times:

Do[x, {i, 20000}] // AbsoluteTiming // First

13.315762


And for 50,000 Symbols resolved 20,000 times:

x = Table[Unique["var"], {50000}];
Do[x, {i, 20000}] // AbsoluteTiming // First

6.758387


So as I see it there is no unexpected inefficiency or rising complexity here. The reason that the dv[] example in my original answer is fast is that there is only one Symbol (var) to resolve for the entire list at each iteration.

I believe I can at least in part explain what you are observing. First consider that not all symbolic expressions experience the slow-down that you describe:

ClearAll[n, var, re, un, dv, x]

re[] := RandomReal[1, n]
un[] := Table[Unique[], {n}]
dv[] := Table[var[i], {i, n}]

time = Module[{a = #[]}, Do[a[[i]], {i, n}] // AbsoluteTiming // First] &;

n = 10000;
time /@ {re, un, dv}

{0.010000, 0.130000, 0.}

n = 20000;
time /@ {re, un, dv}

{0.030000, 0.520001, 0.010000}

n = 40000;
time /@ {re, un, dv}

{0.030000, 9.990026, 0.020000}


Obviously the indexed objects form does not suffer from this slowdown, yet it is of course symbolic.

Now consider that with every extraction you make you also initiate an evaluation. I postulate that individually resolving each of the Symbols created by Unique takes a significant amount of time. First please recall that Part wraps extracted sequences in the original head of the expression, then observe:

n = 40000;
x = HoldComplete @@ un[];
Do[x[[{i}]], {i, n}] // AbsoluteTiming // First

Table[x[[{i}]], {i, 3}]  (* illustration *)

0.030000

{HoldComplete[$40036], HoldComplete[$40037], HoldComplete[\$40038]}


The same extractions are taking place yet they are suddenly very fast, even with the additional complexity of each element being wrapped in HoldComplete. The difference is that HoldComplete fully prevents the Symbols from being resolved, even for the purpose of UpValues etc.

In summary it is my belief that the performance issue you describe has nothing to do with Mathematica's implementation of arrays and part extraction but is instead determined by the resolution of Symbol definitions.

• This is correct. Related:. There was also a related thread on MSE in the past couple of years, but I cannot locate it at the moment. – Daniel Lichtblau Sep 9 '14 at 16:28
• Actually this was the other related link. Older than I had thought. – Daniel Lichtblau Sep 9 '14 at 16:32
• Mr. Wizard, thanks for comments and references. Unfortunately, I have to disagree. As comment field is only 500 characters I am posting a detailed comment below. – Konstantin Konstantinov Sep 9 '14 at 18:01

This looks like a bug, but presumably a workaround is to use a hash table. That is, in @Szabolcs's code, change the first line to:

Table[arr[i] = RandomInteger[10, 5].x^Range[5], {i, 1, 1000000}];

The second line seems to work reasonably quickly.

• Hi, Igor. Thanks for comments. We are not taking about integers but symbols. With numbers it is just fine. I just posted an update with the relevant piece of the code inline the question. – Konstantin Konstantinov Sep 9 '14 at 4:22
• @KonstantinKonstantinov The elements in the example are symbolic expressions, it's the array INDICES that are integers (but don't have to be, I just made them integers for consistency with Szabolcs. – Igor Rivin Sep 9 '14 at 13:39
• It is unfortunate that a workaround might be needed. But yes, this is a viable one (it had a down vote, I know not why). Another approach could be fashioned along lines noted here – Daniel Lichtblau Sep 9 '14 at 17:03
• Hi, Igor. I just posted some more notes above in response to Mr. Wizard comments. Basically if you initialize arr[i] with numbers then the list will contain just numbers and not arr[i]. However, if arr[i] contains some symbolic expressions (which are not numbers) then we are back to square one and the issue comes back. – Konstantin Konstantinov Sep 9 '14 at 18:38

... so, here is why I disagree with you, Mr. Wizard.

If you use an indexed object, which is actually a function definition (like varFunc[ii] instead of ToExpression["var" <> ToString[ii]]) then there is no value associated with varFunc[ii]. Yes, it is true that a long list can then be accessed fast. Unfortunately in that form it is practically useless because varFunc[ii] does not have any values for any ii. However, as soon as we assign values to varFunc[ii], for example by modifying CreationTiming as follows:

CreationTiming[nnn_?IntegerQ, useVariables_?BooleanQ,useRepetitions_?BooleanQ] :=
Module[{timing, noOfRepetitions, repetitionCount},
noOfRepetitions = If[useRepetitions, Ceiling[MaxLength/nnn], 1];

Do[varFunc[ii] := ToExpression["var" <> ToString[ii]], {ii, 1,nnn}];

timing =10^6*Mean[
Table[
If[useVariables,
Timing[vars := Table[varFunc[ii], {ii, 1, nnn}];][[1]],
Timing[vars = Table[RandomReal[], {ii, 1, nnn}];][[1]]
], {repetitionCount, 1, noOfRepetitions}
]
];

Return[timing];
];


then the issue is back again. Using delayed assignment in

 varFunc[ii] := ToExpression["var" <> ToString[ii]]


or in

vars := Table[varFunc[ii], {ii, 1, nnn}]


does not make any difference. The only time when delayed assignment does make a difference if we write:

Do[(x := vars[[ii]]), {ii, 1, nnn}]


Do[(x = vars[[ii]]), {ii, 1, nnn}]


in AccessTiming. But then we are not actually evaluating the value of x! However, as soon as we attempt to extract the value from x, like

Do[(x := vars[[ii]]; y=x), {ii, 1, nnn}]


then the issue comes back. have a bad feeling that the same will be with HoldComplete but I have not tested that thoroughly yet

The workaround, which I posted as a link above does not have that issue even though it is probably not the most efficient Mathematica code. So I disagree that individually resolving all the Symbols created by Unique (or via ToExpression, which I prefer) is very slow. Rather it is implementation of Part, which by the way did not have access (get) issue in version 7 and only had set issue.

• I think the way to use an indexed object in this setting is as follows. In the create... routine, have this definition Do[varFunc[ii] = ToExpression["var" <> ToString[ii]], {ii, 1,nnn}]; and later If[useVariables, Timing[vars = varFunc][[1]].... In the access... function, use Do[(x = vars[ii]), {ii, 1, nnn}];]`. Notice the single brackets in that last. – Daniel Lichtblau Sep 9 '14 at 19:08
• Daniel, thanks. It works! I will now mark the issue as resolved. – Konstantin Konstantinov Sep 9 '14 at 20:09
• @DanielLichtblau. Thanks again for your help. Your answer is the most applicable to the case but I can't mark it as answer because it is in the comments. If you could, please, put it s a separate answer I will then just mark it right away. Thanks. – Konstantin Konstantinov Sep 24 '14 at 17:13
• I'm glad it was useful. But I'm afraid I no longer can figure out enough detail to put this together into a coherent piece of code. For what it's worth, I did manage to put together a bug report for an undesirably slow case (one showing the wrong complexity, more or less). – Daniel Lichtblau Sep 24 '14 at 18:17