# Speed up 4D matrix/array generation

I have to fill a 4D array, whose entries are $\mathrm{sinc}\left[j(a-b)^2+j(c-d)^2-\phi\right]$ for a fixed value of $\phi$ (normally -15) and a fixed value of $j$ (normally about 0.00005). The way I'm trying to speed this up is by first creating a "lookup matrix".

lookup = ParallelTable[Sinc[0.00005 (i^2 + j^2) - 15.],
{i, -2 lim, 2 lim, step},
{j, -2 lim, 2 lim, step}] // Chop


And then use it in this way:

matrix = ParallelTable[lookup[[(a - b) + dim, (c - d) + dim]], range]


where range stands for {a, -lim/step, lim/step, 1}, {b, -lim/step, lim/step, 1}, ... (hence matrix is 4-dimensional) and dim = 2lim/step + 1

In this way (I normally choose lim/step to be an integer) up to values of lim/step $\simeq$ 25 my six-year-old MacBook can do it, but for values higher than that it just freezes up.

The same entries in lookup are read several hundreds if not thousands of times. It seems that there should be a more clever way of doing this. The slow step is the creation of matrix, not the creation of lookup. Am I missing an obvious (or not so obvious) way of optimising this?

-
What does it do with just plain Table? – rcollyer Feb 18 '13 at 17:06
I will try now. – Ziofil Feb 18 '13 at 17:08
It's roughly 1.5x slower. – Ziofil Feb 18 '13 at 17:11
compilation so probably the way to go. look for Compile in the help. – Ajasja Feb 18 '13 at 17:13
I will try it now, but I read that Table and ParallelTable automatically compile their input, so there shouldn't be a performance gain. Besides, the creation of lookup is not the problem, it's reading its values over and over to create matrix. – Ziofil Feb 18 '13 at 17:15

You can reduce a lot of the computations by exploiting the symmetry in the problem. Observe the following example:

Notice that the 4D matrix is actually a Toeplitz matrix where each element is a 2D matrix, which themselves are Toeplitz matrices. So you really need to only compute the first row in the top line of 2D matrices (here, 9 elements of total 81) and let ToeplitzMatrix do rest of the work.

To get the first "row", we need to fix two variables at the lower limit and let the other two vary across the entire range. Then we convert each to a ToeplitzMatrix:

l = ToeplitzMatrix /@ With[{b = -1, c = -1},
Table[f[j (a - b)^2 + j (c - d)^2 - k], {a, -1, 1}, {d, -1, 1}]]


Now we have the first row of the final 4D matrix. To convert this into a Toeplitz matrix, we'll have to slightly redefine ToeplitzMatrix to accept a list of lists as input:

Block[{ToeplitzMatrix},
ToeplitzMatrix[{x__List}] := With[{list = {x} /. List -> \[FormalCapitalL]},
ToeplitzMatrix[List @@ list]];
ToeplitzMatrix[l]
] /. \[FormalCapitalL] -> List


You'll see that the above is the same as the naïve approach with Table. You can replace f in the above with your function and the appropriate limits.

-
This is interesting, but how well does it actually perform? – whuber Feb 18 '13 at 19:11
@whuber For the OP's lim/step of 25, it is slightly faster — ~1.1 s for naïve Table vs. ~0.8 s with mine. However, when you increase the dimensions to lim/step = 50, it is much faster — ~285 s for Table vs ~13.8 s for mine. Although the OP uses Sinc, which is not a computationally expensive function, I can imagine this approach providing a much bigger speed gain when the function f is expensive. I didn't test either approach with parallelization, since I don't have free kernels right now. – R. M. Feb 18 '13 at 19:23
And once again, thinking outside the box saves the day! ^^ – Ziofil Feb 18 '13 at 19:28

You might indeed be missing something.

The natural way to create such arrays is with Array (or ParallelArray) and the natural way to create such a lookup table (which is a great idea) is by "memoizing" the functional argument of Array.

Notice that you could change $j$ and $\phi$ so that the indexes $a,b,c,d$ range from $1$ up to some integer; you are concerned about the case where the upper limit is $50$ or larger. Anticipating this, let's begin with the core calculation, memoized:

g[m_, n_, j_] := g[m, n, j] =  Sinc[m j - n];


Here, m collects the coefficients of $j$ and n plays the role of $\phi$. The value of g is first sought among stored values g[m, n, j]; if not found among them, the value is computed and stored: there's your lookup table, automatically produced and used. Let's also encapsulate the calculation of the array entries in a separate function:

f[a_, b_, c_, d_, {\[Phi]_, j_}] := g[(a - b)^2 + (c - d)^2, \[Phi], j];


Here is an example of its use by ParallelArray to create $100^4$ entries:

a = ParallelArray[f[##, {-15, .00005}] &, {100, 100, 100, 100}]; // AbsoluteTiming


{150.0035797, Null}

(Four cores.) Additional RAM usage was 3.5 GB to hold the result.

### Edit

Sinc is quick to compute, so little is to be gained by optimizing its calculation. However, storing the values of g in an array speeds up the calculation by an order of magnitude:

a = Block[{d1 = 100, d2 = 100, d3 = 100, d4 = 100, j = .00005, n = -15, g},
g = Sinc[# j - n] & /@ Range[0, Max[d1, d2]^2 + Max[d3, d4]^2];
ParallelArray[g[[(#1 - #2)^2 + (#3 - #4)^2 + 1]] &, {d1, d2, d3, d4}]
]; // AbsoluteTiming


{9.6095496, Null}

Compile actually improves things, speeding up the code about another factor of four. Here is a general-purpose solution in which the values of g are cached first and then passed to a compiled version of Array:

 f = Compile[{{g, _Real, 1}, {d, _Integer, 1}},
Table[g[[(i - j)^2 + (k - l)^2 + 1]],
{i, 1, d[[1]]}, {j, 1, d[[2]]}, {k, 1, d[[3]]}, {l, 1, d[[4]]}]
];


All f does is to create the correct tensor structure out of the entries in the array g which have been passed to it. Other structures would be amenable to a similar treatment.

I timed the pre-computation and array creation together:

d = {d1, d2, d3, d4} = {100, 100, 100, 100};
{j0, n} = {.00005, -15};
(g = Table[Sinc[i j0 - n], {i, 0, Max[d1, d2]^2 + Max[d3, d4]^2}] // N;
a = f[g, d]); // AbsoluteTiming


$\{10.9056238,\text{Null}\}$

This was executed on one kernel and used only 800 MB RAM.

-
Thanks whuber, I always overlook memoization! However, the method you describe in your answer is not faster than the method I'm using at the moment. Perhaps you meant using a combination of the two? – Ziofil Feb 18 '13 at 18:10
With the phrase "just freezes up" your question appears to focus on the failure (or unscalability) of your code when the matrix gets big; that's why I provided such a large example. (The case where lim/step is about $25$ corresponds to array dimensions of about $51$ each; using $100$ therefore creates an array almost $16$ times larger.) I was making no effort to provide anything faster (but worked a little to avoid making it any slower). – whuber Feb 18 '13 at 18:14
OK, I'm surprised to see how poorly memoization works compared to indexing into an array. Of course the indexing will be fast--it is made possible by viewing the array as a function of its integer indexes rather than as a function of non-integral steps--but one would think that the hashing needed to look up something like g[m] would be almost as fast as the direct indexing in g[[m]]. – whuber Feb 18 '13 at 19:10
That is surprising indeed! I will have to review some of my older code in view of this. Thanks – Ziofil Feb 18 '13 at 19:14
Another surprise is that compilation actually speeds up Table, enough to be worth consideration in many instances. – whuber Feb 18 '13 at 19:53