Issue with very large lists in Mathematica

Question

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.

Answer 1

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.

Original answer

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.

Answer 2

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.

Answer 3

... 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}]

instead of

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.