# Replacement of Do loops

Here is an example code (to calculate the so called structure factor, please see the comment below the question) for a small data set of 2d coordinates. In reality I have a factor of up to 100 times more coordinates.

Code with AbsoluteTiming:

coordinates = Get@"https://pastebin.com/raw/wFxJ1m4U";

{lX, lZ} = {100, 100};

kX = 2 Pi/lX*(Range[lX] - lX/2) // N;
kZ = 2 Pi/lZ*(Range[lZ] - lZ/2) // N;

cx = Flatten[coordinates[[All, 1]]];
cz = Flatten[coordinates[[All, 2]]];

rX = Flatten@Outer[Differences[{##}] &, cx, cx];
rZ = Flatten@Outer[Differences[{##}] &, cz, cz];

nP = Length[coordinates];

s = Array[0 &, {lZ, lX}];

Do[
Do[
s[[j, i]] = 1/nP*Total[Cos[kX[[i]]*rX + kZ[[j]]*rZ]]
, {i, 1, lX}
]
, {j, 1, lZ}
]; // AbsoluteTiming

{31.3555, Null}


How can I speed up this code?

A huge problem occurs for many coordinates: for about 2400 coordinates the calculation takes about 3 days!

Here are 2403 coordinates: https://pastebin.com/raw/0t1MBAgX . Please use for this test {lX, lZ} = {1600, 1200}

Here is the code for this large data set: https://pastebin.com/raw/MH5Pb4h0

Can the double loop be compiled?

• Can you add information about what this code is supposed to do to your question? – Carl Lange May 8 '19 at 8:17
• @Carl Lange: It calculates the structures factor of a crystaline structure (S(k) algorithm from a paper by Henrich et al., arXiv:1001.3342v1, equation 2.6, pdf file: arxiv.org/pdf/1001.3342.pdf) – lio May 8 '19 at 8:38
• On my computer your first (non-parallel) code only takes 14.4 seconds. Replacing the last double loop with s = Outer[Total[Cos[#1 + #2]] &, KroneckerProduct[kZ, rZ], KroneckerProduct[kX, rX], 1]/nP cuts this down to 8.6 seconds. I have no idea where your excessive runtime of 3714.41 seconds comes from. – Roman May 8 '19 at 8:51
• I agree with @Roman, on my computer your code took about 14s as well. – Turgon May 8 '19 at 8:52
• Have you run your own code in a fresh kernel? What is your Mathematica version? Have you tried replacing the cosine sum by a Fourier or FourierDCT call? – Roman May 8 '19 at 9:01

Table works a little bit better than Do:

s = Table[Total[Cos[kX[[i]]*rX + kZ[[j]]*rZ]], {j, 1, lZ}, {i, 1, lX}]/nP;


Much more speed can be gained by first rewriting the cosines as imaginary exponentials, and then rewriting exponentials of sums as products of exponentials. As a result, rX and rZ aren't needed, and the calculation takes 0.03 seconds, or just over one second with the large data set:

s = Abs[Exp[I*KroneckerProduct[kZ, cz]].Exp[I*KroneckerProduct[cx, kX]]]^2/nP; //
AbsoluteTiming // First
(* 1.22367 *)


A little more speed can be gained with this trick of using an internal function to calculate Abs[...]^2:

s = InternalAbsSquare[
Exp[I*KroneckerProduct[kZ, cz]].Exp[I*KroneckerProduct[cx, kX]]]/nP; //
AbsoluteTiming // First
(* 1.09222 *)


A bit more speed could be gained by replacing the Kronecker products by a NestList, as the values in kX and kZ are uniformly spaced and thus the many exponentials can be replaced by repeated multiplications of a single exponential. This is numerically less stable though, and in any case at this point 50% of the time are spent in the Dot product which is unavoidable, so the best remaining speedup you could hope for is another factor of two.

Lesson of the day: don't think about compiling before you've tried rewriting the algorithm.

## mathematical derivation

Starting from Eq. (2.6) of https://arxiv.org/pdf/1001.3342.pdf

$$S_{\vec{q}}=\frac{1}{N}\sum_{i j} \cos[\vec{q}\cdot(\vec{r}_i-\vec{r}_j)] = \frac{1}{N}\Re\sum_{i j} e^{\text{i} \vec{q}\cdot(\vec{r}_i-\vec{r}_j)} = \frac{1}{N}\Re\sum_{i j} e^{\text{i} \vec{q}\cdot\vec{r}_i}e^{-\text{i}\vec{q}\cdot\vec{r}_j}\\ = \frac{1}{N}\Re\left[\sum_i e^{\text{i} \vec{q}\cdot\vec{r}_i}\right]\left[\sum_j e^{-\text{i} \vec{q}\cdot\vec{r}_j}\right] = \frac{1}{N}\left|\sum_i e^{\text{i} \vec{q}\cdot\vec{r}_i}\right|^2 = \frac{1}{N}\left|\sum_i e^{\text{i} q_xr_{i,x}}e^{\text{i} q_zr_{i,z}}\right|^2$$

From there it's a matter of calculating all the possible combinations of $$e^{\text{i} q_xr_{i,x}}$$ and $$e^{\text{i} q_zr_{i,z}}$$ with all desired $$\vec{q}$$-values (in kX and kZ), and then summing over them. These operations are most efficiently done with vector operations like KroneckerProduct, listable Exp, Dot-products, and finally a listable Abs or AbsSquare.

• Can you check if this result is equal to my s. I named your result as s2 and my as s. And if I check s==s2 it gives False. Please see also (s[[#]] == s2[[#]]) & /@ Range[lZ]. Nevertheless this difference is neglectable. – lio May 8 '19 at 12:10
• Again, please make sure you use a fresh kernel to run calculations. – Roman May 8 '19 at 12:14
• @b3m2a1 see my updated solution: if you rewrite these cosines as exponentials, then you only need to KroneckerProduct with cx instead of rX, and the result easily fits into memory. In effect, rX contains cx twice, and this double-use is replaced by an Abs[...]^2 call. – Roman May 10 '19 at 12:31
• You're very welcome. I hope this allows you to advance in your research. – Roman May 15 '19 at 7:48
• @cphys or simply s-s2 // Abs // Max` – Roman May 15 '19 at 20:46