Multi-dimensional array lookup vs flat array lookup with compute

Question

During Monte-Carlo (MC) simulations on a latice of dimensions $(L_1,L_2,\ldots,L_D)$, all the site varibales are often stored as a $D$ dimensional array array. To lookup a site varibale, one just needs the $D$ dim. coord coord.

For e.g. on a 2 D lattice, a 2 D array may be made where each element is some site variable.(the dimensionalty or nature of the site variable itself is irrelevant). To lookup a var at coord={x,y} one uses, Extract[array,coord]

Often the simulation is performed in sweeps. Each sweep comprise of $N=L_1L_2\ldots L_D$ random lookups on the array.

This is a lot of lookups, especially when a simulation may have a large no of sweeps.

Since memory is ultimately stored as a sequential flat array, it may be better to perform lookups on a flat array instead.

flatArray=Flatten@array

Now a lookup is simpler

flatArray[position]

As before $N$ random positions are generated per sweep.

1. Is this correct? Are flat array lookups faster than nested ones? (for the test case take $D=3$ on a $48\times 48\times48$ lattice)

However, for the flatened array, the lattice structure is lost. To perform calculations which need that structure, one needs the site's coordinate in the lattice. Even though this can be computed every time, one might as well do a lookup instead.

allCords=Flatten[Array[{##}&,{L1,L2,...,LD}],D-1]
toCoord[pos_]:=allCoords[[pos]]

2. Is it correct to assume that coord compute is slower than coord lookup? Keep in mind that the compute is for a general $D$.

3. With this requirement of needing the lattice coord too (via compute or lookup), is the answer to 1. still the same?

Assume in the above that memory is cheap but time expensive.

Answer 1

First let me say that it might be worthwhile to use Part instead of Extract. But this might require to refactor your code and might lead to set-backs at other places. If your random indices come along as index triples, then Extract is quite a good choice. But if you can generate the linear indices directly, then Part might be a better choice:

n = 48;
m = 1000000;
A = RandomReal[{-1, 1}, {n, n, n}];
a = Flatten[A];

r1 = Extract[A, RandomInteger[{1, n}, {m, 3}]]; // RepeatedTiming // First
r2 = Part[a, RandomInteger[{1, n^3}, {m}]]; // RepeatedTiming // First

0.027

0.011

It might even be worthwhile to skip the generation of the indices by using RandomChoice:

r3 = RandomChoice[a, m]; // RepeatedTiming // First

0.0092

However, this would make it impossible to access the lattice coordinates, so this is probably not what you are looking for.

At 1.: Mathematica's packed arrays are indeed flat lists with the nested lookup implemented precisely as you suggest. So there should basically be no difference between a flat memory call with Part in which you compute the linear index from the coordinates by hand or if you do a nested memory call (and let Mathematica do the computation of the linear index from the coordinates). Of course, in practice, there can however be slight differences in performance (e.g., Mathematica performs also some security checks which take longer or it does some clever caching that might speed up things a bit). Typically, when I Compile the code and use Compile`GetElement instead of Part for memory calls, then using linear or nested indexing makes little difference.

At 2.: As you already seem to know, lookups in memory tend to be more time-intensive than computations in machine integers and machine floats (doubles), because of memory lags and limited bandwidth. (The development of memory did not quite keep track with the development of CPUs in the last couple of years.) So sometimes, it is worthwhile to recompute something with data that is already there instead of having to wait for the result of a memory call. In general, one just has to try several alternative implementations to find the one that performs best on the class of problems that one wants to solve.

At 3.: If you really want to go that route, I would suggest to use Quotient and Mod to compute lattice coordinates from linear indices. (One can also use IntegerDigits (with MixedRadix, if the array's dimensions are not all equal), but it seems to be less performant.) The full answer depends quite a lot on how large the lookup tables are (e.g., whether they would lead to frequent cache misses).

As final remark, let me say that Mathematica is a many-purpose, high-level language. It is not meant for and not capable of producing best possibly optimized high performance code. But Mathematica makes up for it by accelerating prototyping quite a lot. So if you have to solve a problem that is memory bound and if you have to do it as fast as you can, at some point, you should consider to port your algorithm to languages that are closer to the machine, e.g., C, C++, or FORTRAN.